# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>

# include "ppc_array.h"
# include "ppc_mesh.h"
# include "ppc_mesh_spec.h"
# include "ppc_twb_quad.h"
# include "ppc_xmalloc.h"

int                  main ( int argc, char **argv );
struct problem_spec *annulus ( int annulus_points );
double               annulus_f ( double x, double y );
void                 annulus_free ( struct problem_spec *spec );
void                 do_demo ( struct problem_spec *spec, double a, 
                       struct TWB_qdat *qdat, int n, char *gv_filename );
void                 draw_triangles_mono ( struct mesh *mesh, FILE *fp );
void                 draw_triangles_zhue ( struct mesh *mesh, FILE *fp );
void                 eval_f ( struct mesh *mesh, 
                       double (*f)(double x, double y) );
void                 get_range ( struct mesh *mesh, double *min, double *max );
float                huefunc ( float s );
void                 plot_with_geomview_mono ( struct mesh *mesh, char *filename );
void                 plot_with_geomview_zhue ( struct mesh *mesh, char *filename );
void                 show_usage ( char *progname );
struct problem_spec *square ( );
double               square_f ( double x, double y );
struct problem_spec *three_holes ( int n );
double               three_holes_f ( double x, double y );
void                 three_holes_free ( struct problem_spec *spec );
void                 timestamp ( );
struct problem_spec *triangle_with_hole();
double               triangle_with_hole_f ( double x, double y );

/******************************************************************************/

int main ( int argc, char** argv )

/******************************************************************************/
/*
  Purpose:
 
    ppc_twb_quad_test() tests ppc_twb_quad().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    John Burkardt
*/
{
  int d;
  double a;
  int annulus_points;
  char *endptr;
  int n;
  struct TWB_qdat *qdat;
  struct problem_spec *spec;
  int three_holes_points;

  timestamp ( );

  printf ( "\n" );
  printf ( "ppc_twb_quad_test():\n" );
  printf ( "  C version.\n" );
  printf ( "  Test ppc_twb_quad()\n" );

  if ( argc != 3 )
  {
    show_usage ( argv[0] );
    return EXIT_FAILURE;
  }

  d = strtol ( argv[1], &endptr, 10 );
  if ( *endptr != '\0' || d < 1 )
  {
    show_usage ( argv[0] );
    return EXIT_FAILURE;
  }

  a = strtod ( argv[2], &endptr );
  if ( *endptr != '\0' || a <= 0.0 )
  {
    show_usage ( argv[0] );
    return EXIT_FAILURE;
  }
/*
  Get the quadrature rule.
*/
  qdat = twb_qdat ( &d, &n );
  printf ( "\n" );
  printf ( "  Quadrature strength = %d:\n", d );
  printf ( "  Number of quadrature points = %d\n", n );
/*
  Estimate integral over the annulus.
*/
  annulus_points = 24;
  spec = annulus ( annulus_points );
  do_demo ( spec, a, qdat, n, "annulus.gv" );
  annulus_free ( spec );
/*
  Estimate integral over the square.
*/
  spec = square ( );
  do_demo ( spec, a, qdat, n, "square.gv" );
/*
  Estimate integral over the three holes region.
*/
  three_holes_points = 12;
  spec = three_holes ( three_holes_points );
  do_demo ( spec, a, qdat, n, "three_holes.gv" );
  three_holes_free ( spec );
/*
  Estimate integral over the triangle with hole.
*/
  spec = triangle_with_hole ( );
  do_demo ( spec, a, qdat, n, "triangle_with_hole.gv" );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "ppc_twb_quad_minimal_test():\n" );
  printf ( "  Normal end of execution.\n" );

  timestamp ( );

  return 0;
}
/******************************************************************************/

struct problem_spec *annulus ( int annulus_points )

