*     D06DBF Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. External Subroutines ..
      EXTERNAL         EX1, EX2
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D06DBF Example Program Results'
      CALL EX1
      CALL EX2
      STOP
      END
*
      SUBROUTINE EX1
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NBEDMX, NVMAX, NELTMAX, LIWORK, NTRANS, LRWORK
      PARAMETER        (NBEDMX=200,NVMAX=900,NELTMAX=2*(NVMAX-1),
     +                 LIWORK=4*NVMAX+2*NELTMAX+2*NBEDMX,NTRANS=1,
     +                 LRWORK=12*NTRANS)
*     .. Local Scalars ..
      DOUBLE PRECISION EPS
      INTEGER          I, IFAIL, IMAX, ITRACE, ITRANS, JMAX, JTRANS, K,
     +                 KTRANS, NEDGE1, NEDGE2, NEDGE3, NELT1, NELT2,
     +                 NELT3, NV1, NV2, NV3, REFTK
      CHARACTER        PMESH
*     .. Local Arrays ..
      DOUBLE PRECISION COOR1(2,NVMAX), COOR2(2,NVMAX), COOR3(2,NVMAX),
     +                 RWORK(LRWORK), TRANS(6,NTRANS)
      INTEGER          CONN1(3,NELTMAX), CONN2(3,NELTMAX),
     +                 CONN3(3,NELTMAX), EDGE1(3,NBEDMX),
     +                 EDGE2(3,NBEDMX), EDGE3(3,NBEDMX), ITYPE(NTRANS),
     +                 IWORK(LIWORK), REFT1(NELTMAX), REFT2(NELTMAX),
     +                 REFT3(NELTMAX)
*     .. External Subroutines ..
      EXTERNAL         D06DAF, D06DBF
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Example 1'
      WRITE (NOUT,*)
      IMAX = 20
      JMAX = IMAX
*
*     Skip heading in data file
*
      READ (NIN,*)
      READ (NIN,*)
*
*     Read the mesh : coordinates and connectivity of the 1st domain
*
      READ (NIN,*) NV1, NELT1
      NV2 = NV1
      NELT2 = NELT1
      DO 20 I = 1, NV1
         READ (NIN,*) COOR1(1,I), COOR1(2,I)
   20 CONTINUE
*
      DO 40 K = 1, NELT1
         READ (NIN,*) CONN1(1,K), CONN1(2,K), CONN1(3,K), REFTK
         REFT1(K) = 1
         REFT2(K) = 2
   40 CONTINUE
*
      READ (NIN,*) PMESH
*
*     the Edges of the boundary
*
      NEDGE1 = 0
      DO 60 I = 1, IMAX - 1
         NEDGE1 = NEDGE1 + 1
         EDGE1(1,NEDGE1) = I
         EDGE1(2,NEDGE1) = I + 1
         EDGE1(3,NEDGE1) = 1
   60 CONTINUE
*
      DO 80 I = 1, JMAX - 1
         NEDGE1 = NEDGE1 + 1
         EDGE1(1,NEDGE1) = I*IMAX
         EDGE1(2,NEDGE1) = (I+1)*IMAX
         EDGE1(3,NEDGE1) = 1
   80 CONTINUE
*
      DO 100 I = 1, IMAX - 1
         NEDGE1 = NEDGE1 + 1
         EDGE1(1,NEDGE1) = IMAX*JMAX - I + 1
         EDGE1(2,NEDGE1) = IMAX*JMAX - I
         EDGE1(3,NEDGE1) = 1
  100 CONTINUE
*
      DO 120 I = 1, JMAX - 1
         NEDGE1 = NEDGE1 + 1
         EDGE1(1,NEDGE1) = (JMAX-I)*IMAX + 1
         EDGE1(2,NEDGE1) = (JMAX-I-1)*IMAX + 1
         EDGE1(3,NEDGE1) = 1
  120 CONTINUE
*
      NEDGE2 = NEDGE1
*
      DO 220 KTRANS = 1, 2
*
*     Translation of the 1st domain to obtain the 2nd domain
*     KTRANS = 1 leading to a domains overlapping
*     KTRANS = 2 leading to a domains partition
*
         IF (KTRANS.EQ.1) THEN
            ITRANS = IMAX - 5
            JTRANS = JMAX - 3
         ELSE
            ITRANS = IMAX - 1
            JTRANS = 0
         END IF
*
         ITYPE(1) = 1
         TRANS(1,1) = DBLE(ITRANS)/DBLE(IMAX-1)
         TRANS(2,1) = DBLE(JTRANS)/DBLE(JMAX-1)
         ITRACE = 0
         IFAIL = 0
