/* nag_mesh_dim2_gen_inc (d06aac) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.3, 2022.
 */

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

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  Integer i1, i, itrace, k, nedge, nelt, nv, nvb, nvb1, nvb2, nvb3, nvmax,
      reftk;
  double coef, power, pi2, r, theta, x0, y0, theta_i;
  /* Arrays */
  char pmesh[2];
  double *bspace = 0, *coor = 0;
  Integer *conn = 0, *edge = 0;
  /* Nag Types */
  Nag_Boolean smooth;
  NagError fail;

  INIT_FAIL(fail);

  exit_status = 0;

  printf("nag_mesh_dim2_gen_inc (d06aac) Example Program Results\n\n");

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

  /* Reading of the geometry sizes */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nvb1,
        &nvb2, &nvb3, &nvmax);

  nvb = nvb1 + nvb2 + nvb3;
  nedge = nvb;

  if (nvb > nvmax) {
    printf("Problem with the array dimensions\n");
    printf(" nvb nvmax %6" NAG_IFMT "%6" NAG_IFMT "\n", nvb, nvmax);
    printf(" Please increase the value of nvmax\n");
    exit_status = -1;
    goto END;
  }

  /* Allocate memory */

  if (!(bspace = NAG_ALLOC(nvb, double)) ||
      !(coor = NAG_ALLOC(2 * nvmax, double)) ||
      !(conn = NAG_ALLOC(3 * (2 * nvmax - 1), Integer)) ||
      !(edge = NAG_ALLOC(3 * nedge, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Coordinates of the boundary mesh vertices and boundary edges */
  i1 = 0;
  pi2 = 2.0 * X01AAC;
  /* Outer circle */
  theta = pi2 / ((double)nvb1);
  r = 1.0;
  x0 = 0.0;
  y0 = 0.0;
  theta_i = 0.0;
  for (i = 0; i < nvb1; ++i) {
    theta_i = theta_i + theta;
    coor[i1] = x0 + r * cos(theta_i);
    coor[i1 + 1] = y0 + r * sin(theta_i);
    i1 = i1 + 2;
  }
  /* Larger inner circle */
  theta = pi2 / ((double)nvb2);
  r = 0.49;
  x0 = -0.5;
  y0 = 0.0;
  theta_i = 0.0;
  for (i = 0; i < nvb2; ++i) {
    theta_i = theta_i + theta;
    coor[i1] = x0 + r * cos(theta_i);
    coor[i1 + 1] = y0 + r * sin(theta_i);
    i1 = i1 + 2;
  }
  /* Smaller inner circle */
  theta = pi2 / ((double)nvb3);
  r = 0.15;
  x0 = -0.5;
  y0 = 0.65;
  theta_i = 0.0;
  for (i = 0; i < nvb3; ++i) {
    theta_i = theta_i + theta;
    coor[i1] = x0 + r * cos(theta_i);
    coor[i1 + 1] = y0 + r * sin(theta_i);
    i1 = i1 + 2;
  }

  /* Boundary edges */
  i1 = 0;
  for (i = 0; i < nedge; ++i) {
    edge[i1] = i + 1;
    edge[i1 + 1] = i + 2;
    edge[i1 + 2] = 1;
    i1 = i1 + 3;
  }
  /* Joins */
  edge[3 * nvb1 - 2] = 1;
  edge[3 * nvb1 + 3 * nvb2 - 2] = nvb1 + 1;
  edge[3 * nvb - 2] = nvb1 + nvb2 + 1;
  scanf(" ' %1s '%*[^\n]", pmesh);

  /* Initialize mesh control parameters */

  for (i = 0; i < nvb; ++i)
    bspace[i] = 0.05;

  smooth = Nag_TRUE;
  itrace = 0;
  coef = 0.75;
  power = 0.25;

  /* Call to the mesh generator */

  /* nag_mesh_dim2_gen_inc (d06aac).
   * Generates a two-dimensional mesh using a simple
   * incremental method
   */
  nag_mesh_dim2_gen_inc(nvb, nvmax, nedge, edge, &nv, &nelt, coor, conn, bspace,
                        smooth, coef, power, itrace, 0, &fail);

  if (fail.code == NE_NOERROR) {
    if (pmesh[0] == '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 = 0; i < nv; ++i)
        printf("  %15.6e  %15.6e\n", coor[2 * i], coor[2 * i + 1]);

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

END:
  NAG_FREE(bspace);
  NAG_FREE(coor);
  NAG_FREE(conn);
  NAG_FREE(edge);

  return exit_status;
}