!   D06DBF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2017.
    Module d06dbfe_mod

!     D06DBF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: fbnd
!     .. Parameters ..
      Integer, Parameter, Public       :: meshout = 7, nin = 5, nout = 6,      &
                                          nvb1 = 19
      Logical, Parameter, Public       :: pmesh = .False.
    Contains
      Function fbnd(i,x,y,ruser,iuser)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: fbnd
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x, y
        Integer, Intent (In)           :: i
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: radius2, x0, y0
!       .. Executable Statements ..
        fbnd = 0.0_nag_wp

        Select Case (i)
        Case (1)

!         inner circle

          x0 = 0.0_nag_wp
          y0 = 0.0_nag_wp
          radius2 = 1.0_nag_wp
          fbnd = (x-x0)**2 + (y-y0)**2 - radius2
        Case (2)

!         outer circle

          x0 = 0.0_nag_wp
          y0 = 0.0_nag_wp
          radius2 = 5.0_nag_wp
          fbnd = (x-x0)**2 + (y-y0)**2 - radius2
        End Select

        Return

      End Function fbnd
    End Module d06dbfe_mod
    Program d06dbfe

!     D06DBF Example Main Program

!     .. Use Statements ..
      Use d06dbfe_mod, Only: nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Executable Statements ..
      Write (nout,*) 'D06DBF Example Program Results'

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use d06dbfe_mod, Only: meshout, nout, nvb1, pmesh
        Use nag_library, Only: d06daf, d06dbf, nag_wp, x01aaf
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: ddx, ddy, dx, eps, pi2, pi4, r_i,    &
                                          r_j
        Integer                        :: i, ifail, imax, itrace, itrans, j,   &
                                          jmax, jtrans, k, ktrans, l, liwork,  &
                                          lrwork, nedge1, nedge2, nedge3,      &
                                          nelt1, nelt2, nelt3, ntrans, nv1,    &
                                          nv2, nv3
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: coor1(:,:), coor2(:,:), coor3(:,:), &
                                          rwork(:), trans(:,:)
        Integer, Allocatable           :: conn1(:,:), conn2(:,:), conn3(:,:),  &
                                          edge1(:,:), edge2(:,:), edge3(:,:),  &
                                          itype(:), iwork(:), reft1(:),        &
                                          reft2(:), reft3(:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: real, sin
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*) 'Example 1'
        Write (nout,*)

        imax = nvb1 + 1
        jmax = imax
        nedge1 = 2*(imax-1) + 2*(jmax-1)
        nedge2 = nedge1
        nedge3 = nedge1 + nedge2
        ntrans = 1
        lrwork = 12*ntrans

!       Allocate for mesh : coordinates and connectivity of the 1st domain,
!                           2nd translated domain and restitched domain.

        nv1 = (nvb1+1)**2
        nelt1 = 2*nvb1*nvb1
        nv2 = nv1
        nv3 = nv1 + nv2
        nelt2 = nelt1
        nelt3 = nelt1 + nelt2
        liwork = 2*nv1 + 3*nv2 + nelt1 + nelt2 + nedge1 + nedge2 + 1024
        Allocate (coor1(2,nv1),coor2(2,nv2),coor3(2,nv3),conn1(3,nelt1),       &
          conn2(3,nelt2),conn3(3,nelt3),reft1(nelt1),reft2(nelt2),             &
          reft3(nelt3),edge1(3,nedge1),edge2(3,nedge2),edge3(3,nedge3),        &
          itype(ntrans),trans(6,ntrans),rwork(lrwork),iwork(liwork))

!       Set up interior mesh as small perturbations on regular grid
!       with regularity on the boundary of square on [0,1]x[0,1].
        dx = 1.0_nag_wp/real(nvb1,kind=nag_wp)
        pi2 = x01aaf(dx) + x01aaf(dx)
        pi4 = pi2 + pi2
        k = 0
        Do j = 1, jmax
          ddy = real(j-1,kind=nag_wp)*dx
          Do i = 1, imax
            k = k + 1
            ddx = real(i-1,kind=nag_wp)*dx
            coor1(1,k) = ddx + dx*0.05_nag_wp*sin(pi4*ddx)*sin(pi4*ddy)
            coor1(2,k) = ddy + dx*0.05_nag_wp*sin(pi2*ddx)*sin(pi2*ddy)
          End Do
        End Do