/******************************************************************************/
/*
  Purpose:

    annulus() defines the annulus region.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    int annulus_points: the number of points used to approximate the circles.
*/
{
  double a = 0.325;
  double b = 2.0 * a;
  int i;
  double PI = 4.0 * atan ( 1.0 );
  struct problem_spec *spec = xmalloc ( sizeof ( *spec ) );
  double t;

  make_vector ( spec->points, 2 * annulus_points );
  make_vector ( spec->segments, 2 * annulus_points );
  make_vector ( spec->holes, 1 );

  for ( i = 0; i < annulus_points; i++ )
  {
    t = 2.0 * i * PI / annulus_points;

    spec->points[i].point_no = i;
    spec->points[i].x = a * cos ( t );
    spec->points[i].y = a * sin ( t );
    spec->points[i].bc = FEM_BC_DIRICHLET;

    spec->points[i+annulus_points].point_no = i + annulus_points;
    spec->points[i+annulus_points].x = b * cos ( t );
    spec->points[i+annulus_points].y = b * sin ( t );
    spec->points[i+annulus_points].bc = FEM_BC_DIRICHLET;
  }

  for ( i = 0; i < annulus_points; i++ )
  {
    spec->segments[i].segment_no = i;
    spec->segments[i].point_no_1 = i;
    spec->segments[i].point_no_2 = i + 1;
    spec->segments[i].bc = FEM_BC_DIRICHLET;

    spec->segments[i+annulus_points].segment_no = i + annulus_points;
    spec->segments[i+annulus_points].point_no_1 = i + annulus_points;
    spec->segments[i+annulus_points].point_no_2 = i + annulus_points + 1;
    spec->segments[i+annulus_points].bc = FEM_BC_DIRICHLET;
  }

  spec->segments[annulus_points-1].point_no_2 = 
    spec->segments[annulus_points-1].point_no_2 - annulus_points;
  spec->segments[2*annulus_points-1].point_no_2 = 
    spec->segments[2*annulus_points-1].point_no_2 - annulus_points;

  spec->holes[0].x = 0.0;
  spec->holes[0].y = 0.0;

  spec->npoints = 2 * annulus_points;
  spec->nsegments = 2 * annulus_points;
  spec->nholes = 1;
  spec->f = annulus_f;
  spec->g = NULL;
  spec->h = NULL;
  spec->eta = NULL;
  spec->u_exact = NULL;

  printf ( "\n" );
  printf ( "  Domain is annulus approximated by a %d-gon\n", annulus_points );
  printf ( "  Inner radius = %g, outer radius = %g\n", a, b );

  return spec;
}
/******************************************************************************/

double annulus_f ( double x, double y )

/******************************************************************************/
/*
  Purpose:

    annulus_f() is the integrand for the annulus problem.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Input:

    double x, y: the evaluation point.

  Output:

    double value: the value of f(x,y).
*/
{
  double theta;
  double value;

  theta = atan2 ( y, x );

  value = 0.20 * cos ( 3.0 * theta );

  return value;
}
/******************************************************************************/

void annulus_free ( struct problem_spec *spec )

/******************************************************************************/
/*
  Purpose:

    annulus_free() frees memory that was used for the annulus problem.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct problem_spec *spec: the structure whose memory is to be freed.
*/
{
  if ( spec != NULL )
  {
    free_vector ( spec->points );
    free_vector ( spec->segments );
    free_vector ( spec->holes );
    free ( spec );
  }
  return;
}
/******************************************************************************/

void do_demo ( struct problem_spec *spec, double a, struct TWB_qdat *qdat, 
  int n, char *gv_filename )

