Program e01eafe
! E01EAF Example Program Text
! Mark 28.3 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