!       Triangulate using skew-diagonals on grid squares
        k = 0
        l = 0
        Do i = 1, nvb1
          Do j = 1, nvb1
            l = l + 1
            k = k + 1
            conn1(1,k) = l
            conn1(2,k) = l + 1
            conn1(3,k) = l + nvb1 + 2
            k = k + 1
            conn1(1,k) = l
            conn1(2,k) = l + nvb1 + 2
            conn1(3,k) = l + nvb1 + 1
          End Do
          l = l + 1
        End Do

        reft1(1:nelt1) = 1
        reft2(1:nelt2) = 2

!       Define the edges of the boundary
        Do i = 1, nvb1
          edge1(1,i) = i
          edge1(2,i) = i + 1
        End Do
        Do i = 1, nvb1
          edge1(1,nvb1+i) = i*imax
          edge1(2,nvb1+i) = (i+1)*imax
        End Do
        Do i = 1, nvb1
          edge1(1,2*nvb1+i) = imax*jmax - i + 1
          edge1(2,2*nvb1+i) = imax*jmax - i
        End Do
        Do i = 1, nvb1
          edge1(1,3*nvb1+i) = (jmax-i)*imax + 1
          edge1(2,3*nvb1+i) = (jmax-i-1)*imax + 1
        End Do
        edge1(3,1:nedge1) = 1

        If (pmesh) Then
!         Print interior mesh of single square
          Do i = 1, nelt1
            Write (meshout,99997) coor1(1,conn1(1,i)), coor1(2,conn1(1,i))
            Write (meshout,99997) coor1(1,conn1(2,i)), coor1(2,conn1(2,i))
            Write (meshout,99997) coor1(1,conn1(3,i)), coor1(2,conn1(3,i))
            Write (meshout,99997) coor1(1,conn1(1,i)), coor1(2,conn1(1,i))
            Write (meshout,*)
          End Do
          Write (meshout,*)
        End If

        Do ktrans = 1, 2

!         Translation of the 1st domain to obtain the 2nd domain
!         KTRANS = 1 leading to a domains (4x2) overlapping
!         KTRANS = 2 leading to a domains partition

          If (ktrans==1) Then
            itrans = nvb1 - 4
            jtrans = nvb1 - 2
          Else
            itrans = nvb1
            jtrans = 0
          End If

          itype(1:ntrans) = (/1/)
          r_i = real(itrans,kind=nag_wp)/real(nvb1,kind=nag_wp)
          r_j = real(jtrans,kind=nag_wp)/real(nvb1,kind=nag_wp)
          trans(1,1:ntrans) = (/r_i/)
          trans(2,1:ntrans) = (/r_j/)
          itrace = 0

          ifail = 0
          Call d06daf(nv2,nedge2,nelt2,ntrans,itype,trans,coor1,edge1,conn1,   &
            coor2,edge2,conn2,itrace,rwork,lrwork,ifail)

          edge2(3,1:nedge2) = 2

!         Call to the restitching driver

          itrace = 0
          eps = 1.E-2_nag_wp

          ifail = 0
          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) Then

!           Output the overlapping or partitioned mesh

            Write (nout,99998) nv3, nelt3, nedge3

            Do i = 1, nelt3
              Write (meshout,99997) coor3(1,conn3(1,i)), coor3(2,conn3(1,i))
              Write (meshout,99997) coor3(1,conn3(2,i)), coor3(2,conn3(2,i))
              Write (meshout,99997) coor3(1,conn3(3,i)), coor3(2,conn3(3,i))
              Write (meshout,99997) coor3(1,conn3(1,i)), coor3(2,conn3(1,i))
              Write (meshout,*)
            End Do
            Write (meshout,*)

            Do k = 1, nelt3
              Write (nout,99996) conn3(1:3,k), reft3(k)
            End Do

            Do k = 1, nedge3
              Write (nout,99998) edge3(1:3,k)
            End Do
          Else

            If (ktrans==1) Then
              Write (nout,*)                                                   &
                'The restitched mesh characteristics in the overlapping case'
            Else
              Write (nout,*)                                                   &
                'The restitched mesh characteristics in the partition case'
            End If

            Write (nout,99999) 'nv    =', nv3
            Write (nout,99999) 'nelt  =', nelt3
            Write (nout,99999) 'nedge =', nedge3
          End If
        End Do

