!   D06DBF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.
    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
!     .. Parameters ..
      Integer, Parameter                   :: nin = 5, nout = 6
    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 nag_library, Only: d06daf, d06dbf, nag_wp
        Use d06dbfe_mod, Only: nin
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: eps
        Integer                              :: i, ifail, imax, itrace,        &
                                                itrans, jmax, jtrans, k,       &
                                                ktrans, liwork, lrwork,        &
                                                nedge1, nedge2, nedge3, nelt1, &
                                                nelt2, nelt3, ntrans, nv1,     &
                                                nv2, nv3
        Character (1)                        :: pmesh
!       .. 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
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*) 'Example 1'
        Write (nout,*)

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

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

!       Read the mesh : coordinates and connectivity of the 1st domain

        Read (nin,*) nv1, nelt1
        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))

        Do i = 1, nv1
          Read (nin,*) coor1(1:2,i)
        End Do

        Do k = 1, nelt1
          Read (nin,*) conn1(1:3,k)
        End Do

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

        Read (nin,*) pmesh

!       the Edges of the boundary

        Do i = 1, imax - 1
          edge1(1,i) = i
          edge1(2,i) = i + 1
        End Do

        Do i = 1, jmax - 1
          edge1(1,imax-1+i) = i*imax
          edge1(2,imax-1+i) = (i+1)*imax
        End Do

        Do i = 1, imax - 1
          edge1(1,imax-1+jmax-1+i) = imax*jmax - i + 1
          edge1(2,imax-1+jmax-1+i) = imax*jmax - i
        End Do

        Do i = 1, jmax - 1
          edge1(1,2*(imax-1)+jmax-1+i) = (jmax-i)*imax + 1
          edge1(2,2*(imax-1)+jmax-1+i) = (jmax-i-1)*imax + 1
        End Do

        edge1(3,1:nedge1) = 1

        Do 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==1) Then
            itrans = imax - 5
            jtrans = jmax - 3
          Else
            itrans = imax - 1
            jtrans = 0
          End If

          itype(1:ntrans) = (/1/)
          trans(1,1:ntrans) = (/real(itrans,kind=nag_wp)/real(imax-1,kind= &
            nag_wp)/)
          trans(2,1:ntrans) = (/real(jtrans,kind=nag_wp)/real(jmax-1,kind= &
            nag_wp)/)
          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)

          Select Case (pmesh)
          Case ('N')

            If (ktrans==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
          Case ('Y')

!           Output the mesh

            Write (nout,99998) nv3, nelt3, nedge3

            Do i = 1, nv3
              Write (nout,99997) coor3(1:2,i)
            End Do

            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

          Case Default
            Write (nout,*) 'Problem with the printing option Y or N'
          End Select

        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 nag_library, Only: d06abf, d06acf, d06baf, d06caf, d06daf, d06dbf, &
                               f16dnf, nag_wp
        Use d06dbfe_mod, Only: fbnd, nin
!       .. 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
        Character (1)                        :: pmesh
!       .. 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 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, 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)

        sdcrus = 0

        Do i = 1, nlines

          If (lined(4,i)<0) Then
            sdcrus = sdcrus + lined(1,i) - 2
          End If

        End Do

        liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus

!       Get max(LINED(1,:)) for computing LRWORK

        Call f16dnf(nlines,lined,4,maxind,maxval)

        lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines

!       The number of connected components to the boundary
!       and their informations

        Read (nin,*) ncomp
        Allocate (nlcomp(ncomp),crus(2,sdcrus),rwork(lrwork),iwork(liwork))

        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)

!       mesh it using Delaunay-Voronoi method

!       Initialise mesh control parameters

        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)

!       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
        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

        Do i = 1, nlines

          If (lined(4,i)<0) Then
            sdcrus = sdcrus + lined(1,i) - 2
          End If

        End Do

        liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus

!       Get max(LINED(1,:)) for computing LRWORK

        Call f16dnf(nlines,lined,4,maxind,maxval)

        lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines

!       The number of connected components to the boundary
!       and their informations

        Read (nin,*) ncomp
        Allocate (nlcomp(ncomp),crus(2,sdcrus),rwork(lrwork),iwork(liwork))

        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)

!       mesh it using the advancing front method

!       Initialise 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)

!       Rotation of the 2nd domain mesh to produce
!       the 3rd mesh domain

        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)

        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)

        Read (nin,*) pmesh

        Select Case (pmesh)
        Case ('N')
          Write (nout,*) 'The restitched mesh characteristics'
          Write (nout,99999) 'NV    =', nv5
          Write (nout,99999) 'NELT  =', nelt5
          Write (nout,99999) 'NEDGE =', nedge5
        Case ('Y')

!         Output the mesh

          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

        Case Default
          Write (nout,*) 'Problem with the printing option Y or N'
        End Select

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
