C*************************************************************** C Test program for subroutine LSEM3 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/LSEM3_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 LSEM3. The program computes C the solution to a Laplace problem exterior to two square plates by C the shell panel method. One plate is assigned the potential one C and the other is at zero potential. The problem may be thought of as C a pair of plates of a capacitor with one plate at zero volts, the C other at 1 volt. C C Background C ---------- C C We wish to solve the Laplace equation C C __ 2 C \/ {\phi} = 0 C C C For the exterior problem, the domain lies exterior to the shells C or plates. The shell condition may be Dirichlet, Robin or C Neumann. It is assumed to have the following general form C C a(q) {\delta}(q) + b(q) {\nu}(q) = f(q) C C A(q) {\Phi}(q) + B(q) V(q) = F(q) C C where {\delta}(q)={\phi+}(q) - {\phi-}(q), C {\nu}(q) = {\nu+}(q) - {\nu-}(q) C {\Phi}(q) = ({\phi+}(q)+{\phi-}(q))/2 C V(q) = ({v+}(q) + {v-}(q))/2 C and {\phi}(q) is the velocity potential at the point q on S, v(q) C is the derivative of {\phi} with respect to the outward normal to S C at q and {\alpha}, {\beta} and f are real-valued functions defined C on S. C C Subroutine LSEM3 accepts a description of the shells and the position C of the exterior points where the solution ({\phi}) C is sought, the shell condition and returns the solution ({\delta}, C {\Phi}, {\nu} and V) on {\Pi} and the value of {\phi} at the exterior C points. C C The test problems C ----------------- C C In this test the domain is the exterior to two square plates of side C 0.1m and 0.01m apart. The solution to the problem with a Dirichlet C shell condition (a=1, b=0, f=0) on both plates and (A=1, B=0, C F=0) on the lower plate and (A=1, B=0, F=1) on the upper plate. C C Each plate is described by a set of 32 planar triangular panels C of equal size. The shell solution points are the centres of the C panels. The solution is sought at the exterior points C (0.05,0.05,z), the axis through the plates with z=0.001,0.003, C 0.005, 0.007, 0.009. C---------------------------------------------------------------------- C The PARAMETER statement C ----------------------- C There are four components in the PARAMETER statement. C integer MAXNS : 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 LEBEM3: Subroutine for solving the exterior Laplace C equation (file LEBEM3.FOR contains the subroutine LEBEM3) C subroutine H3LC: Returns the individual discrete Laplace integral C operators. (file H3LC.FOR contains H3LC and subordinate routines) C subroutine GLS: Solves a general linear system of equations. C (file GLS.FOR contains CGSL and subordinate routines) C file GEOM3D.FOR contains the set of relevant geometric subroutines C The program PROGRAM LSEM3T 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=100) C Limit on the number of vertices (equal to the number of panels) INTEGER MAXNV PARAMETER (MAXNV=100) C Limit on the number of points exterior to the plate, where C potential properties are sought INTEGER MAXNPE PARAMETER (MAXNPE=10) C Constants C Real scalars: 0, 1, 2, pi REAL*8 ZERO,ONE,TWO,THREE,FOUR,PI C Geometrical description of the plate(s) C Number of panels and counter INTEGER NPI,IPI 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,3) C The three nodes that define each panel on the boundaries INTEGER PIELV(MAXNPI,3) C The points exterior to the plate(s) where the potential C properties 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 C Coordinates of the exterior points REAL*8 PEXT(MAXNPE,3) C Shell Conditions C a(p)+b(p)=f(p) C function a at the central points REAL*8 PIA(MAXNPI) C function b at the central points REAL*8 PIB(MAXNPI) C function f at the central points REAL*8 PIF(MAXNPI) C A(p)+B(p)=F(p) C function A at the central points REAL*8 PIAA(MAXNPI) C function B at the central points REAL*8 PIBB(MAXNPI) C function F at the central points REAL*8 PIFF(MAXNPI) C Incident field C The incident potential {\phi} at the central points REAL*8 PIIPHI(MAXNPI) C The derivative of the incident potential {\phi} at the C central points REAL*8 PIIVEL(MAXNPI) C The incident potential {\phi} at the exterior points REAL*8 PEIPHI(MAXNPE) C Validation and control parameters LOGICAL LSOL LOGICAL LVALID REAL*8 EGEOM C Solution REAL*8 PHIDIF(MAXNPI) REAL*8 PHIAV(MAXNPI) REAL*8 VELDIF(MAXNPI) REAL*8 VELAV(MAXNPI) REAL*8 PEPHI(MAXNPE) 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,z coordinates INTEGER ICOORD C Counter through the exterior points INTEGER IPE C Capacitor C Side lenght REAL*8 SQSIDE C Distance between plates REAL*8 GAP C INITIALISATION C -------------- C Set constants ZERO=0.0D0 ONE=1.0D0 TWO=2.0D0 THREE=3.0D0 FOUR=4.0D0 PI=4.0D0*ATAN(ONE) C Describe the nodes and panels that make up the plate C :The two 1x1 plates are each divided into 32 panels, C : (NPI=64). VERTEX and PIELV are defined anti-clockwise around C : the plate when viewed from above so that the normal to the C : plate is assumed to be upward C :Set up nodes C : Set NPI, the number of panels NPI=64 C : Set NV, the number of vertices (equal to the number of panels) NV=50 C : Set coordinates of the nodes C : Set up VERTEX, the coordinates of the vertices of the panels DATA ((VERTEX(IV,ICOORD),ICOORD=1,3),IV=1,50) C Plate 1 * / 0.000D0, 0.000D0, 0.000D0, * 0.250D0, 0.000D0, 0.000D0, * 0.500D0, 0.000D0, 0.000D0, * 0.750D0, 0.000D0, 0.000D0, * 1.000D0, 0.000D0, 0.000D0, * 0.000D0, 0.250D0, 0.000D0, * 0.250D0, 0.250D0, 0.000D0, * 0.500D0, 0.250D0, 0.000D0, * 0.750D0, 0.250D0, 0.000D0, * 1.000D0, 0.250D0, 0.000D0, * 0.000D0, 0.500D0, 0.000D0, * 0.250D0, 0.500D0, 0.000D0, * 0.500D0, 0.500D0, 0.000D0, * 0.750D0, 0.500D0, 0.000D0, * 1.000D0, 0.500D0, 0.000D0, * 0.000D0, 0.750D0, 0.000D0, * 0.250D0, 0.750D0, 0.000D0, * 0.500D0, 0.750D0, 0.000D0, * 0.750D0, 0.750D0, 0.000D0, * 1.000D0, 0.750D0, 0.000D0, * 0.000D0, 1.000D0, 0.000D0, * 0.250D0, 1.000D0, 0.000D0, * 0.500D0, 1.000D0, 0.000D0, * 0.750D0, 1.000D0, 0.000D0, * 1.000D0, 1.000D0, 0.000D0, C Plate 2 * 0.000D0, 0.000D0, 1.0000D0, * 0.250D0, 0.000D0, 1.0000D0, * 0.500D0, 0.000D0, 1.0000D0, * 0.750D0, 0.000D0, 1.0000D0, * 1.000D0, 0.000D0, 1.0000D0, * 0.000D0, 0.250D0, 1.0000D0, * 0.250D0, 0.250D0, 1.0000D0, * 0.500D0, 0.250D0, 1.0000D0, * 0.750D0, 0.250D0, 1.0000D0, * 1.000D0, 0.250D0, 1.0000D0, * 0.000D0, 0.500D0, 1.0000D0, * 0.250D0, 0.500D0, 1.0000D0, * 0.500D0, 0.500D0, 1.0000D0, * 0.750D0, 0.500D0, 1.0000D0, * 1.000D0, 0.500D0, 1.0000D0, * 0.000D0, 0.750D0, 1.0000D0, * 0.250D0, 0.750D0, 1.0000D0, * 0.500D0, 0.750D0, 1.0000D0, * 0.750D0, 0.750D0, 1.0000D0, * 1.000D0, 0.750D0, 1.0000D0, * 0.000D0, 1.000D0, 1.0000D0, * 0.250D0, 1.000D0, 1.0000D0, * 0.500D0, 1.000D0, 1.0000D0, * 0.750D0, 1.000D0, 1.0000D0, * 1.000D0, 1.000D0, 1.0000D0/ C Here they are set up as 0.1mx0.1m plates 0.01m apart. SQSIDE=0.1D0 GAP=0.01D0 DO IPI=1,NV VERTEX(IPI,1)=VERTEX(IPI,1)*SQSIDE VERTEX(IPI,2)=VERTEX(IPI,2)*SQSIDE VERTEX(IPI,3)=VERTEX(IPI,3)*GAP END DO 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 : plate into the potential domain. DATA ((PIELV(IPI,ICOORD),ICOORD=1,3),IPI=1,64) C Plate 1 * / 1, 2, 6, 2, 7, 6, 2, 3, 8, 2, 8, 7, * 3, 4, 8, 4, 9, 8, 4, 5,10, 4,10, 9, * 6, 7,12, 6,12,11, 7, 8,12, 8,13,12, * 8, 9,14, 8,14,13, 9,10,14, 10,15,14, * 11,12,16, 12,17,16, 12,13,18, 12,18,17, * 13,14,18, 14,19,18, 14,15,20, 14,20,19, * 16,17,22, 16,22,21, 17,18,22, 18,23,22, * 18,19,24, 18,24,23, 19,20,24, 20,25,24, C Plate 2 * 26,27,31, 27,32,31, 27,28,33, 27,33,32, * 28,29,33, 29,34,33, 29,30,35, 29,35,34, * 31,32,37, 31,37,36, 32,33,37, 33,38,37, * 33,34,39, 33,39,38, 34,35,39, 35,40,39, * 36,37,41, 37,42,41, 37,38,43, 37,43,42, * 38,39,43, 39,44,43, 39,40,45, 39,45,44, * 41,42,47, 41,47,46, 42,43,47, 43,48,47, * 43,44,49, 43,49,48, 44,45,49, 45,50,49 / C Set the potential on the plates DO 310 IPI=1,NPI PIA(IPI)=1.0D0 PIB(IPI)=0.0D0 PIF(IPI)=0.0D0 PIAA(IPI)=1.0D0 PIBB(IPI)=0.0D0 C Plate 1 IF (IPI.LE.32) PIFF(IPI)=0.0D0 C Plate 2 IF (IPI.GE.33) PIFF(IPI)=1.0D0 310 CONTINUE DO 320 IPI=1,NPI PIIPHI(IPI)=0.0D0 PIIVEL(IPI)=0.0D0 320 CONTINUE DO 330 IPE=1,NPE PEIPHI(IPI)=0.0D0 330 CONTINUE C Set NPE NPE=5 C Set coordinates of chosen exterior points DO 300 IPE=1,NPE PEXT(IPE,1)=SQSIDE/2.0D0 PEXT(IPE,2)=SQSIDE/2.0D0 PEXT(IPE,3)=GAP*DFLOAT(2*IPE-1)/DFLOAT(2*NPE) 300 CONTINUE C :Switch for particular solution LSOL=.TRUE. C :Switch on the validation of LSEM3 LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 CALL LSEM3(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,PEPHI, * WKSPC1,WKSPC2,WKSPC3,WKSPC4,WKSPC5, * WKSPC6,WKSPC7,WKSPC8,WKSPC9,WKSPC0) 999 FORMAT(3F14.4) OPEN(10,FILE='LSEM3.OUT') WRITE(10,*) 'LSEM3: Test Problem Results' WRITE(10,*) ' z computed pot exact pot' DO 600 IPI=1,NPE WRITE(10,999) PEXT(IPI,3),PEPHI(IPI),PEXT(IPI,3)*100 600 CONTINUE END