/* nag_mesh_dim2_gen_delaunay (d06abc) Example Program.
 *
 * Copyright 2021 Numerical Algorithms Group.
 *
 * Mark 27.2, 2021.
 */

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

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  double c, d, dnvint, r, s, theta, theta_i, tmp;
  Integer i, i1, ic, itrace, j, k, nedge, nelt, npropa, nrae;
  Integer nv, nvb, nvint, nvmax, reftk, nearest, nv_near, nelt_near;
  /* Arrays */
  double t[2];
  char pmesh[2];
  double *coor = 0, *weight = 0;
  Integer *conn = 0, *edge = 0;
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);

  exit_status = 0;

  printf("nag_mesh_dim2_gen_delaunay (d06abc) Example Program Results\n\n");

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

  /* Reading of the geometry */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nrae, &nvint,
        &nvmax);
  scanf(" ' %1s '%*[^\n]", pmesh);

  nvb = nvint + 2 * nrae;
  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 (!(coor = NAG_ALLOC(2 * nvmax, double)) ||
      !(weight = NAG_ALLOC(nvint, double)) ||
      !(conn = NAG_ALLOC(3 * (2 * nvmax + 5), 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 */

  /* Circular outer boundary, radius 3 and centre (1,0) */
  theta = 2.0 * X01AAC / ((double)nvint);
  r = 3.0;
  t[0] = 1.0;
  theta_i = 0.0;
  for (i = 0; i < nvint; ++i) {
    coor[2 * i] = r * cos(theta_i) + t[0];
    coor[2 * i + 1] = r * sin(theta_i);
    theta_i = theta_i + theta;
  }

  /* Read data for aerofoil RAE 2822 */
  for (i = nvint; i < nvint + nrae; ++i) {
    scanf("%" NAG_IFMT "%lf%lf%*[^\n] ", &i1, &coor[2 * i], &coor[2 * i + 1]);
  }

  /* Transform RAE 2822 for secondary foil */
  theta = X01AAC / 12.0;
  c = cos(theta);
  s = sin(theta);
  ic = 2 * (nvint + nrae);
  i1 = 2 * nvint;
  /* Copy and rotate coordinates by theta = pi/12 */
  for (i = 0; i < 2 * nrae; ++i) {
    coor[ic + i] = coor[i1 + i];
  }
  i1 = ic;
  for (i = 0; i < nrae; ++i) {
    tmp = coor[i1];
    coor[i1] = s * coor[i1 + 1] + c * tmp;
    coor[i1 + 1] = c * coor[i1 + 1] - s * tmp;
    i1 = i1 + 2;
  }

  /* Reduce by 0.4 and translate to distance 0.25 from intercept at (0.75,0) */
  d = 0.4;
  t[0] = 0.75 + 0.25 * c;
  t[1] = -0.25 * s;
  for (i = 0; i < nrae; ++i) {
    coor[ic + 2 * i] = d * coor[ic + 2 * i] + t[0];
    coor[ic + 2 * i + 1] = d * coor[ic + 2 * i + 1] + t[1];
  }

  /*  Boundary edges */
  i1 = 0;
  for (i = 0; i < nedge; ++i) {
    edge[i1] = i + 1;
    edge[i1 + 1] = i + 2;
    edge[i1 + 2] = 0;
    i1 = i1 + 3;
  }
  /* Tie up end of three boundary edges */
  edge[3 * nvint - 2] = 1;
  edge[3 * (nvint + nrae) - 2] = nvint + 1;
  edge[3 * nedge - 2] = nvint + nrae + 1;

  /* Initialize mesh control parameters */
  itrace = 0;

  /* Generation of interior vertices on the RAE airfoils wake */
  dnvint = 2.5 / (double)(nvint + 1);

  i1 = 2 * nvb;
  for (i = 0; i < nvint; ++i) {
    coor[i1] = (double)(i + 1) * dnvint + 1.38;
    coor[i1 + 1] = -0.27 * coor[i1] + 0.2;
    i1 = i1 + 2;
    weight[i] = 0.01;
  }

  /* Loop on the propagation coef */
  nearest = 250;
  for (j = 0; j < 4; ++j) {
    switch (j) {
    case 0:
      npropa = -5;
      break;
    case 1:
      npropa = -1;
      break;
    case 2:
      npropa = 1;
      break;
    default:
      npropa = 5;
    }

    /* Call to the 2D Delaunay-Voronoi mesh generator */

    /* nag_mesh_dim2_gen_delaunay (d06abc).
     * Generates a two-dimensional mesh using a Delaunay-Voronoi
     * process
     */
    nag_mesh_dim2_gen_delaunay(nvb, nvint, nvmax, nedge, edge, &nv, &nelt, coor,
                               conn, weight, npropa, itrace, 0, &fail);

    if (fail.code == NE_NOERROR) {
      if (pmesh[0] == 'N') {
        printf(" Mesh characteristics with npropa =%6" NAG_IFMT "\n", npropa);
        nv_near = ((nv + nearest / 2) / nearest) * nearest;
        nelt_near = ((nelt + nearest / 2) / nearest) * nearest;
        printf(" nv   =%10" NAG_IFMT " to the nearest %3" NAG_IFMT "\n",
               nv_near, nearest);
        printf(" nelt =%10" NAG_IFMT " to the nearest %3" NAG_IFMT "\n",
               nelt_near, nearest);
      } else if (pmesh[0] == 'Y') {
        /* Output the mesh in a form suitable for printing */

        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_delaunay (d06abc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }

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

  return exit_status;
}