99999   Format (1X,A,I6)
99998   Format (1X,3I10)
99997   Format (2(2X,E13.6))
99996   Format (1X,4I10)
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use d06dbfe_mod, Only: fbnd, meshout, nin, nout, pmesh
        Use nag_library, Only: d06abf, d06acf, d06baf, d06caf, d06daf, d06dbf, &
                               f16dnf, nag_wp
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: eps
        Integer                        :: i, ifail, itrace, j, k, liwork,      &
                                          lrwork, maxind, maxval, ncomp,       &
                                          nedge1, nedge2, nedge3, nedge4,      &
                                          nedge5, nedmx, nelt1, nelt2, nelt3,  &
                                          nelt4, nelt5, nlines, npropa, nqint, &
                                          ntrans, nv1, nv2, nv3, nv4, nv5,     &
                                          nvb1, nvb2, nvfix, nvint, nvmax,     &
                                          sdcrus
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: coor1(:,:), coor2(:,:), coor3(:,:), &
                                          coor4(:,:), coor5(:,:), coorch(:,:), &
                                          crus(:,:), rate(:), rwork(:),        &
                                          trans(:,:), weight(:)
        Real (Kind=nag_wp)             :: ruser(1)
        Integer, Allocatable           :: conn1(:,:), conn2(:,:), conn3(:,:),  &
                                          conn4(:,:), conn5(:,:), edge1(:,:),  &
                                          edge2(:,:), edge3(:,:), edge4(:,:),  &
                                          edge5(:,:), itype(:), iwork(:),      &
                                          lcomp(:), lined(:,:), nlcomp(:),     &
                                          numfix(:), reft1(:), reft2(:),       &
                                          reft3(:), reft4(:), reft5(:)
        Integer                        :: iuser(1)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*) 'Example 2'
        Write (nout,*)

!       Skip headings in data file
        Read (nin,*)
        Read (nin,*)

!       Build and mesh two domains and rotate/translate second to create
!       third.

!       ----------------------------------------------
!       1st domain: Annulus with straight right edge.
!         First two points are end of straight edge.
!         The mesh for inner circle is defined by four NWSE points with mesh
!             points on the quarter-circle between each pair (fbnd, i=1).
!         The mesh for the incomplete outer circle  is defined by three NWS
!             points and the edge points, mesh points computed using fbnd,
!             i=2.

!       The number of lines (1+4+4).
        Read (nin,*) nlines, nvmax, nedmx
        Allocate (coor1(2,nvmax),edge1(3,nedmx),lined(4,nlines),lcomp(nlines), &
          coorch(2,nlines),rate(nlines))

!       Characteristic points of the boundary geometry.
        Read (nin,*) coorch(1,1:nlines)
        Read (nin,*) coorch(2,1:nlines)

!       The Lines of the boundary mesh
        Read (nin,*)(lined(1:4,j),rate(j),j=1,nlines)

!       Allocate workspace for d06baf.

!       sdcrus = 0 in this case.
        sdcrus = 0
        Do i = 1, nlines
          If (lined(4,i)<0) Then
            sdcrus = sdcrus + lined(1,i) - 2
          End If
        End Do
!       Get max(LINED(1,:)) for computing lrwork
        Call f16dnf(nlines,lined,4,maxind,maxval)
        liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus
        lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines
        Allocate (crus(2,sdcrus),rwork(lrwork),iwork(liwork))

!       The number of connected components (outer circle/edge and inner
!       circle)
        Read (nin,*) ncomp
        Allocate (nlcomp(ncomp))

