C*************************************************************** C Test program for subroutine LEBEMA by Stephen Kirkup C*************************************************************** C C Copyright 2001- Stephen Kirkup C School of Computing Engineering and Physical Sciences C University of Central Lancashire - www.uclan.ac.uk C smkirkup@uclan.ac.uk C http://www.researchgate.net/profile/Stephen_Kirkup C C This open source code can be found at C www.boundary-element-method.com/fortran/LEBEMA_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 LEBEMA. The program computes C the solution to the Laplace equation C C __ 2 C \/ {\phi} = 0 C C exterior to a sphere by the boundary panel method where {\phi} is C known as the potential. C C For the exterior problem, the domain lies exterior to a closed C boundary S. The boundary condition may be Dirichlet, Robin or C Neumann. It is assumed to have the following general form C C {\alpha}(q) {\phi}(q) + {\beta}(q) v(q) = f(q) C C where {\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 LEBEMA accepts a description of the boundary of the domain C and the position of the exterior points where the solution ({\phi}) C is sought, the boundary 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 a sphere of side 1 C (metre). The solution to the problem with a Dirichlet boundary C condition ({\alpha}=1, {\beta}=0) and with a Neumann boundary C condition ({\alpha}=0, beta=1) are sought. C C In the R-z plane, the boundary conditions are specified through C taking the solution to be determined by C C {\phi} = 1/R , C C which gives v=1 on the surface. C C The boundary is described by a set of NS=32 panels of equal size, C so that each side comprises eight panels. The boundary solution C points are the centres of the panels. C The *s show the exterior points at which the solution is sought; C the points (0,2), (1,1), (0,100). C---------------------------------------------------------------------- C The PARAMETER statement C ----------------------- C There are four components in the PARAMETER statement. C integer MAXNS : The limit on the number of boundary 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 LEBEMA: Subroutine for solving the exterior Laplace C equation. (file LEBEMA.FOR contains LEBEMA 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 LEBEMAT 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 MAXNS PARAMETER (MAXNS=32) C Limit on the number of vertices (equal to the number of panels) INTEGER MAXNV PARAMETER (MAXNV=MAXNS) C Limit on the number of points exterior to the boundary, 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 boundary(ies) C Number of panels and counter INTEGER NS,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 SELV(MAXNS,2) C The points exterior to the boundary(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 LEBEMA. C SALPHA(j) is assigned the value of {\alpha} at the centre of the C j-th panel. REAL*8 SALPHA(MAXNS) C SBETA(j) is assigned the value of {\beta} at the centre of the C j-th panel. REAL*8 SBETA(MAXNS) C SF(j) is assigned the value of f at the centre of the j-th panel. REAL*8 SF(MAXNS) C The incident field C The incident potential on the panels of the shell(s) REAL*8 SIPHI(MAXNS) C The derivative of the incident potential on the panels of the shell(s) REAL*8 SIVEL(MAXNS) C The incident potential at the domain points REAL*8 PEIPHI(MAXNS) C Validation and control parameters for LEBEMA 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 boundary. REAL*8 EGEOM C Output from subroutine LEBEMA C The velocity potential (phi - the solution) at the centres of the C panels REAL*8 SPHI(MAXNS) C The normal derivative of the velocity potential at the centres of the C panels REAL*8 SVEL(MAXNS) C The velocity potential (phi - the solution) at exterior points REAL*8 PEPHI(MAXNPE) C Workspace for LEBEMA REAL*8 WKSPC1(MAXNS,MAXNS) REAL*8 WKSPC2(MAXNS,MAXNS) REAL*8 WKSPC3(MAXNPE,MAXNS) REAL*8 WKSPC4(MAXNPE,MAXNS) REAL*8 WKSPC5(MAXNS) REAL*8 WKSPC6(MAXNS) LOGICAL WKSPC7(MAXNS) C Counter through the x,y coordinates INTEGER ICOORD C The coordinates of the centres of the panels REAL*8 SELCNT(MAXNS,2) C Other variables used in specifying the boundary condition REAL*8 Z,R 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 boundary C :The circle that generates the sphere is divided into NS=18 uniform C : panels. VERTEX and SELV are defined anti-clockwise around the C : boundary so that the normal to the boundary is assumed to be C : outward C :Set up nodes C : Set NS, the number of panels NS=18 C : Set NV, the number of vertices (equal to the number of panels) NV=NS+1 C : Set coordinates of the nodes DATA ((VERTEX(IV,ICOORD),ICOORD=1,2),IV=1,19) * / 0.000D0, 1.000D0, * 0.174D0, 0.985D0, * 0.342D0, 0.940D0, * 0.500D0, 0.866D0, * 0.643D0, 0.766D0, * 0.766D0, 0.643D0, * 0.866D0, 0.500D0, * 0.940D0, 0.342D0, * 0.985D0, 0.174D0, * 1.000D0, 0.000D0, * 0.985D0,-0.174D0, * 0.940D0,-0.342D0, * 0.866D0,-0.500D0, * 0.766D0,-0.643D0, * 0.643D0,-0.766D0, * 0.500D0,-0.866D0, * 0.342D0,-0.940D0, * 0.174D0,-0.985D0, * 0.000D0,-1.000D0 / C :Describe the panels that make up the two boundarys C : Set NS, the number of panels NS=18 C : Set nodal indices that describe the panels of the boundarys. C : The indices refer to the nodes in VERTEX. The order of the C : nodes in SELV dictates that the normal is outward from the C : boundary into the domain. DATA ((SELV(IS,ICOORD),ICOORD=1,2),IS=1,18) * / 1, 2, * 2, 3, * 3, 4, * 4, 5, * 5, 6, * 6, 7, * 7, 8, * 8, 9, * 9, 10, * 10, 11, * 11, 12, * 12, 13, * 13, 14, * 14, 15, * 15, 16, * 16, 17, * 17, 18, * 18, 19 / C Set the centres of the panels, the central points DO IS=1,NS SELCNT(IS,1)=(VERTEX(SELV(IS,1),1) * +VERTEX(SELV(IS,2),1))/TWO SELCNT(IS,2)=(VERTEX(SELV(IS,1),2) * +VERTEX(SELV(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, 2.000D0, * 1.000D0, 1.000D0, * 0.000D0, 100.000D0/ C The number of points on the boundary is equal to the number of C panels NSP=NS C :The test problem is an pulsating sphere so that C {\phi}= 1 C :This gives the following v on S C v = -1 DO 250 ISP=1,NSP/2 R=SELCNT(ISP,1) Z=SELCNT(ISP,2) SALPHA(ISP)=1.0D0 SBETA(ISP)=0.0D0 SF(ISP)= 1 250 CONTINUE DO 260 ISP=NSP/2+1,NSP R=SELCNT(ISP,1) Z=SELCNT(ISP,2) SALPHA(ISP)=0.0D0 SBETA(ISP)=1.0D0 SF(ISP)= -1 260 CONTINUE DO 640 ISP=1,NSP SIPHI(ISP)=0.0D0 SIVEL(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 LEBEMA LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 C Set up particular alpha and beta functions for this wavenumber C and type of boundary condition CALL LEBEMA(MAXNV,NV,VERTEX,MAXNS,NS,SELV, * MAXNPE,NPE,PEXT, * SALPHA,SBETA,SF,SIPHI,SIVEL,PEIPHI, * LSOL,LVALID,EGEOM, * SPHI,SVEL,PEPHI, * WKSPC1,WKSPC2,WKSPC3,WKSPC4, * WKSPC5,WKSPC6,WKSPC7) C Output the solutions C Open file for the output data OPEN(UNIT=20,FILE='LEBEMA.OUT') WRITE(20,*) 'LEBEMA: 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),PEPHI(IPE), * 1.0/SQRT(PEXT(IPE,1)**2+PEXT(IPE,2)**2) 800 CONTINUE 999 FORMAT(4F10.4) CLOSE(20) END