```/* nag_interp_dim2_triangulate (e01eac) Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
*
* Mark 27.0, 2019.
*/

#include <stdio.h>
#include <nag.h>

#ifdef __cplusplus
extern "C"
{
#endif
static void NAG_CALL  print_triang(Integer, const double[],
const double[], const Integer[]);
static void NAG_CALL  triang2list(Integer, const Integer[], Integer *,
Integer[], Integer *, Integer[]);
static void NAG_CALL  convex_hull(Integer, const Integer[], Integer[],
Integer *);
#ifdef __cplusplus
}
#endif

int main(void)
{

/* Scalars */
Integer pr_tr = 0;
Integer exit_status, i, j, m, n, nt, nb, nbn;

/* Arrays */
double *f = 0, *pf = 0, *px = 0, *py = 0, *x = 0, *y = 0;
Integer *triang = 0, *tri = 0, *bnd = 0, *bnodes = 0;

/* Nag Types */
NagError fail;

exit_status = 0;
INIT_FAIL(fail);

printf("nag_interp_dim2_triangulate (e01eac) Example Program Results\n\n");

scanf("%*[^\n]");
scanf("%" NAG_IFMT "%*[^\n]", &n);
scanf("%" NAG_IFMT "%*[^\n]", &m);

if (!(x = NAG_ALLOC(n, double)) ||
!(y = NAG_ALLOC(n, double)) ||
!(f = NAG_ALLOC(n, double)) ||
!(triang = NAG_ALLOC(7 * n, Integer)) ||
!(tri = NAG_ALLOC(6 * n - 15, Integer)) ||
!(bnd = NAG_ALLOC(2 * n - 5, Integer)) ||
!(bnodes = NAG_ALLOC(n, Integer)) ||
!(px = NAG_ALLOC(m, double)) ||
!(py = NAG_ALLOC(m, double)) ||
!(pf = NAG_ALLOC(m, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}

/* Read scattered 2d data points and function values. */
for (i = 0; i < n; i++) {
scanf("%lf%lf%lf%*[^\n]", &x[i], &y[i], &f[i]);
}

/* Obtain triangulation of scattered points (x,y) using
* nag_interp_dim2_triangulate (e01eac).
*/
nag_interp_dim2_triangulate(n, x, y, triang, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_interp_dim2_triangulate (e01eac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

/* Convert array triang to an integer array tri[]
* such that tri[3*(k-1)] to tri[3*(k-1)+2] hold node indices for triangle k.
* nt = number of triangles. nb = number of triangles with boundary edge.
* If bnd[i] > 0 Then triangle i+1 is a triangle with boundary edge;
*    bnd[i] = 4      triangle i+1 has two boundary edges,
*    bnd[i] = 1,2,3  node tri[3*i+bnd[i]-1] is an internal vertex.
*/
triang2list(n,triang,&nt,tri,&nb,bnd);

printf("Triangle Information\n");
printf(" %-43s = %3" NAG_IFMT "\n","Number of triangles", nt);
printf(" %-43s = %3" NAG_IFMT "\n","Number of triangles with boundary edges",
nb);

/* Also find the boundary nodes in anti-clockwise order using triang
* That is, find the the convex hull.
*/
convex_hull(n,triang,bnodes,&nbn);

/* Print boundary nodes */
printf("\nBoundary Node Information\n");
printf(" %-43s = %3" NAG_IFMT "\n","Number of boundary nodes", nbn);
printf("\n  node     Coordinates\n");
for (i = 0; i<nbn; ++i) {
j = bnodes[i];
printf("%5" NAG_IFMT "   (%7.4f, %7.4f)\n", j+1, x[j], y[j]);
}
printf("\n");

/* Read points at which interpolated values required. */
for (i = 0; i < m; i++) {
scanf("%lf%lf%*[^\n]", &px[i], &py[i]);
}
/* Use triangulation to perform barycentric interpolation on to the
* set of m points (px,py) using nag_interp_dim2_triang_bary_eval (e01ebc).
*/
nag_interp_dim2_triang_bary_eval(m, n, x, y, f, triang, px, py, pf, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_interp_dim2_triang_bary_eval (e01ebc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* Display results */
printf("Interpolation Results\n");
printf("    %4s   %7s   %19s\n", "px", "py", "Interpolated Value");
for (i = 0; i < m; i++) {
printf("   %7.4f   %7.4f       %7.4f\n", px[i], py[i], pf[i]);
}
if (pr_tr!=0)
print_triang(n,x,y,triang);

END:
NAG_FREE(f);
NAG_FREE(pf);
NAG_FREE(px);
NAG_FREE(py);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(triang);
NAG_FREE(tri);
NAG_FREE(bnd);
NAG_FREE(bnodes);

return exit_status;
}
static void NAG_CALL  print_triang(Integer n, const double x[],
const double y[], const Integer triang[])
{
Integer i_k, j, j_k, k, k2;

printf("\n  Triangulation as a set of line segments\n\n");
j_k = 0;
for (k = 0; k < n; ++k) {
i_k = j_k;
j_k = triang[6*n+k];
for (j = i_k; j < j_k; ++j) {
k2 = triang[j] - 1;
if (k2 > k) {
printf("%7.4f %7.4f\n", x[k], y[k]);
printf("%7.4f %7.4f\n\n", x[k2], y[k2]);
}
}
}
}
static void NAG_CALL  triang2list(Integer n, const Integer triang[],
Integer *nt,  Integer tri[], Integer *nb,
Integer bnd[])
{
Integer    i, i_k, i_ts, j, jj, j_k, k, m;
Integer    t, b;

/* m
*    Initial setting:     0
*    During calculation:  current triangle number
*    On final exit:       total number of unique triangles
*/
m = -1;

/* Define lastindx to be index of last connected node for node k */
#define lastindx(K) triang[6*n + K - 1]
/* Define node to be node number under index I */
#define node(I) triang[I - 1]

j_k = 0;
*nb = 0;
for (k = 1; k <= n; ++k) {
/* 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 = lastindx(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 = k;
t = node(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 (node(j_k) > 0) {
/* 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;
}

/* Now loop over the remaining connecting nodes */
for (j = i_k+1; j <= jj; ++j) {
if (j == j_k + 1) {
/* final triangle; last node is i_k */
t = node(i_k);
} else {
t = node(j);
}
/* 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>k && t>k) {
/* new triangle here */
m = m + 1;
tri[3*m] = t;
tri[3*m+1] = t;
tri[3*m+2] = t;

/* 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 = lastindx(t);
b = lastindx(t);
b = lastindx(t);
b = node(b);
b = node(b);
b = node(b);
i_ts = b + b + b;

if (i_ts == 0) {
/* triangle lies on corner */
bnd[m] = 4;
*nb = *nb + 1;
} else {
for (i = 1; i <= 3; ++i) {
if (b[i-1] == i_ts) {
/* Triangle edge is on boundary
*  Set bnd(m) to index of vertex that is internal
*/
bnd[m] = i;
*nb = *nb + 1;
}
}
}
}

/* 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 = t;
}
}
*nt = m + 1;
#undef lastindx
#undef node
}

static void NAG_CALL  convex_hull(Integer n, const Integer triang[],
Integer bnodes[], Integer *nbn)
{
Integer    i_k, j, jj, j_k, k_j, k, m, found;

/* 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
*/

/* Define lastindx to be index of last connected node for node k */
#define lastindx(K) triang[6*n + K - 1]
/* Define node to be node under index I */
#define node(I) triang[I - 1]

/* 1. Find first boundary node (index lastnode is zero) */
k = 1;
j_k = lastindx(k);
while (node(j_k) > 0 && k < n) {
/* Node k is internal, move to next. */
k = k + 1;
j_k = lastindx(k);
}

/* Either node k is boundary or something has gone wrong */
if (k == n) {
bnodes = -1;
printf("\n could not find first boundary point\n");
return;
}

/* We have first boundary node k */
m = 1;
bnodes[m-1] = k-1;

/* Loop until we come back round to node k */
found = 1;
while (found==1) {
/* 2. Which nodes are connected to current node? */
if (k==1) {
i_k = 1;
} else {
i_k = lastindx(k-1) + 1;
}
j_k = lastindx(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;
j = -2;
for (jj = i_k; jj <= j_k && j != k; ++jj) {
j = node(jj);
k_j = lastindx(j);
if (node(k_j) == 0) {
k = j;
}
}

/* Check that boundary node found */
if (k == 0) {
/* could not find bnode to connect to mode k bnodes(m)
* set flag in bnodes(m+1) and exit
*/
bnodes[m+1] = -1;
found = 0;
} else if (k-1 == bnodes) {
/* We have gone right round the boundary now, so exit */
found = 0;
} else {
/* Update new boundary node and cycle */
m = m + 1;
bnodes[m-1] = k-1;
}
}
*nbn = m;
#undef lastindx
#undef node
}
```