C*************************************************************** C Test program for subroutine HEBEM3 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/HEBEM3_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 C This program is a test for the subroutine HEBEM3. The program computes C the solution to an acoustic/Helmholtz problem exterior to a sphere C by the boundary element method. C C Background C ---------- 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 complex-valued functions defined C on S. C C Subroutine HEBEM3 accepts the wavenumber, a description of the C boundary of the domain and the position of the exterior 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 exterior points. C C The test problems C ----------------- C C In this test the domain is a sphere of radius 1 (metre). The acoustic C medium is air (at 20 celcius and 1 atmosphere, c=344.0 (metres per C second), density {\rho}=1.205 (kilograms per cubic metre) and the C solution to the problem with a Dirichlet boundary condition C ({\alpha}=1, {\beta}=0) and with a Neumann boundary condition C ({\alpha}=0, beta=1) are sought. For both problems the frequency is C 4000Hz (hence specifying k). C Assuming the sphere is centred at the origin the boundary conditions C are specified through taking the solution to be determined by C C {\phi} = exp(jkr)/r , C C which is clearly a fundamental solution of the Helmholtz equation. C C The boundary is described by a set of NS=36 planar triangular elements C of approximately equal size. The boundary solution points are the C centres of the elements. C The solution is sought at the exterior points (0,0), (0,0.5), C (0.25,0.25). 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 elements. 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 HEBEM3: Subroutine for solving the exterior Helmholtz C equation (file EBEM2.FOR contains EBEM2 and subordinate routines) C subroutine H3LC: Returns the individual discrete Helmholtz integral C operators. (file H3LC.FOR contains H3LC and subordinate routines) C subroutine CGLS: Solves a general linear system of equations. C (file CGLS.FOR contains CGSL and subordinate routines) C file GEOM3D.FOR contains the set of relevant geometric subroutines C The program PROGRAM HEBEM3T IMPLICIT NONE C VARIABLE DECLARATION C -------------------- C PARAMETERs for storing the limits on the dimension of arrays C Limit on the number of elements INTEGER MAXNS PARAMETER (MAXNS=40) C Limit on the number of vertices (equal to the number of elements) INTEGER MAXNV PARAMETER (MAXNV=30) C Limit on the number of test problems INTEGER MAXTEST PARAMETER (MAXTEST=5) C Limit on the number of points exterior to the boundary, where C acoustic properties are sought INTEGER MAXNPE PARAMETER (MAXNPE=6) C Constants C Real scalars: 0, 1, 2, pi REAL*8 ZERO,ONE,TWO,THREE,FOUR,PI C Complex scalars: (0,0), (1,0), (0,1) COMPLEX*16 CZERO,CONE,CIMAG C Wavenumber parameter for HEBEM3 COMPLEX*16 K C Geometrical description of the boundary(ies) C Number of elements and counter INTEGER NS,IS C Number of collocation 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,3) C The three nodes that define each element on the boundaries INTEGER SELV(MAXNS,3) C The points exterior to the boundary(ies) where the acoustic 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,IPE C Coordinates of the exterior points REAL*8 PEXT(MAXNPE,3) C Number of test problems and counter INTEGER NTEST,ITEST C Data structures that contain the parameters that define the test C problems C The wavenumber for each test. KVAL(i) is assigned the wavenumber C of the i-th test problem. COMPLEX*16 KVAL(MAXTEST) C The nature of the boundary condition is specified by assigning C values to the data structures SALVAL and SBEVAL. C SALVAL(i,j) is assigned the value of {\alpha} at the center of the C j-th element for the i-th test problem. COMPLEX*16 SALVAL(MAXTEST,MAXNS) C SBEVAL(i,j) is assigned the value of {\beta} at the center of the C j-th element for the i-th test problem. COMPLEX*16 SBEVAL(MAXTEST,MAXNS) C The actual boundary condition is specified by assigning values to C the data structure SFVAL. C SFVAL(i,j) is assigned the value of f at the center of the j-th C element for the i-th test problem. COMPLEX*16 SFVAL(MAXTEST,MAXNS) C The incident potential at the centres of the elements in each test C case COMPLEX*16 SPHIIN(MAXTEST,MAXNS) C The derivative of the incident potential at the centres of the C elements in each test COMPLEX*16 SVELIN(MAXTEST,MAXNS) C The incident potential at the exterior points COMPLEX*16 PPHIIN(MAXTEST,MAXNPE) C Data structures that are used to define each test problem in turn C and are input parameters to HEBEM3. C SALPHA(j) is assigned the value of {\alpha} at the centre of the C j-th element. COMPLEX*16 SALPHA(MAXNS) C SBETA(j) is assigned the value of {\beta} at the centre of the C j-th element. COMPLEX*16 SBETA(MAXNS) C SF(j) is assigned the value of f at the centre of the j-th element. COMPLEX*16 SF(MAXNS) C The incident potential at the centres of the elements in each test C case COMPLEX*16 SFFPHI(MAXNS) C The derivative of the incident potential at the centres of the C elements in each test COMPLEX*16 SFFVEL(MAXNS) C The incident potential at the exterior points COMPLEX*16 PFFPHI(MAXNPE) C Validation and control parameters for HEBEM3 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 The parameter COMPLEX*16 MU C Output from subroutine HEBEM3 C The velocity potential (phi - the solution) at the centres of the C elements COMPLEX*16 SPHI(MAXNS) C The normal derivative of the velocity potential at the centres of the C elements COMPLEX*16 SVEL(MAXNS) C The velocity potential (phi - the solution) at exterior points COMPLEX*16 PEPHI(MAXNPE) C Workspace for HEBEM3 COMPLEX*16 WKSPC1(MAXNS,MAXNS) COMPLEX*16 WKSPC2(MAXNS,MAXNS) COMPLEX*16 WKSPC3(MAXNPE,MAXNS) COMPLEX*16 WKSPC4(MAXNPE,MAXNS) COMPLEX*16 WKSPC5(MAXNS) COMPLEX*16 WKSPC6(MAXNS) LOGICAL WKSPC7(MAXNS) REAL*8 SIZE3,DOT3 C Acoustic properties. These data structures are appended after each C execution of HEBEM3 and contain the numerical solution to the test C problems. C Exact solutions at the exterior points COMPLEX*16 EXPHI(MAXTEST,MAXNPE) C Computed solutions at the exterior points COMPLEX*16 PPHIS(MAXTEST,MAXNPE) C Counter through the x,y coordinates INTEGER ICOORD C The coordinates of the centres of the elements REAL*8 SELCNT(MAXNS,3) C Other variables used in specifying the boundary condition REAL*8 PSOURCE(3),QA(3),QB(3),QC(3),QNORM(3),RR(3),DRBDN REAL*8 R REAL*8 EPS 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) CZERO=CMPLX(ZERO,ZERO) CONE=CMPLX(ONE,ZERO) CIMAG=CMPLX(ZERO,ONE) EPS=1.0E-10 C Describe the nodes and elements that make up the boundary C :The unit sphere, centred at the origin is divided C : into NS=36 uniform elements. 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 elements NS=36 C : Set NV, the number of vertices (equal to the number of elements) NV=NS C : Set coordinates of the nodes C Set up VERTEX, the coordinates of the vertices of the elements NV=20 DATA ((VERTEX(IV,ICOORD),ICOORD=1,3),IV=1,20) * / 0.000D0, 0.000D0, 1.000D0, * 0.000D0, 0.745D0, 0.667D0, * 0.645D0, 0.372D0, 0.667D0, * 0.645D0,-0.372D0, 0.667D0, * 0.000D0,-0.745D0, 0.667D0, * -0.645D0,-0.372D0, 0.667D0, * -0.645D0, 0.372D0, 0.667D0, * 0.500D0, 0.866D0, 0.000D0, * 1.000D0, 0.000D0, 0.000D0, * 0.500D0,-0.866D0, 0.000D0, * -0.500D0,-0.866D0, 0.000D0, * -1.000D0, 0.000D0, 0.000D0, * -0.500D0, 0.866D0, 0.000D0, * 0.000D0, 0.745D0,-0.667D0, * 0.645D0, 0.372D0,-0.667D0, * 0.645D0,-0.372D0,-0.667D0, * 0.000D0,-0.745D0,-0.667D0, * -0.645D0,-0.372D0,-0.667D0, * -0.645D0, 0.372D0,-0.667D0, * 0.000D0, 0.000D0,-1.000D0/ C : Set nodal indices that describe the elements 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 acoustic domain. DATA ((SELV(IS,ICOORD),ICOORD=1,3),IS=1,36) * / 1, 3, 2, 1, 4, 3, 1, 5, 4, 1, 6, 5, * 1, 7, 6, 1, 2, 7, 2, 3, 8, 3, 9, 8, * 3, 4, 9, 4,10, 9, 4, 5,10, 5,11,10, * 5, 6,11, 6,12,11, 6, 7,12, 7,13,12, * 7, 2,13, 2, 8,13, 8,15,14, 8, 9,15, * 9,16,15, 9,10,16, 10,17,16, 10,11,17, * 11,18,17, 11,12,18, 12,19,18, 12,13,19, * 13,14,19, 13, 8,14, 14,15,20, 15,16,20, * 16,17,20, 17,18,20, 18,19,20, 19,14,20/ C Set the centres of the elements, the collocation points DO IS=1,NS SELCNT(IS,1)=(VERTEX(SELV(IS,1),1) * +VERTEX(SELV(IS,2),1)+VERTEX(SELV(IS,3),1))/THREE SELCNT(IS,2)=(VERTEX(SELV(IS,1),2) * +VERTEX(SELV(IS,2),2)+VERTEX(SELV(IS,3),2))/THREE SELCNT(IS,3)=(VERTEX(SELV(IS,1),3) * +VERTEX(SELV(IS,2),3)+VERTEX(SELV(IS,3),3))/THREE END DO C Set the points in the acoustic domain where the acoustic properties C are sought, PINT. C : Let them be the points (0.000,0.000), (0.000,0.500), (0.250,0.250). NPE=4 DATA ((PEXT(IPE,ICOORD),ICOORD=1,3),IPE=1,4) * / 0.000D0, 0.000D0, 2.000D0, * 0.000D0, 0.000D0, 4.000D0, * 0.000D0, 0.000D0, 8.000D0, * 0.000D0, 0.000D0, -2.000D0/ C The number of points on the boundary is equal to the number of C elements NSP=NS C Set up test problems C :Set the number of test problems NTEST=2 C TEST PROBLEM 1 C ============== C : Set the wavenumber in KVAL KVAL(1)=CMPLX(ONE,ZERO) C :Set nature of the boundary condition by prescribing the values of C the boundary functions SALVAL and SBEVAL at the collocation points C :In this case a Dirichlet (phi-valued) boundary condition DO 160 ISP=1,NSP SALVAL(1,ISP)=CONE SBEVAL(1,ISP)=CZERO 160 CONTINUE C :The test problem is devised so that C {\phi}=exp(jkr)/r C where r is the distance from PSOURCE the source point to the C boundary point C :Set PSOURCE PSOURCE(1)=0.0D0 PSOURCE(2)=0.0D0 PSOURCE(3)=0.0D0 C :Set K, the wavenumber K=KVAL(1) DO 170 ISP=1,NSP RR(1)=SELCNT(ISP,1)-PSOURCE(1) RR(2)=SELCNT(ISP,2)-PSOURCE(2) RR(3)=SELCNT(ISP,3)-PSOURCE(3) R=SIZE3(RR) SFVAL(1,ISP)=(COS(K*R)+CIMAG*SIN(K*R))/R 170 CONTINUE DO 180 ISP=1,NSP SPHIIN(1,ISP)=CZERO SVELIN(1,ISP)=CZERO 180 CONTINUE DO 190 IPE=1,NPE PPHIIN(1,IPE)=CZERO 190 CONTINUE DO 195 IPE=1,NPE RR(1)=PEXT(IPE,1)-PSOURCE(1) RR(2)=PEXT(IPE,2)-PSOURCE(2) RR(3)=PEXT(IPE,3)-PSOURCE(3) R=SIZE3(RR) EXPHI(1,IPE)=(COS(K*R)+CIMAG*SIN(K*R))/R 195 CONTINUE C TEST PROBLEM 2 C ============== C : Set the wavenumber in KVAL KVAL(2)=CMPLX(TWO,ZERO) C :Set nature of the boundary condition by prescribing the values of C the boundary functions SALVAL and SBEVAL at the collocation points C :In this case a Dirichlet (phi-valued) boundary condition DO 200 ISP=1,NSP SALVAL(2,ISP)=CZERO SBEVAL(2,ISP)=CONE 200 CONTINUE DO 210 ISP=1,NSP SPHIIN(2,ISP)=CZERO SVELIN(2,ISP)=CZERO 210 CONTINUE DO 220 IPE=1,NPE PPHIIN(2,IPE)=CZERO 220 CONTINUE C :The test problem is devised so that C d{\phi}/dn=exp(jkr)(jkr-1)/r/r C where r is the distance from PSOURCE the source point to the C boundary point C :Set PSOURCE PSOURCE(1)=0.0D0 PSOURCE(2)=0.0D0 PSOURCE(3)=0.0D0 C :Set K, the wavenumber K=KVAL(2) DO 300 ISP=1,NSP RR(1)=SELCNT(ISP,1)-PSOURCE(1) RR(2)=SELCNT(ISP,2)-PSOURCE(2) RR(3)=SELCNT(ISP,3)-PSOURCE(3) R=SIZE3(RR) QA(1)=VERTEX(SELV(ISP,1),1) QA(2)=VERTEX(SELV(ISP,1),2) QA(3)=VERTEX(SELV(ISP,1),3) QB(1)=VERTEX(SELV(ISP,2),1) QB(2)=VERTEX(SELV(ISP,2),2) QB(3)=VERTEX(SELV(ISP,2),3) QC(1)=VERTEX(SELV(ISP,3),1) QC(2)=VERTEX(SELV(ISP,3),2) QC(3)=VERTEX(SELV(ISP,3),3) CALL NORM3(QA,QB,QC,QNORM) DRBDN=DOT3(RR,QNORM)/R SFVAL(2,ISP)=(COS(K*R) + * CIMAG*SIN(K*R))*(-ONE+CIMAG*K*R)/R/R*DRBDN 300 CONTINUE DO 295 IPE=1,NPE RR(1)=PEXT(IPE,1)-PSOURCE(1) RR(2)=PEXT(IPE,2)-PSOURCE(2) RR(3)=PEXT(IPE,3)-PSOURCE(3) R=SIZE3(RR) EXPHI(2,IPE)=(COS(K*R)+CIMAG*SIN(K*R))/R 295 CONTINUE C :Switch for the particular solution LSOL=.TRUE. C :Switch on the validation of HEBEM3 LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 C Loop(ITEST) through the test problems DO 500 ITEST=1,NTEST C Set OMEGA, the angular frequency omega and K, the wavenumber K=KVAL(ITEST) C Set up particular alpha and beta functions for this wavenumber C and type of boundary condition DO 640 ISP=1,NSP SALPHA(ISP)=SALVAL(ITEST,ISP) SBETA(ISP)=SBEVAL(ITEST,ISP) SF(ISP)=SFVAL(ITEST,ISP) SFFPHI(ISP)=SPHIIN(ITEST,ISP) SFFVEL(ISP)=SVELIN(ITEST,ISP) 640 CONTINUE DO 650 IPE=1,NPE PFFPHI(IPE)=PPHIIN(ITEST,IPE) 650 CONTINUE C Set the parameter IF(DIMAG(K).GT.DREAL(K)) THEN MU=CZERO ELSE MU=CIMAG/(DREAL(K)-DIMAG(K)+ONE) END IF CALL HEBEM3(K, * MAXNV,NV,VERTEX,MAXNS,NS,SELV, * MAXNPE,NPE,PEXT, * SALPHA,SBETA,SF,SFFPHI,SFFVEL,PFFPHI, * LSOL,LVALID,EGEOM,MU, * SPHI,SVEL,PEPHI, * WKSPC1,WKSPC2,WKSPC3,WKSPC4, * WKSPC5,WKSPC6,WKSPC7) DO 695 IPE=1,NPE PPHIS(ITEST,IPE)=PEPHI(IPE) 695 CONTINUE C Close loop(ITEST) through the test problems 500 CONTINUE C Output the solutions C Open file for the output data OPEN(UNIT=20,FILE='HEBEM3.OUT') C Formats for output 2810 FORMAT(1X,'Wavenumber = ',F8.2/) 2845 FORMAT(4X,'Potential at the exterior points',/) 2850 FORMAT(5X,'index',7X,'Potential',19X,'Exact Potential'/) 2860 FORMAT(4X,I4,2X,E10.4,'+ ',E10.4,' i ', * E10.4, '+ ',E10.4,' i '/) WRITE(20,*) 'HEBEM3: Computed solution to the Helmholtz problems' WRITE(20,*) C Loop(ITEST) through the test problems. DO 2000 ITEST=1,NTEST C Output the wavenumber WRITE(20,*) WRITE(20,*) WRITE(20,*) 'Test problem ',ITEST WRITE(20,*) WRITE(20,2810) KVAL(ITEST) WRITE(20,*) WRITE(20,2845) WRITE(20,*) WRITE(20,2850) DO 2020 IPE=1,NPE WRITE(20,2860) IPE,DREAL(PPHIS(ITEST,IPE)), * DIMAG(PPHIS(ITEST,IPE)), * DREAL(EXPHI(ITEST,IPE)), * DIMAG(EXPHI(ITEST,IPE)) 2020 CONTINUE WRITE(20,*) WRITE(20,*) 2000 CONTINUE CLOSE(20) END