!   D06CCF Example Program Text
!   Mark 25 Release. NAG Copyright 2014.
    Module d06ccfe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp, x01aaf
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                               :: create_mesh
!     .. Parameters ..
      Integer, Parameter, Public           :: matout = 7, nout = 6, nvb1 = 40, &
                                              nvb2 = 30, nvb3 = 30, nvmax = 260
      Integer, Parameter, Public           :: nvb = nvb1 + nvb2 + nvb3
      Integer, Parameter, Public           :: nedge = nvb
      Logical, Parameter, Public           :: pmesh = .False.

    Contains
      Subroutine create_mesh(edge,coor,nv,nelt,conn)

!       .. Use Statements ..
        Use nag_library, Only: d06aaf
!       .. Scalar Arguments ..
        Integer, Intent (Out)                :: nelt, nv
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: coor(2,nvmax)
        Integer, Intent (Out)                :: conn(3,2*nvmax-2), edge(3,nedge)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: coef, pi2, power, r, theta,    &
                                                theta_i, x0, y0
        Integer                              :: i, ifail, itrace, liwork, lrwork
        Logical                              :: smooth
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable      :: bspace(:), rwork(:)
        Integer, Allocatable                 :: iwork(:)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: cos, max, real, sin
!       .. Executable Statements ..

        lrwork = nvmax
        liwork = 16*nvmax + 2*nedge + max(4*nvmax+2,nedge-14)
        Allocate (bspace(nvb),rwork(lrwork),iwork(liwork))

!       Outer circle
        pi2 = 2.0_nag_wp*x01aaf(theta)
        theta = pi2/real(nvb1,kind=nag_wp)
        r = 1.0_nag_wp
        x0 = 0.0_nag_wp
        y0 = 0.0_nag_wp
        Do i = 1, nvb1
          theta_i = theta*real(i,kind=nag_wp)
          coor(1,i) = x0 + r*cos(theta_i)
          coor(2,i) = y0 + r*sin(theta_i)
        End Do
!       Larger inner circle
        theta = pi2/real(nvb2,kind=nag_wp)
        r = 0.49_nag_wp
        x0 = -0.5_nag_wp
        y0 = 0.0_nag_wp
        Do i = 1, nvb2
          theta_i = theta*real(i,kind=nag_wp)
          coor(1,nvb1+i) = x0 + r*cos(theta_i)
          coor(2,nvb1+i) = y0 + r*sin(theta_i)
        End Do
!       Smaller inner circle
        theta = pi2/real(nvb3,kind=nag_wp)
        r = 0.15_nag_wp
        x0 = -0.5_nag_wp
        y0 = 0.65_nag_wp
        Do i = 1, nvb3
          theta_i = theta*real(i,kind=nag_wp)
          coor(1,nvb1+nvb2+i) = x0 + r*cos(theta_i)
          coor(2,nvb1+nvb2+i) = y0 + r*sin(theta_i)
        End Do

!       Boundary edges
        Do i = 1, nedge
          edge(1,i) = i
          edge(2,i) = i + 1
          edge(3,i) = 1
        End Do
        edge(2,nvb1) = 1
        edge(2,nvb1+nvb2) = nvb1 + 1
        edge(2,nvb) = nvb1 + nvb2 + 1

!       Initialise mesh control parameters
        bspace(1:nvb) = 0.05E0_nag_wp
        smooth = .True.
        itrace = 0
        coef = 0.75E0_nag_wp
        power = 0.25E0_nag_wp

!       Call to the mesh generator
        ifail = 0
        Call d06aaf(nvb,nvmax,nedge,edge,nv,nelt,coor,conn,bspace,smooth,coef, &
          power,itrace,rwork,lrwork,iwork,liwork,ifail)

        Deallocate (bspace,rwork,iwork)

        Return
      End Subroutine create_mesh
    End Module d06ccfe_mod
    Program d06ccfe

!     D06CCF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d06cbf, d06ccf, nag_wp
      Use d06ccfe_mod, Only: create_mesh, matout, nedge, nout, nvmax, pmesh
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Integer                              :: i, ifail, itrace, liwork,        &
                                              lrwork, nelt, nnz, nnzmax, nv
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: coor(:,:), rwork(:)
      Integer, Allocatable                 :: conn(:,:), edge(:,:), icol(:),   &
                                              irow(:), iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: max
!     .. Executable Statements ..
      Write (nout,*) 'D06CCF Example Program Results'
      Flush (nout)

!     Allocate arrays defining mesh
      Allocate (conn(3,2*nvmax-2),edge(3,nedge),coor(2,nvmax))

!     Define boundary mesh and Generate interior mesh
      Call create_mesh(edge,coor,nv,nelt,conn)

      nnzmax = nv**2
      liwork = max(nnzmax,20*nv)
      lrwork = nv
      Allocate (irow(nnzmax),icol(nnzmax),iwork(liwork),rwork(lrwork))

!     Compute the sparsity of the FE matrix
!     from the input geometry

      ifail = 0
      Call d06cbf(nv,nelt,nnzmax,conn,nnz,irow,icol,ifail)

      Write (nout,*)

      If (pmesh) Then

!       Output the sparsity of the mesh
        Write (matout,99998)(irow(i),icol(i),i=1,nnz)

      Else
        Write (nout,*) 'Matrix Sparsity characteristics before renumbering'
        Write (nout,99999) 'nv   =', nv
        Write (nout,99999) 'nnz  =', nnz
        Write (nout,99999) 'nelt =', nelt
      End If
      Flush (nout)

!     Call the renumbering routine and get the new sparsity

      itrace = 1

      ifail = 0
      Call d06ccf(nv,nelt,nedge,nnzmax,nnz,coor,edge,conn,irow,icol,itrace, &
        iwork,liwork,rwork,lrwork,ifail)

      If (pmesh) Then

!       Output the sparsity of the renumbered mesh
        Write (matout,*)
        Write (matout,*)
        Write (matout,99998)(irow(i),icol(i),i=1,nnz)

!       Output the renumbered mesh
        Write (nout,99998) nv, nelt
        Write (nout,99997)(coor(1:2,i),i=1,nv)
        Write (nout,99996)(conn(1:3,i),i=1,nelt)

      Else
        Write (nout,*)
        Write (nout,*) 'Matrix Sparsity characteristics after renumbering'
        Write (nout,99999) 'nv   =', nv
        Write (nout,99999) 'nnz  =', nnz
        Write (nout,99999) 'nelt =', nelt
      End If


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