/******************************************************************************/
/*
  Purpose:

    do_demo() estimates an integral over a triangulated region.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct problem_spec *spec: the structure whose region is to be triangulated.

    double a: the size of the maximal triangle.

    struct TWB_qdat *qdat: the chosen quadrature rule.

    int n: the number of quadrature points.

    char *gv_filename: a name for the geomview() file to be created.
*/
{
  double area_sum;
  int i;
  int j;
  int k;
  double lambda[3];
  struct mesh *mesh;
  double sum;
  double x;
  double y;

  mesh = make_mesh ( spec, a );

  printf ( "  vertices = %d\n", mesh->nnodes );
  printf ( "  edges    = %d\n", mesh->nedges );
  printf ( "  elements = %d\n", mesh->nelems );

  eval_f ( mesh, spec->f );

  plot_with_geomview_mono (mesh, gv_filename );
/*
  Integrate over the region.
*/
  area_sum = 0.0;
  sum = 0.0;
/*
  Choose the i-th triangular element.
*/
  for ( i = 0; i < mesh->nelems; i++ )
  {
    struct elem *ep = &mesh->elems[i];

    area_sum = area_sum + ep->area;
/*
  Compute the location of the J-th quadrature point in the element.
*/
    for ( j = 0; j < n; j++ )
    {
      lambda[0] = qdat[j].lambda1;
      lambda[1] = qdat[j].lambda2;
      lambda[2] = qdat[j].lambda3;

      x = 0.0;
      y = 0.0;
      for ( k = 0; k < 3; k++ )
      {
        x = x + lambda[k] * ep->n[k]->x;
        y = y + lambda[k] * ep->n[k]->y;
      }

      sum = sum + ep->area * qdat[j].weight * spec->f ( x, y );
    }
  }

  sum = sum / TWB_STANDARD_AREA;

  printf ( "  Integral estimate = %g\n", sum );
  printf ( "  Area of region    = %g\n", area_sum );
/*
  Free memory.
*/
  free_mesh ( mesh );

  return;
}
/******************************************************************************/

void draw_triangles_mono ( struct mesh *mesh, FILE *fp )

/******************************************************************************/
/*
  Purpose:

    draw_triangles_mono ( ) plots node grays and triangles.

  Discussion:

    In an OFF file, the triangles may be specified as
    3 n1 n2 n3 R G B
    where n1 n2 n3 are the triangle's vertex numbers and R G B
    (each between 0 and 1) define the triangle's color.  Alternatively,
    we may omit the R G B, as in:
    3 n1 n2 n3
    in which case the triangle will be shown in a default silvery color.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int i;

  fputs ( "OFF\n", fp );
  fprintf ( fp, "%d %d -1\n", mesh->nnodes, mesh->nelems );

  for ( i = 0; i < mesh->nnodes; i++ )
  {
    struct node *np = &mesh->nodes[i];
    fprintf ( fp, "%g %g %g\n", np->x, np->y, np->z );
  }

  for ( i = 0; i < mesh->nelems; i++ )
  {
    struct elem *ep = &mesh->elems[i];

    fprintf ( fp, "3 %d %d %d ",
      ep->n[0]->nodeno,
      ep->n[1]->nodeno,
      ep->n[2]->nodeno);
    fputs ( "1.0 1.0 0.0\n", fp );
  }

  return;
}
/******************************************************************************/

void draw_triangles_zhue ( struct mesh *mesh, FILE *fp )

/******************************************************************************/
/* 
  Purpose:

    draw_triangles_zhue ( ) specifies node colors and triangles.

  Discussion:

    In a COFF file, colors are specified at polygon vertices.
    These are interpolated smoothly to shade the polygon.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
# define REDHUE(s)       ( huefunc((s)-23.0/64) )
# define GREENHUE(s)     ( huefunc((s) -7.0/64) )
# define BLUEHUE(s)      ( huefunc((s) +9.0/64) )

  int i;
  double max;
  double min;
  double s;

  get_range ( mesh, &min, &max );

  fputs ( "COFF\n", fp );
  fprintf ( fp, "%d %d %d\n", mesh->nnodes, mesh->nelems, mesh->nedges );

  for ( i = 0; i < mesh->nnodes; i++ )
  {
    struct node *np = &mesh->nodes[i];
    if ( max - min < 1.0e-7 )
    {
      s = 0.5;
    }
    else
    {
      s = ( np->z - min ) / ( max - min );
      fprintf ( fp, "%g %g %g %g %g %g 1.0\n", np->x, np->y, np->z,
        REDHUE(s), GREENHUE(s), BLUEHUE(s) );
    }
  }

  for ( i = 0; i < mesh->nelems; i++ )
  {
    struct elem *ep = &mesh->elems[i];

    fprintf ( fp, "3 %d %d %d\n",
      ep->n[0]->nodeno,
      ep->n[1]->nodeno,
      ep->n[2]->nodeno );
  }
  return;
}
/******************************************************************************/

