C ************************************************************** C * Test program for subroutine L3ALC 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/L3ALC.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 L3ALC. The program computes C the solution of the Laplace problem exterior to a sphere centred C at the origin and with an axisymmetric solution via the integral C equation (M+I/2)\phi = L v. The boundary condition and solution are C assumed to be independent of theta (ie axisymmetric). C In order to use L3ALC, 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 position of one source point and one sink point which C determine the boundary condition and the exact solution, C (6) 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 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 LIMNPI: The limit on the number of interior points. C External modules related to the test program C ============================================ C Subroutine LINSL: This computes the solution to a linear system of C equations C Subroutine GLRULE: This generates the Gauss-Legendre quadrature rules. C real function FNSQRT: Returns the square root. C External modules related to the package C ======================================= C Subroutine L3ALC: This computes the discrete Laplace integral C operators. PROGRAM IBEM3A C Declaration of variables C Limits on the size of the arrays 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 interior points where a solution is sought INTEGER LIMNPi C Set PARAMETERs PARAMETER (LIMN=16,LIMNGQ=32,LIMNPI=2) PARAMETER (LIMNTQ=LIMN*LIMNGQ) C Constants REAL*8 ZERO,HALF,ONE,TWO,FOUR,EIGHT,PI,ROOT2 REAL*8 CIMAG C Number of elements INTEGER N C Geometrical information INTEGER NPINT REAL*8 P(2),NORMP(2),Q(2),QA(2),QB(2) REAL*8 PINT(LIMNPI,2) REAL*8 SEGANG,ANGL,ANGU,ANGM,RADMID,GLEN,SGLEN,CIRMID,TDIV REAL*8 RAD,ANGP,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 L3ALC LOGICAL LVALID,LFAIL REAL*8 EGEOM,EQRULE LOGICAL LL,LM,LMT,LN C Storage of discrete integral operators in L3ALC REAL*8 DISL,DISM,DISMT,DISN C Matrices corresponding to the Laplace integral operators REAL*8 M(LIMN,LIMN),L(LIMN,LIMN) C Composite matrices used in the Burton and Miller equation REAL*8 A(LIMN,LIMN),B(LIMN,LIMN) C Exact Dirichlet boundary condition, exact solution and computed c solution REAL*8 PHI(LIMN),VEL(LIMN),APPVEL(LIMN) C Solution at the exterior point REAL*8 PHIPT C Temporary variables used in computing the solutions REAL*8 APHI(LIMN),SUM C Counters and indices INTEGER*4 I,IDIV,IQ,J,IPE C Work space parameter for L3ALC REAL*8 WKSPCE(2*LIMNTQ+LIMNGQ) C Set up constants ZERO=0.0D0 HALF=0.5D0 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=.FALSE. EGEOM=1D-08 EQRULE=1D-08 OPEN(UNIT=10,FILE='L3ALC.ERR',STATUS='UNKNOWN') OPEN(UNIT=11,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 exterior points (R,z coordinates), where an approximation to the C solution phi is sought NPINT=1 PINT(1,1)=0.0D0 PINT(1,2)=0.5D0 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 (NPINT.GT.LIMNPI) THEN WRITE(*,*) 'ERROR - NPINT must be less than or equal to LIMNPI' LERROR=.TRUE. END IF DO 50 IPE=1,NPINT IF (PINT(IPE,1)*PINT(IPE,1)+PINT(IPE,2)*PINT(IPE,2).GT.RAD*RAD) * THEN WRITE(*,*) 'ERROR - PINT must be interior 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 L3ALC' WRITE(*,*) '=========================' WRITE(*,*) 'Test problem is that of a sphere centred at the' WRITE(*,*) 'origin with a Dirichlet boundary condition' WRITE(*,*) 'The radius of the sphere is ',RAD WRITE(*,*) 'Number of boundary elements is',N 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 L3ALC 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) LL=.TRUE. LM=.TRUE. LMT=.FALSE. LN=.FALSE. C Call of L3ALC routine to compute [Lk], [Mk] CALL L3ALC(P,NORMP,QA,QB,LPONEL, * LIMNGQ,NGQ,AGQ,WGQ,LIMNTQ,NTQ,ATQ,WTQ, * LVALID,EGEOM,EQRULE,LFAIL, * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN, * WKSPCE) C Storage of results in matrices M(I,J)=DISM L(I,J)=DISL C End loop(J) through elements 110 CONTINUE C End loop(I) through collocation points 100 CONTINUE C Set up the Neumann boundary condition and the exact solution C Here \phi = z which has the solution v = z. ANGU=0.0D0 DO 160 J=1,N ANGL=ANGU ANGU=ANGL+SEGANG ANGM=(ANGU+ANGL)/TWO Q(1)=RAD*SIN(ANGM) Q(2)=RAD*COS(ANGM) PHI(J)=Q(2) VEL(J)=Q(2) 160 CONTINUE C Construct composite matrices A and B for the linear system DO 170 I=1,N DO 180 J=1,N A(I,J)=M(I,J) B(I,J)=L(I,J) 180 CONTINUE A(I,I)=A(I,I)+HALF 170 CONTINUE C Solve linear system A PHI = B APPVEL to give APPVEL, the C approximation to VEL. DO 190 I=1,N APHI(I)=ZERO DO 200 J=1,N APHI(I)=APHI(I)+A(I,J)*PHI(J) 200 CONTINUE 190 CONTINUE C Call the subroutine for solving the linear system of equations CALL LINSL(LIMN,N,B,APPVEL,APHI) C Output of solution on boundary 999 FORMAT(2F12.6,' ',2F12.6) WRITE(*,*) 'EXACT, APPROX VEL ON S' DO 210 I=1,N WRITE(*,999) VEL(I),APPVEL(I) 210 CONTINUE C Loop(IPE) through the points in the interior DO 300 IPE=1,NPINT P(1)=PINT(IPE,1) P(2)=PINT(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 C Here \phi = z PHIPT=P(2) 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. LL=.TRUE. LM=.TRUE. LMT=.FALSE. LN=.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 L3ALC routine to compute [Lk], [Mk] CALL L3ALC(P,VECP,QA,QB,LPONEL, * LIMNGQ,NGQ,AGQ,WGQ,LIMNTQ,NTQ,ATQ,WTQ, * LVALID,EGEOM,EQRULE,LFAIL, * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN, * WKSPCE) C Accumulate SUM SUM=SUM-DISM*PHI(J)+DISL*APPVEL(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 CLOSE(10) END C Subroutine LINSL C ================ C This subroutine implements a Gaussian elimination procedure with C partial pivotting in order to solve the linear system of equations C A X = B. AMAT(MAXN,MAXN) stores the NxN maxtrix on input but is altered C on output. BVEC(MAXN) stores the N vector on input and XVEC(MAXN) the C solution on output. SUBROUTINE LINSL(MAXN,N,AMAT,XVEC,BVEC) IMPLICIT REAL*8 (A-H,O-Z) INTEGER MAXN,N REAL*8 AMAT(MAXN,MAXN),BVEC(MAXN),XVEC(MAXN) REAL*8 TEMP,CSUM,RATIO REAL*8 EPS,AA,COLMAX INTEGER II,IIMAX,J,I EPS=1.0D-20 DO 10 I=1,N-1 COLMAX=0 DO 20 II=I,N AA=ABS(AMAT(II,I)) IF (AA.GT.COLMAX) THEN COLMAX=AA IIMAX=II ENDIF 20 CONTINUE DO 30 J=I,N TEMP=AMAT(I,J) AMAT(I,J)=AMAT(IIMAX,J) AMAT(IIMAX,J)=TEMP 30 CONTINUE DO 40 II=I+1,N IF (ABS(AMAT(I,I)).LT.EPS) THEN WRITE(*,*) 'ERROR(LINSL) - Input matrix A is singular' GOTO 998 END IF RATIO=AMAT(II,I)/AMAT(I,I) DO 50 J=I+1,N AMAT(II,J)=AMAT(II,J)-RATIO*AMAT(I,J) 50 CONTINUE BVEC(II)=BVEC(II)-RATIO*BVEC(I) 40 CONTINUE 10 CONTINUE XVEC(N)=BVEC(N)/AMAT(N,N) DO 60 I=N-1,1,-1 CSUM=BVEC(I) DO 70 J=I+1,N CSUM=CSUM-AMAT(I,J)*XVEC(J) 70 CONTINUE XVEC(I)=CSUM/AMAT(I,I) 60 CONTINUE 998 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 condions C required of its corresponding value in the NAG routine 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 real function FNSQRT: Returns the square root. REAL*8 FUNCTION FNSQRT(X) REAL*8 X FNSQRT=SQRT(X) END