PROGRAM d06cafe

!      D06CAF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. 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 Functions ..
       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 initialiser 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))

!      Initialise the seed

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

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

       DEALLOCATE (state)
       ALLOCATE (state(lstate))

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

20     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