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

NAG CL Interface Introduction
Example description
/* nag_mesh_dim2_smooth_bary (d06cac) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.5, 2022.
 *
 *
 */

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

#define EDGE(I, J) edge[3 * ((J)-1) + (I)-1]
#define CONN(I, J) conn[3 * ((J)-1) + (I)-1]
#define COOR(I, J) coor[2 * ((J)-1) + (I)-1]

int main(void) {
  /* Integer scalar and array declarations */
  const Integer nvfix = 0;
  Integer exit_status = 0;
  Integer i, imax, imaxm1, ind, itrace, j, jmax, jmaxm1, k, me1, me2, me3,
      nedge, nelt, nqint, nv, reftk, lstate;
  Integer *conn = 0, *edge = 0, *numfix = 0, *state = 0;

  /* Character scalar and array declarations */
  char pmesh[2];

  /* NAG structures */
  NagError fail;

  /* Double scalar and array declarations */
  double delta, hx, hy, pi, dpi, r, rad, sk, theta, x1, x2, x3, y1, y2, y3;
  double one_draw[1];
  double *coor = 0;

  /* Choose the base generator */
  Nag_BaseRNG genid = Nag_Basic;
  Integer subid = 0;

  /* Set the seed */
  Integer seed[] = {1762543};
  Integer lseed = 1;

  /* Initialize the error structure */
  INIT_FAIL(fail);

  /* Get the length of the state array */
  lstate = -1;
  nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf(" nag_mesh_dim2_smooth_bary (d06cac) Example Program Results\n\n");
  fflush(stdout);

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Read imax and jmax, the number of vertices */
  /* in the x and y directions respectively. */
  scanf("%" NAG_IFMT "", &imax);
  scanf("%" NAG_IFMT "", &jmax);
  scanf("%*[^\n] ");

  /* Read distortion percentage and calculate radius */
  /* of distortion neighbourhood so that cross-over */
  /* can only occur at 100% or greater. */
  scanf("%lf", &delta);
  scanf("%*[^\n] ");

  nv = imax * jmax;
  imaxm1 = imax - 1;
  jmaxm1 = jmax - 1;
  nelt = 2 * imaxm1 * jmaxm1;
  nedge = 2 * (imaxm1 + jmaxm1);

  /* Allocate memory */
  if (!(coor = NAG_ALLOC(2 * nv, double)) ||
      !(conn = NAG_ALLOC(3 * nelt, Integer)) ||
      !(edge = NAG_ALLOC(3 * nedge, Integer)) ||
      !(state = NAG_ALLOC(lstate, Integer)) ||
      !(numfix = NAG_ALLOC(1, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  scanf(" ' %1s '", pmesh);
  scanf("%*[^\n] ");

  hx = 1.0 / (double)imaxm1;
  hy = 1.0 / (double)jmaxm1;
  rad = 0.01 * delta * (hx > hy ? hy : hx) / 2.0;
  pi = 4.0 * atan(1.0);

  /* Initialize the generator to a repeatable sequence */
  nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Generate a simple uniform mesh and then distort it */
  /* randomly within the distortion neighbourhood of each */
  /* node. */
  ind = 0;
  for (j = 1; j <= jmax; ++j) {
    for (i = 1; i <= imax; ++i) {
      /* Generate one uniform variate between 0 and rad */
      nag_rand_dist_uniform(1, 0.0, rad, state, one_draw, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_rand_dist_uniform (g05sqc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
      r = one_draw[0];

      dpi = 2.0 * pi;

      /* Generate one uniform variate between 0 and dpi */
      nag_rand_dist_uniform(1, 0.0, dpi, state, one_draw, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_rand_dist_uniform (g05sqc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
      theta = one_draw[0];

      if (i == 1 || i == imax || j == 1 || j == jmax)
        r = 0.0;
      k = (j - 1) * imax + i;
      COOR(1, k) = (i - 1) * hx + r * cos(theta);
      COOR(2, k) = (j - 1) * hy + r * sin(theta);
      if (i < imax && j < jmax) {
        ++ind;
        CONN(1, ind) = k;
        CONN(2, ind) = k + 1;
        CONN(3, ind) = k + imax + 1;
        ++ind;
        CONN(1, ind) = k;
        CONN(2, ind) = k + imax + 1;
        CONN(3, ind) = k + imax;
      }
    }
  }

  if (pmesh[0] == 'N') {
    printf(" The complete distorted mesh characteristics\n");
    printf(" nv   =%6" NAG_IFMT "\n", nv);
    printf(" nelt =%6" NAG_IFMT "\n\n", nelt);
  } else if (pmesh[0] == 'Y') {
    /* Output the mesh to view it using the NAG Graphics Library */

    printf(" %10" NAG_IFMT "%10" NAG_IFMT "\n", nv, nelt);

    for (i = 1; i <= nv; ++i)
      printf("  %15.6e  %15.6e  \n", COOR(1, i), COOR(2, i));
  } else {
    printf("Problem with the printing option Y or N\n");
  }

  reftk = 0;

  for (k = 1; k <= nelt; ++k) {
    me1 = CONN(1, k);
    me2 = CONN(2, k);
    me3 = CONN(3, k);

    x1 = COOR(1, me1);
    x2 = COOR(1, me2);
    x3 = COOR(1, me3);
    y1 = COOR(2, me1);
    y2 = COOR(2, me2);
    y3 = COOR(2, me3);

    sk = 0.5 * ((x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1));
    if (sk < 0.0) {
      printf("Error the surface of the element is negative\n");
      printf(" k = %6" NAG_IFMT "\n", k);
      printf(" sk = %15.6e\n", sk);
      exit_status = -1;
      goto END;
    }

    if (pmesh[0] == 'Y')
      printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "\n",
             CONN(1, k), CONN(2, k), CONN(3, k), reftk);
  }

  /* Boundary edges */

  ind = 0;
  for (i = 1; i <= imaxm1; ++i) {
    ++ind;
    EDGE(1, ind) = i;
    EDGE(2, ind) = i + 1;
    EDGE(3, ind) = 0;
  }

  for (i = 1; i <= jmaxm1; ++i) {
    ++ind;
    EDGE(1, ind) = i * imax;
    EDGE(2, ind) = (i + 1) * imax;
    EDGE(3, ind) = 0;
  }

  for (i = 1; i <= (imax - 1); ++i) {
    ++ind;
    EDGE(1, ind) = imax * jmax - i + 1;
    EDGE(2, ind) = imax * jmax - i;
    EDGE(3, ind) = 0;
  }

  for (i = 1; i <= jmaxm1; ++i) {
    ++ind;
    EDGE(1, ind) = (jmax - i) * imax + 1;
    EDGE(2, ind) = (jmax - i - 1) * imax + 1;
    EDGE(3, ind) = 0;
  }

  itrace = 1;
  nqint = 10;

  /* Call the smoothing routine */

  /* nag_mesh_dim2_smooth_bary (d06cac).
   * Uses a barycentering technique to smooth a given mesh
   */
  fflush(stdout);
  nag_mesh_dim2_smooth_bary(nv, nelt, nedge, coor, edge, conn, nvfix, numfix,
                            itrace, 0, nqint, &fail);
  if (fail.code == NE_NOERROR) {
    if (pmesh[0] == 'N') {
      printf("\n The complete smoothed mesh characteristics\n");
      printf(" nv   =%6" NAG_IFMT "\n", nv);
      printf(" nelt =%6" NAG_IFMT "\n", nelt);
    } else if (pmesh[0] == 'Y') {
      /* Output the mesh to view it using the NAG Graphics Library */

      printf(" %10" NAG_IFMT "%10" NAG_IFMT "\n", nv, nelt);

      for (i = 1; i <= nv; ++i)
        printf("  %15.6e  %15.6e  \n", COOR(1, i), COOR(2, i));

      reftk = 0;

      for (k = 1; k <= nelt; ++k)
        printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT
               "\n",
               CONN(1, k), CONN(2, k), CONN(3, k), reftk);
    } else {
      printf("Problem with the printing option Y or N\n");
      exit_status = -1;
      goto END;
    }
  } else {
    printf("Error from nag_mesh_dim2_smooth_bary (d06cac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

END:
  NAG_FREE(coor);
  NAG_FREE(conn);
  NAG_FREE(edge);
  NAG_FREE(numfix);
  NAG_FREE(state);

  return exit_status;
}