C*************************************************************** C * Test program for subroutine H2LC 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/H2LC_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 C This program is a test for the subroutine H2LC. The program computes C the solution to the Helmholtz problem exterior to a circle via the C integral equation of Burton & Miller. In order to use H2LC, the C circle is approximated by a regular polygon with each side being one C element. The boundary functions are approximated by a constant at the C 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 circle, C (3) the number of elements (sides on the approximating polygon), 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 C at the selected exterior points. The program also give the exact C solution at the same points so that computed and exact solutions C may be compared. C The PARAMETER statement C ======================= C There are four components in the PARAMETER statement. C integer LIMK : The limit on the number of wavenumbers. C integer LIMN : The limit on the number of elements. C integer LIMNQ : The limit on the number of quadrature points allowed C in the quadrature rules. C integer LIMNPE : The limit on the number of exterior points. C External modules related to the test program C ============================================ C Subroutine FNHANK: This computes Hankel functions of the first kind C and of order zero and one. C Subroutine CLINSL: This computes the solution to a linear system of C equations. C Subroutine GLRULE: This generates Gauss-Legendre quadrature rules. C real function FNLOG: Returns the natural logarithm. C real function FNSQRT: Returns the square root. C External modules related to the package C ======================================= C Subroutine H2LC: This computes the discrete Helmholtz integral C operators. PROGRAM EBEM2 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 quadrature points in any quadrature rule INTEGER LIMNQ C Limit on the number of exterior points where a solution is sought INTEGER LIMNPE C Set PARAMETERs PARAMETER (LIMK=4,LIMN=64,LIMNQ=16,LIMNPE=4) C Constants REAL*8 ZERO,ONE,TWO,FOUR,PI,TWOPI 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),VECP(2) REAL*8 P1(2),P2(2),RR1(2),RR2(2),NORMQ(2),PEXT(LIMNPE,2) REAL*8 SEGANG,ANGL,ANGU,RAD,ANGP REAL*8 ANG,R1,R2,DRDNQ1,DRDNQ2 LOGICAL LPONEL C Error flag used 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(LIMNQ),WQ(LIMNQ) C Particular quadrature rule C Number of points INTEGER NQT C Points and weights REAL*8 AQT(LIMNQ),WQT(LIMNQ) C Validation and control variables in H2LC LOGICAL LVALID,LFAIL REAL*8 EK,EGEOM,EQRULE LOGICAL LLK,LMK,LMKT,LNK C Storage of discrete integral operators in H2LC 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 variable used in computing the solutions COMPLEX*16 BVEL(LIMN),SUM C Counters and indices INTEGER I,J,IQ,IK,IPE C Variables used in determining the test problem COMPLEX*16 H1(0:1),H2(0:1) COMPLEX*16 KR1,KR2 C Set up constants ZERO=0.0D0 ONE=1.0D0 TWO=2.0D0 FOUR=4.0D0 CIMAG=CMPLX(ZERO,ONE) PI=FOUR*ATAN(ONE) TWOPI=TWO*PI C Set up control and tolerances and open file for storing the possible C error messages LVALID=.TRUE. EK=1D-08 EGEOM=1D-08 EQRULE=1D-08 OPEN(UNIT=10,FILE='H2LC.ERR',STATUS='UNKNOWN') C Set up standard Gaussian Quadrature rule NQ=8 CALL GLRULE(NQ,ZERO,ONE,WQ,AQ) C Set radius of circle (RAD) RAD=1.0D0 C Set the number of elements (N) N=64 C Set the number and values of wavenumber required (NOK and KVAL) NOK=2 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),(1.0D0,0.0D0), * (0.0D0,0.0D0),(0.0D0,0.0D0)/ C Set interior points that determine the Neumann boundary condition and C exact solution. P1 is a unit source point and 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, where an approximation the the solution phi is C 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.LIMNQ) THEN WRITE(*,*) 'ERROR - NQ must be less than or equal to LIMNQ' 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 (LERROR) STOP C Output description of test the problem WRITE(*,*) 'TEST FOR SUBROUTINE H2LC' WRITE(*,*) '=========================' WRITE(*,*) 'Test problem is that of a circle centred at the' WRITE(*,*) 'origin with a Neumann boundary condition' WRITE(*,*) 'The radius of the circle 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 each point makes with the C y-axis at the origin. SEGANG=TWOPI/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 H2LC as being C outward from the boundary of the circle. 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 LPONEL IF (I.EQ.J) THEN LPONEL=.TRUE. ELSE LPONEL=.FALSE. END IF C Select the quadrature rule and store in NQT,AQT,WQT. If LPONEL is C false then use the existing quadrature rule NQ,AQ,WQ. If LPONEL C is true then split the NQ,AQ,WQ quadrature rule into two at the C centre. This deals with the dicontinuity at the centre. IF (LPONEL) THEN NQT=2*NQ DO 120 IQ=1,NQ WQT(IQ)=WQ(IQ)/TWO WQT(NQ+IQ)=WQ(IQ)/TWO AQT(IQ)=AQ(IQ)/TWO AQT(NQ+IQ)=(ONE+AQ(IQ))/TWO 120 CONTINUE ELSE NQT=NQ DO 130 IQ=1,NQ WQT(IQ)=WQ(IQ) AQT(IQ)=AQ(IQ) 130 CONTINUE ENDIF 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 H2LC routine to compute [Lk], [Mk], [Mkt], [Nk] CALL H2LC(K,P,NORMP,QA,QB,LPONEL, * LIMNQ,NQT,AQT,WQT, * LVALID,EK,EGEOM,EQRULE,LFAIL, * LLK,LMK,LMKT,LNK,DISLK,DISMK,DISMKT,DISNK) 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 solution C on the boundary (PHI) ANG=-SEGANG/TWO DO 140 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)) DRDNQ1=-(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)) DRDNQ2=-(RR2(1)*NORMQ(1)+RR2(2)*NORMQ(2))/R2 IF (DBLE(K).GT.EK) THEN KR1=K*R1 CALL FNHANK(KR1,H1) KR2=K*R2 CALL FNHANK(KR2,H2) PHI(I)=CIMAG*(H1(0)-H2(0))/FOUR VEL(I)=-CIMAG*K*(DRDNQ1*H1(1)-DRDNQ2*H2(1))/FOUR ELSE PHI(I)=-(LOG(R1)-LOG(R2))/TWOPI VEL(I)=-(DRDNQ1/R1-DRDNQ2/R2)/TWOPI END IF 140 CONTINUE C Construct composite matrices AK and BK of the Burton and Miller C integral equation DO 150 I=1,N DO 160 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) 160 CONTINUE AK(I,I)=AK(I,I)-ALPHA(IK)/TWO BK(I,I)=BK(I,I)+BETA(IK)/TWO 150 CONTINUE C Solve linear system AK APPPHI = BK VEL to give APPPHI, the C approximation to PHI. DO 170 I=1,N BVEL(I)=ZERO DO 180 J=1,N BVEL(I)=BVEL(I)+BK(I,J)*VEL(J) 180 CONTINUE 170 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 190 I=1,N WRITE(*,999) DBLE(PHI(I)),DIMAG(PHI(I)), * DREAL(APPPHI(I)),DIMAG(APPPHI(I)) 190 CONTINUE C Loop(IPE) through the points in the exterior DO 400 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)-PEXT(IPE,1) RR1(2)=P1(2)-PEXT(IPE,2) RR2(1)=P2(1)-PEXT(IPE,1) RR2(2)=P2(2)-PEXT(IPE,2) R1=SQRT(RR1(1)*RR1(1)+RR1(2)*RR1(2)) R2=SQRT(RR2(1)*RR2(1)+RR2(2)*RR2(2)) IF (DBLE(K).GT.EK) THEN KR1=K*R1 KR2=K*R2 CALL FNHANK(KR1,H1) CALL FNHANK(KR2,H2) PHIPT=CIMAG*(H1(0)-H2(0))/FOUR ELSE PHIPT=-(LOG(R1)-LOG(R2))/TWOPI END IF 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 200 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 NQT=NQ DO 210 IQ=1,NQ WQT(IQ)=WQ(IQ) AQT(IQ)=AQ(IQ) 210 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 H2LC routine to compute [Lk], [Mk] CALL H2LC(K,P,VECP,QA,QB,LPONEL, * LIMNQ,NQT,AQT,WQT, * LVALID,EK,EGEOM,EQRULE,LFAIL, * LLK,LMK,LMKT,LNK,DISLK,DISMK,DISMKT,DISNK) C Accumulate SUM SUM=SUM+DISMK*APPPHI(J)-DISLK*VEL(J) C End loop(J) through elements 200 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 exterior points 400 CONTINUE C End loop(IK) through wavenumbers 250 CONTINUE CLOSE(10) END C Subroutine FNHANK: This computes Hankel functions of the first kind C and of order zero and one. C Input parameter C complex Z : The argument of the Hankel function (Im(Z)=0). C Output parameter C H(0:1) : The Hankel functions of order zero and one respectively. C External routines C S17AEF: From the NAG library. C S17AFF: From the NAG library. C S17ACF: From the NAG library. C S17ADF: From the NAG library. SUBROUTINE FNHANK(Z,H) COMPLEX*16 Z,J0,J1,Y0,Y1,H(0:1) REAL*8 X,REH0,REH1,IMH0,IMH1 REAL*8 S17AEF,S17AFF,S17ACF,S17ADF INTEGER IFAIL IF (ABS(DIMAG(Z)).GT.1.0D-6) THEN WRITE(*,*) 'ERROR(FNHANK) - Im(Z) must be zero' STOP END IF IFAIL=0 X=DREAL(Z) J0=S17AEF(X,IFAIL) J1=S17AFF(X,IFAIL) Y0=S17ACF(X,IFAIL) Y1=S17ADF(X,IFAIL) REH0=DREAL(J0)-DIMAG(Y0) IMH0=DIMAG(J0)+DREAL(Y0) REH1=DREAL(J1)-DIMAG(Y1) IMH1=DIMAG(J1)+DREAL(Y1) H(0)=CMPLX(REH0,IMH0) H(1)=CMPLX(REH1,IMH1) 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 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 modules C D01BBF : From the NAG library. C D01BAZ : From the NAG library. SUBROUTINE GLRULE(N,A,B,WTS,PTS) INTEGER N REAL*8 A,B,WTS(N),PTS(N) INTEGER IFAIL,ITYPE EXTERNAL D01BAZ ITYPE=1 IFAIL=0 CALL D01BBF(D01BAZ,A,B,ITYPE,N,WTS,PTS,IFAIL) END C real function FNLOG: Returns the natural logarithm of X. REAL*8 FUNCTION FNLOG(X) REAL*8 X FNLOG=LOG(X) END C real function FNSQRT: Returns the square root of X. REAL*8 FUNCTION FNSQRT(X) REAL*8 X FNSQRT=SQRT(X) END