Program d06aafe

!     D06AAF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: d06aaf, nag_wp, x01aaf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: meshout = 7, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: coef, pi2, power, r, theta, theta_i, &
                                          x0, y0
      Integer                          :: i, ifail, itrace, liwork, lrwork,    &
                                          nedge, nelt, nv, nvb, nvb1, nvb2,    &
                                          nvb3, nvmax
      Logical                          :: smooth
      Character (1)                    :: pmesh
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: bspace(:), coor(:,:), rwork(:)
      Integer, Allocatable             :: conn(:,:), edge(:,:), iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: cos, max, real, sin
!     .. Executable Statements ..
      Write (nout,*) 'D06AAF Example Program Results'

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

!     Reading of the geometry
!     Coordinates of the boundary mesh vertices and
!     edges references.

      Read (nin,*) nvb1, nvb2, nvb3, nvmax
      nvb = nvb1 + nvb2 + nvb3
      nedge = nvb

      lrwork = nvmax
      liwork = 16*nvmax + 2*nedge + max(4*nvmax+2,nedge-14)
      Allocate (bspace(nvb),coor(2,nvmax),rwork(lrwork),conn(3,2*(nvmax-       &
        1)),edge(3,nedge),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

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

      Write (nout,*)
      Read (nin,*) pmesh

      Select Case (pmesh)
      Case ('N')
        Write (nout,99999) 'NV   =', nv
        Write (nout,99999) 'NELT =', nelt
      Case ('Y')

!       Output the mesh in a form suitable for printing

        Write (meshout,*) '# D06ABF Example Program Mesh results'
        Do i = 1, nelt
          Write (meshout,99998) coor(1,conn(1,i)), coor(2,conn(1,i))
          Write (meshout,99998) coor(1,conn(2,i)), coor(2,conn(2,i))
          Write (meshout,99998) coor(1,conn(3,i)), coor(2,conn(3,i))
          Write (meshout,99998) coor(1,conn(1,i)), coor(2,conn(1,i))
          Write (meshout,*)
        End Do
        Write (meshout,*)
      Case Default
        Write (nout,*) 'Problem with the printing option Y or N'
      End Select

99999 Format (1X,A,I6)
99998 Format (2(2X,E13.6))
    End Program d06aafe