*     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
