# include # include # include # include # include # include "ppc_array.h" # include "ppc_poisson.h" # include "ppc_xmalloc.h" int main ( int argc, char **argv ); struct problem_spec *annulus ( int n ); double annulus_f ( double x, double y ); void do_demo ( struct problem_spec *spec, double a, int d, char *gv_filename ); void draw_triangles_mono ( struct mesh *mesh, FILE *fp ); void draw_triangles_zhue ( struct mesh *mesh, FILE *fp ); void free_annulus ( struct problem_spec *spec ); void free_three_holes ( struct problem_spec *spec ); 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 ); double square_u_exact ( double x, double y ); struct problem_spec *three_holes ( int n ); double three_holes_f ( double x, double y ); 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_poisson_test() tests ppc_poisson(). Licensing: This code is distributed under the MIT license. Modified: 01 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 argc: the command line count. char **argv: the command line arguments. */ { double a; int d; char *endptr; int n; struct problem_spec *spec; timestamp ( ); printf ( "\n" ); printf ( "ppc_poisson_test ( ):\n" ); printf ( " C version.\n" ); printf ( " Test ppc_poisson().\n" ); /* Expect two arguments, d: the quadrature strength; a: the maximal triangle area. */ if ( argc != 3 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } d = strtol ( argv[1], &endptr, 10 ); if ( *endptr != '\0' || d <= 0 ) { 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; } printf ( "\n" ); printf ( " User requested quadrature strength %d\n", d ); printf ( " User requested maximal triangle area %g\n", a ); /* Annulus. By default, use 24 segments to draw the circular boundaries. */ if ( true ) { n = 24; spec = annulus ( n ); do_demo ( spec, a, d, "annulus.gv" ); free_annulus ( spec ); } /* Square. Exact solution is known. */ if ( true ) { spec = square ( ); do_demo ( spec, a, d, "square.gv" ); } /* Three holes. */ if ( true ) { n = 12; spec = three_holes ( n ); do_demo ( spec, a, d, "three_holes.gv" ); free_three_holes ( spec ); } /* Triangle with hole. */ if ( true ) { spec = triangle_with_hole ( ); do_demo ( spec, a, d, "triangle_with_hole.gv" ); } /* Terminate. */ printf ( "\n" ); printf ( "ppc_poisson_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 = 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", n ); printf ( " Inner radius = %g, outer radius = %g\n", a, b ); return spec; } /******************************************************************************/ double annulus_f ( double x, double y ) /******************************************************************************/ /* Purpose: annulus_f() defines the force term for the annulus region. 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 Input: double x, y: the evaluation point. Output: double value: the value of f(x,y). */ { double t; double value; t = atan2 ( y, x ); value = 20.0 * cos ( 3.0 * t ); return value; } /******************************************************************************/ void do_demo ( struct problem_spec *spec, double a, int d, char *gv_filename ) /******************************************************************************/ /* Purpose: do_demo() demonstrates the meshing of four regions. Licensing: This code is distributed under the MIT license. Modified: 30 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. int d: the quadrature rule strength. char *gv_filename: a name for the Geomview file to be created. */ { struct mesh *mesh; mesh = make_mesh ( spec, a ); printf ( " vertices = %d\n", mesh->nnodes ); printf ( " edges = %d\n", mesh->nedges ); printf ( " elements = %d\n", mesh->nelems ); poisson_solve ( spec, mesh, d ); /* If an exact solution was given, compute and report the errors. */ if ( spec->u_exact != NULL ) { struct errors errors = eval_errors ( spec, mesh, d ); printf ( "Linfty error = %g\n", errors.Linfty ); printf ( "L2norm error = %g\n", errors.L2norm ); printf ( "energy error = %g\n", errors.energy ); } plot_with_geomview_zhue ( mesh, gv_filename ); 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 Input: struct mesh *mesh: a pointer to the mesh structure. FILE *fp: the file to which the OFF image is written. */ { 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 Input: struct mesh *mesh: a pointer to the mesh structure. FILE *fp: the file to which the OFF image is written. */ { # 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 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 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 Input: struct mesh *mesh: a pointer to the mesh structure. Output: double *min, *max: the minimum and maximum of the scalar value for which a contour map (color or gray scaled) is to be created. */ { int i; for ( i = 0; i < mesh->nnodes; i++ ) { if ( i == 0 ) { *min = mesh->nodes[i].z; *max = mesh->nodes[i].z; } else { 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 REDUE, 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 Input: float s: Output: float value: */ { 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: 30 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 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: 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, 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 = square_u_exact, }; printf ( "\n" ); printf ( " This domain is a square.\n" ); return &spec; } /******************************************************************************/ double square_f ( double x, double y ) /******************************************************************************/ /* Purpose: square_f() defines the force term for the square region. 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 Input: double x, y: the evaluation point. Output: double value: the value of f(x,y). */ { double value; value = 32.0 * ( x * ( 1.0 - x ) + y * ( 1.0 - y ) ); return value; } /******************************************************************************/ double square_u_exact ( double x, double y ) /******************************************************************************/ /* Purpose: square_u_exact() defines the exact solution for the square region. 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 Input: double x, y: the evaluation point. Output: double value: the value of the exact solution at (x,y). */ { double value; value = 16.0 * x * ( 1.0 - x ) * y * ( 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: 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 = 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() defines the force term for the three_holes region. 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 Input: double x, y: the evaluation point. Output: double value: the value of f(x,y). */ { double value; value = 100.0 * ( x * x + y * y ); return value; } /******************************************************************************/ 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 = triangle_with_hole_f, .g = NULL, .h = NULL, .eta = NULL, .u_exact = NULL, }; printf ( "\n" ); printf ( " The domain is a triangle with a hole.\n" ); return &spec; } /******************************************************************************/ double triangle_with_hole_f ( double x, double y ) /******************************************************************************/ /* Purpose: triangle_with_hole_f() defines the force term for the triangle_with_hole region. 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 Input: double x, y: the evaluation point. Output: double value: the value of f(x,y). */ { double value; value = 30.0 * x * y * ( 2.0 - x ) * ( 2.0 - y ); return value; }