Program d06cafe

!     D06CAF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: d06caf, g05kff, g05sqf, nag_wp, x01aaf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: delta, hx, hy, pi2, r, rad, sk,      &
                                          theta, x1, x2, x3, y1, y2, y3
      Integer                          :: genid, i, ifail, imax, ind, itrace,  &
                                          j, jmax, k, liwork, lrwork, lseed,   &
                                          lstate, me1, me2, me3, nedge, nelt,  &
                                          nqint, nv, nvfix, reftk, subid
      Character (1)                    :: pmesh
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: coor(:,:), rwork(:), variates(:)
      Integer, Allocatable             :: conn(:,:), edge(:,:), iwork(:),      &
                                          numfix(:), seed(:), state(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: cos, min, real, sin
!     .. Executable Statements ..
      Write (nout,*) 'D06CAF Example Program Results'
      Flush (nout)

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

!     Read IMAX and JMAX, the number of vertices
!     in the x and y directions respectively.

      Read (nin,*) imax, jmax
      nv = imax*jmax
      nelt = 2*(imax-1)*(jmax-1)
      nedge = 2*(imax-1) + 2*(jmax-1)
      liwork = 8*nelt + 2*nv
      lrwork = 2*nv + nelt

!     The array VARIATES will be used when distorting the mesh

      Allocate (variates(2*nv),coor(2,nv),conn(3,nelt),edge(3,nedge),          &
        iwork(liwork),rwork(lrwork))

!     Read distortion percentage and calculate radius
!     of distortion neighbourhood so that cross-over
!     can only occur at 100% or greater.

      Read (nin,*) delta

      hx = 1.0E0_nag_wp/real(imax-1,kind=nag_wp)
      hy = 1.0E0_nag_wp/real(jmax-1,kind=nag_wp)
      rad = 0.005E0_nag_wp*delta*min(hx,hy)
      pi2 = 2.0E0_nag_wp*x01aaf(pi2)

!     GENID identifies the base generator

      genid = 1
      subid = 1

!     For GENID = 1 only one seed is required
!     The initializer is first called in query mode to get the value of
!     LSTATE for the chosen base generator

      lseed = 1
      lstate = -1
      Allocate (seed(lseed),state(lstate))

!     Initialize the seed

      seed(1:lseed) = (/1762541/)

      ifail = 0
      Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

      Deallocate (state)
      Allocate (state(lstate))

!     Initialize the generator to a repeatable sequence

      ifail = 0
      Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!     Generate two sets of uniform random variates

      ifail = 0
      Call g05sqf(nv,0.0E0_nag_wp,rad,state,variates,ifail)

      ifail = 0
      Call g05sqf(nv,0.0E0_nag_wp,pi2,state,variates(nv+1),ifail)

!     Generate a simple uniform mesh and then distort it
!     randomly within the distortion neighbourhood of each
!     node.

      k = 0
      ind = 0

      Do j = 1, jmax

        Do i = 1, imax
          k = k + 1
          r = variates(k)
          theta = variates(nv+k)

          If (i==1 .Or. i==imax .Or. j==1 .Or. j==jmax) Then
            r = 0.E0_nag_wp
          End If

          coor(1,k) = real(i-1,kind=nag_wp)*hx + r*cos(theta)
          coor(2,k) = real(j-1,kind=nag_wp)*hy + r*sin(theta)

          If (i<imax .And. j<jmax) Then
            ind = ind + 1
            conn(1,ind) = k
            conn(2,ind) = k + 1
            conn(3,ind) = k + imax + 1
            ind = ind + 1
            conn(1,ind) = k
            conn(2,ind) = k + imax + 1
            conn(3,ind) = k + imax
          End If

        End Do

      End Do

      Read (nin,*) pmesh

      Write (nout,*)

      Select Case (pmesh)
      Case ('N')
        Write (nout,*) 'The complete distorted mesh characteristics'
        Write (nout,99999) 'Number of vertices =', nv
        Write (nout,99999) 'Number of elements =', nelt
      Case ('Y')

!       Output the mesh

        Write (nout,99998) nv, nelt

        Do i = 1, nv
          Write (nout,99997) coor(1,i), coor(2,i)
        End Do

      Case Default
        Write (nout,*) 'Problem: the printing option must be Y or N'
        Go To 100
      End Select

      reftk = 0

      Do k = 1, nelt
        me1 = conn(1,k)
        me2 = conn(2,k)
        me3 = conn(3,k)

        x1 = coor(1,me1)
        x2 = coor(1,me2)
        x3 = coor(1,me3)
        y1 = coor(2,me1)
        y2 = coor(2,me2)
        y3 = coor(2,me3)

        sk = ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))/2.E0_nag_wp

        If (sk<0.E0_nag_wp) Then
          Write (nout,*) 'Error: the surface of the element is negative'
          Write (nout,99999) 'element number  = ', k
          Write (nout,99995) 'element surface = ', sk
          Go To 100
        End If

        If (pmesh=='Y') Then
          Write (nout,99996) conn(1,k), conn(2,k), conn(3,k), reftk
        End If

      End Do
      Flush (nout)

!     Boundary edges

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

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

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

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

      edge(3,1:nedge) = 0

      nvfix = 0
      Allocate (numfix(nvfix))

      itrace = 1
      nqint = 10

!     Call the smoothing routine

      ifail = 0
      Call d06caf(nv,nelt,nedge,coor,edge,conn,nvfix,numfix,itrace,nqint,      &
        iwork,liwork,rwork,lrwork,ifail)

      Select Case (pmesh)
      Case ('N')
        Write (nout,*)
        Write (nout,*) 'The complete smoothed mesh characteristics'
        Write (nout,99999) 'Number of vertices =', nv
        Write (nout,99999) 'Number of elements =', nelt
      Case ('Y')

!       Output the mesh

        Write (nout,99998) nv, nelt

        Do i = 1, nv
          Write (nout,99997) coor(1,i), coor(2,i)
        End Do

        reftk = 0

        Do k = 1, nelt
          Write (nout,99996) conn(1,k), conn(2,k), conn(3,k), reftk
        End Do

      End Select

100   Continue

99999 Format (1X,A,I6)
99998 Format (1X,2I10)
99997 Format (2(2X,E13.6))
99996 Format (1X,4I10)
99995 Format (1X,A,E13.6)
    End Program d06cafe