*
         CALL D06DAF(NV2,NEDGE2,NELT2,NTRANS,ITYPE,TRANS,COOR1,EDGE1,
     +               CONN1,COOR2,EDGE2,CONN2,ITRACE,RWORK,LRWORK,IFAIL)
*
         DO 140 I = 1, NEDGE2
            EDGE2(3,I) = 2
  140    CONTINUE
*
*     Call to the restitching driver
*
         ITRACE = 0
         IFAIL = 0
         EPS = 1.D-2
*
         CALL D06DBF(EPS,NV1,NELT1,NEDGE1,COOR1,EDGE1,CONN1,REFT1,NV2,
     +               NELT2,NEDGE2,COOR2,EDGE2,CONN2,REFT2,NV3,NELT3,
     +               NEDGE3,COOR3,EDGE3,CONN3,REFT3,ITRACE,IWORK,LIWORK,
     +               IFAIL)
*
         IF (PMESH.EQ.'N') THEN
            IF (KTRANS.EQ.1) THEN
               WRITE (NOUT,*) 'The restitched mesh characteristics'
               WRITE (NOUT,*) 'in the overlapping case'
            ELSE
               WRITE (NOUT,*) 'in the partition case'
            END IF
            WRITE (NOUT,99999) 'NV    =', NV3
            WRITE (NOUT,99999) 'NELT  =', NELT3
            WRITE (NOUT,99999) 'NEDGE =', NEDGE3
         ELSE IF (PMESH.EQ.'Y') THEN
*
*     Output the mesh to view it using the NAG Graphics Library
*
            WRITE (NOUT,99998) NV3, NELT3, NEDGE3
            DO 160 I = 1, NV3
               WRITE (NOUT,99997) COOR3(1,I), COOR3(2,I)
  160       CONTINUE
*
            DO 180 K = 1, NELT3
               WRITE (NOUT,99996) CONN3(1,K), CONN3(2,K), CONN3(3,K),
     +           REFT3(K)
  180       CONTINUE
*
            DO 200 K = 1, NEDGE3
               WRITE (NOUT,99998) EDGE3(1,K), EDGE3(2,K), EDGE3(3,K)
  200       CONTINUE
         ELSE
            WRITE (NOUT,*) 'Problem with the printing option Y or N'
            STOP
         END IF
  220 CONTINUE
*
99999 FORMAT (1X,A,I6)
99998 FORMAT (1X,3I10)
99997 FORMAT (2(2X,E12.6))
99996 FORMAT (1X,4I10)
      END
*
      SUBROUTINE EX2
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NEDMX, NVMAX, NLINESX, NUS, NCOMPX, NVIMX,
     +                 MAXCAN, LRWORK, LIWORK
      PARAMETER        (NEDMX=200,NVMAX=700,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 EPS
      INTEGER          I, IFAIL, ITRACE, J, K, NCOMP, NEDGE1, NEDGE2,
     +                 NEDGE3, NEDGE4, NEDGE5, NELT1, NELT2, NELT3,
     +                 NELT4, NELT5, NLINES, NPROPA, NQINT, NTRANS, NV1,
     +                 NV2, NV3, NV4, NV5, NVB1, NVB2, NVFIX, NVINT
      CHARACTER        PMESH
*     .. Local Arrays ..
      DOUBLE PRECISION COOR1(2,NVMAX), COOR2(2,NVMAX), COOR3(2,NVMAX),
     +                 COOR4(2,NVMAX), COOR5(2,NVMAX),
     +                 COORCH(2,NLINESX), COORUS(2,NUS), RATE(NLINESX),
     +                 RUSER(1), RWORK(LRWORK), TRANS(6,1),
     +                 WEIGHT(NVIMX)
      INTEGER          CONN1(3,2*NVMAX+5), CONN2(3,2*NVMAX+5),
     +                 CONN3(3,2*NVMAX+5), CONN4(3,2*NVMAX+5),
     +                 CONN5(3,2*NVMAX+5), EDGE1(3,NEDMX),
     +                 EDGE2(3,NEDMX), EDGE3(3,NEDMX), EDGE4(3,NEDMX),
     +                 EDGE5(3,NEDMX), ITYPE(1), IUSER(1),
     +                 IWORK(LIWORK), LCOMP(NLINESX), LINE(4,NLINESX),
     +                 NLCOMP(NCOMPX), NUMFIX(1), REFT1(2*NVMAX+5),
     +                 REFT2(2*NVMAX+5), REFT3(2*NVMAX+5),
     +                 REFT4(2*NVMAX+5), REFT5(2*NVMAX+5)
*     .. External Functions ..
      DOUBLE PRECISION FBND
      EXTERNAL         FBND
*     .. External Subroutines ..
      EXTERNAL         D06ABF, D06ACF, D06BAF, D06CAF, D06DAF, D06DBF
*     .. Intrinsic Functions ..
      INTRINSIC        ABS
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Example 2'
      WRITE (NOUT,*)
*
*     Skip heading in data file
*
      READ (NIN,*)
*
*     Build the mesh of the 1st domain
*
*     Initialise boundary mesh inputs:
*     the number of line and of the characteristic points of
*     the boundary mesh
*
*
      READ (NIN,*) NLINES
*
*     Characteristic points of the boundary geometry
*
      READ (NIN,*) (COORCH(1,J),J=1,NLINES)
      READ (NIN,*) (COORCH(2,J),J=1,NLINES)
*
*     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
*
      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,NVB1,COOR1,NEDGE1,EDGE1,ITRACE,
     +            RUSER,IUSER,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
*
*     mesh it 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)
*
      DO 40 K = 1, NELT1
         REFT1(K) = 1
   40 CONTINUE
