C*************************************************************** C Test program for subroutine LSEMA by Stephen Kirkup C*************************************************************** C C Copyright 2001- Stephen Kirkup C John Tyndall Nuclear Research Institute C School of Computing Engineering and Physical Sciences C University of Central Lancashire - www.uclan.ac.uk C Westlakes Campus C Samuel Lindow Building C West Lakes Science and Technology Park C Whitehaven C Cumbria CA24 3JY C United Kingdom C smkirkup@uclan.ac.uk C C This open source code can be found at C www.boundary-element-method.com/fortran/LSEMA_T.FOR C C Issued under the GNU General Public License 2007, see gpl.txt C C Part of the the author's open source BEM packages. C All codes and manuals can be downloaded from C www.boundary-element-method.com C C*************************************************************** C This program is a test for the subroutine LSEMA. The program computes C the solution to the Laplace equation C C __ 2 C \/ {\phi} = 0 C C exterior to a (set of) thin shell(s). C C Subroutine LSEMA accepts a description of the shell of the domain C and the position of the exterior points where the solution ({\phi}) C is sought, the shell condition and returns the solution ({\phi} C and v) on S and the value of {\phi} at the exterior points. C C C The test problems C ----------------- C C In this test the domain is the exterior to two circular plates of C radius 1 (metre) and 0.1 apart. The solution to the problem with a C Dirichlet shell condition ({a}=1, {b}=0) with C {\phi}=f=1.0 on the upper plate and {\phi}=f=0.0 on the lower plate. C C We would expect a simple gradient in potential between the two C plates and this is investigated in this test program. C C The plates are described by a set of NPI=20 panels of equal length C along the generator. The shell solution points are the centres of the C generator of the panels. C The solution is sought the points (0,0,0.025), (0,0,0.05), C (0,0,0.075). C---------------------------------------------------------------------- C The PARAMETER statement C ----------------------- C There are four components in the PARAMETER statement. C integer MAXNPI : The limit on the number of shell panels. C integer MAXNV : The limit on the number of vertices. C integer MAXNPE : The limit on the number of exterior points. C External modules related to the package C --------------------------------------- C subroutine LSEMA: Subroutine for solving the exterior Laplace C equation. (file LSEMA.FOR contains LSEMA and subordinate routines) C subroutine L3ALC: Returns the individual discrete Laplace integral C operators. (file L3ALC.FOR contains L3ALC and subordinate routines) C subroutine GLS: Solves a general linear system of equations. C (file GLS.FOR contains CGSL and subordinate routines) C file GEOM2D.FOR contains the set of relevant geometric subroutines C The program PROGRAM LSEMAT IMPLICIT NONE C VARIABLE DECLARATION C -------------------- C PARAMETERs for storing the limits on the dimension of arrays C Limit on the number of panels INTEGER MAXNPI PARAMETER (MAXNPI=32) C Limit on the number of vertices (equal to the number of panels) INTEGER MAXNV PARAMETER (MAXNV=MAXNPI) C Limit on the number of points exterior to the shell, where C the solutions are sought INTEGER MAXNPE PARAMETER (MAXNPE=6) C Constants C Real scalars: 0, 1, 2, pi REAL*8 ZERO,ONE,TWO,FOUR,PI C Geometrical description of the shell(ies) C Number of panels and counter INTEGER NPI,IS C Number of central points (on S) and counter INTEGER NSP,ISP C Number of vetices and counter INTEGER NV,IV C Index of nodal coordinate for defining boundaries (standard unit is C metres) REAL*8 VERTEX(MAXNV,2) C The two nodes that define each panel on the boundaries INTEGER PIELV(MAXNPI,2) C The points exterior to the shell(ies) where the solutions C are sought and the directional vectors at those points. C [Only necessary if an exterior solution is sought.] C Number of exterior points and counter INTEGER NPE,IPE C Coordinates of the exterior points REAL*8 PEXT(MAXNPE,2) C Data structures that are used to define each test problem in turn C and are input parameters to LSEMA. C PIA(j) is assigned the value of {a} at the centre of the C j-th panel. REAL*8 PIA(MAXNPI) C PIB(j) is assigned the value of {b} at the centre of the C j-th panel. REAL*8 PIB(MAXNPI) C PIF(j) is assigned the value of f at the centre of the j-th panel. REAL*8 PIF(MAXNPI) C PIA(j) is assigned the value of A at the centre of the C j-th panel. REAL*8 PIAA(MAXNPI) C PIB(j) is assigned the value of B at the centre of the C j-th panel. REAL*8 PIBB(MAXNPI) C PIF(j) is assigned the value of F at the centre of the j-th panel. REAL*8 PIFF(MAXNPI) C The incident field C The incident potential on the shell(s) REAL*8 PIIPHI(MAXNPI) C The derivative of the incident potential on the shell(s) REAL*8 PIIVEL(MAXNPI) C The incident potential at the exterior points REAL*8 PEIPHI(MAXNPI) C Validation and control parameters for LSEMA C Switch for particular solution LOGICAL LSOL C Validation switch LOGICAL LVALID C The maximum absolute error in the parameters that describe the C geometry of the shell. REAL*8 EGEOM C Output from subroutine LSEMA C The average potential at the centres of the panels REAL*8 PHIAV(MAXNPI) C The difference in potential at the centres of the panels REAL*8 PHIDIF(MAXNPI) C The average derivative of the potential at the centres of the C panels REAL*8 VELAV(MAXNPI) C The difference in the derivative of the potential at the centres C of the panels REAL*8 VELDIF(MAXNPI) C Workspace for LSEMA C Working space REAL*8 WKSPC1(2*MAXNPI,2*MAXNPI) REAL*8 WKSPC2(2*MAXNPI,2*MAXNPI) REAL*8 WKSPC3(2*MAXNPI) REAL*8 WKSPC4(2*MAXNPI) REAL*8 WKSPC5(2*MAXNPI) REAL*8 WKSPC6(2*MAXNPI) REAL*8 WKSPC7(2*MAXNPI) REAL*8 WKSPC8(2*MAXNPI) REAL*8 WKSPC9(2*MAXNPI) LOGICAL WKSPC0(2*MAXNPI) C Counter through the x,y coordinates INTEGER ICOORD C The coordinates of the centres of the panels REAL*8 SELCNT(MAXNPI,2) REAL*8 EPS C INITIALISATION C -------------- C Set constants ZERO=0.0D0 ONE=1.0D0 TWO=2.0D0 FOUR=4.0D0 PI=4.0D0*ATAN(ONE) EPS=1.0E-10 C Describe the nodes and panels that make up the shell C :The plates are circles, each divided into 10 panels. VERTEX and C : PIELV are defined out from the centres of the plates so that the C : normal to the plates is assumed to be upward (+z direction). C :Set up nodes C : Set NPI, the number of panels NPI=20 C : Set NV, the number of vertices (equal to the number of panels) NV=22 C : Set coordinates of the nodes DATA ((VERTEX(IV,ICOORD),ICOORD=1,2),IV=1,22) * / 0.000D0, 0.000D0, * 0.100D0, 0.000D0, * 0.200D0, 0.000D0, * 0.300D0, 0.000D0, * 0.400D0, 0.000D0, * 0.500D0, 0.000D0, * 0.600D0, 0.000D0, * 0.700D0, 0.000D0, * 0.800D0, 0.000D0, * 0.900D0, 0.000D0, * 1.000D0, 0.000D0, * 0.000D0, 0.100D0, * 0.100D0, 0.100D0, * 0.200D0, 0.100D0, * 0.300D0, 0.100D0, * 0.400D0, 0.100D0, * 0.500D0, 0.100D0, * 0.600D0, 0.100D0, * 0.700D0, 0.100D0, * 0.800D0, 0.100D0, * 0.900D0, 0.100D0, * 1.000D0, 0.100D0/ C :Describe the panels that make up the two shells C : Set NPI, the number of panels NPI=20 C : Set nodal indices that describe the panels of the shells. C : The indices refer to the nodes in VERTEX. The order of the C : nodes in PIELV dictates that the normal is outward from the C : shell into the domain. DATA ((PIELV(IS,ICOORD),ICOORD=1,2),IS=1,20) * / 1, 2, * 2, 3, * 3, 4, * 4, 5, * 5, 6, * 6, 7, * 7, 8, * 8, 9, * 9, 10, * 10, 11, * 12, 13, * 13, 14, * 14, 15, * 15, 16, * 16, 17, * 17, 18, * 18, 19, * 19, 20, * 20, 21, * 21, 22 / C Set the centres of the panels, the central points DO IS=1,NPI SELCNT(IS,1)=(VERTEX(PIELV(IS,1),1) * +VERTEX(PIELV(IS,2),1))/TWO SELCNT(IS,2)=(VERTEX(PIELV(IS,1),2) * +VERTEX(PIELV(IS,2),2))/TWO END DO C Set the points in the domain where the solutions are sought, PEXT. NPE=3 DATA ((PEXT(IPE,ICOORD),ICOORD=1,2),IPE=1,3) * / 0.000D0, 0.025D0, * 0.000D0, 0.050D0, * 0.000D0, 0.075D0/ C The number of points on the shell is equal to the number of C panels NSP=NPI C :In the test problem the plate 1 has C {\phi}= 1 C : and plate 2 has C {\phi}= 0 DO 250 ISP=1,NSP/2 PIA(ISP)=1.0D0 PIB(ISP)=0.0D0 PIF(ISP)= 0.0D0 PIAA(ISP)=1.0D0 PIBB(ISP)=0.0D0 PIFF(ISP)= 0.0D0 250 CONTINUE DO 260 ISP=NSP/2+1,NSP PIA(ISP)=1.0D0 PIB(ISP)=0.0D0 PIF(ISP)= 0.0D0 PIAA(ISP)=1.0D0 PIBB(ISP)=0.0D0 PIFF(ISP)= 1.0D0 260 CONTINUE C There is no incident field DO 640 ISP=1,NSP PIIPHI(ISP)=0.0D0 PIIVEL(ISP)=0.0D0 640 CONTINUE DO 650 IPE=1,NPE PEIPHI(IPE)=0.0D0 650 CONTINUE C Set up validation and control parameters C :Switch for particular solution LSOL=.TRUE. C :Switch on the validation of LSEMA LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 C Set up particular a,b,f,A,B,F functions for this wavenumber C and type of shell condition CALL LSEMA(MAXNV,NV,VERTEX,MAXNPI,NPI,PIELV, * MAXNPE,NPE,PEXT, * PIA,PIB,PIF,PIAA,PIBB,PIFF, * PIIPHI,PIIVEL,PEIPHI, * LSOL,LVALID,EGEOM, * PHIDIF,PHIAV,VELDIF,VELAV,PEIPHI, * WKSPC1,WKSPC2,WKSPC3,WKSPC4,WKSPC5, * WKSPC6,WKSPC7,WKSPC8,WKSPC9,WKSPC0) C Output the solutions C Open file for the output data OPEN(UNIT=20,FILE='LSEMA.OUT') WRITE(20,*) 'LSEMA: Computed solutions in the domain' WRITE(20,*) WRITE(20,*) ' P(1) P(2) Computed Exact' DO 800 IPE=1,NPE WRITE(20,999) PEXT(IPE,1),PEXT(IPE,2),PEIPHI(IPE), * PEXT(IPE,2)/10.0D0 800 CONTINUE 999 FORMAT(4F10.4) CLOSE(20) END