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

NAG CL Interface Introduction
Example description
/* nag_interp_dim2_triangulate (e01eac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.0, 2023.
 */

#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
}