* D06BAF Example Program Text * Mark 20 Release. NAG Copyright 2001. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NEDMX, NVMAX, NLINESX, NUS, NCOMPX, NVIMX, + MAXCAN, LRWORK, LIWORK PARAMETER (NEDMX=300,NVMAX=1000,NLINESX=50,NUS=100, + NCOMPX=5,NVIMX=20,MAXCAN=10000, + LRWORK=12*NVMAX+3*MAXCAN+15, + LIWORK=8*NEDMX+55*NVMAX+MAXCAN+78) * .. Local Scalars .. DOUBLE PRECISION X0, XA, XB, XMAX, XMIN, Y0, YMAX, YMIN INTEGER I, IFAIL, ITRACE, J, K, NCOMP, NEDGE, NELT, + NLINES, NPROPA, NV, NVB, NVINT, REFTK CHARACTER PMESH * .. Local Arrays .. DOUBLE PRECISION COOR(2,NVMAX), COORCH(2,NLINESX), COORUS(2,NUS), + RATE(NLINESX), RUSER(4), RWORK(LRWORK), + WEIGHT(NVIMX) INTEGER CONN(3,2*NVMAX+5), EDGE(3,NEDMX), IUSER(1), + IWORK(LIWORK), LCOMP(NLINESX), LINE(4,NLINESX), + NLCOMP(NCOMPX) * .. External Functions .. DOUBLE PRECISION FBND EXTERNAL FBND * .. External Subroutines .. EXTERNAL D06ABF, D06ACF, D06BAF * .. Intrinsic Functions .. * INTRINSIC ABS * .. Executable Statements .. WRITE (NOUT,*) 'D06BAF Example Program Results' WRITE (NOUT,*) * * Skip heading in data file * READ (NIN,*) * * Initialise boundary mesh inputs: * the number of line and of the characteristic points of * the boundary mesh * READ (NIN,*) NLINES * * The ellipse boundary which envelops the NAg Logo * the N, the A and the G * READ (NIN,*) (COORCH(1,J),J=1,NLINES) READ (NIN,*) (COORCH(2,J),J=1,NLINES) * READ (NIN,*) (COORUS(1,J),J=1,4) READ (NIN,*) (COORUS(2,J),J=1,4) * * The Lines of the boundary mesh * READ (NIN,*) ((LINE(I,J),I=1,4),RATE(J),J=1,NLINES) * * The number of connected components to the boundary * and their informations * READ (NIN,*) NCOMP J = 1 DO 20 I = 1, NCOMP READ (NIN,*) NLCOMP(I) * READ (NIN,*) (LCOMP(K),K=J,J+ABS(NLCOMP(I))-1) J = J + ABS(NLCOMP(I)) 20 CONTINUE * READ (NIN,*) PMESH * * Data passed to the user-supplied function * XMIN = COORCH(1,4) XMAX = COORCH(1,2) YMIN = COORCH(2,1) YMAX = COORCH(2,3) * XA = (XMAX-XMIN)/2.D0 XB = (YMAX-YMIN)/2.D0 * X0 = (XMIN+XMAX)/2.D0 Y0 = (YMIN+YMAX)/2.D0 * RUSER(1) = XA RUSER(2) = XB RUSER(3) = X0 RUSER(4) = Y0 IUSER(1) = 0 * ITRACE = -1 * * Call to the boundary mesh generator * IFAIL = 0 * CALL D06BAF(NLINES,COORCH,LINE,FBND,COORUS,NUS,RATE,NCOMP,NLCOMP, + LCOMP,NVMAX,NEDMX,NVB,COOR,NEDGE,EDGE,ITRACE,RUSER, + IUSER,RWORK,LRWORK,IWORK,LIWORK,IFAIL) * IF (PMESH.EQ.'N') THEN WRITE (NOUT,*) 'Boundary mesh characteristics' WRITE (NOUT,99999) 'NVB =', NVB WRITE (NOUT,99999) 'NEDGE =', NEDGE ELSE IF (PMESH.EQ.'Y') THEN * * Output the mesh to view it using the NAG Graphics Library * WRITE (NOUT,99998) NVB, NEDGE * DO 40 I = 1, NVB WRITE (NOUT,99997) I, COOR(1,I), COOR(2,I) 40 CONTINUE * DO 60 I = 1, NEDGE WRITE (NOUT,99996) I, EDGE(1,I), EDGE(2,I), EDGE(3,I) 60 CONTINUE ELSE WRITE (NOUT,*) 'Problem with the printing option Y or N' STOP END IF * * Initialise mesh control parameters * ITRACE = 0 NPROPA = 1 NVINT = 0 IFAIL = 0 * * Call to the 2D Delaunay-Voronoi mesh generator * CALL D06ABF(NVB,NVINT,NVMAX,NEDGE,EDGE,NV,NELT,COOR,CONN,WEIGHT, + NPROPA,ITRACE,RWORK,LRWORK,IWORK,LIWORK,IFAIL) * IF (PMESH.EQ.'N') THEN WRITE (NOUT,*) + 'Complete mesh (via the 2D Delaunay-Voronoi mesh' WRITE (NOUT,*) 'generator) characteristics' WRITE (NOUT,99999) 'NV =', NV WRITE (NOUT,99999) 'NELT =', NELT ELSE IF (PMESH.EQ.'Y') THEN * * Output the mesh to view it using the NAG Graphics Library * WRITE (NOUT,99998) NV, NELT DO 80 I = 1, NV WRITE (NOUT,99995) COOR(1,I), COOR(2,I) 80 CONTINUE * REFTK = 0 DO 100 K = 1, NELT WRITE (NOUT,99994) CONN(1,K), CONN(2,K), CONN(3,K), REFTK 100 CONTINUE END IF * * Call to the 2D Advancing front mesh generator * IFAIL = 0 * CALL D06ACF(NVB,NVINT,NVMAX,NEDGE,EDGE,NV,NELT,COOR,CONN,WEIGHT, + ITRACE,RWORK,LRWORK,IWORK,LIWORK,IFAIL) * IF (PMESH.EQ.'N') THEN WRITE (NOUT,*) 'Complete mesh (via the 2D Advancing front mesh' WRITE (NOUT,*) 'generator) characteristics' WRITE (NOUT,99999) 'NV =', NV WRITE (NOUT,99999) 'NELT =', NELT ELSE IF (PMESH.EQ.'Y') THEN * * Output the mesh to view it using the NAG Graphics Library * WRITE (NOUT,99998) NV, NELT DO 120 I = 1, NV WRITE (NOUT,99995) COOR(1,I), COOR(2,I) 120 CONTINUE * REFTK = 0 DO 140 K = 1, NELT WRITE (NOUT,99994) CONN(1,K), CONN(2,K), CONN(3,K), REFTK 140 CONTINUE END IF * STOP * 99999 FORMAT (1X,A,I6) 99998 FORMAT (1X,2I10) 99997 FORMAT (2X,I4,2(2X,E12.6)) 99996 FORMAT (1X,4I4) 99995 FORMAT (2(2X,E12.6)) 99994 FORMAT (1X,4I10) END * DOUBLE PRECISION FUNCTION FBND(I,X,Y,RUSER,IUSER) * .. Scalar Arguments .. DOUBLE PRECISION X, Y INTEGER I * .. Array Arguments .. DOUBLE PRECISION RUSER(*) INTEGER IUSER(*) * .. Local Scalars .. * DOUBLE PRECISION RADIUS2, X0, XA, XB, Y0 * .. Executable Statements .. XA = RUSER(1) XB = RUSER(2) X0 = RUSER(3) Y0 = RUSER(4) * FBND = 0.D0 IF (I.EQ.1) THEN * * line 1,2,3, and 4: ellipse centred in (X0,Y0) with * XA and XB as coefficients * FBND = ((X-X0)/XA)**2 + ((Y-Y0)/XB)**2 - 1.D0 ELSE IF (I.EQ.2) THEN * * line 24, 27, 33 and 38 are a circle centred in (X0,Y0) * with radius SQRT(RADIUS2) * X0 = 20.5D0 Y0 = 4.D0 RADIUS2 = 4.25D0 FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2 ELSE IF (I.EQ.3) THEN X0 = 17.D0 Y0 = 8.5D0 RADIUS2 = 5.D0 FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2 ELSE IF (I.EQ.4) THEN X0 = 17.D0 Y0 = 8.5D0 RADIUS2 = 5.D0 FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2 ELSE IF (I.EQ.5) THEN X0 = 19.5D0 Y0 = 4.D0 RADIUS2 = 1.25D0 FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2 END IF * RETURN END