/* nag_interp_dim2_triangulate (e01eac) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.3, 2022.
*/
#include <nag.h>
#include <stdio.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");
/* Skip heading in data file and read array lengths */
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 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[3], b[3];
/* 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[0] = k;
t[1] = 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[2] = node(i_k);
} else {
t[2] = 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[1] > k && t[2] > k) {
/* new triangle here */
m = m + 1;
tri[3 * m] = t[0];
tri[3 * m + 1] = t[1];
tri[3 * m + 2] = t[2];
/* 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[0] = lastindx(t[0]);
b[1] = lastindx(t[1]);
b[2] = lastindx(t[2]);
b[0] = node(b[0]);
b[1] = node(b[1]);
b[2] = node(b[2]);
i_ts = b[0] + b[1] + b[2];
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[1] = t[2];
}
}
*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[0] = -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[0]) {
/* 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
}