/* nag_mesh2d_smooth (d06cac) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 *
 * Mark 8 revised, 2004
 *
 */

#include <stdio.h>
#include <nag.h>
#include <nagx04.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nagd06.h>
#include <nagg05.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(int argc, char *argv[])
{
  FILE          *fpin, *fpout;
  /* 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];
  char          *outfile = 0;

  /* 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;

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

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  (void) nag_example_file_io(argc, argv, "-nag_write", &outfile);

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

  fprintf(fpout, " nag_mesh2d_smooth (d06cac) Example Program Results\n\n");

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

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

  /* Read distortion percentage and calculate radius */
  /* of distortion neighbourhood so that cross-over */
  /* can only occur at 100% or greater. */
  fscanf(fpin, "%lf", &delta);
  fscanf(fpin, "%*[^\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)))
    {
      fprintf(fpout, "Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  fscanf(fpin, " ' %1s '", pmesh);
  fscanf(fpin, "%*[^\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);

  /* Initialise the generator to a repeatable sequence */
  nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_rand_init_repeatable (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_uniform(1, 0.0, rad, state, one_draw, &fail);
          if (fail.code != NE_NOERROR)
            {
              fprintf(fpout, "Error from nag_rand_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_uniform(1, 0.0, dpi, state, one_draw, &fail);
          if (fail.code != NE_NOERROR)
            {
              fprintf(fpout, "Error from nag_rand_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')
    {
      fprintf(fpout, " The complete distorted mesh characteristics\n");
      fprintf(fpout, " nv   =%6ld\n", nv);
      fprintf(fpout, " nelt =%6ld\n\n", nelt);
    }
  else if (pmesh[0] == 'Y')
    {
      /* Output the mesh to view it using the NAG Graphics Library */

      fprintf(fpout, " %10ld%10ld\n", nv, nelt);

      for (i = 1; i <= nv; ++i)
        fprintf(fpout, "  %15.6e  %15.6e  \n", COOR(1, i), COOR(2, i));
    }
  else
    {
      fprintf(fpout, "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)
        {
          fprintf(fpout, "Error the surface of the element is negative\n");
          fprintf(fpout, " k = %6ld\n", k);
          fprintf(fpout, " sk = %15.6e\n", sk);
          exit_status = -1;
          goto END;
        }

      if (pmesh[0] == 'Y')
        fprintf(fpout, " %10ld%10ld%10ld%10ld\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_mesh2d_smooth (d06cac).
   * Uses a barycentering technique to smooth a given mesh
   */
  if (outfile) fclose(fpout);
  nag_mesh2d_smooth(nv, nelt, nedge, coor, edge, conn, nvfix, numfix, itrace,
                    outfile, nqint, &fail);
  if (outfile && !(fpout = fopen(outfile, "a")))
    {
      exit_status = 2;
      goto END;
    }
  if (fail.code == NE_NOERROR)
    {
      if (pmesh[0] == 'N')
        {
          fprintf(fpout, "\n The complete smoothed mesh characteristics\n");
          fprintf(fpout, " nv   =%6ld\n", nv);
          fprintf(fpout, " nelt =%6ld\n", nelt);
        }
      else if (pmesh[0] == 'Y')
        {
          /* Output the mesh to view it using the NAG Graphics Library */

          fprintf(fpout, " %10ld%10ld\n", nv, nelt);

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

          reftk = 0;

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

 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (coor) NAG_FREE(coor);
  if (conn) NAG_FREE(conn);
  if (edge) NAG_FREE(edge);
  if (numfix) NAG_FREE(numfix);
  if (state) NAG_FREE(state);

  return exit_status;
}