* Tutorial Example Program Text * Mark 20 Release. NAG Copyright 2001 C .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NEDMX, NVMAX, NELTMX, NLINEX, NUS, NCOMPX, NVIMX, + MAXCAN, NNZMAX, NVFXMX, LRWORK, LIWORK PARAMETER (NEDMX=1000,NVMAX=6000,NELTMX=2*NVMAX+5, + NLINEX=50,NUS=100,NCOMPX=5,NVIMX=20,MAXCAN=10000, + NNZMAX=50*NVMAX,NVFXMX=20, + LRWORK=12*NVMAX+3*MAXCAN+15, + LIWORK=8*NEDMX+55*NVMAX+MAXCAN+78) DOUBLE PRECISION C12, C24, HALF, MHALF, ZERO PARAMETER (C12=1.D0/12.D0,C24=1.D0/24.D0,HALF=0.5D0, + MHALF=-HALF,ZERO=0.D0) C .. Local Arrays .. DOUBLE PRECISION AMAT(NNZMAX), CHKMAT(NVMAX), COOR(2,NVMAX), + COOR1(2,NVMAX), COOR2(2,NVMAX), COOR3(2,NVMAX), + COOR4(2,NVMAX), COORCH(2,NLINEX), COORUS(2,NUS), + DXDX(3,3), DXDY(3,3), DYDX(3,3), DYDY(3,3), + F(NVMAX), PHI(3,3), RATE(NLINEX), RHS(NVMAX), + RUSER(1), RWORK(LRWORK), TRANS(6,1), UDIR(NVMAX), + UNEU(NVMAX), USOL(NVMAX), UTRUE(NVMAX), + WEIGHT(NVIMX) INTEGER CONN(3,NELTMX), CONN1(3,NELTMX), CONN2(3,NELTMX), + CONN3(3,NELTMX), CONN4(3,NELTMX), EDGE(3,NEDMX), + EDGE1(3,NEDMX), EDGE2(3,NEDMX), EDGE3(3,NEDMX), + EDGE4(3,NEDMX), ICOL(NNZMAX), IROW(NNZMAX), + IROW1(NVMAX+1), ITYPE(1), IUSER(1), + IWORK(LIWORK), LCOMP(NLINEX), LINE(4,NLINEX), + NLCOMP(NCOMPX), NUMFIX(NVFXMX), REFT(NELTMX), + REFT1(NELTMX), REFT2(NELTMX), REFT3(NELTMX), + REFT4(NELTMX) C .. Local Scalars .. DOUBLE PRECISION ALPHA, ALPHAK, BETA, BETAK, DL12, DTOL, EPS, FK1, + FK2, FK3, GAMMAK, H1, H1NRM, H1NRM0, H2, H3, + HSIZE, JFK, L2NRM, L2NRM0, RNORM, TOL, XK1, XK2, + XK3, XL1, XL2, YK1, YK2, YK3, YL1, YL2 INTEGER I, IFAIL, IJ, IK, ITN, ITRACE, J, JK, K, K1, K2, + K3, L, L1, L2, LFILL, LL1, LL2, M, MAXITN, NCOMP, + NEDGE, NEDGE1, NEDGE2, NEDGE3, NEDGE4, NELT, + NELT1, NELT2, NELT3, NELT4, NLINES, NNZ, NNZC, + NPIVM, NPROPA, NQINT, NTRANS, NV, NV1, NV2, NV3, + NV4, NVB1, NVB2, NVFIX, NVINT CHARACTER CHECK, MILU, PSTRAT, SOLUT, TRANS1 CHARACTER*2 TYPE CHARACTER*8 METHOD C .. External Subroutines .. EXTERNAL D06ABF, D06ACF, D06BAF, D06CAF, D06CCF, D06DAF, + D06DBF, F06DBF, F11DAF, F11DBF, F11DCF C .. External Functions .. DOUBLE PRECISION FBND, X02AJF EXTERNAL FBND, X02AJF C .. Intrinsic Functions .. INTRINSIC ABS, DABS, MAX, SQRT C .. Data statements .. DATA PHI/C12, C24, C24, C24, C12, C24, C24, C24, C12/ DATA DXDX/HALF, MHALF, ZERO, MHALF, HALF, ZERO, ZERO, + ZERO, ZERO/ DATA DYDY/HALF, ZERO, MHALF, ZERO, ZERO, ZERO, MHALF, + ZERO, HALF/ DATA DXDY/HALF, MHALF, ZERO, ZERO, ZERO, ZERO, MHALF, + HALF, ZERO/ DATA DYDX/HALF, ZERO, MHALF, MHALF, ZERO, HALF, ZERO, + ZERO, ZERO/ C .. Executable Statements .. WRITE (NOUT,FMT='(A)') 'Tutorial Example Results ' * * Build the mesh of the 1st domain * * Initialise boundary mesh inputs: * the number of line and characteristic points on * the boundary mesh * READ (NIN,FMT=*) * READ (NIN,FMT=*) NLINES * * Characteristic points on the boundary geometry * READ (NIN,FMT=*) (COORCH(1,J),J=1,NLINES) READ (NIN,FMT=*) (COORCH(2,J),J=1,NLINES) * * The lines of the boundary mesh * READ (NIN,FMT=*) ((LINE(I,J),I=1,4),RATE(J),J=1,NLINES) * * The number of connected components on the boundary * and their specification * READ (NIN,FMT=*) NCOMP J = 1 DO 20 I = 1, NCOMP READ (NIN,FMT=*) NLCOMP(I) * READ (NIN,FMT=*) (LCOMP(K),K=J,J+ABS(NLCOMP(I))-1) J = J + ABS(NLCOMP(I)) 20 CONTINUE * * Call to the 2D boundary mesh generator * ITRACE = 0 IFAIL = 0 * CALL D06BAF(NLINES,COORCH,LINE,FBND,COORUS,NUS,RATE,NCOMP,NLCOMP, + LCOMP,NVMAX,NEDMX,NVB1,COOR1,NEDGE1,EDGE1,ITRACE, + RUSER,IUSER,RWORK,LRWORK,IWORK,LIWORK,IFAIL) * * Mesh using Delaunay-Voronoi method * * Initialise mesh control parameters * ITRACE = 0 NPROPA = 1 NVINT = 0 IFAIL = 0 * * Call to the 2D Delaunay-Voronoi mesh generator * CALL D06ABF(NVB1,NVINT,NVMAX,NEDGE1,EDGE1,NV1,NELT1,COOR1,CONN1, + WEIGHT,NPROPA,ITRACE,RWORK,LRWORK,IWORK,LIWORK,IFAIL) * * build the mesh of the 2nd domain * * Initialise boundary mesh inputs: * the number of line and characteristic points on * the boundary mesh * READ (NIN,FMT=*) NLINES * * Characteristic points on the boundary geometry * READ (NIN,FMT=*) (COORCH(1,J),J=1,NLINES) READ (NIN,FMT=*) (COORCH(2,J),J=1,NLINES) * * The lines of the boundary mesh * READ (NIN,FMT=*) ((LINE(I,J),I=1,4),RATE(J),J=1,NLINES) * * The number of connected components on the boundary * and their specification * READ (NIN,FMT=*) NCOMP J = 1 DO 40 I = 1, NCOMP READ (NIN,FMT=*) NLCOMP(I) * READ (NIN,FMT=*) (LCOMP(K),K=J,J+ABS(NLCOMP(I))-1) J = J + ABS(NLCOMP(I)) 40 CONTINUE * * Solution method 'D' direct method or 'I' iterative method * READ (NIN,FMT=*) SOLUT * * The iterative method to be used RGMRES CGS BICGSTAB or TFQMR * READ (NIN,FMT=*) METHOD * * Type of problem which will be solved P1 or P2 * READ (NIN,FMT=*) TYPE * ITRACE = 0 * * Call to the 2D boundary mesh generator * IFAIL = 0 * CALL D06BAF(NLINES,COORCH,LINE,FBND,COORUS,NUS,RATE,NCOMP,NLCOMP, + LCOMP,NVMAX,NEDMX,NVB2,COOR2,NEDGE2,EDGE2,ITRACE, + RUSER,IUSER,RWORK,LRWORK,IWORK,LIWORK,IFAIL) * * Mesh using the advancing front method * * Initialise mesh control parameters * ITRACE = 0 NVINT = 0 IFAIL = 0 * * Call to the 2D advancing front mesh generator * CALL D06ACF(NVB2,NVINT,NVMAX,NEDGE2,EDGE2,NV2,NELT2,COOR2,CONN2, + WEIGHT,ITRACE,RWORK,LRWORK,IWORK,LIWORK,IFAIL) * * Rotation of the 2nd domain mesh to produce * the 3rd mesh domain * NTRANS = 1 ITYPE(1) = 3 TRANS(1,1) = 6.D0 TRANS(2,1) = -1.D0 TRANS(3,1) = -90.D0 ITRACE = 0 IFAIL = 0 * CALL D06DAF(NV2,NEDGE2,NELT2,NTRANS,ITYPE,TRANS,COOR2,EDGE2,CONN2, + COOR3,EDGE3,CONN3,ITRACE,RWORK,LRWORK,IFAIL) * NV3 = NV2 NELT3 = NELT2 NEDGE3 = NEDGE2 * * Tag triangles and boundary edges for all 3 domains * DO 60 K = 1, NELT1 REFT1(K) = 1 60 CONTINUE * DO 80 K = 1, NELT2 REFT2(K) = 2 80 CONTINUE * DO 100 K = 1, NELT3 REFT3(K) = 3 100 CONTINUE * DO 120 K = 1, NEDGE1 EDGE1(3,K) = 1 120 CONTINUE * DO 140 K = 1, NEDGE2 EDGE2(3,K) = 2 140 CONTINUE * DO 160 K = 1, NEDGE3 EDGE3(3,K) = 3 160 CONTINUE * * Restitching of the mesh 1 and 2 to form the mesh 4 * EPS = 1.D-3 ITRACE = 0 IFAIL = 0 * CALL D06DBF(EPS,NV1,NELT1,NEDGE1,COOR1,EDGE1,CONN1,REFT1,NV2, + NELT2,NEDGE2,COOR2,EDGE2,CONN2,REFT2,NV4,NELT4,NEDGE4, + COOR4,EDGE4,CONN4,REFT4,ITRACE,IWORK,LIWORK,IFAIL) * * Restitching of the mesh 3 and 4 to form the final mesh * ITRACE = 0 IFAIL = 0 * CALL D06DBF(EPS,NV4,NELT4,NEDGE4,COOR4,EDGE4,CONN4,REFT4,NV3, + NELT3,NEDGE3,COOR3,EDGE3,CONN3,REFT3,NV,NELT,NEDGE, + COOR,EDGE,CONN,REFT,ITRACE,IWORK,LIWORK,IFAIL) * * Call the smoothing routine * NVFIX = 0 NQINT = 10 * CALL D06CAF(NV,NELT,NEDGE,COOR,EDGE,CONN,NVFIX,NUMFIX,ITRACE, + NQINT,IWORK,LIWORK,RWORK,LRWORK,IFAIL) * * Call the renumbering routine and get the sparsity * IFAIL = 0 ITRACE = 0 CALL D06CCF(NV,NELT,NEDGE,NNZMAX,NNZ,COOR,EDGE,CONN,IROW,ICOL, + ITRACE,IWORK,LIWORK,RWORK,LRWORK,IFAIL) * * Special tag for the Neumann boundary * EPS = 1.D-5 DO 180 K = 1, NEDGE IF ((DABS(COOR(1,EDGE(1,K))**2+COOR(2,EDGE(1,K))**2-1.D0) + .LT.EPS) .AND. (DABS(COOR(1,EDGE(2,K))**2+COOR(2,EDGE(2,K)) + **2-1.D0).LT.EPS)) THEN EDGE(3,K) = 5 END IF 180 CONTINUE * * Find number of elements in each row. * CALL F06DBF(NV,0,IWORK,1) * DO 200 I = 1, NNZ IWORK(IROW(I)) = IWORK(IROW(I)) + 1 AMAT(I) = 0.D0 200 CONTINUE * * Set up pointer to start of each row and * set data for the Helmholtz problem * ALPHA = 1.D0 BETA = 1.D0 IROW1(1) = 1 * DO 220 I = 1, NV IROW1(I+1) = IROW1(I) + IWORK(I) RHS(I) = 0.D0 CHKMAT(I) = 0.D0 * IF (TYPE.EQ.'P2') THEN UTRUE(I) = COOR(1,I)**2 + COOR(2,I)**2 F(I) = BETA*UTRUE(I) - 4.D0*ALPHA ELSE IF (TYPE.EQ.'P1') THEN UTRUE(I) = 1.D0 + COOR(1,I) + COOR(2,I) F(I) = BETA*UTRUE(I) END IF * UDIR(I) = 0.D0 UNEU(I) = 0.D0 USOL(I) = 0.D0 220 CONTINUE * * Boundary conditions * DO 260 L = 1, NEDGE L1 = EDGE(1,L) L2 = EDGE(2,L) * IF (EDGE(3,L).EQ.1 .OR. EDGE(3,L).EQ.2 .OR. EDGE(3,L).EQ.3) + THEN UDIR(L1) = UTRUE(L1) UDIR(L2) = UTRUE(L2) * USOL(L1) = UDIR(L1) USOL(L2) = UDIR(L2) ELSE IF (EDGE(3,L).EQ.5) THEN IF (TYPE.EQ.'P2') THEN UNEU(L1) = -2.D0 UNEU(L2) = -2.D0 ELSE IF (TYPE.EQ.'P1') THEN UNEU(L1) = -COOR(1,L1) - COOR(2,L1) UNEU(L2) = -COOR(1,L2) - COOR(2,L2) END IF END IF 260 CONTINUE * * Compute the FE matrix and the right-hand side of the linear system * HSIZE = -1.D0 DO 340 K = 1, NELT K1 = CONN(1,K) K2 = CONN(2,K) K3 = CONN(3,K) * XK1 = COOR(1,K1) YK1 = COOR(2,K1) * XK2 = COOR(1,K2) YK2 = COOR(2,K2) * XK3 = COOR(1,K3) YK3 = COOR(2,K3) * H1 = SQRT((XK2-XK3)**2+(YK2-YK3)**2) H2 = SQRT((XK1-XK3)**2+(YK1-YK3)**2) H3 = SQRT((XK2-XK1)**2+(YK2-YK1)**2) * HSIZE = MAX(H1,HSIZE) HSIZE = MAX(H2,HSIZE) HSIZE = MAX(H3,HSIZE) * FK1 = F(K1) FK2 = F(K2) FK3 = F(K3) * JFK = (XK2-XK1)*(YK3-YK1) - (XK3-XK1)*(YK2-YK1) ALPHAK = (YK3-YK1)**2 + (XK3-XK1)**2 BETAK = (YK2-YK1)**2 + (XK2-XK1)**2 GAMMAK = (YK3-YK1)*(YK1-YK2) + (XK1-XK3)*(XK2-XK1) * DO 320 I = 1, 3 IK = CONN(I,K) CHKMAT(IK) = CHKMAT(IK) + BETA*JFK/6.D0 * * Element contribution to the FE matrix * DO 300 J = 1, 3 JK = CONN(J,K) IFAIL = 0 * L2NRM = 1.D0/JFK*(ALPHAK*DXDX(I,J)+GAMMAK*(DXDY(I,J) + +DYDX(I,J))+BETAK*DYDY(I,J)) L2NRM0 = JFK*PHI(I,J) * DO 280 IJ = IROW1(IK), IROW1(IK+1) - 1 IF (ICOL(IJ).EQ.JK) THEN AMAT(IJ) = AMAT(IJ) + ALPHA*L2NRM + BETA*L2NRM0 END IF 280 CONTINUE 300 CONTINUE * * Right-hand side contributions * RHS(IK) = RHS(IK) + JFK*(FK1*PHI(I,1)+FK2*PHI(I,2) + +FK3*PHI(I,3)) 320 CONTINUE 340 CONTINUE * * Check of the FE matrix * EPS = 1.D+1*X02AJF() * DO 380 I = 1, NV L2NRM = 0.D0 DO 360 J = IROW1(I), IROW1(I+1) - 1 L2NRM = L2NRM + AMAT(J) 360 CONTINUE * IF (ABS(CHKMAT(I)-L2NRM).GT.EPS) THEN WRITE (*,FMT=*) ' Prob with FE Matrix ', I, IROW1(I), + IROW1(I+1) - 1, L2NRM, CHKMAT(I) END IF 380 CONTINUE * * Set the boundary conditions * IFAIL = 0 DO 440 L = 1, NEDGE L1 = EDGE(1,L) L2 = EDGE(2,L) * XL1 = COOR(1,L1) YL1 = COOR(2,L1) * XL2 = COOR(1,L2) YL2 = COOR(2,L2) * DL12 = SQRT((XL2-XL1)**2+(YL2-YL1)**2) * * Dirichlet * IF (EDGE(3,L).EQ.1 .OR. EDGE(3,L).EQ.2 .OR. EDGE(3,L).EQ.3) + THEN RHS(L1) = UDIR(L1) RHS(L2) = UDIR(L2) DO 400 LL1 = IROW1(L1), IROW1(L1+1) - 1 AMAT(LL1) = 0.D0 * IF (ICOL(LL1).EQ.L1) THEN AMAT(LL1) = 1.D0 END IF 400 CONTINUE * DO 420 LL2 = IROW1(L2), IROW1(L2+1) - 1 AMAT(LL2) = 0.D0 * IF (ICOL(LL2).EQ.L2) THEN AMAT(LL2) = 1.D0 END IF 420 CONTINUE * * Neumann * ELSE IF (EDGE(3,L).EQ.5) THEN RHS(L1) = RHS(L1) + 0.5D0*ALPHA*DL12*UNEU(L1) RHS(L2) = RHS(L2) + 0.5D0*ALPHA*DL12*UNEU(L2) END IF 440 CONTINUE * IF (SOLUT.EQ.'D') THEN * * Direct solution * Compute the LU factorization of the FE matrix * LFILL = -1 DTOL = 0.D0 PSTRAT = 'C' MILU = 'N' IFAIL = 0 * CALL F11DAF(NV,NNZ,AMAT,NNZMAX,IROW,ICOL,LFILL,DTOL,PSTRAT, + MILU,IWORK,IWORK(NV+1),IROW1,IWORK(2*NV+1),NNZC, + NPIVM,IWORK(3*NV+2),(LIWORK-3*NV-1),IFAIL) * * Solve the system AMAT*USOL = RHS * TRANS1 = 'N' CHECK = 'C' IFAIL = 0 * CALL F11DBF(TRANS1,NV,AMAT,NNZMAX,IROW,ICOL,IWORK,IWORK(NV+1), + IROW1,IWORK(2*NV+1),CHECK,RHS,USOL,IFAIL) * WRITE (NOUT,FMT='(A,I10)') 'Dimension of the system = ', NV WRITE (NOUT,FMT='(A,A,I10)') 'Number of triangles of the ', + 'mesh = ', NELT WRITE (NOUT,FMT='(A,I10)') 'Number of edges = ', NEDGE ELSE IF (SOLUT.EQ.'I') THEN * * Iterative solution * Compute the incomplete LU factorisation * LFILL = 1 DTOL = 0.D0 PSTRAT = 'C' MILU = 'N' IFAIL = 0 * CALL F11DAF(NV,NNZ,AMAT,NNZMAX,IROW,ICOL,LFILL,DTOL,PSTRAT, + MILU,IWORK(1),IWORK(NV+1),IROW1,IWORK(2*NV+1),NNZC, + NPIVM,IWORK(3*NV+2),(LIWORK-3*NV-1),IFAIL) * * Solve the system AMAT*USOL = RHS * TOL = 1.D-12 MAXITN = NV M = 10 IFAIL = -1 * CALL F11DCF(METHOD,NV,NNZ,AMAT,NNZMAX,IROW,ICOL,IWORK(1), + IWORK(NV+1),IROW1,IWORK(2*NV+1),RHS,M,TOL,MAXITN, + USOL,RNORM,ITN,RWORK,LRWORK,IFAIL) * WRITE (NOUT,FMT='(A,I10)') 'Dimension of the system = ', NV WRITE (NOUT,FMT='(A,A,I10)') 'Number of triangles of the ', + 'mesh = ', NELT WRITE (NOUT,FMT='(A,I10)') 'Number of edges = ', NEDGE WRITE (NOUT,FMT='(A,I10,A)') 'Converged in', ITN, ' iterations' WRITE (NOUT,FMT='(A,1P,D16.3)') 'Final residual norm =', RNORM WRITE (NOUT,FMT=*) END IF * * L2 and H1 Error estimation * L2NRM = 0.D0 L2NRM0 = 0.D0 H1NRM = 0.D0 H1NRM0 = 0.D0 DO 460 K = 1, NELT K1 = CONN(1,K) K2 = CONN(2,K) K3 = CONN(3,K) * XK1 = COOR(1,K1) YK1 = COOR(2,K1) * XK2 = COOR(1,K2) YK2 = COOR(2,K2) * XK3 = COOR(1,K3) YK3 = COOR(2,K3) * JFK = (XK2-XK1)*(YK3-YK1) - (XK3-XK1)*(YK2-YK1) ALPHAK = (YK3-YK1)**2 + (XK3-XK1)**2 BETAK = (YK2-YK1)**2 + (XK2-XK1)**2 GAMMAK = (YK3-YK1)*(YK1-YK2) + (XK1-XK3)*(XK2-XK1) * FK1 = USOL(K1) FK2 = USOL(K2) FK3 = USOL(K3) * XK1 = UTRUE(K1) XK2 = UTRUE(K2) XK3 = UTRUE(K3) * L2NRM0 = L2NRM0 + JFK*(XK1*XK1+XK2*XK2+XK3*XK3)/6.D0 H1NRM0 = H1NRM0 + 1.D0/(2.D0*JFK)*(ALPHAK*(XK2-XK1)*(XK2-XK1) + +2.D0*GAMMAK*(XK2-XK1)*(XK3-XK1)+BETAK*(XK3-XK1) + *(XK3-XK1)) * XK1 = XK1 - FK1 XK2 = XK2 - FK2 XK3 = XK3 - FK3 * L2NRM = L2NRM + JFK*(XK1*XK1+XK2*XK2+XK3*XK3)/6.D0 H1NRM = H1NRM + 1.D0/(2.D0*JFK)*(ALPHAK*(XK2-XK1)*(XK2-XK1) + +2.D0*GAMMAK*(XK2-XK1)*(XK3-XK1)+BETAK*(XK3-XK1) + *(XK3-XK1)) 460 CONTINUE * L2NRM = SQRT(L2NRM) L2NRM0 = SQRT(L2NRM0) H1NRM = SQRT(H1NRM) H1NRM0 = SQRT(H1NRM0) ALPHAK = L2NRM/L2NRM0 BETAK = H1NRM/H1NRM0 * WRITE (NOUT,FMT=*) 'L2 Error = ', L2NRM WRITE (NOUT,FMT=*) 'Relative L2 Error = ', ALPHAK WRITE (NOUT,FMT=*) 'H1 Error = ', H1NRM WRITE (NOUT,FMT=*) 'Relative H1 Error = ', BETAK WRITE (NOUT,FMT=*) ' Mesh step = ', HSIZE * STOP END * DOUBLE PRECISION FUNCTION FBND(I,X,Y,RUSER,IUSER) C .. Scalar Arguments .. DOUBLE PRECISION X, Y INTEGER I C .. Array Arguments .. DOUBLE PRECISION RUSER(*) INTEGER IUSER(*) C .. Local Scalars .. DOUBLE PRECISION R2, X0, Y0 C .. Executable Statements .. FBND = 0.D0 IF (I.EQ.1) THEN * * inner circle * X0 = 0.D0 Y0 = 0.D0 R2 = 1.D0 FBND = (X-X0)**2 + (Y-Y0)**2 - R2 ELSE IF (I.EQ.2) THEN * * outer circle * X0 = 0.D0 Y0 = 0.D0 R2 = 5.D0 FBND = (X-X0)**2 + (Y-Y0)**2 - R2 END IF * RETURN END