/* nag_mesh2d_join (d06dbc) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 * Mark 7b revised, 2004.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd06.h>
static double fbnd(Integer , double , double , Nag_Comm *);

extern int ex1(void), ex2(void);

#define EDGE1(I,J) edge1[3*((J)-1)+(I)-1]
#define EDGE2(I,J) edge2[3*((J)-1)+(I)-1]
#define EDGE3(I,J) edge3[3*((J)-1)+(I)-1]
#define EDGE4(I,J) edge4[3*((J)-1)+(I)-1]
#define EDGE5(I,J) edge5[3*((J)-1)+(I)-1]
#define CONN1(I,J) conn1[3*((J)-1)+(I)-1]
#define CONN2(I,J) conn2[3*((J)-1)+(I)-1]
#define CONN3(I,J) conn3[3*((J)-1)+(I)-1]
#define CONN4(I,J) conn4[3*((J)-1)+(I)-1]
#define CONN5(I,J) conn5[3*((J)-1)+(I)-1]
#define COOR1(I,J) coor1[2*((J)-1)+(I)-1]
#define COOR2(I,J) coor2[2*((J)-1)+(I)-1]
#define COOR3(I,J) coor3[2*((J)-1)+(I)-1]
#define COOR4(I,J) coor4[2*((J)-1)+(I)-1]
#define COOR5(I,J) coor5[2*((J)-1)+(I)-1]
#define TRANS(I,J) trans[6*((J)-1)+(I)-1]
#define LINED(I,J) lined[4*((J)-1)+(I)-1]
#define COORCH(I,J) coorch[2*(J-1)+I-1]
#define COORUS(I,J) coorus[2*(J-1)+I-1]

int main(void)
{
  Vprintf("nag_mesh2d_join (d06dbc) Example Program Results\n");

  ex1();
  ex2();

  return 0;
}


int ex1(void)
{
  const Integer nvmax=900, nedmx=200, neltmx=2*nvmax+5, ntrans=1, mode=0;
  double eps;
  Integer exit_status, i, imax, itrace, itrans, jmax, jtrans, k, ktrans,
    nedge1, nedge2, nedge3, nelt1, nelt2, nelt3, nv1, nv2, nv3, reftk;
  Integer imaxm1, jmaxm1, ind;
  char pmesh[2];
  double *coor1=0, *coor2=0, *coor3=0, *trans=0;
  Integer *conn1=0, *conn2=0, *conn3=0, *edge1=0,
    *edge2=0, *edge3=0, *itype=0, *reft1=0, *reft2=0, *reft3=0;
  NagError fail;

  INIT_FAIL(fail);
  exit_status = 0;

  Vprintf("\nExample 1\n\n");

  /* Skip heading in data file */

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

  /* Read the mesh: coordinates and connectivity of the 1st domain */

  Vscanf("%ld", &nv1);
  Vscanf("%ld", &nelt1);
  Vscanf("%*[^\n] ");

  nv2 = nv1;
  nelt2 = nelt1;

  imax = 20;
  jmax = imax;
  imaxm1 = imax - 1;
  jmaxm1 = jmax - 1;
  nedge1 = 2*(imaxm1 + jmaxm1);
  nedge2 = nedge1;

  /* Allocate memory */

  if ( !(coor1 = NAG_ALLOC(2*nv1, double)) ||
       !(coor2 = NAG_ALLOC(2*nv2, double)) ||
       !(coor3 = NAG_ALLOC(2*nvmax, double)) ||
       !(trans = NAG_ALLOC(6*ntrans, double)) ||
       !(conn1 = NAG_ALLOC(3*nelt1, Integer)) ||
       !(conn2 = NAG_ALLOC(3*nelt2, Integer)) ||
       !(conn3 = NAG_ALLOC(3*neltmx, Integer)) ||
       !(edge1 = NAG_ALLOC(3*nedge1, Integer)) ||
       !(edge2 = NAG_ALLOC(3*nedge2, Integer)) ||
       !(edge3 = NAG_ALLOC(3*nedmx, Integer)) ||
       !(itype = NAG_ALLOC(ntrans, Integer)) ||
       !(reft1 = NAG_ALLOC(nelt1, Integer)) ||
       !(reft2 = NAG_ALLOC(nelt2, Integer)) ||
       !(reft3 = NAG_ALLOC(neltmx, Integer)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  for (i = 1; i <= nv1; ++i)
    {
      Vscanf("%lf", &COOR1(1,i));
      Vscanf("%lf", &COOR1(2,i));
      Vscanf("%*[^\n] ");
    }

  for (k = 1; k <= nelt1; ++k)
    {
      Vscanf("%ld", &CONN1(1,k));
      Vscanf("%ld", &CONN1(2,k));
      Vscanf("%ld", &CONN1(3,k));
      Vscanf("%ld", &reftk);
      Vscanf("%*[^\n] ");

      reft1[k-1] = 1;
      reft2[k-1] = 2;
    }

  Vscanf(" ' %1s '", pmesh);
  Vscanf("%*[^\n] ");

  /* the edges of the boundary */

  ind = 0;

  for (i = 1; i <= imaxm1; ++i)
    {
      ++ind;
      EDGE1(1,ind) = i;
      EDGE1(2,ind) = i + 1;
      EDGE1(3,ind) = 1;
    }

  for (i = 1; i <= jmaxm1; ++i)
    {
      ++ind;
      EDGE1(1,ind) = i*imax;
      EDGE1(2,ind) = (i+1)*imax;
      EDGE1(3,ind) = 1;
    }

  for (i = 1; i <= imaxm1; ++i)
    {
      ++ind;
      EDGE1(1,ind) = imax*jmax - i + 1;
      EDGE1(2,ind) = imax*jmax - i;
      EDGE1(3,ind) = 1;
    }

  for (i = 1; i <= jmaxm1; ++i)
    {
      ++ind;
      EDGE1(1,ind) = (jmax - i)*imax + 1;
      EDGE1(2,ind) = (jmax - i - 1)*imax + 1;
      EDGE1(3,ind) = 1;
    }

  for (ktrans = 0; ktrans < 2; ++ktrans)
    {
      /* Translation of the 1st domain to obtain the 2nd domain */
      /* KTRANS = 0 leading to a domain overlapping */
      /* KTRANS = 1 leading to a domain partition */

      if (ktrans == 0)
        {
          itrans = imax - 5;
          jtrans = jmax - 3;
        }
      else
        {
          itrans = imax - 1;
          jtrans = 0;
        }

      itype[0] = 1;
      TRANS(1, 1) = (double)itrans/(imax - 1.0);
      TRANS(2, 1) = (double)jtrans/(jmax - 1.0);
      itrace = 0;

      /* nag_mesh2d_trans (d06dac).
       * Generates a mesh resulting from an affine transformation
       * of a given mesh
       */
      nag_mesh2d_trans(mode, nv2, nedge2, nelt2, ntrans, itype, trans, coor1,
                       edge1, conn1, coor2, edge2, conn2, itrace, 0, &fail);

      if (fail.code == NE_NOERROR)
        {
          for (i = 1; i <= nedge2; ++i) EDGE2(3, i) = 2;

          /* Call to the restitching driver */

          itrace = 0;
          eps = 0.01;

          /* nag_mesh2d_join (d06dbc).
           * Joins together two given adjacent (possibly overlapping)
           * meshes
           */
          nag_mesh2d_join(eps, nv1, nelt1, nedge1, coor1, edge1, conn1, reft1,
                          nv2, nelt2, nedge2, coor2, edge2, conn2, reft2, &nv3,
                          &nelt3, &nedge3, coor3, edge3, conn3, reft3, itrace,
                          0, &fail);

          if (fail.code == NE_NOERROR)
            {
              if (pmesh[0] == 'N')
                {
                  if (ktrans == 0)
                    {
                      Vprintf("The restitched mesh characteristics\n");
                      Vprintf("in the overlapping case\n");
                    }
                  else
                    {
                      Vprintf("in the partition case\n");
                    }
                  Vprintf(" nv    =%6ld\n", nv3);
                  Vprintf(" nelt  =%6ld\n", nelt3);
                  Vprintf(" nedge =%6ld\n", nedge3);
                }
              else if (pmesh[0] == 'Y')
                {

                  /* Output the mesh to view it using the
                     NAG Graphics Library */

                  Vprintf(" %10ld%10ld%10ld\n", nv3,
                          nelt3, nedge3);

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

                  for (k = 1; k <= nelt3; ++k)
                    Vprintf(" %10ld%10ld%10ld"
                            "%10ld\n", CONN3(1,k), CONN3(2,k),
                            CONN3(3,k), reft3[k - 1]);

                  for (k = 1; k <= nedge3; ++k)
                    Vprintf(" %10ld%10ld%10ld\n",
                            EDGE3(1,k), EDGE3(2,k), EDGE3(3,k));
                }
              else
                {
                  Vprintf("Problem with the printing option Y or N\n");
                  exit_status = -1;
                  goto END;
                }
            }
          else
            {
              Vprintf("Error from nag_mesh2d_join (d06dbc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }
        }
      else
        {
          Vprintf("Error from nag_mesh2d_trans (d06dac).\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }
    }

 END:

  if (coor1) NAG_FREE(coor1);
  if (coor2) NAG_FREE(coor2);
  if (coor3) NAG_FREE(coor3);
  if (trans) NAG_FREE(trans);
  if (conn1) NAG_FREE(conn1);
  if (conn2) NAG_FREE(conn2);
  if (conn3) NAG_FREE(conn3);
  if (edge1) NAG_FREE(edge1);
  if (edge2) NAG_FREE(edge2);
  if (edge3) NAG_FREE(edge3);
  if (itype) NAG_FREE(itype);
  if (reft1) NAG_FREE(reft1);
  if (reft2) NAG_FREE(reft2);
  if (reft3) NAG_FREE(reft3);

  return exit_status;
}



int ex2(void)
{
  const Integer nvmax=900, nedmx=200, neltmx=2*nvmax+5,
    ntrans=1, nus=0, nvint=0, nvfix=0;
  double eps;
  Integer exit_status, i, itrace, j, k, l, ncomp, nedge1,
    nedge2, nedge3, nedge4, nedge5, nelt1, nelt2, nelt3,
    nelt4, nelt5, nlines, npropa, nqint, nv1, nv2, nv3, nv4,
    nv5, nvb1, nvb2, mode;
  char pmesh[2];
  double *coor1=0, *coor2=0, *coor3=0, *coor4=0, *coor5=0,
    *coorch=0, *coorus=0, *rate=0, *trans=0, *weight=0;
  Integer *conn1=0, *conn2=0, *conn3=0, *conn4=0, *conn5=0,
    *edge1=0, *edge2=0, *edge3=0, *edge4=0, *edge5=0, *itype=0,
    *lcomp=0, *lined=0, *nlcomp=0, *numfix=0, *reft1=0, *reft2=0,
    *reft3=0, *reft4=0, *reft5=0;
  NagError fail;
  Nag_Comm comm;

  INIT_FAIL(fail);
  exit_status = 0;

  Vprintf("\nExample 2\n\n");

  /* Skip heading in data file */

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

  /* Build the mesh of the 1st domain */
  /* Initialise boundary mesh inputs: */
  /* the number of line and of the characteristic points of */
  /* the boundary mesh */

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

  /* Allocate memory */

  if ( !(coor1 = NAG_ALLOC(2*nvmax, double)) ||
       !(coor2 = NAG_ALLOC(2*nvmax, double)) ||
       !(coor3 = NAG_ALLOC(2*nvmax, double)) ||
       !(coor4 = NAG_ALLOC(2*nvmax, double)) ||
       !(coor5 = NAG_ALLOC(2*nvmax, double)) ||
       !(coorch = NAG_ALLOC(2*nlines, double)) ||
       !(coorus = NAG_ALLOC(1, double)) ||
       !(rate = NAG_ALLOC(nlines, double)) ||
       !(trans = NAG_ALLOC(6*ntrans, double)) ||
       !(weight = NAG_ALLOC(1, double)) ||
       !(conn1 = NAG_ALLOC(3*neltmx, Integer)) ||
       !(conn2 = NAG_ALLOC(3*neltmx, Integer)) ||
       !(conn3 = NAG_ALLOC(3*neltmx, Integer)) ||
       !(conn4 = NAG_ALLOC(3*neltmx, Integer)) ||
       !(conn5 = NAG_ALLOC(3*neltmx, Integer)) ||
       !(edge1 = NAG_ALLOC(3*nedmx, Integer)) ||
       !(edge2 = NAG_ALLOC(3*nedmx, Integer)) ||
       !(edge3 = NAG_ALLOC(3*nedmx, Integer)) ||
       !(edge4 = NAG_ALLOC(3*nedmx, Integer)) ||
       !(edge5 = NAG_ALLOC(3*nedmx, Integer)) ||
       !(itype = NAG_ALLOC(ntrans, Integer)) ||
       !(lcomp = NAG_ALLOC(nlines, Integer)) ||
       !(lined = NAG_ALLOC(4*nlines, Integer)) ||
       !(numfix = NAG_ALLOC(1, Integer)) ||
       !(reft1 = NAG_ALLOC(2*nvmax+5, Integer)) ||
       !(reft2 = NAG_ALLOC(2*nvmax+5, Integer)) ||
       !(reft3 = NAG_ALLOC(2*nvmax+5, Integer)) ||
       !(reft4 = NAG_ALLOC(2*nvmax+5, Integer)) ||
       !(reft5 = NAG_ALLOC(2*nvmax+5, Integer)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Characteristic points of the boundary geometry */

  for (j = 1; j <= nlines; ++j)
    {
      Vscanf("%lf", &COORCH(1,j));
    }
  Vscanf("%*[^\n] ");

  for (j = 1; j <= nlines; ++j)
    {
      Vscanf("%lf", &COORCH(2,j));
    }
  Vscanf("%*[^\n] ");

  /* The lines of the boundary mesh */

  for (j = 1; j <= nlines; ++j)
    {
      for (i = 1; i <= 4; ++i) Vscanf("%ld", &LINED(i,j));
      Vscanf("%lf", &rate[j - 1]);
    }
  Vscanf("%*[^\n] ");

  /* The number of connected components */
  /* on the boundary and their data */

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

  /* Allocate memory */

  if (  !(nlcomp = NAG_ALLOC(ncomp, Integer)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  j = 1;

  for (i = 1; i <= ncomp; ++i)
    {
      Vscanf("%ld", &nlcomp[i - 1]);
      Vscanf("%*[^\n] ");
      l = j + abs(nlcomp[i - 1]) - 1;
      for (k = j; k <= l; ++k) Vscanf("%ld", &lcomp[k - 1]);
      Vscanf("%*[^\n] ");
      j += abs(nlcomp[i - 1]);
    }

  itrace = 0;

  /* Call to the 2D boundary mesh generator */

  /* nag_mesh2d_bound (d06bac).
   * Generates a boundary mesh
   */
  nag_mesh2d_bound(nlines, coorch, lined, fbnd, coorus, nus, rate, ncomp,
                   nlcomp, lcomp, nvmax, nedmx, &nvb1, coor1, &nedge1,
                   edge1, itrace, 0, &comm, &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_mesh2d_bound (d06bac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Generate mesh using Delaunay-Voronoi method */

  /* Initialise mesh control parameters */

  itrace = 0;
  npropa = 1;

  /* nag_mesh2d_delaunay (d06abc).
   * Generates a two-dimensional mesh using a Delaunay-Voronoi
   * process
   */
  nag_mesh2d_delaunay(nvb1, nvint, nvmax, nedge1, edge1, &nv1, &nelt1, coor1,
                      conn1, weight, npropa, itrace, 0, &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_mesh2d_delaunay (d06abc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  for (k = 1; k <= nelt1; ++k) reft1[k - 1] = 1;

  /* Call the smoothing routine */

  nqint = 10;
  /* nag_mesh2d_smooth (d06cac).
   * Uses a barycentering technique to smooth a given mesh
   */
  nag_mesh2d_smooth(nv1, nelt1, nedge1, coor1, edge1, conn1, nvfix, numfix,
                    itrace, 0, nqint, &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_mesh2d_smooth (d06cac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Build the mesh of the 2nd domain */

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

  /* Characteristic points of the boundary geometry */

  for (j = 1; j <= nlines; ++j) Vscanf("%lf", &COORCH(1,j));
  Vscanf("%*[^\n] ");

  for (j = 1; j <= nlines; ++j) Vscanf("%lf", &COORCH(2,j));
  Vscanf("%*[^\n] ");

  /* The lines of the boundary mesh */

  for (j = 1; j <= nlines; ++j)
    {
      for (i = 1; i <= 4; ++i) Vscanf("%ld", &LINED(i,j));
      Vscanf("%lf", &rate[j - 1]);
    }
  Vscanf("%*[^\n] ");

  /* The number of connected components */
  /* to the boundary and their data */

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

  j = 1;
  for (i = 1; i <= ncomp; ++i)
    {
      Vscanf("%ld", &nlcomp[i - 1]);
      Vscanf("%*[^\n] ");

      for (k = j; k <= j+abs(nlcomp[i-1])-1; ++k)
        Vscanf("%ld", &lcomp[k - 1]);
      Vscanf("%*[^\n] ");
      j += abs(nlcomp[i-1]);
    }

  Vscanf(" ' %1s '", pmesh);
  Vscanf("%*[^\n] ");

  itrace = 0;

  /* Call to the 2D boundary mesh generator */

  /* nag_mesh2d_bound (d06bac), see above. */
  nag_mesh2d_bound(nlines, coorch, lined, fbnd, coorus, nus, rate, ncomp,
                   nlcomp, lcomp, nvmax, nedmx, &nvb2, coor2, &nedge2, edge2,
                   itrace, 0, &comm, &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_mesh2d_bound (d06bac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Generate mesh using the advancing front method */

  itrace = 0;

  /* nag_mesh2d_front (d06acc).
   * Generates a two-dimensional mesh using an Advancing-front
   * method
   */
  nag_mesh2d_front(nvb2, nvint, nvmax, nedge2, edge2, &nv2, &nelt2, coor2,
                   conn2, weight, itrace, 0, &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_mesh2d_front (d06acc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  for (k = 1; k <= nelt2; ++k) reft2[k - 1] = 2;

  /* Rotation of the 2nd domain mesh */
  /* to produce the 3rd mesh domain */

  itype[0] = 3;
  TRANS(1, 1) = 6.0;
  TRANS(2, 1) = -1.0;
  TRANS(3, 1) = -90.0;
  itrace = 0;
  mode = 0;

  /* nag_mesh2d_trans (d06dac), see above. */
  nag_mesh2d_trans(mode, nv2, nedge2, nelt2, ntrans, itype, trans, coor2,
                   edge2, conn2, coor3, edge3, conn3, itrace, 0, &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_mesh2d_trans (d06dac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  nv3 = nv2;
  nelt3 = nelt2;
  nedge3 = nedge2;

  for (k = 1; k <= nelt3; ++k) reft3[k - 1] = 3;

  /* Restitching meshes 1 and 2 to form mesh 4 */

  eps = 0.001;
  itrace = 0;

  /* nag_mesh2d_join (d06dbc), see above. */
  nag_mesh2d_join(eps, nv1, nelt1, nedge1, coor1, edge1, conn1, reft1, nv2,
                  nelt2, nedge2, coor2, edge2, conn2, reft2, &nv4, &nelt4,
                  &nedge4, coor4, edge4, conn4, reft4, itrace, 0, &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_mesh2d_join (d06dbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Restitching meshes 3 and 4 to form mesh 5 */

  itrace = 0;

  /* nag_mesh2d_join (d06dbc), see above. */
  nag_mesh2d_join(eps, nv4, nelt4, nedge4, coor4, edge4, conn4, reft4, nv3,
                  nelt3, nedge3, coor3, edge3, conn3, reft3, &nv5, &nelt5,
                  &nedge5, coor5, edge5, conn5, reft5, itrace, 0, &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_mesh2d_join (d06dbc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  if (pmesh[0] == 'N')
    {
      Vprintf("The restitched mesh characteristics\n");
      Vprintf(" nv    =%6ld\n", nv5);
      Vprintf(" nelt  =%6ld\n", nelt5);
      Vprintf(" nedge =%6ld\n", nedge5);
    }
  else if (pmesh[0] == 'Y')
    {

      /* Output the mesh to view it using the NAG Graphics Library */

      Vprintf(" %10ld%10ld%10ld\n", nv5, nelt5, nedge5);

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

      for (k = 1; k <= nelt5; ++k)
        Vprintf("%10ld%10ld%10ld%10ld\n",
                CONN5(1,k), CONN5(2,k), CONN5(3,k), reft5[k - 1]);

      for (k = 1; k <= nedge5; ++k)
        Vprintf(" %10ld%10ld%10ld\n",
                EDGE5(1,k), EDGE5(2,k), EDGE5(3,k));
    }
  else
    {
      Vprintf("Problem with the printing option Y or N\n");
    }

 END:

  if (coor1) NAG_FREE(coor1);
  if (coor2) NAG_FREE(coor2);
  if (coor3) NAG_FREE(coor3);
  if (coor4) NAG_FREE(coor4);
  if (coor5) NAG_FREE(coor5);
  if (coorch) NAG_FREE(coorch);
  if (coorus) NAG_FREE(coorus);
  if (rate) NAG_FREE(rate);
  if (trans) NAG_FREE(trans);
  if (weight) NAG_FREE(weight);
  if (conn1) NAG_FREE(conn1);
  if (conn2) NAG_FREE(conn2);
  if (conn3) NAG_FREE(conn3);
  if (conn4) NAG_FREE(conn4);
  if (conn5) NAG_FREE(conn5);
  if (edge1) NAG_FREE(edge1);
  if (edge2) NAG_FREE(edge2);
  if (edge3) NAG_FREE(edge3);
  if (edge4) NAG_FREE(edge4);
  if (edge5) NAG_FREE(edge5);
  if (itype) NAG_FREE(itype);
  if (lcomp) NAG_FREE(lcomp);
  if (lined) NAG_FREE(lined);
  if (nlcomp) NAG_FREE(nlcomp);
  if (numfix) NAG_FREE(numfix);
  if (reft1) NAG_FREE(reft1);
  if (reft2) NAG_FREE(reft2);
  if (reft3) NAG_FREE(reft3);
  if (reft4) NAG_FREE(reft4);
  if (reft5) NAG_FREE(reft5);

  return exit_status;
}
static double fbnd(Integer i, double x, double y, Nag_Comm *pcomm)
{
  double radius2, x0, y0, ret_val;

  ret_val = 0.0;
  switch (i)
    {
    case 1:

      /* inner circle */

      x0 = 0.0;
      y0 = 0.0;
      radius2 = 1.0;
      ret_val = (x-x0)*(x-x0) + (y-y0)*(y-y0) - radius2;
      break;

    case 2:

      /* outer circle */

      x0 = 0.0;
      y0 = 0.0;
      radius2 = 5.0;
      ret_val = (x-x0)*(x-x0) + (y-y0)*(y-y0) - radius2;
      break;
    }

  return ret_val;
}