void eval_f ( struct mesh *mesh, double (*f)(double x, double y) )

/******************************************************************************/
/*
  Purpose:

    eval_f() evaluates the function at every mesh node and stores in z.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int i;
  double x;
  double y;

  for ( i = 0; i < mesh->nnodes; i++ )
  {
    x = mesh->nodes[i].x;
    y = mesh->nodes[i].y;
    mesh->nodes[i].z = f ( x, y );
  }
  return;
}
/******************************************************************************/

void get_range ( struct mesh *mesh, double *min, double *max )

/******************************************************************************/
/*
  Purpose:

    get_range() computes the min and max of the z values of the nodes of a mesh.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int i;

  *min = mesh->nodes[0].z;
  *max = mesh->nodes[0].z;

  for ( i = 1; i < mesh->nnodes; i++ )
  {
    double z = mesh->nodes[i].z;

    if ( z < *min )
    {
      *min = z;
    }
    else if ( *max < z )
    {
      *max = z;
    }
  }
  return;
}
/******************************************************************************/

float huefunc ( float s )

/******************************************************************************/
/*
  Purpose:

    huefunc() defines the RGB color scale.

  Discussion:

    This defines a continuous, piecewise-linear function defined on
    the entire real line.  Its graph is trapezoid-shaped and connects
    the points (0,0) (16/64,1), (33/64,1), (49/64,0).  The function
    is zero outside the support of the trapezoid.
    In combination with REDHUE, GREENHUE, BLUEHUE macros defined above,
    produces an RGB scale used for our ZHUE rendering which coincides
    with Matlab's JET colormap.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  if ( s < 0.0 )
  {
    return 0.0;
  }
  else if ( s < 16.0 / 64.0 )
  {
    return 4.0 * s;
  }
  else if ( s < 33.0 / 64.0 )
  {
    return 1.0;
  }
  else if ( s < 49.0 / 64.0 )
  {
    return 1.0 - 4.0 * ( s - 33.0 / 64.0 );
  }
  else
  {
    return 0.0;
  }
}
/******************************************************************************/

void plot_with_geomview_mono ( struct mesh *mesh, char *filename )

/******************************************************************************/
/*
  Purpose:

    plot_with_geomview_mono ( ): grayscale plot of a scalar over a mesh.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  FILE *fp;

  fp = fopen ( filename, "w" );

  if ( fp == NULL )
  {
    fprintf ( stderr, "unable to open file %s for writing\n", filename );
    return;
  }

  fputs ( "{ appearance { +edge }\n", fp );
  draw_triangles_mono ( mesh, fp );
  fputs ( "}\n", fp );
  fclose ( fp );
  printf ( "geomview output written to file %s\n", filename );

  return;
}
/******************************************************************************/

void plot_with_geomview_zhue ( struct mesh *mesh, char *filename )

/******************************************************************************/
/*
  Purpose:

    plot_with_geomview_zhue ( ) creates a color plot of a scalar over a mesh.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  FILE *fp;

  fp = fopen ( filename, "w" );

  if ( fp == NULL )
  {
    fprintf ( stderr, "unable to open file %s for writing\n", filename );
    return;
  }

  fputs ( "{ appearance { +edge +csmooth }\n", fp );
  draw_triangles_zhue ( mesh, fp );
  fputs ( "}\n", fp );
  fclose ( fp );
  printf ( "geomview output written to file %s\n", filename );

  return;
}
/******************************************************************************/

