NAG Library Manual, Mark 28.5
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Program e01eafe

!     E01EAF Example Program Text
!     Mark 28.5 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: e01eaf, e01ebf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
      Logical, Parameter               :: pr_tr = .False.
!     .. Local Scalars ..
      Integer                          :: i, ifail, j, m, n, nb, nbn, nt
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: f(:), pf(:), px(:), py(:), x(:),     &
                                          y(:)
      Integer, Allocatable             :: bnd(:), bnodes(:), tri(:,:),         &
                                          triang(:)
!     .. Executable Statements ..

      Write (nout,*) 'E01EAF Example Program Results'

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

      Read (nin,*) n
      Allocate (x(n),y(n),f(n),triang(7*n),tri(3,2*n-5),bnd(2*n-5))
      Read (nin,*)(x(i),y(i),f(i),i=1,n)

!     Triangulate data
      ifail = 0
      Call e01eaf(n,x,y,triang,ifail)

!     Convert array triang to an integer array tri(3,2*n-5)
!     such that tri(1:3,k) holds the node indices for the k-th triangle.
!     nt = number of triangles. nb = number of triangles with boundary edge.
!     If bnd(i) > 0 Then tri(:,i) is a triangle with boundary edge;
!       bnd(i) = 4     means that tri(:,i) has two boundary edges,
!       bnd(i) = 1,2,3 means that tri(bnd(i),i) is internal vertex.

      Call triang2list(n,triang,nt,tri,nb,bnd)

      Write (nout,*)
      Write (nout,*) ' Triangle Information'
      Write (nout,99999) '  Number of triangles                     = ', nt
      Write (nout,99999) '  Number of triangles with boundary edges = ', nb
99999 Format (1X,A,I3)

!     Also find the boundary nodes in anti-clockwise order using triang
!     That is, find the convex hull.
      Allocate (bnodes(n))
      Call convex_hull(n,triang,bnodes,nbn)

!     Print boundary nodes
      Write (nout,*)
      Write (nout,*) ' Boundary Node Information'
      Write (nout,99999) '  Number of boundary nodes                = ', nbn
      Write (nout,*)
      Write (nout,*) '  node     Coordinates'
      Do i = 1, nbn
        j = bnodes(i)
        Write (nout,99998) j, x(j), y(j)
      End Do
99998 Format (1X,I5,3X,'(',F7.4,', ',F7.4,')')

      Read (nin,*) m
      Allocate (px(m),py(m),pf(m))
      Read (nin,*)(px(i),py(i),i=1,m)

!     Interpolate data
      ifail = 0
      Call e01ebf(m,n,x,y,f,triang,px,py,pf,ifail)

!     Display results
      Write (nout,*)
      Write (nout,*) ' Interpolation Results'
      Write (nout,99997) 'px', 'py', 'Interpolated Value'
      Write (nout,99996)(px(i),py(i),pf(i),i=1,m)

      If (pr_tr) Then
        Call print_triang
      End If

99997 Format (2X,A4,4X,A4,4X,A19)
99996 Format (1X,F7.4,1X,F7.4,8X,F7.4)
    Contains
      Subroutine print_triang

!       .. Implicit None Statement ..
        Implicit None
!       .. Local Scalars ..
        Integer                        :: i_k, j, j_k, k
!       .. Executable Statements ..

!       Print a sequence of unique line segments for plotting triangulation
        Write (nout,*)
        Write (nout,*) '  Triangulation as a set of line segments'
        Write (nout,*)
        j_k = 0
        Do k = 1, n
          i_k = j_k + 1
          j_k = triang(6*n+k)
          Do j = i_k, j_k
            If (triang(j)>k) Then
              Write (nout,99999) x(k), y(k)
              Write (nout,99999) x(triang(j)), y(triang(j))
              Write (nout,*)
            End If
          End Do
        End Do
        Return
99999   Format (1X,F7.4,1X,F7.4)
      End Subroutine print_triang
      Subroutine triang2list(n,triang,nt,tri,nb,bnd)

!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: n
        Integer, Intent (Out)          :: nb, nt
!       .. Array Arguments ..
        Integer, Intent (Out)          :: bnd(2*n-5), tri(3,2*n-5)
        Integer, Intent (In)           :: triang(7*n)
!       .. Local Scalars ..
        Integer                        :: i, i_k, i_ts, j, jj, j_k, k, m
!       .. Local Arrays ..
        Integer                        :: b(3), t(3)
!       .. Executable Statements ..

!       m
!          Initial setting:     0
!          During calculation:  current triangle number
!          On final exit:       total number of unique triangles
        m = 0
        j_k = 0
        nb = 0
        Do k = 1, n
