C*************************************************************** C Test program for subroutine L2LC 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/L2LC_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 L2LC. The program computes C the solution to the Laplace problem interior In order to use L2LC, C the circle is approximated by a regular polygon with each side being C one element. The boundary functions are approximated by a constant C 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 circle, C (3) the number of elements (sides on the approximating polygon), C (4) the position of one source point and one sink point which C determine the boundary condition and the exact solution, C (5) the points in the interior 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 interior 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 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 LIMNPI : The limit on the number of interior points. C External modules related to the test program C ============================================ 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 L2LC: This computes the discrete Laplace integral C operators. PROGRAM IBEM2 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 quadrature points in any quadrature rule INTEGER LIMNQ C Limit on the number of interior points where a solution is sought INTEGER LIMNPI C Set PARAMETERs PARAMETER (LIMN=64,LIMNQ=16,LIMNPI=4) C Constants REAL*8 ZERO,HALF,ONE,TWO,FOUR,PI,TWOPI REAL*8 CIMAG C Number of elements INTEGER N C Geometrical information INTEGER NPEXT REAL*8 P(2),NORMP(2),Q(2),QA(2),QB(2),VECP(2) REAL*8 NORMQ(2),PINT(LIMNPI,2) REAL*8 SEGANG,ANGL,ANGU,RAD,ANGP REAL*8 ANG 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 L2LC LOGICAL LVALID,LFAIL REAL*8 EGEOM,EQRULE LOGICAL LL,LM,LMT,LN C Storage of discrete integral operators in L2LC REAL*8 DISL,DISM,DISMT,DISN C Matrices corresponding to the Laplace integral operators REAL*8 L(LIMN,LIMN),M(LIMN,LIMN) C Composite matrices 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 interior point REAL*8 PHIPT C Temporary variable used in computing the solutions REAL*8 APHI(LIMN),SUM C Counters and indices INTEGER I,J,IQ,IPI C Set up constants ZERO=0.0D0 HALF=0.5D0 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. EGEOM=1D-08 EQRULE=1D-08 OPEN(UNIT=10,FILE='L2LC.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 interior points, where an approximation the the solution phi is C sought NPEXT=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.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 (NPEXT.GT.LIMNPI) THEN WRITE(*,*) 'ERROR - NPEXT must be less than or equal to LIMNPI' LERROR=.TRUE. END IF DO 50 IPI=1,NPEXT IF (PINT(IPI,1)*PINT(IPI,1)+PINT(IPI,2)*PINT(IPI,2).GT.RAD*RAD) * THEN WRITE(*,*) 'ERROR - PINT must be interior to the circle' LERROR=.TRUE. END IF 50 CONTINUE IF (LERROR) STOP C Output description of test the problem WRITE(*,*) 'TEST FOR SUBROUTINE L2LC' WRITE(*,*) '=========================' WRITE(*,*) 'Test problem is that of a circle centred at the' WRITE(*,*) 'origin with a Dirichlet boundary condition' WRITE(*,*) 'The radius of the circle 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 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 L2LC 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 Only integral operators L and M are required LL=.TRUE. LM=.TRUE. LMT=.FALSE. LN=.FALSE. C Call of L2LC routine to compute [Lk], [Mk] CALL L2LC(P,NORMP,QA,QB,LPONEL, * LIMNQ,NQT,AQT,WQT, * LVALID,EGEOM,EQRULE,LFAIL, * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN) C Storage of results in matrices L(I,J)=DISL M(I,J)=DISM C End loop(J) through elements 110 CONTINUE C End loop(I) through collocation points 100 CONTINUE C Set up test Dirichlet boundary condition (PHI) and the exact solution C on the boundary (VEL). Here \phi = y which has the exact solution C v=y 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) PHI(I)=Q(2) VEL(I)=Q(2) 140 CONTINUE C Construct composite matrices A and B for the linear system DO 150 I=1,N DO 160 J=1,N A(I,J)=M(I,J) B(I,J)=L(I,J) 160 CONTINUE A(I,I)=A(I,I)+HALF 150 CONTINUE C Solve linear system B APPVEL = A PHI to give APPVEL, the C approximation to PHI. DO 170 I=1,N APHI(I)=ZERO DO 180 J=1,N APHI(I)=APHI(I)+A(I,J)*PHI(J) 180 CONTINUE 170 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 PHI ON S APPROX PHI ON S' DO 190 I=1,N WRITE(*,999) VEL(I), APPVEL(I) 190 CONTINUE C Loop(IPI) through the points in the interior DO 400 IPI=1,NPEXT P(1)=PINT(IPI,1) P(2)=PINT(IPI,2) C Output the coordinates of the point in the interior WRITE(*,*) 'Point in the interior is ',P C Compute the exact solution PHIPT at the interior point C Exact solution is \phi = y PHIPT=P(2) C Initialise SUM, this ultimately stores the value of the summation C that approximates to phi at the interior 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. LL=.TRUE. LM=.TRUE. LMT=.FALSE. LN=.FALSE. C Since P is in the interior then LPONEL is always false LPONEL=.FALSE. C Value of VECP is not relevant C Call of L2LC routine to compute [Lk], [Mk] CALL L2LC(P,VECP,QA,QB,LPONEL, * LIMNQ,NQT,AQT,WQT, * LVALID,EGEOM,EQRULE,LFAIL, * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN) C Accumulate SUM SUM=SUM-DISM*PHI(J)+DISL*APPVEL(J) C End loop(J) through elements 200 CONTINUE C Ouput the solution at the interior point WRITE(*,*) 'EXACT PHI AT POINT = ',PHIPT WRITE(*,*) 'APPROX PHI AT POINT = ',SUM C End loop(IPI) through interior points 400 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 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 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