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

# include "ppc_array.h"
# include "ppc_mesh.h"
# include "ppc_problem_spec.h"
# include "ppc_xmalloc.h"

int main ( );
struct problem_spec *annulus ( int n );
void do_demo ( struct problem_spec *spec, double a, char *eps_filename );
void free_annulus ( struct problem_spec *spec );
void free_three_holes ( struct problem_spec *spec );
void show_usage ( char *progname );
struct problem_spec *square ( );
struct problem_spec *three_holes ( int n );
void timestamp ( );
struct problem_spec *triangle_with_hole();

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

int main ( int argc, char **argv )

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

    ppc_mesh_test() tests ppc_mesh().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    21 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
*/
{
  double a;
  char *endptr;
  int n;
  struct problem_spec *spec;

  timestamp ( );
  printf ( "\n" );
  printf ( "ppc_mesh_test ( ):\n" );
  printf ( "  C version.\n" );
  printf ( "  Test ppc_mesh().\n" );

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

  a = strtod ( argv[1], &endptr );

  if ( *endptr != '\0' || a <= 0.0 )
  {
    show_usage ( argv[0] );
    return EXIT_FAILURE;
  }

  n = 24;
  spec = annulus ( n );
  do_demo ( spec, a, "annulus.eps" );
  printf ( "\n" );

  do_demo ( square ( ), a, "square.eps" );
  printf ( "\n" );

  n = 12;
  spec = three_holes ( n );
  do_demo ( spec, a, "three_holes.eps" );
  printf ( "\n" );

  do_demo ( triangle_with_hole ( ), a, "triangle_with_hole.eps" );
  printf ( "\n" );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "ppc_mesh_test ( ):\n" );
  printf ( "  Normal end of execution.\n" );
  timestamp ( );

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

struct problem_spec *annulus ( int n )

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

    annulus() defines the annulus region.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    21 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

  Input:

    int n: 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 * n );
  make_vector ( spec->segments, 2 * n );
  make_vector ( spec->holes, 1 );

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

    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+n].point_no = i + n;
    spec->points[i+n].x = b * cos ( t );
    spec->points[i+n].y = b * sin ( t );
    spec->points[i+n].bc = FEM_BC_DIRICHLET;
  }

  for ( i = 0; i < n; 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+n].segment_no = i + n;
    spec->segments[i+n].point_no_1 = i + n;
    spec->segments[i+n].point_no_2 = i + n + 1;
    spec->segments[i+n].bc = FEM_BC_DIRICHLET;
  }

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

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

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

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

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

void do_demo ( struct problem_spec *spec, double a, char *eps_filename )

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

    do_demo() demonstates the meshing of four regions.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    21 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

  Input:

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

    double a: the size of the maximal triangle.

    char *eps_filename: a name for the EPS file to be created.
*/
{
  struct mesh *mesh = make_mesh ( spec, a );

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

  mesh_to_eps ( mesh, eps_filename );

  free_mesh ( mesh );

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

void free_annulus ( struct problem_spec *spec )

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

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

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    21 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

  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 free_three_holes ( struct problem_spec *spec )

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

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

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    21 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

  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 show_usage ( char *progname )

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

    show_usage() prints usage instructions.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    21 May 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 a\n", progname );
  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:

    21 May 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, -1.0,  -1.0,  FEM_BC_DIRICHLET },
    { 1,  1.0,  -1.0,  FEM_BC_DIRICHLET },
    { 2,  1.0,   1.0,  FEM_BC_DIRICHLET },
    { 3, -1.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         = NULL,
    .g         = NULL,
    .h         = NULL,
    .eta       = NULL,
    .u_exact   = NULL,
  };

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

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

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:

    21 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

  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 = NULL;
  spec->g = NULL;
  spec->h = NULL;
  spec->eta = NULL;
  spec->u_exact = NULL;

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

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

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:

    21 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
*/
{
  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         = NULL,
    .g         = NULL,
    .h         = NULL,
    .eta       = NULL,
    .u_exact   = NULL,
  };

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

  return &spec;
}