C ************************************************************** C * Test program for subroutine L3LC 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/l3lc_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 subroutine L3LC. The program computes the C solution to the Laplace problem interior to a cube.In order to use C L3LC, the cube is represented by 24 uniform triangles. The boundary C functions are approximated by a constant on 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 geometry of the surface (of the approximating polyhedron), C (3) the set of wavenumbers, C (4) the position of one source point which determines the boundary C 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 at C the selected interior point. 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 five components in the PARAMETER statement. C integer LIMN : The limit on the number of elements. C integer LIMNV : The limit on the number of vertices. 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 LINSL: This computes the solution to a linear system of C equations. C real function FNSQRT: Returns the square root. C External modules related to the package C ======================================= C Subroutine L3LC: This computes the discrete Laplace integral C operators. PROGRAM IBEM3 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 vertices INTEGER LIMNV 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=24,LIMNV=20,LIMNQ=30,LIMNPI=4) C Constants REAL*8 ZERO,ONE,TWO,FOUR,EIGHT,THIRD,HALF,THREE,PI,ROOT2 C Number of elements INTEGER*4 N C Geometrical information REAL*8 VERTPT(LIMNV,3),ELNORM(LIMN,3) INTEGER ELVERT(LIMN,3),NOVERT,NPINT REAL*8 P(3),NORMP(3),Q(3),QA(3),QB(3),QC(3),NORMQ(3) REAL*8 PINT(LIMNPI,3),VECP(3) 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 (x and y coordinates) and weights REAL*8 XQ(7),YQ(7),WQ(7) C Particular quadrature rule C Number of points INTEGER NQT C Points (x and y coordinates) and weights REAL*8 XQT(LIMNQ),YQT(LIMNQ),WQT(LIMNQ) C Validation and control variables in L3LC LOGICAL LVALID,LFAIL REAL*8 EGEOM,EQRULE LOGICAL LL,LM,LMT,LN C Storage of discrete integral operators on L3LC REAL*8 DISL,DISM,DISMT,DISN C Matrices corresponding to the Laplace integral operators REAL*8 L(LIMN,LIMN),M(LIMN,LIMN) 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*4 I,IQ,J,I1,I2,I3 C Counters INTEGER IPI C Set up constants ZERO=0.0D0 THIRD=1.0D0/3.0D0 HALF=0.5D0 ONE=1.0D0 TWO=2.0D0 THREE=3.0D0 FOUR=4.0D0 EIGHT=8.0D0 PI=FOUR*ATAN(ONE) ROOT2=SQRT(TWO) 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='L3LC.ERR',STATUS='UNKNOWN') C Set up standard Gaussian Quadrature rule. Stores the values C for the 7 point Gaussian quadrature over the standard triangle. C Taken from "Some criteria for numerically integrated matrices and C quadrature formulas for triangles" by M.E.Laursen and M.Gellert in C International Journal for numerical methods in engineering Vol 12 C 67-76 1978. NQ=7 DATA WQ/ 0.225000000000000D0, * 0.125939180544827D0, * 0.125939180544827D0, * 0.125939180544827D0, * 0.132394152788506D0, * 0.132394152788506D0, * 0.132394152788506D0/ DATA XQ/ 0.333333333333333D0, * 0.797426985353087D0, * 0.101286507323456D0, * 0.101286507323456D0, * 0.470142064105115D0, * 0.470142064105115D0, * 0.059715871789770D0/ DATA YQ/ 0.333333333333333D0, * 0.101286507323456D0, * 0.797426985353087D0, * 0.101286507323456D0, * 0.470142064105115D0, * 0.059715871789770D0, * 0.470142064105115D0/ C Set up the geometry of the structure by completing the data C structures VERTPT and ELVERT. C Set up VERTPT, the coordinates of the vertices of the elements NOVERT=14 DATA ((VERTPT(I,J),J=1,3),I=1,14) * /-1.0D0,-1.0D0, 1.0D0, * -1.0D0, 1.0D0, 1.0D0, * -1.0D0, 1.0D0,-1.0D0, * -1.0D0,-1.0D0,-1.0D0, * 1.0D0,-1.0D0, 1.0D0, * 1.0D0, 1.0D0, 1.0D0, * 1.0D0, 1.0D0,-1.0D0, * 1.0D0,-1.0D0,-1.0D0, * 0.0D0, 0.0D0, 1.0D0, * 0.0D0, 1.0D0, 0.0D0, * -1.0D0, 0.0D0, 0.0D0, * 0.0D0, 0.0D0,-1.0D0, * 1.0D0, 0.0D0, 0.0D0, * 0.0D0,-1.0D0, 0.0D0/ C Set the number of elements (N) N=24 C Set up ELEMVERT, the list of elements and the indices of the C vertices that describe the elements. Since the integral equation C that is applied assumes the normal to be outward, the vertices are C enumerated anticlockwise about the element when viewed from a C point near to the element but interior to the surface. This is C so that the normals to the elements are implicitly defined as C being outward within the subroutine L3LC. DATA ((ELVERT(I,J),J=1,3),I=1,24) / * 9, 1, 5, 9, 5, 6, 9, 6, 2, 9, 2, 1, * 10, 2, 6, 10, 6, 7, 10, 7, 3, 10, 3, 2, * 11, 2, 3, 11, 1, 2, 11, 4, 1, 11, 3, 4, * 12, 3, 7, 12, 4, 3, 12, 8, 4, 12, 7, 8, * 13, 7, 6, 13, 8, 7, 13, 5, 8, 13, 6, 5, * 14, 4, 8, 14, 1, 4, 14, 5, 1, 14, 8, 5 / C The outward normals to the elements are stated explicitly. Equally, C they could have been computed from the description of each element DATA ((ELNORM(I,J),J=1,3),I=1,12) / * 0.0D0, 0.0D0, 1.0D0, * 0.0D0, 0.0D0, 1.0D0, * 0.0D0, 0.0D0, 1.0D0, * 0.0D0, 0.0D0, 1.0D0, * 0.0D0, 1.0D0, 0.0D0, * 0.0D0, 1.0D0, 0.0D0, * 0.0D0, 1.0D0, 0.0D0, * 0.0D0, 1.0D0, 0.0D0, * -1.0D0, 0.0D0, 0.0D0, * -1.0D0, 0.0D0, 0.0D0, * -1.0D0, 0.0D0, 0.0D0, * -1.0D0, 0.0D0, 0.0D0/ DATA ((ELNORM(I,J),J=1,3),I=13,24) / * 0.0D0, 0.0D0,-1.0D0, * 0.0D0, 0.0D0,-1.0D0, * 0.0D0, 0.0D0,-1.0D0, * 0.0D0, 0.0D0,-1.0D0, * 1.0D0, 0.0D0, 0.0D0, * 1.0D0, 0.0D0, 0.0D0, * 1.0D0, 0.0D0, 0.0D0, * 1.0D0, 0.0D0, 0.0D0, * 0.0D0,-1.0D0, 0.0D0, * 0.0D0,-1.0D0, 0.0D0, * 0.0D0,-1.0D0, 0.0D0, * 0.0D0,-1.0D0, 0.0D0/ C Set interior points, where an approximation the the solution phi is C sought NPINT=1 PINT(1,1)=0.0D0 PINT(1,2)=0.0D0 PINT(1,3)=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 (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 IPI=1,NPINT IF (ABS(PINT(IPI,1)).GE.ONE.AND.ABS(PINT(IPI,2)).GE.ONE. * AND.ABS(PINT(IPI,3)).GE.ONE) THEN WRITE(*,*) 'ERROR - PINT must be interior to the circle' LERROR=.TRUE. END IF 50 CONTINUE IF (LERROR) STOP C Output description of the test problem WRITE(*,*) 'TEST FOR SUBROUTINE L3LC' WRITE(*,*) '========================' WRITE(*,*) 'Test problem is that of a cube with sides' WRITE(*,*) 'of length 2.0 and centred at the origin with' WRITE(*,*) 'a Dirichlet boundary condition. Solution is via' WRITE(*,*) 'the integral equation (Mk+I/2) \phi=Lk v' WRITE(*,*) 'Number of elements is',N C Loop(I) through each collocation point DO 100 I=1,N C Set P(1..3) (the collocation point) and NORMP(1..3) the normal at C P(1..3) I1=ELVERT(I,1) I2=ELVERT(I,2) I3=ELVERT(I,3) P(1)=(VERTPT(I1,1)+VERTPT(I2,1)+VERTPT(I3,1))/THREE P(2)=(VERTPT(I1,2)+VERTPT(I2,2)+VERTPT(I3,2))/THREE P(3)=(VERTPT(I1,3)+VERTPT(I2,3)+VERTPT(I3,3))/THREE NORMP(1)=ELNORM(I,1) NORMP(2)=ELNORM(I,2) NORMP(3)=ELNORM(I,3) C Loop(J) through each element DO 110 J=1,N C Set QA(1..3), QB(1..3) and QC(1..3), the vertices of the element. C The vertices are enumerated in anti-clockwise direction (when C viewed from a point near to the element and interior to the C surface). Hence the normal will be implicitly defined as outward C from the boundary in subroutine L3LC. QA(1)=VERTPT(ELVERT(J,1),1) QA(2)=VERTPT(ELVERT(J,1),2) QA(3)=VERTPT(ELVERT(J,1),3) QB(1)=VERTPT(ELVERT(J,2),1) QB(2)=VERTPT(ELVERT(J,2),2) QB(3)=VERTPT(ELVERT(J,2),3) QC(1)=VERTPT(ELVERT(J,3),1) QC(2)=VERTPT(ELVERT(J,3),2) QC(3)=VERTPT(ELVERT(J,3),3) C Set LPONEL IF (I.EQ.J) THEN LPONEL=.TRUE. ELSE LPONEL=.FALSE. END IF C Set up quadrature rule data. If LPONEL is false then use the standard C Gaussian quadrature rule above. If LPONEL is true then then a C quadrature rule with 3 times as many points is used, this is made C up from three standard quadrature rules with the quadrature points C translated to the three triangles that each have the cetroid and two C of the original vertices as its vertices. IF (LPONEL) THEN NQT=21 DO 120 IQ=1,7 XQT(IQ)=XQ(IQ)*THIRD+YQ(IQ) YQT(IQ)=XQ(IQ)*THIRD WQT(IQ)=WQ(IQ)/THREE XQT(IQ+7)=XQ(IQ)*THIRD YQT(IQ+7)=XQ(IQ)*THIRD+YQ(IQ) WQT(IQ+7)=WQ(IQ)/THREE XQT(IQ+14)=THIRD*(ONE+TWO*XQ(IQ)-YQ(IQ)) YQT(IQ+14)=THIRD*(ONE-XQ(IQ)+TWO*YQ(IQ)) WQT(IQ+14)=WQ(IQ)/THREE 120 CONTINUE ELSE NQT=7 DO 130 IQ=1,7 XQT(IQ)=XQ(IQ) YQT(IQ)=YQ(IQ) WQT(IQ)=WQ(IQ) 130 CONTINUE END IF C Only integral operators L and M are required LL=.TRUE. LM=.TRUE. LMT=.FALSE. LN=.FALSE. C Value of VECP is irrelevant C Call of L3LC subroutine to compute [Lk], [Mk] CALL L3LC(P,NORMP,QA,QB,QC,LPONEL, * LIMNQ,NQT,XQT,YQT,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) C Here the boundary condition is \phi = z which has the exact solution C vel = n_3 (the z-component of the unit normal to the surface) DO 140 I=1,N I1=ELVERT(I,1) I2=ELVERT(I,2) I3=ELVERT(I,3) Q(1)=(VERTPT(I1,1)+VERTPT(I2,1)+VERTPT(I3,1))/THREE Q(2)=(VERTPT(I1,2)+VERTPT(I2,2)+VERTPT(I3,2))/THREE Q(3)=(VERTPT(I1,3)+VERTPT(I2,3)+VERTPT(I3,3))/THREE NORMQ(1)=ELNORM(I,1) NORMQ(2)=ELNORM(I,2) NORMQ(3)=ELNORM(I,3) PHI(I)=Q(3) VEL(I)=ELNORM(I,3) 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 A PHI = B APPVEL 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 300 IPI=1,NPINT P(1)=PINT(IPI,1) P(2)=PINT(IPI,2) P(3)=PINT(IPI,3) 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 Here \phi = z PHIPT=P(3) C Initialise SUM, this ultimately stores the value of the summation C that approximates to phi at the interior point SUM=ZERO C Loop(J) through the elements DO 200 J=1,N C Set QA(1..3), QB(1..3) and QC(1..3) the vertices of the element. QA(1)=VERTPT(ELVERT(J,1),1) QA(2)=VERTPT(ELVERT(J,1),2) QA(3)=VERTPT(ELVERT(J,1),3) QB(1)=VERTPT(ELVERT(J,2),1) QB(2)=VERTPT(ELVERT(J,2),2) QB(3)=VERTPT(ELVERT(J,2),3) QC(1)=VERTPT(ELVERT(J,3),1) QC(2)=VERTPT(ELVERT(J,3),2) QC(3)=VERTPT(ELVERT(J,3),3) C Set up quadrature rule NQT=7 DO 210 IQ=1,7 WQT(IQ)=WQ(IQ) XQT(IQ)=XQ(IQ) YQT(IQ)=YQ(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 L3LC routine to compute [Lk], [Mk] CALL L3LC(P,VECP,QA,QB,QC,LPONEL, * LIMNQ,NQT,XQT,YQT,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 the interior 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 REAL*8 FUNCTION FNSQRT(X) REAL*8 X FNSQRT=SQRT(X) END