# include # include # include # include # include # include "ppc_array.h" # include "ppc_bvp.h" # include "ppc_xmalloc.h" int main ( int argc, char **argv ); void do_demo ( struct problem_spec *spec, double a, int twb_d, int gauss_d, char *gv_filename ); void draw_triangles_mono ( struct mesh *mesh, FILE *fp ); void draw_triangles_zhue ( struct mesh *mesh, FILE *fp ); 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 *square1 ( ); double square1_f ( double x, double y ); double square1_g ( double x, double y ); double square1_h ( double x, double y ); double square1_eta ( double x, double y ); double square1_u_exact ( double x, double y ); struct problem_spec *square2 ( ); double square2_f ( double x, double y ); double square2_g ( double x, double y ); double square2_h ( double x, double y ); double square2_eta ( double x, double y ); double square2_u_exact ( double x, double y ); struct problem_spec *square3 ( ); double square3_f ( double x, double y ); double square3_g ( double x, double y ); double square3_h ( double x, double y ); double square3_eta ( double x, double y ); double square3_u_exact ( double x, double y ); struct problem_spec *square4 ( ); double square4_f ( double x, double y ); double square4_g ( double x, double y ); double square4_h ( double x, double y ); double square4_eta ( double x, double y ); void timestamp ( ); /******************************************************************************/ int main ( int argc, char **argv ) /******************************************************************************/ /* Purpose: ppc_bvp_test() tests ppc_bvp(). Licensing: This code is distributed under the MIT license. Modified: 14 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; char *endptr; int gauss_d; struct problem_spec *spec; int twb_d; timestamp ( ); printf ( "\n" ); printf ( "ppc_bvp_test ( ):\n" ); printf ( " C version.\n" ); printf ( " Test ppc_bvp().\n" ); /* Expect three arguments, twb_d: the TWB quadrature strength; a: the maximal triangle area. */ if ( argc != 4 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } twb_d = strtol ( argv[1], &endptr, 10 ); if ( *endptr != '\0' || twb_d <= 0 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } gauss_d = strtol ( argv[2], &endptr, 10 ); if ( *endptr != '\0' || gauss_d <= 0 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } a = strtod ( argv[3], &endptr ); if ( *endptr != '\0' || a <= 0.0 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } printf ( "\n" ); printf ( "User input:\n" ); printf ( " TWB quadrature strength %d\n", twb_d ); printf ( " Gauss quadrature strength %d\n", gauss_d ); printf ( " maximal triangle area %g\n", a ); /* Square1. Exact solution is known. */ if ( true ) { spec = square1 ( ); do_demo ( spec, a, twb_d, gauss_d, "square1.gv" ); } /* Square2. Exact solution is known. */ if ( true ) { spec = square2 ( ); do_demo ( spec, a, twb_d, gauss_d, "square2.gv" ); } /* Square3. Exact solution is known. */ if ( true ) { spec = square3 ( ); do_demo ( spec, a, twb_d, gauss_d, "square3.gv" ); } /* Square4. Exact solution not known. */ if ( true ) { spec = square4 ( ); do_demo ( spec, a, twb_d, gauss_d, "square4.gv" ); } /* Terminate. */ printf ( "\n" ); printf ( "ppc_bvp_test ( ):\n" ); printf ( " Normal end of execution.\n" ); timestamp ( ); return EXIT_SUCCESS; } /******************************************************************************/ void do_demo ( struct problem_spec *spec, double a, int twb_d, int gauss_d, char *gv_filename ) /******************************************************************************/ /* Purpose: do_demo() sets up and solves a boundary value problem. Licensing: This code is distributed under the MIT license. Modified: 13 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. int twb_d: the TWB quadrature rule strength. int gauss_d: the Gauss 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 ); bvp_solve ( spec, mesh, twb_d, gauss_d ); /* If an exact solution was given, compute and report the errors. */ if ( spec->u_exact != NULL ) { struct errors errors = eval_errors ( spec, mesh, twb_d ); printf ( " Linfty error = %g\n", errors.Linfty ); printf ( " L2norm error = %g\n", errors.L2norm ); } 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 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 twb_d gauss_d a\n", progname ); printf ( " twb_d = triangle quadrature rule strength.\n" ); printf ( " gauss_d = interval quadrature rule strength.\n" ); printf ( " a = maximal triangle area.\n" ); return; } /******************************************************************************/ struct problem_spec *square1 ( ) /******************************************************************************/ /* Purpose: square1() defines the square1 region. Licensing: This code is distributed under the MIT license. Modified: 13 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 spec = { .points = points, .npoints = sizeof ( points ) / sizeof ( points[0] ), .segments = segments, .nsegments = sizeof ( segments ) / sizeof (segments[0] ), .holes = NULL, .nholes = 0, .f = square1_f, .g = square1_g, .h = square1_h, .eta = square1_eta, .u_exact = square1_u_exact, }; printf ( "\n" ); printf ( " Domain is square1.\n" ); return &spec; } /******************************************************************************/ double square1_f ( double x, double y ) /******************************************************************************/ /* Purpose: square1_f() defines the force term for the square1 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of f(x,y). */ { double value; value = - 32.0 * x * ( 1.0 - x ) + 8.0 * ( 1.0 - y ) * ( 1.0 - 4.0 * y ); return value; } /******************************************************************************/ double square1_g ( double x, double y ) /******************************************************************************/ /* Purpose: square1_g() defines the Dirichlet boundary term for the square1 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of g(x,y). */ { double value; if ( y == 0.0 ) { value = 4.0 * x * ( 1.0 - x ); } else { value = 0.0; } return value; } /******************************************************************************/ double square1_h ( double x, double y ) /******************************************************************************/ /* Purpose: square1_h() defines the Neumann boundary term for the square1 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of h(x,y). */ { double value; value = 0.0; return value; } /******************************************************************************/ double square1_eta ( double x, double y ) /******************************************************************************/ /* Purpose: square1_eta() defines the eta function for the square1 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of eta(x,y). */ { double value; value = 1.0; return value; } /******************************************************************************/ double square1_u_exact ( double x, double y ) /******************************************************************************/ /* Purpose: square1_u_exact() defines the exact solution for the square1 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of the exact solution at (x,y). */ { double value; value = 4.0 * x * ( 1.0 - x ) * ( 1.0 - y ) * ( 1.0 - 4.0 * y ); return value; } /******************************************************************************/ struct problem_spec *square2 ( ) /******************************************************************************/ /* Purpose: square2() defines the square2 region. Licensing: This code is distributed under the MIT license. Modified: 12 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_NEUMANN }, { 1, 1, 2, FEM_BC_DIRICHLET }, { 2, 2, 3, FEM_BC_DIRICHLET }, { 3, 3, 0, FEM_BC_DIRICHLET }, }; static struct problem_spec spec = { .points = points, .npoints = sizeof ( points ) / sizeof ( points[0] ), .segments = segments, .nsegments = sizeof ( segments ) / sizeof (segments[0] ), .holes = NULL, .nholes = 0, .f = square2_f, .g = square2_g, .h = square2_h, .eta = square2_eta, .u_exact = square2_u_exact, }; printf ( "\n" ); printf ( " Domain is square2.\n" ); return &spec; } /******************************************************************************/ double square2_f ( double x, double y ) /******************************************************************************/ /* Purpose: square2_f() defines the force term for the square2 region. Licensing: This code is distributed under the MIT license. Modified: 12 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: double x, y: the evaluation point. Output: double value: the value of f(x,y). */ { double value; value = - 32.0 * x * ( 1.0 - x ) + 8.0 * ( 1.0 - y ) * ( 1.0 - 4.0 * y ); return value; } /******************************************************************************/ double square2_g ( double x, double y ) /******************************************************************************/ /* Purpose: square2_g() defines the Dirichlet boundary term for the square2 region. Licensing: This code is distributed under the MIT license. Modified: 12 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: double x, y: the evaluation point. Output: double value: the value of g(x,y). */ { double value; value = 0.0; return value; } /******************************************************************************/ double square2_h ( double x, double y ) /******************************************************************************/ /* Purpose: square2_h() defines the Neumann boundary term for the square2 region. Licensing: This code is distributed under the MIT license. Modified: 12 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: double x, y: the evaluation point. Output: double value: the value of h(x,y). */ { double value; value = 20.0 * x * ( 1.0 - x ); return value; } /******************************************************************************/ double square2_eta ( double x, double y ) /******************************************************************************/ /* Purpose: square2_eta() defines the eta function for the square2 region. Licensing: This code is distributed under the MIT license. Modified: 12 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: double x, y: the evaluation point. Output: double value: the value of eta(x,y). */ { double value; value = 1.0; return value; } /******************************************************************************/ double square2_u_exact ( double x, double y ) /******************************************************************************/ /* Purpose: square2_u_exact() defines the exact solution for the square2 region. Licensing: This code is distributed under the MIT license. Modified: 12 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: double x, y: the evaluation point. Output: double value: the value of the exact solution at (x,y). */ { double value; value = 4.0 * x * ( 1.0 - x ) * ( 1.0 - y ) * ( 1.0 - 4.0 * y ); return value; } /******************************************************************************/ struct problem_spec *square3 ( ) /******************************************************************************/ /* Purpose: square3() defines the square3 region. Licensing: This code is distributed under the MIT license. Modified: 13 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 spec = { .points = points, .npoints = sizeof ( points ) / sizeof ( points[0] ), .segments = segments, .nsegments = sizeof ( segments ) / sizeof (segments[0] ), .holes = NULL, .nholes = 0, .f = square3_f, .g = square3_g, .h = square3_h, .eta = square3_eta, .u_exact = square3_u_exact, }; printf ( "\n" ); printf ( " Domain is square3.\n" ); return &spec; } /******************************************************************************/ double square3_f ( double x, double y ) /******************************************************************************/ /* Purpose: square3_f() defines the force term for the square3 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of f(x,y). */ { double value; value = 16.0 * pow ( x - 0.125, 2 ) + 16.0 * pow ( y - 0.125, 2 ) + 15.0 / 2.0; return value; } /******************************************************************************/ double square3_g ( double x, double y ) /******************************************************************************/ /* Purpose: square3_g() defines the Dirichlet boundary term for the square3 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of g(x,y). */ { double value; value = 1.0 - 2.0 * pow ( x - 0.5, 2 ) - 2.0 * pow ( y - 0.5, 2 ); return value; } /******************************************************************************/ double square3_h ( double x, double y ) /******************************************************************************/ /* Purpose: square3_h() defines the Neumann boundary term for the square3 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of h(x,y). */ { double value; value = 0.0; return value; } /******************************************************************************/ double square3_eta ( double x, double y ) /******************************************************************************/ /* Purpose: square3_eta() defines the eta function for the square3 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of eta(x,y). */ { double value; value = 1.0 + x * x + y * y; return value; } /******************************************************************************/ double square3_u_exact ( double x, double y ) /******************************************************************************/ /* Purpose: square3_u_exact() defines the exact solution for the square3 region. Licensing: This code is distributed under the MIT license. Modified: 13 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: double x, y: the evaluation point. Output: double value: the value of the exact solution at (x,y). */ { double value; value = 1.0 - 2.0 * pow ( x - 0.5, 2 ) - 2.0 * pow ( y - 0.5, 2 ); return value; } /******************************************************************************/ struct problem_spec *square4 ( ) /******************************************************************************/ /* Purpose: square4() defines the square4 region. Licensing: This code is distributed under the MIT license. Modified: 14 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, 0.5, 0.0, FEM_BC_DIRICHLET }, { 2, 1.0, 0.0, FEM_BC_NEUMANN }, { 3, 1.0, 1.0, FEM_BC_DIRICHLET }, { 4, 0.0, 1.0, FEM_BC_DIRICHLET }, }; static struct problem_spec_segment segments[] = { { 0, 0, 1, FEM_BC_DIRICHLET }, { 1, 1, 2, FEM_BC_NEUMANN }, { 2, 2, 3, FEM_BC_NEUMANN }, { 3, 3, 4, FEM_BC_DIRICHLET }, { 4, 4, 0, FEM_BC_NEUMANN }, }; static struct problem_spec spec = { .points = points, .npoints = sizeof ( points ) / sizeof ( points[0] ), .segments = segments, .nsegments = sizeof ( segments ) / sizeof (segments[0] ), .holes = NULL, .nholes = 0, .f = square4_f, .g = square4_g, .h = square4_h, .eta = square4_eta, .u_exact = NULL, }; printf ( "\n" ); printf ( " Domain is square4.\n" ); return &spec; } /******************************************************************************/ double square4_f ( double x, double y ) /******************************************************************************/ /* Purpose: square4_f() defines the force term for the square4 region. Discussion: I could not get geomview to display the solution when I set f(x,y)=0. Setting f(x,y)=1, I got a plausible result. Licensing: This code is distributed under the MIT license. Modified: 14 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: double x, y: the evaluation point. Output: double value: the value of f(x,y). */ { double value; value = 1.0; return value; } /******************************************************************************/ double square4_g ( double x, double y ) /******************************************************************************/ /* Purpose: square4_g() defines the Dirichlet boundary term for the square4 region. Licensing: This code is distributed under the MIT license. Modified: 14 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: double x, y: the evaluation point. Output: double value: the value of g(x,y). */ { double value; if ( y == 0.0 ) { value = 1.0; } else { value = 0.0; } return value; } /******************************************************************************/ double square4_h ( double x, double y ) /******************************************************************************/ /* Purpose: square4_h() defines the Neumann boundary term for the square4 region. Licensing: This code is distributed under the MIT license. Modified: 14 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: double x, y: the evaluation point. Output: double value: the value of h(x,y). */ { double value; value = 0.0; return value; } /******************************************************************************/ double square4_eta ( double x, double y ) /******************************************************************************/ /* Purpose: square4_eta() defines the eta function for the square4 region. Licensing: This code is distributed under the MIT license. Modified: 14 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: double x, y: the evaluation point. Output: double value: the value of eta(x,y). */ { double value; value = 1.0; 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 }