*
*     Call the smoothing routine
*
      NVFIX = 0
      NQINT = 10
      IFAIL = 0
*
      CALL D06CAF(NV1,NELT1,NEDGE1,COOR1,EDGE1,CONN1,NVFIX,NUMFIX,
     +            ITRACE,NQINT,IWORK,LIWORK,RWORK,LRWORK,IFAIL)
*
*     build the mesh of the 2nd domain
*
*     Initialise boundary mesh inputs:
*     the number of line and of the characteristic points of
*     the boundary mesh
*
      READ (NIN,*) NLINES
*
*     Characteristic points of the boundary geometry
*
      READ (NIN,*) (COORCH(1,J),J=1,NLINES)
      READ (NIN,*) (COORCH(2,J),J=1,NLINES)
*
*     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 60 I = 1, NCOMP
         READ (NIN,*) NLCOMP(I)
*
         READ (NIN,*) (LCOMP(K),K=J,J+ABS(NLCOMP(I))-1)
         J = J + ABS(NLCOMP(I))
   60 CONTINUE
*
      READ (NIN,*) PMESH
*
      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 it 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)
*
      DO 80 K = 1, NELT2
         REFT2(K) = 2
   80 CONTINUE
*
*     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
*
      DO 100 K = 1, NELT3
         REFT3(K) = 3
  100 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 mesh 5
*
      ITRACE = 0
      IFAIL = 0
*
      CALL D06DBF(EPS,NV4,NELT4,NEDGE4,COOR4,EDGE4,CONN4,REFT4,NV3,
     +            NELT3,NEDGE3,COOR3,EDGE3,CONN3,REFT3,NV5,NELT5,NEDGE5,
     +            COOR5,EDGE5,CONN5,REFT5,ITRACE,IWORK,LIWORK,IFAIL)
*
      IF (PMESH.EQ.'N') THEN
         WRITE (NOUT,*) 'The restitched mesh characteristics'
         WRITE (NOUT,99999) 'NV    =', NV5
         WRITE (NOUT,99999) 'NELT  =', NELT5
         WRITE (NOUT,99999) 'NEDGE =', NEDGE5
      ELSE IF (PMESH.EQ.'Y') THEN
*
*     Output the mesh to view it using the NAG Graphics Library
*
         WRITE (NOUT,99998) NV5, NELT5, NEDGE5
         DO 120 I = 1, NV5
            WRITE (NOUT,99997) COOR5(1,I), COOR5(2,I)
  120    CONTINUE
*
         DO 140 K = 1, NELT5
            WRITE (NOUT,99996) CONN5(1,K), CONN5(2,K), CONN5(3,K),
     +        REFT5(K)
  140    CONTINUE
*
         DO 160 K = 1, NEDGE5
            WRITE (NOUT,99998) EDGE5(1,K), EDGE5(2,K), EDGE5(3,K)
  160    CONTINUE
      ELSE
         WRITE (NOUT,*) 'Problem with the printing option Y or N'
         STOP
      END IF
*
99999 FORMAT (1X,A,I6)
99998 FORMAT (1X,3I10)
99997 FORMAT (2(2X,E12.6))
99996 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, Y0
*     .. Executable Statements ..
      FBND = 0.D0
      IF (I.EQ.1) THEN
*
*     inner circle
*
         X0 = 0.D0
         Y0 = 0.D0
         RADIUS2 = 1.D0
         FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2
      ELSE IF (I.EQ.2) THEN
*
*     outer circle
*
         X0 = 0.D0
         Y0 = 0.D0
         RADIUS2 = 5.D0
         FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2
      END IF
*
      RETURN
      END