!         Which nodes are connected to node k?

!         Node k is connected to nodes i_k to j_k
          i_k = j_k + 1
          j_k = triang(6*n+k)

!         Let n_k (= j_k - i_k + 1) be the number of nodes

!         first connection setup, first two vertices of first triangle:
!                                 node k and node i_k
          t(1) = k
          t(2) = triang(i_k)

!         For the remaining connected nodes, we need to know whether node k
!         is internal or on the boundary.
!         An internal node has n_k associated triangles;
!         A boundary node has n_k-1 associated triangles.
          If (triang(j_k)>0) Then
!           Node k is internal
!           Connected node order is anti-clockwise;
!           so for the last triangle connected to node k has nodes
!           k, j_k and i_k.
!           We need to loop round final connected point to first point
            jj = j_k + 1
          Else
!           Node k is a boundary points;
!           There are n_k-1 triangles and no need to look at zero marker.
            jj = j_k - 1
          End If

!         Now loop over the remaining connecting nodes
          Do j = i_k + 1, jj
            If (j==j_k+1) Then
!             final triangle; last node is i_k
              t(3) = triang(i_k)
            Else
              t(3) = triang(j)
            End If
!           Each triangle will be visited a number of times, but since we
!           are going through in ascending node number k, a new triangle
!           is identified by having connecting node numbers greater than k.
            If (t(2)>k .And. t(3)>k) Then
!             new triangle here
              m = m + 1
              tri(1:3,m) = t(1:3)

!             First assume that no edge of triangle is on boundary
              bnd(m) = 0

!             If two if the following terms are zero then
!                the triangle has an edge on the boundary
              b(1:3) = triang(6*n+t(1:3))
              b(1:3) = triang(b(1:3))
              i_ts = b(1) + b(2) + b(3)
              If (i_ts==0) Then
!               triangle lies on corner
                bnd(m) = 4
                nb = nb + 1
              Else
                Do i = 1, 3
                  If (b(i)==i_ts) Then
!                   Triangle edge is on boundary
!                   Set bnd(m) to index of vertex that is internal
                    bnd(m) = i
                    nb = nb + 1
                  End If
                End Do
              End If
            End If
!           shuffle down (round, anti-clockwise) for next connected node
!           The last point in current triangle connected to node k becomes
!           the second point in next triangle.
            t(2) = t(3)
          End Do
        End Do
        nt = m
        Return
      End Subroutine triang2list
      Subroutine convex_hull(n,triang,bnodes,nbn)

!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: n
        Integer, Intent (Out)          :: nbn
!       .. Array Arguments ..
        Integer, Intent (Out)          :: bnodes(n)
        Integer, Intent (In)           :: triang(7*n)
!       .. Local Scalars ..
        Integer                        :: i_k, j, jj, j_k, k, k_j, m
        Logical                        :: found
!       .. Executable Statements ..

!       Algorithm:
!          1. Find first node in sequence that is boundary node: m = 1
!          2. Get list of nodes connecting to current node
!          3. Find first connecting node that is boundary node and
!             make this the current node: m = m + 1
!          4. Repeat steps 2. and 3. until current mode is first node


!       1. Find first boundary node
        k = 1
        j_k = triang(6*n+k)
        Do While (triang(j_k)>0 .And. k<n)
!         Node k is internal, move to next.
          k = k + 1
          j_k = triang(6*n+k)
        End Do

!       Either node k is boundary or something has gone wrong
        If (k==n) Then
          bnodes(1:nt) = -1
          Return
        End If

!       We have first boundary node k
        m = 1
        bnodes(m) = k

!       Loop until we come back round to node k
        found = .True.
        Do While (found)

!         2. Which nodes are connected to current node?

          If (k==1) Then
            i_k = 1
          Else
            i_k = triang(6*n+k-1) + 1
          End If
          j_k = triang(6*n+k) - 1
!         Node k is connected to nodes triang(i_k) to triang(j_k)

!         3. Find next connecting node that is boundary node

          k = 0
          Do jj = i_k, j_k
            j = triang(jj)
            k_j = triang(6*n+j)
            If (triang(k_j)==0) Then
              k = j
              Exit
            End If
          End Do

!         Check that boundary node found
          If (k==0) Then
!           could not find bnode to connect to mode k bnodes(m)
!           set flag in bnodes(m+1) and exit
            bnodes(m+1) = -1
            found = .False.
          Else If (k==bnodes(1)) Then
!           We have gone right round the boundary now, so exit
            found = .False.
          Else
!           Update new boundary node and cycle
            m = m + 1
            bnodes(m) = k
          End If
        End Do
        nbn = m
        Return
      End Subroutine convex_hull

    End Program e01eafe