!       Read the lines comprising each connected component
        j = 1
        Do i = 1, ncomp
          Read (nin,*) nlcomp(i)
          k = j + abs(nlcomp(i)) - 1
          Read (nin,*) lcomp(j:k)
          j = k + 1
        End Do

        itrace = 0
!       Call to the 2D boundary mesh generator
        ifail = 0
        Call d06baf(nlines,coorch,lined,fbnd,crus,sdcrus,rate,ncomp,nlcomp,    &
          lcomp,nvmax,nedmx,nvb1,coor1,nedge1,edge1,itrace,ruser,iuser,rwork,  &
          lrwork,iwork,liwork,ifail)

        Deallocate (rwork,iwork)

!       Generate mesh using Delaunay-Voronoi method

!       Initialize mesh control parameters and allocate workspace.
        itrace = 0
        npropa = 1
        nvint = 0
        lrwork = 12*nvmax + 15
        liwork = 6*nedge1 + 32*nvmax + 2*nvb1 + 78
        Allocate (weight(nvint),rwork(lrwork),iwork(liwork),                   &
          conn1(3,2*nvmax+5))

!       Call to the 2D Delaunay-Voronoi mesh generator
        ifail = 0
        Call d06abf(nvb1,nvint,nvmax,nedge1,edge1,nv1,nelt1,coor1,conn1,       &
          weight,npropa,itrace,rwork,lrwork,iwork,liwork,ifail)

        Deallocate (rwork,iwork)

!       Call the smoothing routine
        nvfix = 0
        nqint = 10
        lrwork = 2*nv1 + nelt1
        liwork = 8*nelt1 + 2*nv1
        Allocate (numfix(nvfix),rwork(lrwork),iwork(liwork))
        ifail = 0
        Call d06caf(nv1,nelt1,nedge1,coor1,edge1,conn1,nvfix,numfix,itrace,    &
          nqint,iwork,liwork,rwork,lrwork,ifail)

        Deallocate (rwork,iwork,coorch,lined,lcomp,rate,nlcomp,crus)

!       ----------------------------------------------
!       2nd domain: a rectangle (4 by 2) abutting straight edge of 1st domain.

        Read (nin,*) nlines
        Allocate (lined(4,nlines),lcomp(nlines),coorch(2,nlines),rate(nlines), &
          coor2(2,nvmax),edge2(3,nedmx))

!       Characteristic points of the boundary geometry
        Read (nin,*) coorch(1,1:nlines)
        Read (nin,*) coorch(2,1:nlines)

!       The Lines of the boundary mesh
        Read (nin,*)(lined(1:4,j),rate(j),j=1,nlines)

!       sdcrus = 0 again here.
        sdcrus = 0
        Do i = 1, nlines
          If (lined(4,i)<0) Then
            sdcrus = sdcrus + lined(1,i) - 2
          End If
        End Do

!       Generate initial mesh using Delaunay-Voronoi method

        liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus
        Call f16dnf(nlines,lined,4,maxind,maxval)
        lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines

        Allocate (crus(2,sdcrus),rwork(lrwork),iwork(liwork))

!       The number of connected components (1 rectangle)
        Read (nin,*) ncomp
        Allocate (nlcomp(ncomp))

!       Four lines of rectangle
        j = 1
        Do i = 1, ncomp
          Read (nin,*) nlcomp(i)
          k = j + abs(nlcomp(i)) - 1
          Read (nin,*) lcomp(j:k)
          j = k + 1
        End Do

        itrace = 0
!       Call to the 2D boundary mesh generator
        ifail = 0
        Call d06baf(nlines,coorch,lined,fbnd,crus,sdcrus,rate,ncomp,nlcomp,    &
          lcomp,nvmax,nedmx,nvb2,coor2,nedge2,edge2,itrace,ruser,iuser,rwork,  &
          lrwork,iwork,liwork,ifail)

        Deallocate (rwork,iwork,weight)

!       Remesh 2nd domain using the advancing front method

!       Initialize mesh control parameters
        itrace = 0
        nvint = 0
        lrwork = 12*nvmax + 30015
        liwork = 8*nedge2 + 53*nvmax + 2*nvb2 + 10078
        Allocate (weight(nvint),rwork(lrwork),iwork(liwork),                   &
          conn2(3,2*nvmax+5))

