/* nag_mesh_dim2_join (d06dbc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.2, 2024.
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static double NAG_CALL fbnd(Integer, double, double, Nag_Comm *);
#ifdef __cplusplus
}
#endif

static int ex1(void);
static int ex2(void);

int main(void) {
  Integer exit_status_ex1 = 0;
  Integer exit_status_ex2 = 0;

  printf("nag_mesh_dim2_join (d06dbc) Example Program Results\n");
  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();

  return (exit_status_ex1 == 0 && exit_status_ex2 == 0) ? 0 : 1;
}

int ex1(void) {
  const Integer nvmax = 900, nedmx = 200, neltmx = 2 * nvmax + 5, ntrans = 1,
                mode = 0;
  double eps, dx, xint, pi2, pi4, ddx, ddy;
  Integer exit_status, i, imax, itrace, itrans, j, jmax, jtrans, k, ktrans, l;
  Integer nedge1, nedge2, nedge3, nelt1, nelt2, nelt3, nv1, nv2, nv3;
  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;
  Integer *edge3 = 0, *itype = 0, *reft1 = 0, *reft2 = 0, *reft3 = 0;
  NagError fail;

  INIT_FAIL(fail);
  exit_status = 0;

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

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

  /* Set sizes for the mesh: coordinates and connectivity of the 1st domain */
  nv1 = 400;
  nelt1 = 722;

  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))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Set up interior mesh as small perturbations on regular grid
   * with regularity on the boundary of square on [0,1]x[0,1].
   */
  dx = 1.0 / ((double)(imaxm1));
  xint = 1.0 / ((double)(imax));
  pi2 = X01AAC + X01AAC;
  pi4 = pi2 + pi2;
  k = 0;
  ddy = 0.0;
  for (j = 0; j < jmax; ++j) {
    ddx = 0.0;
    for (i = 0; i < imax; ++i) {
      coor1[k] = ddx + dx * xint * sin(pi4 * ddx) * sin(pi4 * ddy);
      coor1[k + 1] = ddy + dx * xint * sin(pi2 * ddx) * sin(pi2 * ddy);
      k = k + 2;
      ddx = ddx + dx;
    }
    ddy = ddy + dx;
  }

  /* Triangulate using skew-diagonals on grid squares */
  k = 0;
  l = 0;
  for (i = 0; i < imaxm1; ++i) {
    for (j = 0; j < jmaxm1; ++j) {
      l = l + 1;
      conn1[k] = l;
      conn1[k + 1] = l + 1;
      conn1[k + 2] = l + imaxm1 + 2;
      k = k + 3;
      conn1[k] = l;
      conn1[k + 1] = l + imaxm1 + 2;
      conn1[k + 2] = l + imaxm1 + 1;
      k = k + 3;
    }
    l = l + 1;
  }

  for (i = 0; i < nelt1; ++i) {
    reft1[i] = 1;
    reft2[i] = 2;
  }

  pmesh[0] = 'N';

  /* the edges of the boundary */
  ind = 0;

  for (i = 1; i <= imaxm1; ++i) {
    edge1[ind] = i;
    edge1[ind + 1] = i + 1;
    edge1[ind + 2] = 1;
    ind = ind + 3;
  }

  for (i = 1; i <= jmaxm1; ++i) {
    edge1[ind] = i * imax;
    edge1[ind + 1] = (i + 1) * imax;
    edge1[ind + 2] = 1;
    ind = ind + 3;
  }

  for (i = 1; i <= imaxm1; ++i) {
    edge1[ind] = imax * jmax - i + 1;
    edge1[ind + 1] = imax * jmax - i;
    edge1[ind + 2] = 1;
    ind = ind + 3;
  }

  for (i = 1; i <= jmaxm1; ++i) {
    edge1[ind] = (jmax - i) * imax + 1;
    edge1[ind + 1] = (jmax - i - 1) * imax + 1;
    edge1[ind + 2] = 1;
    ind = ind + 3;
  }

  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[0] = (double)itrans / (imax - 1.0);
    trans[1] = (double)jtrans / (jmax - 1.0);
    itrace = 0;

    /* nag_mesh_dim2_transform_affine (d06dac).
     * Generates a mesh resulting from an affine transformation
     * of a given mesh
     */
    nag_mesh_dim2_transform_affine(mode, nv2, nedge2, nelt2, ntrans, itype,
                                   trans, coor1, edge1, conn1, coor2, edge2,
                                   conn2, itrace, 0, &fail);
    if (fail.code == NE_NOERROR) {
      for (i = 0; i < nedge2; ++i)
        edge2[3 * i + 2] = 2;

      /* Call to the restitching driver */

      itrace = 0;
      eps = 0.01;

      /* nag_mesh_dim2_join (d06dbc).
       * Joins together two given adjacent (possibly overlapping)
       * meshes
       */
      nag_mesh_dim2_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) {
            printf("The restitched mesh characteristics\n");
            printf("in the overlapping case\n");
          } else {
            printf("in the partition case\n");
          }
          printf(" nv    =%6" NAG_IFMT "\n", nv3);
          printf(" nelt  =%6" NAG_IFMT "\n", nelt3);
          printf(" nedge =%6" NAG_IFMT "\n", nedge3);
        } else if (pmesh[0] == 'Y') {
          /* Output the mesh to view it using the
             NAG Graphics Library */

          printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "\n", nv3, nelt3,
                 nedge3);

          for (i = 0; i < nv3; ++i)
            printf("  %15.6e  %15.6e\n", coor3[2 * i], coor3[2 * i + 1]);

          for (k = 0; k < nelt3; ++k)
            printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT ""
                   "%10" NAG_IFMT "\n",
                   conn3[3 * k], conn3[3 * k + 1], conn3[3 * k + 2], reft3[k]);

          for (k = 0; k < nedge3; ++k)
            printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "\n",
                   edge3[3 * k], edge3[3 * k + 1], edge3[3 * k + 2]);
        } else {
          printf("Problem with the printing option Y or N\n");
          exit_status = -1;
          goto END;
        }
      } else {
        printf("Error from nag_mesh_dim2_join (d06dbc).\n%s\n", fail.message);
        exit_status = 1;
        goto END;
      }
    } else {
      printf("Error from nag_mesh_dim2_transform_affine (d06dac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }

END:
  NAG_FREE(coor1);
  NAG_FREE(coor2);
  NAG_FREE(coor3);
  NAG_FREE(trans);
  NAG_FREE(conn1);
  NAG_FREE(conn2);
  NAG_FREE(conn3);
  NAG_FREE(edge1);
  NAG_FREE(edge2);
  NAG_FREE(edge3);
  NAG_FREE(itype);
  NAG_FREE(reft1);
  NAG_FREE(reft2);
  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;
  static double ruser[1] = {-1.0};
  double eps;
  Integer exit_status, i, itrace, j, k, l, ncomp, nedge1, nedge2, nedge3;
  Integer nedge4, nedge5, nelt1, nelt2, nelt3, nelt4, nelt5, nlines;
  Integer npropa, nqint, nv1, nv2, nv3, nv4, nv5, nvb1, nvb2, mode;
  char pmesh[2];
  double *coor1 = 0, *coor2 = 0, *coor3 = 0, *coor4 = 0, *coor5 = 0;
  double *coorch = 0, *coorus = 0, *rate = 0, *trans = 0, *weight = 0;
  Integer *conn1 = 0, *conn2 = 0, *conn3 = 0, *conn4 = 0, *conn5 = 0;
  Integer *edge1 = 0, *edge2 = 0, *edge3 = 0, *edge4 = 0, *edge5 = 0;
  Integer *itype = 0, *lcomp = 0, *lined = 0, *nlcomp = 0, *numfix = 0;
  Integer *reft1 = 0, *reft2 = 0, *reft3 = 0, *reft4 = 0, *reft5 = 0;
  NagError fail;
  Nag_Comm comm;

  INIT_FAIL(fail);

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  exit_status = 0;

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

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

  /* Build the mesh of the 1st domain */
  /* Initialize boundary mesh inputs: */
  /* the number of line and of the characteristic points of */
  /* the boundary mesh */
  scanf("%" NAG_IFMT "", &nlines);
  scanf("%*[^\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))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Characteristic points of the boundary geometry */

  for (j = 0; j < nlines; ++j) {
    scanf("%lf", &coorch[2 * j]);
  }
  scanf("%*[^\n] ");

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

  /* The lines of the boundary mesh */

  for (j = 0; j < nlines; ++j) {
    for (i = 0; i < 4; ++i)
      scanf("%" NAG_IFMT "", &lined[4 * j + i]);
    scanf("%lf", &rate[j]);
  }
  scanf("%*[^\n] ");

  /* The number of connected components */
  /* on the boundary and their data */
  scanf("%" NAG_IFMT "", &ncomp);
  scanf("%*[^\n] ");

  /* Allocate memory */

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

  j = 0;

  for (i = 0; i < ncomp; ++i) {
    scanf("%" NAG_IFMT "", &nlcomp[i]);
    scanf("%*[^\n] ");
    l = j + NAG_IABS(nlcomp[i]);
    for (k = j; k < l; ++k)
      scanf("%" NAG_IFMT "", &lcomp[k]);
    scanf("%*[^\n] ");
    j += NAG_IABS(nlcomp[i]);
  }

  itrace = 0;

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

  /* nag_mesh_dim2_gen_boundary (d06bac).
   * Generates a boundary mesh
   */
  nag_mesh_dim2_gen_boundary(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) {
    printf("Error from nag_mesh_dim2_gen_boundary (d06bac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Generate mesh using Delaunay-Voronoi method */

  /* Initialize mesh control parameters */

  itrace = 0;
  npropa = 1;

  /* nag_mesh_dim2_gen_delaunay (d06abc).
   * Generates a two-dimensional mesh using a Delaunay-Voronoi
   * process
   */
  nag_mesh_dim2_gen_delaunay(nvb1, nvint, nvmax, nedge1, edge1, &nv1, &nelt1,
                             coor1, conn1, weight, npropa, itrace, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_mesh_dim2_gen_delaunay (d06abc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

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

  /* Call the smoothing routine */

  nqint = 10;
  /* nag_mesh_dim2_smooth_bary (d06cac).
   * Uses a barycentering technique to smooth a given mesh
   */
  nag_mesh_dim2_smooth_bary(nv1, nelt1, nedge1, coor1, edge1, conn1, nvfix,
                            numfix, itrace, 0, nqint, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_mesh_dim2_smooth_bary (d06cac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Build the mesh of the 2nd domain */
  scanf("%" NAG_IFMT "", &nlines);
  scanf("%*[^\n] ");

  /* Characteristic points of the boundary geometry */
  for (j = 0; j < nlines; ++j)
    scanf("%lf", &coorch[2 * j]);
  scanf("%*[^\n] ");
  for (j = 0; j < nlines; ++j)
    scanf("%lf", &coorch[2 * j + 1]);
  scanf("%*[^\n] ");

  /* The lines of the boundary mesh */

  for (j = 0; j < nlines; ++j) {
    for (i = 0; i < 4; ++i)
      scanf("%" NAG_IFMT "", &lined[4 * j + i]);
    scanf("%lf", &rate[j]);
  }
  scanf("%*[^\n] ");

  /* The number of connected components */
  /* to the boundary and their data */
  scanf("%" NAG_IFMT "", &ncomp);
  scanf("%*[^\n] ");

  j = 0;
  for (i = 0; i < ncomp; ++i) {
    scanf("%" NAG_IFMT "", &nlcomp[i]);
    scanf("%*[^\n] ");

    for (k = j; k <= j + NAG_IABS(nlcomp[i]); ++k)
      scanf("%" NAG_IFMT "", &lcomp[k]);
    scanf("%*[^\n] ");
    j += NAG_IABS(nlcomp[i]);
  }
  scanf(" ' %1s '", pmesh);
  scanf("%*[^\n] ");

  itrace = 0;

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

  /* nag_mesh_dim2_gen_boundary (d06bac), see above. */
  nag_mesh_dim2_gen_boundary(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) {
    printf("Error from nag_mesh_dim2_gen_boundary (d06bac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Generate mesh using the advancing front method */

  itrace = 0;

  /* nag_mesh_dim2_gen_front (d06acc).
   * Generates a two-dimensional mesh using an Advancing-front
   * method
   */
  nag_mesh_dim2_gen_front(nvb2, nvint, nvmax, nedge2, edge2, &nv2, &nelt2,
                          coor2, conn2, weight, itrace, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_mesh_dim2_gen_front (d06acc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

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

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

  itype[0] = 3;
  trans[0] = 6.0;
  trans[1] = -1.0;
  trans[2] = -90.0;
  itrace = 0;
  mode = 0;

  /* nag_mesh_dim2_transform_affine (d06dac), see above. */
  nag_mesh_dim2_transform_affine(mode, nv2, nedge2, nelt2, ntrans, itype, trans,
                                 coor2, edge2, conn2, coor3, edge3, conn3,
                                 itrace, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_mesh_dim2_transform_affine (d06dac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

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

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

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

  eps = 0.001;
  itrace = 0;

  /* nag_mesh_dim2_join (d06dbc), see above. */
  nag_mesh_dim2_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) {
    printf("Error from nag_mesh_dim2_join (d06dbc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

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

  itrace = 0;

  /* nag_mesh_dim2_join (d06dbc), see above. */
  nag_mesh_dim2_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) {
    printf("Error from nag_mesh_dim2_join (d06dbc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  if (pmesh[0] == 'N') {
    printf("The restitched mesh characteristics\n");
    printf(" nv    =%6" NAG_IFMT "\n", nv5);
    printf(" nelt  =%6" NAG_IFMT "\n", nelt5);
    printf(" nedge =%6" NAG_IFMT "\n", nedge5);
  } else if (pmesh[0] == 'Y') {
    /* Output the mesh to view it using the NAG Graphics Library */

    printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "\n", nv5, nelt5,
           nedge5);

    for (i = 0; i < nv5; ++i)
      printf("  %15.6e  %15.6e\n", coor5[2 * i], coor5[2 * i + 1]);

    for (k = 0; k < nelt5; ++k)
      printf("%10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "\n",
             conn5[3 * k], conn5[3 * k + 1], conn5[3 * k + 2], reft5[k]);

    for (k = 0; k < nedge5; ++k)
      printf(" %10" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "\n", edge5[3 * k],
             edge5[3 * k + 1], edge5[3 * k + 2]);
  } else {
    printf("Problem with the printing option Y or N\n");
  }

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

  return exit_status;
}

double NAG_CALL fbnd(Integer i, double x, double y, Nag_Comm *pcomm) {
  double radius2, x0, y0, ret_val;

  if (pcomm->user[0] == -1.0) {
    printf("(User-supplied callback fbnd, first invocation.)\n");
    pcomm->user[0] = 0.0;
  }

  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;
}