void show_usage ( char *progname )

/******************************************************************************/
/*
  Purpose:

    show_usage() prints usage instructions.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 June 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    char *progname: the name of the program being executed.
*/
{
  printf ( "\n" );
  printf ( "show_usage():\n" );
  printf ( "  %s d  a\n", progname );
  printf ( "  d = quadrature rule strength.\n" );
  printf ( "  a = maximal triangle area.\n" );

  return;
}
/******************************************************************************/

struct problem_spec *square ( )

/******************************************************************************/
/*
  Purpose:

    square() defines the square region.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  static struct problem_spec_point points[] =
  {
    { 0,  0.0,   0.0,  FEM_BC_DIRICHLET },
    { 1,  1.0,   0.0,  FEM_BC_DIRICHLET },
    { 2,  1.0,   1.0,  FEM_BC_DIRICHLET },
    { 3,  0.0,   1.0,  FEM_BC_DIRICHLET },
  };

  static struct problem_spec_segment segments[] =
  {
    { 0, 0, 1, FEM_BC_DIRICHLET },
    { 1, 1, 2, FEM_BC_DIRICHLET },
    { 2, 2, 3, FEM_BC_DIRICHLET },
    { 3, 3, 0, FEM_BC_DIRICHLET },
  };

  static struct problem_spec_hole holes[] =
  {
  };

  static struct problem_spec spec = 
  {
    .points    = points,
    .segments  = segments,
    .holes     = holes,
    .npoints   = sizeof ( points ) / sizeof ( points[0] ),
    .nsegments = sizeof ( segments ) / sizeof (segments[0] ),
    .nholes    = 0,
    .f         = square_f,
    .g         = NULL,
    .h         = NULL,
    .eta       = NULL,
    .u_exact   = NULL,
  };

  printf ( "\n" );
  printf ( "  This domain is a square.\n" );

  return &spec;
}
/******************************************************************************/

double square_f ( double x, double y )

/******************************************************************************/
/*
  Purpose:

    square_f() is the integrand for the square problem.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Input:

    double x, y: the evaluation point.

  Output:

    double value: the value of f(x,y).
*/
{
  double value;

  value = 36.0 * x * y * ( 1.0 - x ) * ( 1.0 - y );
  
  return value;
}
/******************************************************************************/

struct problem_spec *three_holes ( int n )