!       Call to the 2D Advancing front mesh generator
        ifail = 0
        Call d06acf(nvb2,nvint,nvmax,nedge2,edge2,nv2,nelt2,coor2,conn2,       &
          weight,itrace,rwork,lrwork,iwork,liwork,ifail)

        Deallocate (rwork,iwork)

!       ----------------------------------------------
!       3rd domain: rotation and translation of the 2nd domain mesh

        ntrans = 1
        lrwork = 12*ntrans
        Allocate (rwork(lrwork),itype(ntrans),trans(6,ntrans),coor3(2,nv2),    &
          edge3(3,nedge2),conn3(3,nelt2))

        itype(1:ntrans) = (/3/)
        trans(1,1:ntrans) = (/6.0_nag_wp/)
        trans(2,1:ntrans) = (/-1.0_nag_wp/)
        trans(3,1:ntrans) = (/-90.0_nag_wp/)
        itrace = 0

        ifail = 0
        Call d06daf(nv2,nedge2,nelt2,ntrans,itype,trans,coor2,edge2,conn2,     &
          coor3,edge3,conn3,itrace,rwork,lrwork,ifail)

        Deallocate (rwork)

!       ----------------------------------------------
!       Combine Meshes

        nv3 = nv2
        nelt3 = nelt2
        nedge3 = nedge2

        nv4 = nv1 + nv2
        nelt4 = nelt1 + nelt2
        nedge4 = nedge1 + nedge2
        liwork = 2*nv1 + 3*nv2 + nelt1 + nelt2 + nedge1 + nedge2 + 1024
        Allocate (iwork(liwork),coor4(2,nv4),edge4(3,nedge4),conn4(3,nelt4),   &
          reft4(nelt4),reft1(nelt1),reft2(nelt2))

!       Restitching of the mesh 1 and 2 to form the mesh 4
        reft1(1:nelt1) = 1
        reft2(1:nelt2) = 2
        eps = 1.E-3_nag_wp
        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)

        Deallocate (iwork)

        nv5 = nv4 + nv3
        nelt5 = nelt4 + nelt3
        nedge5 = nedge4 + nedge3
        liwork = 2*nv4 + 3*nv3 + nelt4 + nelt3 + nedge4 + nedge3 + 1024
        Allocate (iwork(liwork),coor5(2,nv5),edge5(3,nedge5),conn5(3,nelt5),   &
          reft5(nelt5),reft3(nelt3))

!       Restitching of the mesh 3 and 4 to form the mesh 5
        reft3(1:nelt3) = 3
        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) Then
!         Output the mesh
          Do i = 1, nelt5
            Write (meshout,99997) coor5(1,conn5(1,i)), coor5(2,conn5(1,i))
            Write (meshout,99997) coor5(1,conn5(2,i)), coor5(2,conn5(2,i))
            Write (meshout,99997) coor5(1,conn5(3,i)), coor5(2,conn5(3,i))
            Write (meshout,99997) coor5(1,conn5(1,i)), coor5(2,conn5(1,i))
            Write (meshout,*)
          End Do
          Write (meshout,*)

          Write (nout,99998) nv5, nelt5, nedge5
          Do i = 1, nv5
            Write (nout,99997) coor5(1:2,i)
          End Do
          Do k = 1, nelt5
            Write (nout,99996) conn5(1,k), conn5(2,k), conn5(3,k), reft5(k)
          End Do
          Do k = 1, nedge5
            Write (nout,99998) edge5(1:3,k)
          End Do
        Else
          Write (nout,*) 'The restitched mesh characteristics'
          Write (nout,99999) 'nv    =', nv5
          Write (nout,99999) 'nelt  =', nelt5
          Write (nout,99999) 'nedge =', nedge5
        End If

99999   Format (1X,A,I6)
99998   Format (1X,3I10)
99997   Format (2(2X,E13.6))
99996   Format (1X,4I10)
      End Subroutine ex2
    End Program d06dbfe