C*************************************************************** C Test program for subroutine LIBEM2 by Stephen Kirkup C*************************************************************** C C Version 2 (Sept 2015) (tidy up of previous version Jul 2015) C (previous version 2001) Stephen Kirkup C School of Engineering 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/LIBEM2_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 LIBEM2. The program computes C the solution to an Laplace problem interior to a square C by the boundary panel method. C C Background C ---------- C C The test problem solves the Laplace equation C C __ 2 C \/ {\phi} = 0 C C where {\phi} is known as the potential. C C For the interior problem, the domain lies interior 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 complex-values functions defined C on S. C C Subroutine LIBEM2 accepts a description of the C boundary of the domain and the position of the interior points C where the solution ({\phi}) is sought, the boundary condition and C returns the solution ({\phi} and v) on S and the value of {\phi} C at the interior points. C C The test problems C ----------------- C C In this test the domain is a square of side 0.1 (metre) with a C Dirichlet boundary condition ({\alpha}=1, {\beta}=0). C Assuming the square has vertices (0,0), (0,0.1), (0.1,0.1) and (0.1,0) C in the x-y plane, the boundary conditions are specified through C taking the solution to be determined by C C 2 2 C {\phi} = x - y (in the x-y plane). C C The form of the test problems are illustrated in the following C diagram. C C ^ C y | C | ^n C | C 0.1 -------------------------------------- C |B C| C | | C | * * | C | 2 2 | C <- | {\phi} = x - y | -> C n | | n C | * | C | | C | | C | | C | * * | C | | C |A D| C 0.0 -------------------------------------- C | n C 0.0 \ / 0.1 ---> x C C Note that on side CD v = d {\phi}/dx, on AB v = - d {\phi}/dx, C on BC v = d {\phi}/dy and on DA v = - d{\phi}/dy: the sign C is negative if the normal is in the negative x or y direction. 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 interior points at which the solution is sought; C the points (0.025,0.025), (0.075,0.025), (0.025,0.075), C (0.075,0.075), and (0.05,0.05). 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 MAXNPI : The limit on the number of interior points. C External modules related to the package C --------------------------------------- C subroutine LIBEM2: Subroutine for solving the interior Laplace C equation. (file LIBEM2.FOR contains LIBEM2 and subordinate routines) C subroutine L2LC: Returns the individual discrete Laplace integral C operators. (file L2LC.FOR contains L2LC and subordinate routines) C subroutine GLS: Solves a general linear system of equations. C (file GLS.FOR contains GLS and subordinate routines) C file GEOM2D.FOR contains the set of relevant geometric subroutines C The program PROGRAM LIBEM2T 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 interior to the boundary, where C properties are sought INTEGER MAXNPI PARAMETER (MAXNPI=6) C Constants C Real scalars: 0, 1, 2, pi REAL*8 ZERO,ONE,TWO,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 interior to the boundary(ies) where the C properties are sought and the directional vectors at those points. C [Only necessary if an interior solution is sought.] C Number of interior points and counter INTEGER NPI,IPI C Coordinates of the interior points REAL*8 PINT(MAXNPI,2) C Data structures that are used to define each test problem in turn C and are input parameters to LIBEM2. 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 SIPHI(j) is the value of the free-field potential at the centre C of the j-th panel. REAL*8 SIPHI(MAXNS) C PDIPHI(j) is the value of the free-field potential at the interior C points. REAL*8 PDIPHI(MAXNPI) C Validation and control parameters for LIBEM2 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 LIBEM2 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 interior points REAL*8 PDPHI(MAXNPI) C Workspace for LIBEM2 REAL*8 L_SS(MAXNS,MAXNS) REAL*8 M_SS(MAXNS,MAXNS) REAL*8 L_DS(MAXNPI,MAXNS) REAL*8 M_DS(MAXNPI,MAXNS) REAL*8 C(MAXNS) REAL*8 WKSPC1(MAXNS) REAL*8 WKSPC2(MAXNS) REAL*8 WKSPC3(MAXNS) INTEGER*4 PERM(MAXNS) LOGICAL XORY(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 The normals to the panels REAL*8 NORMAL(MAXNS,2) C Functions defining the exact {\phi} and d{\phi}/dn REAL*8 PHI,DPHIDN C Temporary variables REAL*8 Q(2),QNORM(2),QA(2),QB(2) C Other variables REAL*8 EPS C INITIALISATION C -------------- C Set constants ZERO=0.0D0 ONE=1.0D0 TWO=2.0D0 PI=4.0D0*ATAN(ONE) EPS=1.0E-10 C Describe the nodes and panels that make up the boundary C :The square with vertices (0,0), (0.1,0), (0.1,0.1), (0,0.1) is divided C : into NS=32 uniform panels. VERTEX and SELV are defined C : anti-clockwise around the boundary so that the normal to the C : boundary is assumed to be outward C :Set up nodes C : Set NS, the number of panels NS=32 C : Set NV, the number of vertices (equal to the number of panels) NV=NS C : Set coordinates of the nodes DATA ((VERTEX(IV,ICOORD),ICOORD=1,2),IV=1,32) * / 0.0000000000D0, 0.00000000000D0, * 0.0000000000D0, 0.1250000000D0, * 0.0000000000D0, 0.2500000000D0, * 0.0000000000D0, 0.3750000000D0, * 0.0000000000D0, 0.5000000000D0, * 0.0000000000D0, 0.6250000000D0, * 0.0000000000D0, 0.7500000000D0, * 0.0000000000D0, 0.8750000000D0, * 0.0000000000D0, 1.0000000000D0, * 0.1250000000D0, 1.0000000000D0, * 0.2500000000D0, 1.0000000000D0, * 0.3750000000D0, 1.0000000000D0, * 0.5000000000D0, 1.0000000000D0, * 0.6250000000D0, 1.0000000000D0, * 0.7500000000D0, 1.0000000000D0, * 0.8750000000D0, 1.0000000000D0, * 1.0000000000D0, 1.0000000000D0, * 1.0000000000D0, 0.8750000000D0, * 1.0000000000D0, 0.7500000000D0, * 1.0000000000D0, 0.6250000000D0, * 1.0000000000D0, 0.5000000000D0, * 1.0000000000D0, 0.3750000000D0, * 1.0000000000D0, 0.2500000000D0, * 1.0000000000D0, 0.1250000000D0, * 1.00000000000D0, 0.0000000000D0, * 0.8750000000D0, 0.0000000000D0, * 0.7500000000D0, 0.0000000000D0, * 0.6250000000D0, 0.0000000000D0, * 0.5000000000D0, 0.0000000000D0, * 0.3750000000D0, 0.0000000000D0, * 0.2500000000D0, 0.0000000000D0, * 0.1250000000D0, 0.0000000000D0 / C :Describe the panels that make up the two boundarys C : Set NS, the number of panels NS=32 C : Set nodal indices that describe the panels of the boundary(s). 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,32) * / 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, * 19, 20, * 20, 21, * 21, 22, * 22, 23, * 23, 24, * 24, 25, * 25, 26, * 26, 27, * 27, 28, * 28, 29, * 29, 30, * 30, 31, * 31, 32, * 32, 1 / C Set the centres of the panels, the central points DO 100 IS=1,NS QA(1)=VERTEX(SELV(IS,1),1) QA(2)=VERTEX(SELV(IS,1),2) QB(1)=VERTEX(SELV(IS,2),1) QB(2)=VERTEX(SELV(IS,2),2) SELCNT(IS,1)=(QA(1)+QB(1))/TWO SELCNT(IS,2)=(QA(2)+QB(2))/TWO CALL NORM2(QA,QB,QNORM) NORMAL(IS,1)=QNORM(1) NORMAL(IS,2)=QNORM(2) 100 CONTINUE C Set the points in the domain where the properties C are sought, PINT. C : Let them be the points (0.025,0.025), (0.075,0.025), (0.025,0.075), C : (0.075,0.075) and (0.5,0.5). NPI=5 DATA ((PINT(IPI,ICOORD),ICOORD=1,2),IPI=1,5) * / 0.250D0, 0.250D0, * 0.750D0, 0.250D0, * 0.250D0, 0.750D0, * 0.750D0, 0.750D0, * 0.500D0, 0.500D0 / C The number of points on the boundary is equal to the number of C panels NSP=NS C Set up test problems C Set up validation and control parameters C :Switch for particular solution LSOL=.TRUE. C :Switch on the validation of LIBEM2 LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 C Set up particular alpha and beta functions for the boundary C condition C :Set nature of the boundary condition by prescribing the values of C the boundary functions SALVAL and SBEVAL at the central points DO 640 ISP=1,NSP/2 C :In this case a Dirichlet (phi-valued) boundary condition SALPHA(ISP)=ONE SBETA(ISP)=ZERO Q(1)=SELCNT(ISP,1) Q(2)=SELCNT(ISP,2) SF(ISP)=PHI(Q) SIPHI(ISP)=ZERO 640 CONTINUE DO 645 ISP=NSP/2+1,NSP C :In this case a Neumann boundary condition SALPHA(ISP)=ZERO SBETA(ISP)=ONE Q(1)=SELCNT(ISP,1) Q(2)=SELCNT(ISP,2) QNORM(1)=NORMAL(ISP,1) QNORM(2)=NORMAL(ISP,2) SF(ISP)=DPHIDN(Q,QNORM) SIPHI(ISP)=ZERO C(ISP)=SIPHI(ISP) 645 CONTINUE DO 650 IPI=1,NPI PDIPHI(IPI)=ZERO 650 CONTINUE CALL LIBEM2(MAXNV,NV,VERTEX,MAXNS,NS,SELV, * MAXNPI,NPI,PINT, * SALPHA,SBETA,SF,SIPHI,PDIPHI, * LSOL,LVALID,EGEOM, * SPHI,SVEL,PDPHI, * L_SS,M_SS,L_DS,M_DS, * C,PERM,XORY,WKSPC1,WKSPC2,WKSPC3) C Output the solutions C Open file for the output data OPEN(UNIT=20,FILE='LIBEM2.OUT') C Formats for output WRITE(20,*) 'LIBEM2: Solution to the test problem LIBEM2_T' WRITE(20,*) C Loop(ISP) through the points on the boundary WRITE(20,*) ' x y phi_app phi_ex' DO 210 IPI=1,NPI Q(1)=PINT(IPI,1) Q(2)=PINT(IPI,2) WRITE(20,999) Q(1),Q(2),PDPHI(IPI),PHI(Q) 210 CONTINUE C Close file CLOSE(20) 999 FORMAT(2F10.3,' ',2F10.5) END C PHI = x**2 - y**2 REAL*8 FUNCTION PHI(Q) REAL*8 Q(2) PHI=2*(Q(1)**2-Q(2)**2) END REAL*8 FUNCTION DPHIDN(Q,QNORM) REAL*8 Q(2),QNORM(2) DPHIDN=2*(2*Q(1)*QNORM(1)-2*Q(2)*QNORM(2)) END