C*************************************************************** C * Test program for subroutine H3ALC by Stephen Kirkup C*************************************************************** C C Copyright 1998- Stephen Kirkup C School of Computing Engineering and Physical Sciences 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/H3ALC_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 H3ALC. The program computes C the solution of the Helmholtz problem exterior to a sphere centred C at the origin and with an axisymmetric solution via the integral C equation of Burton & Miller. The boundary condition and solution are C assumed to be independent of theta (ie axisymmetric). C In order to use H3ALC, the sphere is approximated by a set of conical C elements. Each element is decribed by a straight line on the generator C (R-z plane) swept through 2 pi (in the theta direction). The boundary C functions are approximated by a constant at the centre of each element. C Simple changes in the program allow the set up of the following data: C (1) the choice of quadrature rule, C (2) the radius of the sphere, C (3) the number of elements, C (4) the set of wavenumbers, C (5) the parameter(s) in the Burton & Miller integral equation, C (6) the position of one source point and one sink point which C determine the boundary condition and the exact solution, C (7) the points in the exterior where the solution is sought. C During execution the program gives the solution at the collocation C points (the points at the centre of each element) and the solution at C the selected exterior points. The program also give the exact C solution at the same points so that computed and exact solutions may C be compared. C The PARAMETER statement C ======================= C There are five components in the PARAMETER statements. C integer LIMK : The limit on the number of wavenumbers. C integer LIMN : The limit on the number of elements. C integer LIMNGQ: The limit on the number of quadrature points allowed C in the quadrature rules in the generator direction. C integer LIMNTQ: The limit on the number of quadrature points allowed C in the quadrature rules in the theta direction. C integer LIMNPE: The limit on the number of exterior points. C External modules related to the test program C ============================================ C Subroutine CLINSL: This computes the solution to a linear system of C equations C Subroutine GLRULE: This generates the Gauss-Legendre quadrature rules. C complex function FNEXP: Returns the exponential. C real function FNSQRT: Returns the square root. C External modules related to the package C ======================================= C Subroutine H3ALC: This computes the discrete Helmholtz integral C operators. PROGRAM EBEM3A C Declaration of variables C Limits on the size of the arrays C Limit on the number of wavenumbers INTEGER LIMK C Limit on the number of elements INTEGER LIMN C Limit on the number of generator-direction quadrature points in any C quadrature rule INTEGER LIMNGQ C Limit on the number of theta-direction quadrature points in any C quadrature rule INTEGER LIMNTQ C Limit on the number of exterior points where a solution is sought INTEGER LIMNPE C Set PARAMETERs PARAMETER (LIMK=4,LIMN=16,LIMNGQ=32,LIMNPE=2) PARAMETER (LIMNTQ=LIMN*LIMNGQ) C Constants REAL*8 ZERO,ONE,TWO,FOUR,EIGHT,PI,ROOT2 COMPLEX*16 CIMAG C Number of elements INTEGER N C Wavenumber information INTEGER NOK COMPLEX*16 K,KVAL(LIMK) C Geometrical information INTEGER NPEXT REAL*8 P(2),NORMP(2),Q(2),QA(2),QB(2),NORMQ(2) REAL*8 P1(2),P2(2),RR1(2),RR2(2),R1,R2,PEXT(LIMNPE,2) REAL*8 SEGANG,ANGL,ANGU,ANGM,RADMID,GLEN,SGLEN,CIRMID,TDIV REAL*8 ANG,RAD,ANGP,DR1DNQ,DR2DNQ,VECP(2) LOGICAL LPONEL C Error flag for checking the set-up information LOGICAL LERROR C Quadrature rules C General quadrature rule C Number of points INTEGER NQ C Points and weights REAL*8 AQ(LIMNGQ),WQ(LIMNGQ) C Generator-direction quadrature rule C Number of points INTEGER NGQ C Points and weights REAL*8 AGQ(LIMNGQ),WGQ(LIMNGQ) C Number subdivisions of composite rule INTEGER NDIV C Theta-direction quadrature rule C Number of points INTEGER NTQ C Points and weights REAL*8 ATQ(LIMNTQ),WTQ(LIMNTQ) C Validation and control variables in H3ALC LOGICAL LVALID,LFAIL REAL*8 EK,EGEOM,EQRULE LOGICAL LLK,LMK,LMKT,LNK C Storage of discrete integral operators in H3ALC COMPLEX*16 DISLK,DISMK,DISMKT,DISNK C Matrices corresponding to the Helmholtz integral operators COMPLEX*16 LK(LIMN,LIMN),MK(LIMN,LIMN) COMPLEX*16 MKT(LIMN,LIMN),NK(LIMN,LIMN) C Parameters in the Burton and Miller integral equation for each C wavenumber COMPLEX*16 ALPHA(LIMK),BETA(LIMK) C Composite matrices used in the Burton and Miller equation COMPLEX*16 AK(LIMN,LIMN),BK(LIMN,LIMN) C Exact Neumann boundary condition, exact solution and computed c solution COMPLEX*16 VEL(LIMN),PHI(LIMN),APPPHI(LIMN) C Solution at the exterior point COMPLEX*16 PHIPT C Temporary variables used in computing the solutions COMPLEX*16 BVEL(LIMN),SUM C Counters and indices INTEGER*4 I,IDIV,IQ,J,IPE,IK C Variables used in determining of the test problem COMPLEX*16 KR1,KR2,IKR1,IKR2,E1,E2 C Work space parameter for H3ALC REAL*8 WKSPCE(2*LIMNTQ+LIMNGQ) C Set up constants ZERO=0.0D0 ONE=1.0D0 TWO=2.0D0 FOUR=4.0D0 EIGHT=8.0D0 CIMAG=CMPLX(ZERO,ONE) PI=FOUR*ATAN(ONE) ROOT2=SQRT(TWO) C Set up control and tolerances and open file for the possible error C messages LVALID=.TRUE. EK=1D-08 EGEOM=1D-08 EQRULE=1D-08 OPEN(UNIT=10,FILE='H3ALC.ERR',STATUS='UNKNOWN') C Set up standard Gaussian Quadrature rule. This quadrature rule is C used directly along the generator for most integrals. In the C case when the point P lies at the centre of the element a composite C rule consisting of two sets of the standard rule is used on the C generator. For the theta direction rule a composite rule is formed C from the standard rule so that the density of quadrature points in C the theta direction is approximately the same as the density of C points in the generator direction. NQ=8 CALL GLRULE(NQ,ZERO,ONE,WQ,AQ) C Set radius of the sphere (RAD) RAD=1.0D0 C Set number of elements (N) N=16 C Set the number and values of wavenumber required (NOK and KVAL) NOK=4 DATA KVAL/(0.0D0,0.0D0),(1.0D0,0.0D0), * (0.0D0,1.0D0),(1.0D0,1.0D0)/ C Set the parameters of the Burton and Miller equation corresponding to C the wavenumbers (ALPHA and BETA) DATA ALPHA/(1.0D0,0.0D0),(1.0D0,0.0D0), * (1.0D0,0.0D0),(1.0D0,0.0D0)/ DATA BETA/(0.0D0,0.0D0),(0.0D0,1.0D0), * (0.0D0,0.0D0),(0.0D0,0.0D0)/ C Set interior points (R,z coordinates) that determine the Neumann C boundary condition and exact solution. P1 is a unit source point and C P2 is a unit sink. P1(1)=0.0D0 P1(2)=0.5D0 P2(1)=0.0D0 P2(2)=-0.5D0 C Set exterior points (R,z coordinates), where an approximation to the C solution phi is sought NPEXT=1 PEXT(1,1)=0.0D0 PEXT(1,2)=2.0D0 C Check that the set-up data is satisfactory LERROR=.FALSE. IF (NQ.GT.LIMNGQ) THEN WRITE(*,*) 'ERROR - NQ must be less than or equal to LIMNGQ' LERROR=.TRUE. END IF IF (RAD.LE.ZERO) THEN WRITE(*,*) 'ERROR - RAD must be greater than 0' LERROR=.TRUE. END IF IF (N.GT.LIMN) THEN WRITE(*,*) 'ERROR - N must be less than or equal to LIMN' LERROR=.TRUE. END IF IF (NOK.GT.LIMK) THEN WRITE(*,*) 'ERROR - NOK must be less than or equal to LIMK' LERROR=.TRUE. END IF IF (P1(1)*P1(1)+P1(2)*P1(2).GE.RAD*RAD) THEN WRITE(*,*) 'ERROR - P1 must be interior to the circle' LERROR=.TRUE. END IF IF (P2(1)*P2(1)+P2(2)*P2(2).GE.RAD*RAD) THEN WRITE(*,*) 'ERROR - P2 must be interior to the circle' LERROR=.TRUE. END IF IF (NPEXT.GT.LIMNPE) THEN WRITE(*,*) 'ERROR - NPEXT must be less than or equal to LIMNPE' LERROR=.TRUE. END IF DO 50 IPE=1,NPEXT IF (PEXT(IPE,1)*PEXT(IPE,1)+PEXT(IPE,2)*PEXT(IPE,2).LT.RAD*RAD) * THEN WRITE(*,*) 'ERROR - PEXT must be exterior to the circle' LERROR=.TRUE. END IF 50 CONTINUE IF (NQ*N.GT.LIMNTQ) THEN WRITE(*,*) 'ERROR - LIMNTQ should be at least NQ*N' LERROR=.TRUE. END IF IF (LERROR) STOP C Output description of test the problem WRITE(*,*) 'TEST FOR SUBROUTINE H3ALC' WRITE(*,*) '=========================' WRITE(*,*) 'Test problem is that of a sphere centred at the' WRITE(*,*) 'origin with a Neumann boundary condition' WRITE(*,*) 'The radius of the sphere is ',RAD WRITE(*,*) 'Solution is via the integral equation of Burton' WRITE(*,*) 'and Miller.' WRITE(*,*) 'Number of boundary elements is',N C Loop(IK) through each wavenumber DO 250 IK=1,NOK C Set wavenumber K K=KVAL(IK) C Output the value of the wavenumber WRITE(*,*) 'Wavenumber = ',K C The point P is defined by setting the angle to -SEGANG/2 then adding C SEGANG systematically to give the angle that each point makes with C the z-axis at the origin. SEGANG=PI/DBLE(N) ANGP=-SEGANG/TWO C Loop(I) through each collocation point DO 100 I=1,N C Set P(1..2) (the collocation point) and NORMP(1..2) the normal at C P(1..2). ANGP=ANGP+SEGANG P(1)=RAD*COS(SEGANG/TWO)*SIN(ANGP) P(2)=RAD*COS(SEGANG/TWO)*COS(ANGP) NORMP(1)=SIN(ANGP) NORMP(2)=COS(ANGP) ANGU=0.0D0 C Loop(J) through the elements DO 110 J=1,N C Set QA(1..2) and QB(1..2) the vertices of the element. Since the C integral that is applied assumes the normal to be outward, C the vertices rotate clockwise around the boundary. This is C so that the normal is implicitly defined in H3ALC as being C outward from the boundary of the sphere. ANGL=ANGU ANGU=ANGL+SEGANG ANGM=(ANGU+ANGL)/TWO QA(1)=RAD*SIN(ANGL) QA(2)=RAD*COS(ANGL) QB(1)=RAD*SIN(ANGU) QB(2)=RAD*COS(ANGU) C Set LPONEL IF (I.EQ.J) THEN LPONEL=.TRUE. ELSE LPONEL=.FALSE. END IF C Select the generator-direction quadrature rule and store in NQT, C AQT,WQT. If LPONEL is false then use the existing quadrature rule C NQ,AQ,WQ. If LPONEL is true then split the NQ,AQ,WQ quadrature C rule into two at the centre. This deals with the dicontinuity at C the centre. IF (LPONEL) THEN NGQ=2*NQ DO 120 IQ=1,NQ WGQ(NQ+1-IQ)=WQ(IQ)/TWO WGQ(2*NQ+1-IQ)=WQ(IQ)/TWO AGQ(NQ+1-IQ)=AQ(IQ)/TWO AGQ(2*NQ+1-IQ)=(ONE+AQ(IQ))/TWO 120 CONTINUE ELSE NGQ=NQ DO 130 IQ=1,NQ WGQ(NQ+1-IQ)=WQ(IQ) AGQ(NQ+1-IQ)=AQ(IQ) 130 CONTINUE ENDIF C Quadrature rule in the theta direction is constructed out of individual C Gauss rules so that the length of each is approximately equal to the C length of the element at the generator. RADMID=(QA(1)+QB(1))/TWO SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+ * (QA(2)-QB(2))*(QA(2)-QB(2)) GLEN=SQRT(SGLEN) CIRMID=PI*RADMID NDIV=1+CIRMID/GLEN TDIV=ONE/DBLE(NDIV) NTQ=NDIV*NQ DO 140 IDIV=1,NDIV DO 150 IQ=1,NQ WTQ((IDIV-1)*NQ+IQ)=WQ(IQ)/DBLE(NDIV) ATQ((IDIV-1)*NQ+IQ)=AQ(IQ)/DBLE(NDIV)+TDIV*DBLE(IDIV-1) 150 CONTINUE 140 CONTINUE C For primary stage of method all integral operators are required C in general (for the Burton and Miller equation) LLK=.TRUE. LMK=.TRUE. LMKT=.TRUE. LNK=.TRUE. C Call of H3ALC routine to compute [Lk], [Mk], [Mkt], [Nk] CALL H3ALC(K,P,NORMP,QA,QB,LPONEL, * LIMNGQ,NGQ,AGQ,WGQ,LIMNTQ,NTQ,ATQ,WTQ, * LVALID,EK,EGEOM,EQRULE,LFAIL, * LLK,LMK,LMKT,LNK,DISLK,DISMK,DISMKT,DISNK, * WKSPCE) C Storage of results in matrices LK(I,J)=DISLK MK(I,J)=DISMK MKT(I,J)=DISMKT NK(I,J)=DISNK C End loop(J) through elements 110 CONTINUE C End loop(I) through collocation points 100 CONTINUE C Output the position of the source point and sink point that C prescribe the boundary condition and exact solution WRITE(*,*) 'Neumann boundary condition is that produced by' WRITE(*,*) 'a point source at',P1 WRITE(*,*) 'and a point sink at',P2 C Set up test Neumann boundary condition (VEL) and the exact C solution on the boundary (PHI) ANG=-SEGANG/TWO DO 160 I=1,N ANG=ANG+SEGANG Q(1)=RAD*SIN(ANG) Q(2)=RAD*COS(ANG) NORMQ(1)=SIN(ANG) NORMQ(2)=COS(ANG) RR1(1)=P1(1)-Q(1) RR1(2)=P1(2)-Q(2) R1=SQRT(RR1(1)*RR1(1)+RR1(2)*RR1(2)) DR1DNQ=-(RR1(1)*NORMQ(1)+RR1(2)*NORMQ(2))/R1 RR2(1)=P2(1)-Q(1) RR2(2)=P2(2)-Q(2) R2=SQRT(RR2(1)*RR2(1)+RR2(2)*RR2(2)) DR2DNQ=-(RR2(1)*NORMQ(1)+RR2(2)*NORMQ(2))/R2 KR1=K*R1 IKR1=CIMAG*KR1 KR2=K*R2 IKR2=CIMAG*KR2 E1=EXP(IKR1) E2=EXP(IKR2) PHI(I)=E1/R1-E2/R2 VEL(I)=E1*(IKR1-ONE)*DR1DNQ/R1/R1-E2*(IKR2-ONE)*DR2DNQ/R2/R2 160 CONTINUE C Construct composite matrices AK and BK of the Burton and Miller C integral equation DO 170 I=1,N DO 180 J=1,N AK(I,J)=ALPHA(IK)*MK(I,J)+BETA(IK)*NK(I,J) BK(I,J)=ALPHA(IK)*LK(I,J)+BETA(IK)*MKT(I,J) 180 CONTINUE AK(I,I)=AK(I,I)-ALPHA(IK)/TWO BK(I,I)=BK(I,I)+BETA(IK)/TWO 170 CONTINUE C Solve linear system AK APPPHI = BK VEL to give APPPHI, the C approximation to PHI. DO 190 I=1,N BVEL(I)=ZERO DO 200 J=1,N BVEL(I)=BVEL(I)+BK(I,J)*VEL(J) 200 CONTINUE 190 CONTINUE C Call the subroutine for solving the linear system of equations CALL CLINSL(LIMN,N,AK,APPPHI,BVEL) C Output of solution on boundary 999 FORMAT(2F12.6,' ',2F12.6) WRITE(*,*) ' EXACT PHI ON S APPROX PHI ON S' DO 210 I=1,N WRITE(*,999) DBLE(PHI(I)),DIMAG(PHI(I)), * DREAL(APPPHI(I)),DIMAG(APPPHI(I)) 210 CONTINUE C Loop(IPE) through the points in the exterior DO 300 IPE=1,NPEXT P(1)=PEXT(IPE,1) P(2)=PEXT(IPE,2) C Output the coordinates of the point in the exterior WRITE(*,*) 'Point in the exterior is ',P C Compute the exact solution PHIPT at the exterior point RR1(1)=P1(1)-P(1) RR1(2)=P1(2)-P(2) R1=SQRT(RR1(1)*RR1(1)+RR1(2)*RR1(2)) RR2(1)=P2(1)-P(1) RR2(2)=P2(2)-P(2) R2=SQRT(RR2(1)*RR2(1)+RR2(2)*RR2(2)) PHIPT=EXP(CIMAG*K*R1)/R1-EXP(CIMAG*K*R2)/R2 C Initialise SUM, this ultimately stores the value of the summation C that approximates to phi at the exterior point SUM=ZERO ANGU=0.0D0 C Loop(J) through the elements DO 320 J=1,N C Set QA(1..2) and QB(1..2) the vertices of the element. ANGL=ANGU ANGU=ANGL+SEGANG QA(1)=RAD*SIN(ANGL) QA(2)=RAD*COS(ANGL) QB(1)=RAD*SIN(ANGU) QB(2)=RAD*COS(ANGU) C Set up the quadrature rule C Generator direction NGQ=NQ DO 330 IQ=1,NQ WGQ(IQ)=WQ(NQ+1-IQ) AGQ(IQ)=AQ(NQ+1-IQ) 330 CONTINUE C Quadrature rule in the theta direction is constructed out of C individual Gauss rules so that the length of each is C approximately equal to the length of the element at the C generator. RADMID=(QA(1)+QB(1))/TWO SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+ * (QA(2)-QB(2))*(QA(2)-QB(2)) GLEN=SQRT(SGLEN) CIRMID=PI*RADMID NDIV=1+CIRMID/GLEN TDIV=ONE/DBLE(NDIV) NTQ=NDIV*NQ DO 340 IDIV=1,NDIV DO 350 IQ=1,NQ WTQ((IDIV-1)*NQ+IQ)=WQ(IQ)/DBLE(NDIV) ATQ((IDIV-1)*NQ+IQ)=AQ(IQ)/DBLE(NDIV)+TDIV*DBLE(IDIV-1) 350 CONTINUE 340 CONTINUE C For secondary stage of method Lk and Mk are required. LLK=.TRUE. LMK=.TRUE. LMKT=.FALSE. LNK=.FALSE. C Since P is in the exterior then LPONEL is always false LPONEL=.FALSE. C Value of VECP is not relevant C Call of H3ALC routine to compute [Lk], [Mk] CALL H3ALC(K,P,VECP,QA,QB,LPONEL, * LIMNGQ,NGQ,AGQ,WGQ,LIMNTQ,NTQ,ATQ,WTQ, * LVALID,EK,EGEOM,EQRULE,LFAIL, * LLK,LMK,LMKT,LNK,DISLK,DISMK,DISMKT,DISNK, * WKSPCE) C Accumulate SUM SUM=SUM+DISMK*APPPHI(J)-DISLK*VEL(J) C End loop(J) through elements 320 CONTINUE C Ouput the solution at the exterior point WRITE(*,*) 'EXACT PHI AT POINT = ',PHIPT WRITE(*,*) 'APPROX PHI AT POINT = ',SUM C End loop(IPE) through the exterior points 300 CONTINUE C End loop(IK) through wavenumbers 250 CONTINUE CLOSE(10) END C Subroutine CLINSL C ================= C This subroutine calculates the solution to the linear system C of complex equations Ax=b. C Input parameters C integer MAXN : The maximum dimension of the matrix. Its value must not C be changed between calls of CLINSL. C integer N : The dimension of the matrix. C complex*16 A(MAXN,MAXN) : The matrix `A'. C complex*16 B(MAXN) : The vector `b'. C Output parameters C complex*16 X(MAXN) : The vector `x'. SUBROUTINE CLINSL(MAXN,N,A,X,B) INTEGER MAXN,N COMPLEX*16 A(MAXN,MAXN),X(MAXN),B(MAXN) COMPLEX*16 TEMP,CSUM,RATIO REAL*8 COLMAX,AA INTEGER I,J,II,IIMAX IF (MAXN.LT.N) THEN WRITE(6,*) 'ERROR(CLINSL) - Must have MAXN>=N' STOP ENDIF DO 10 I=1,N-1 COLMAX=0 DO 20 II=I,N AA=ABS(A(II,I)) IF (AA.GT.COLMAX) THEN COLMAX=AA IIMAX=II ENDIF 20 CONTINUE DO 30 J=I,N TEMP=A(I,J) A(I,J)=A(IIMAX,J) A(IIMAX,J)=TEMP 30 CONTINUE TEMP=B(I) B(I)=B(IIMAX) B(IIMAX)=TEMP DO 40 II=I+1,N RATIO=A(II,I)/A(I,I) DO 50 J=I+1,N A(II,J)=A(II,J)-RATIO*A(I,J) 50 CONTINUE B(II)=B(II)-RATIO*B(I) 40 CONTINUE 10 CONTINUE X(N)=B(N)/A(N,N) DO 60 I=N-1,1,-1 CSUM=B(I) DO 70 J=I+1,N CSUM=CSUM-A(I,J)*X(J) 70 CONTINUE X(I)=CSUM/A(I,I) 60 CONTINUE END C Subroutine GLRULE: This generates the Gauss-Legendre quadrature rules. C This subroutine supplies the weights and abscissae for an N point C Gaussian quadrature rule on the interval [A,B]. C Input parameters C integer N : The number of Gaussian points. N must satisfy the C conditions required of its corresponding value in the NAG routine C D01BBF. C real*8 A : The lower limit of the quadrature rule. C real*8 B : The upper limit of the quadrature rule. C Output parameters C real*8 WTS(N) : The Gaussian quadrature weights. C real*8 PTS(N) : The Gaussian quadraure points. C External routines C D01BBF : From the NAG library. C D01BAZ : From the NAG library. SUBROUTINE GLRULE(N,A,B,WTS,PTS) INTEGER*4 N REAL*8 A,B,WTS(N),PTS(N) INTEGER*4 IFAIL,ITYPE EXTERNAL D01BAZ ITYPE=1 IFAIL=0 CALL D01BBF(D01BAZ,A,B,ITYPE,N,WTS,PTS,IFAIL) END C complex function FNEXP: Returns the exponential. COMPLEX*16 FUNCTION FNEXP(Z) COMPLEX*16 Z FNEXP=EXP(Z) END C real function FNSQRT: Returns the square root. REAL*8 FUNCTION FNSQRT(X) REAL*8 X FNSQRT=SQRT(X) END