/* 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 <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(void)
{
  const Integer nvfix=0;
  const double zero=0.0;
  double delta, hx, hy, pi, dpi, r, rad, sk, theta, x1, x2, x3, y1, y2, y3;
  Integer exit_status, i, imax, imaxm1, ind, itrace, j, jmax, jmaxm1, k,
    me1, me2, me3, nedge, nelt, nqint, nv, reftk;
  char pmesh[2];
  double *coor=0;
  Integer *conn=0, *edge=0, *numfix=0;
  /* nag_rngs_init_repeatable (g05kbc) requires iseed[0]=<num> and igen=0 for
   * it to mimic nag_random_init_repeatable(<num>).
   */
  Integer igen = 0, iseed[] = {0, 0, 0, 0};
  double one_draw[1];
  NagError fail;


  INIT_FAIL(fail);
  exit_status = 0;

  Vprintf(" nag_mesh2d_smooth (d06cac) Example Program Results\n\n");

  /* Skip heading in data file */

  Vscanf("%*[^\n] ");

  /* Read imax and jmax, the number of vertices */
  /* in the x and y directions respectively. */

  Vscanf("%ld", &imax);
  Vscanf("%ld", &jmax);
  Vscanf("%*[^\n] ");

  /* Read distortion percentage and calculate radius */
  /* of distortion neighbourhood so that cross-over */
  /* can only occur at 100% or greater. */

  Vscanf("%lf", &delta);
  Vscanf("%*[^\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)) ||
       !(numfix = NAG_ALLOC(1, Integer)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  Vscanf(" ' %1s '", pmesh);
  Vscanf("%*[^\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);
  /* nag_rngs_init_repeatable (g05kbc).
   * Initialize seeds of a given generator for random number
   * generating functions (that pass seeds explicitly) to give
   * a repeatable sequence
   */
  nag_rngs_init_repeatable(&igen, iseed);

  /* 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)
        {
          /* nag_rngs_uniform (g05lgc).
           * Generates a vector of random numbers from a uniform
           * distribution, seeds and generator number passed
           * explicitly
           */
          nag_rngs_uniform(zero, rad, 1, one_draw, igen, iseed, NAGERR_DEFAULT);
          r = one_draw[0];

          dpi = 2.0*pi;
          /* nag_rngs_uniform (g05lgc), see above. */
          nag_rngs_uniform(zero, dpi, 1, one_draw, igen, iseed, NAGERR_DEFAULT);
          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')
    {
      Vprintf(" The complete distorted mesh characteristics\n");
      Vprintf(" nv   =%6ld\n", nv);
      Vprintf(" nelt =%6ld\n", nelt);
    }
  else if (pmesh[0] == 'Y')
    {
      /* Output the mesh to view it using the NAG Graphics Library */

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

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

      if (pmesh[0] == 'Y')
        Vprintf(" %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
   */
  nag_mesh2d_smooth(nv, nelt, nedge, coor, edge, conn, nvfix, numfix, itrace,
                    0, nqint, &fail);

  if (fail.code == NE_NOERROR)
    {
      if (pmesh[0] == 'N')
        {
          Vprintf(" The complete smoothed mesh characteristics\n");
          Vprintf(" nv   =%6ld\n", nv);
          Vprintf(" nelt =%6ld\n", nelt);
        }
      else if (pmesh[0] == 'Y')
        {
          /* Output the mesh to view it using the NAG Graphics Library */

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

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

          reftk = 0;

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

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

  return exit_status;
}