/* nag_mesh_dim2_renumber (d06ccc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */

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

int main(void) {
  /* Scalars */
  const Integer nvb1 = 40, nvb2 = 30, nvb3 = 30, nvmax = 260, pmesh = 0;
  Integer exit_status = 0;
  Integer i, i1, itrace, nedge, nelt, nnz, nnzmax, nv, nvb, reftk;
  double coef, pi2, power, r, theta, theta_i, x0, y0;
  /* Arrays */
  double *coor = 0, *bspace = 0;
  Integer *conn = 0, *edge = 0, *icol = 0, *irow = 0;
  /* Nag Types */
  Nag_Boolean smooth;
  NagError fail;

  INIT_FAIL(fail);

  printf(" nag_mesh_dim2_renumber (d06ccc) Example Program Results\n\n");
  fflush(stdout);

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

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

  /* Stage 1: Create Mesh using nag_mesh_dim2_gen_inc (d06aac) */

  pi2 = 2.0 * X01AAC;
  i1 = 0;
  /* 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] = 0;
    i1 = i1 + 3;
  }
  /* Tie up end of three boundary edges */
  edge[3 * nvb1 - 2] = 1;
  edge[3 * (nvb1 + nvb2) - 2] = nvb1 + 1;
  edge[3 * nvb - 2] = nvb1 + nvb2 + 1;

  /* 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(nvb, nvmax, nedge, edge, &nv, &nelt, coor, conn, bspace,
                        smooth, coef, power, itrace, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_mesh_dim2_gen_inc (d06aac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Stage 2
   * Compute the sparsity of the FE matrix from the conn as generated above.
   *
   * nag_mesh_dim2_sparsity (d06cbc).
   * Generates a sparsity pattern of a Finite Element matrix
   * associated with a given mesh
   */
  nag_mesh_dim2_sparsity(nv, nelt, nnzmax, conn, &nnz, irow, icol, &fail);

  if (fail.code == NE_NOERROR) {
    if (pmesh == 0) {
      printf(" Matrix Sparsity characteristics before renumbering\n");
      printf(" nv   =%6" NAG_IFMT "\n", nv);
      printf(" nnz  =%6" NAG_IFMT "\n", nnz);
      printf(" nelt =%6" NAG_IFMT "\n", nelt);
    } else {
      /* Output the sparsity of the mesh to view */

      printf(" %10" NAG_IFMT "%10" NAG_IFMT "\n", nv, nnz);
      for (i = 0; i < nnz; ++i)
        printf(" %10" NAG_IFMT "%10" NAG_IFMT "\n", irow[i], icol[i]);
    }
  } else {
    printf("Error from nag_mesh_dim2_sparsity (d06cbc).\n%s\n", fail.message);
    exit_status = 2;
    goto END;
  }

  /* Stage 3: Call the renumbering routine and get the new sparsity */

  itrace = 1;

  /* nag_mesh_dim2_renumber (d06ccc).
   * Renumbers a given mesh using Gibbs method
   */
  fflush(stdout);
  nag_mesh_dim2_renumber(nv, nelt, nedge, nnzmax, &nnz, coor, edge, conn, irow,
                         icol, itrace, 0, &fail);
  if (fail.code == NE_NOERROR) {
    if (pmesh == 0) {
      printf("\n Matrix Sparsity characteristics after renumbering\n");
      printf(" nv   =%6" NAG_IFMT "\n", nv);
      printf(" nnz  =%6" NAG_IFMT "\n", nnz);
      printf(" nelt =%6" NAG_IFMT "\n", nelt);
    } else {
      /* Output the sparsity of the renumbered mesh */
      /* to view it using the NAG Graphics Library */

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

      for (i = 0; i < nnz; ++i)
        printf(" %10" NAG_IFMT "%10" NAG_IFMT "\n", irow[i], icol[i]);

      /* Output the renumbered 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 (i = 0; i < nelt; ++i)
        printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT
               "\n",
               conn[3 * i], conn[3 * i + 1], conn[3 * i + 2], reftk);
    }
  } else {
    printf("Error from nag_mesh_dim2_renumber (d06ccc).\n%s\n", fail.message);
    exit_status = 3;
    goto END;
  }

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

  return exit_status;
}