/******************************************************************************/
/*
  Purpose:

    three_holes() defines the three holes region.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    int n: the number of points used to approximate the circles.
*/
{
  int i;
  double PI = 4.0 * atan ( 1.0 );
  struct problem_spec *spec = xmalloc ( sizeof ( *spec ) );
  double t;
/*
  Points.
*/
  spec->npoints = 3 * n + 6;
  make_vector ( spec->points, spec->npoints );

  for ( i = 0; i < n; i++ )
  {
    t = 2.0 * i * PI / n;

    spec->points[i].point_no = i;
    spec->points[i].x = 0.5 + 0.125 * cos ( t );
    spec->points[i].y = 0.5 + 0.125 * sin ( t );
    spec->points[i].bc = FEM_BC_DIRICHLET;

    spec->points[i+n].point_no = i + n;
    spec->points[i+n].x = 0.5 + 0.125 * cos ( t );
    spec->points[i+n].y = 1.5 + 0.125 * sin ( t );
    spec->points[i+n].bc = FEM_BC_DIRICHLET;

    spec->points[i+2*n].point_no = i + 2 * n;
    spec->points[i+2*n].x = 1.5 + 0.125 * cos ( t );
    spec->points[i+2*n].y = 1.5 + 0.125 * sin ( t );
    spec->points[i+2*n].bc = FEM_BC_DIRICHLET;
  }
  i = 3 * n;
  spec->points[i].point_no = i;
  spec->points[i].x = 0.0;
  spec->points[i].y = 0.0;
  spec->points[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 1;
  spec->points[i].point_no = i;
  spec->points[i].x = 1.0;
  spec->points[i].y = 0.0;
  spec->points[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 2;
  spec->points[i].point_no = i;
  spec->points[i].x = 1.0;
  spec->points[i].y = 1.0;
  spec->points[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 3;
  spec->points[i].point_no = i;
  spec->points[i].x = 2.0;
  spec->points[i].y = 1.0;
  spec->points[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 4;
  spec->points[i].point_no = i;
  spec->points[i].x = 2.0;
  spec->points[i].y = 2.0;
  spec->points[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 5;
  spec->points[i].point_no = i;
  spec->points[i].x = 0.0;
  spec->points[i].y = 2.0;
  spec->points[i].bc = FEM_BC_DIRICHLET;
/*
  Segments.
*/
  spec->nsegments = 3 * n + 6;
  make_vector ( spec->segments, spec->nsegments );

  for ( i = 0; i < n; i++ )
  {
    spec->segments[i].segment_no = i;
    spec->segments[i].point_no_1 = i;
    if ( i == n - 1 )
    {
      spec->segments[i].point_no_2 = 0;
    }
    else
    {
      spec->segments[i].point_no_2 = i + 1;
    }
    spec->segments[i].bc = FEM_BC_DIRICHLET;

    spec->segments[i+n].segment_no = i + n;
    spec->segments[i+n].point_no_1 = i + n;
    if ( i == n - 1 )
    {
      spec->segments[i+n].point_no_2 = n;
    }
    else
    {
      spec->segments[i+n].point_no_2 = i + n + 1;
    }
    spec->segments[i+n].bc = FEM_BC_DIRICHLET;

    spec->segments[i+2*n].segment_no = i + 2 * n;
    spec->segments[i+2*n].point_no_1 = i + 2 * n;
    if ( i == n - 1 )
    {
      spec->segments[i+2*n].point_no_2 = 2 * n;
    }
    else
    {
      spec->segments[i+2*n].point_no_2 = i + 2 * n + 1;
    }
    spec->segments[i+2*n].bc = FEM_BC_DIRICHLET;
  }

  i = 3 * n;
  spec->segments[i].segment_no = i;
  spec->segments[i].point_no_1 = i;
  spec->segments[i].point_no_2 = i + 1;
  spec->segments[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 1;
  spec->segments[i].segment_no = i;
  spec->segments[i].point_no_1 = i;
  spec->segments[i].point_no_2 = i + 1;
  spec->segments[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 2;
  spec->segments[i].segment_no = i;
  spec->segments[i].point_no_1 = i;
  spec->segments[i].point_no_2 = i + 1;
  spec->segments[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 3;
  spec->segments[i].segment_no = i;
  spec->segments[i].point_no_1 = i;
  spec->segments[i].point_no_2 = i + 1;
  spec->segments[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 4;
  spec->segments[i].segment_no = i;
  spec->segments[i].point_no_1 = i;
  spec->segments[i].point_no_2 = i + 1;
  spec->segments[i].bc = FEM_BC_DIRICHLET;

  i = 3 * n + 5;
  spec->segments[i].segment_no = i;
  spec->segments[i].point_no_1 = i;
  spec->segments[i].point_no_2 = 3 * n;
  spec->segments[i].bc = FEM_BC_DIRICHLET;
/*
  Holes.
*/

/*
  spec->nholes = 1;
  make_vector ( spec->holes, spec->nholes );
  spec->holes[0].x = 0.5;
  spec->holes[0].y = 0.5;
*/
  spec->nholes = 3;
  make_vector ( spec->holes, spec->nholes );

  spec->holes[0].x = 0.5;
  spec->holes[0].y = 0.5;

  spec->holes[1].x = 0.5;
  spec->holes[1].y = 1.5;

  spec->holes[2].x = 1.5;
  spec->holes[2].y = 1.5;

  spec->f = three_holes_f;
  spec->g = NULL;
  spec->h = NULL;
  spec->eta = NULL;
  spec->u_exact = NULL;

  printf ( "\n" );
  printf ( "  Domain is an L-shaped region, with three holes\n" );
  printf ( "  Each hole is approximated by a %d-gon\n", n );

  return spec;
}
/******************************************************************************/

double three_holes_f ( double x, double y )

/******************************************************************************/
/*
  Purpose:

    three_holes_f() is the integrand for the three hole problem.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Input:

    double x, y: the evaluation point.

  Output:

    double value: the value of f(x,y).
*/
{
  double value;

  value = x * x + y * y;

  return value;
}
/******************************************************************************/

void three_holes_free ( struct problem_spec *spec )

/******************************************************************************/
/*
  Purpose:

    three_holes_free() frees memory that was used for the three_holes region.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct problem_spec *spec: the structure whose memory is to be freed.
*/
{
  if ( spec != NULL )
  {
    free_vector ( spec->points );
    free_vector ( spec->segments );
    free_vector ( spec->holes );
    free ( spec );
  }
  return;
}
/******************************************************************************/

void timestamp ( )

/******************************************************************************/
/*
  Purpose:
 
    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}
/******************************************************************************/

struct problem_spec *triangle_with_hole ( )

/******************************************************************************/
/*
  Purpose:

    triangle_with_hole() defines a triangular region with rectangular hole.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  static struct problem_spec_point points[] =
  {
    { 0, -1.0,  0.0,  FEM_BC_DIRICHLET },
    { 1,  1.0,  0.0,  FEM_BC_DIRICHLET },
    { 2,  0.0,  2.0,  FEM_BC_DIRICHLET },
    { 3, -0.25, 0.25, FEM_BC_DIRICHLET },
    { 4, +0.25, 0.25, FEM_BC_DIRICHLET },
    { 5, +0.25, 1.0,  FEM_BC_DIRICHLET },
    { 6, -0.25, 1.0,  FEM_BC_DIRICHLET },
  };

  static struct problem_spec_segment segments[] =
  {
    { 0, 0, 1, FEM_BC_DIRICHLET },
    { 1, 1, 2, FEM_BC_DIRICHLET },
    { 2, 2, 0, FEM_BC_DIRICHLET },
    { 3, 3, 4, FEM_BC_DIRICHLET },
    { 4, 4, 5, FEM_BC_DIRICHLET },
    { 5, 5, 6, FEM_BC_DIRICHLET },
    { 6, 6, 3, FEM_BC_DIRICHLET },
  };

  static struct problem_spec_hole holes[] =
  {
    { 0.0, 0.75 }
  };

  static struct problem_spec spec = 
  {
    .points    = points,
    .segments  = segments,
    .holes     = holes,
    .npoints   = sizeof ( points ) / sizeof ( points[0] ),
    .nsegments = sizeof ( segments ) / sizeof (segments[0] ),
    .nholes    = sizeof ( holes ) / sizeof ( holes[0] ),
    .f         = triangle_with_hole_f,
    .g         = NULL,
    .h         = NULL,
    .eta       = NULL,
    .u_exact   = NULL,
  };

  printf ( "\n" );
  printf ( "  This domain is a triangle with a hole.\n" );

  return &spec;
}
/******************************************************************************/

double triangle_with_hole_f ( double x, double y )

/******************************************************************************/
/*
  Purpose:

    triangle_with_hole_f() is the integrand for the triangle with a hole.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Input:

    double x, y: the evaluation point.

  Output:

    double value: the value of f(x,y).
*/
{
  double value;

  value = 0.4 * x * y * ( 2.0 - x ) * ( 2.0 - y );

  return value;
}