# include # include # include # include # include # include "tetrahedron_witherden_rule.h" int main ( ); void tetrahedron_witherden_rule_test01 ( int p ); void tetrahedron_witherden_rule_test02 ( int p ); void tetrahedron_witherden_rule_test03 ( int p_lo, int p_hi ); void tetrahedron_witherden_rule_test04 ( ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: tetrahedron_witherden_rule_test() tests tetrahedron_witherden_rule(). Licensing: This code is distributed under the GNU GPL license. Modified: 11 July 2023 Author: John Burkardt. */ { int p; int p_hi; int p_lo; timestamp ( ); printf ( "\n" ); printf ( "tetrahedron_witherden_rule_test():\n" ); printf ( " C version\n" ); printf ( " Test tetrahedron_witherden_rule().\n" ); p = 5; tetrahedron_witherden_rule_test01 ( p ); p = 5; tetrahedron_witherden_rule_test02 ( p ); p_lo = 0; p_hi = 10; tetrahedron_witherden_rule_test03 ( p_lo, p_hi ); tetrahedron_witherden_rule_test04 ( ); /* Terminate. */ printf ( "\n" ); printf ( "tetrahedron_witherden_rule_test()\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void tetrahedron_witherden_rule_test01 ( int p ) /******************************************************************************/ /* Purpose: tetrahedron_witherden_rule_test01() prints the rule of precision P. Licensing: This code is distributed under the GNU GPL license. Modified: 27 May 2023 Author: John Burkardt. */ { double *a; double *b; double *c; double *d; int i; int n; double *w; double w_sum; printf ( "\n" ); printf ( "tetrahedron_witherden_rule_test01():\n" ); printf ( " Quadrature rule for the tetrahedron,\n" ); printf ( " given in barycentric coordinates.\n" ); printf ( " Precision p = %d\n", p ); /* Retrieve the quadrature rule. */ n = rule_order ( p ); printf ( "\n" ); printf ( " Number of nodes N = %d\n", n ); w = ( double * ) malloc ( n * sizeof ( double ) ); a = ( double * ) malloc ( n * sizeof ( double ) ); b = ( double * ) malloc ( n * sizeof ( double ) ); c = ( double * ) malloc ( n * sizeof ( double ) ); d = ( double * ) malloc ( n * sizeof ( double ) ); tetrahedron_witherden_rule ( p, n, a, b, c, d, w ); /* Print the rule. */ printf ( "\n" ); printf ( " I W A B C D\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { printf ( " %4d %10.6g %10.6g %10.6g %10.6g %10.6g\n", i, w[i], a[i], b[i], c[i], d[i] ); } /* Verify weights sum to 1. */ w_sum = 0.0; for ( i = 0; i < n; i++ ) { w_sum = w_sum + w[i]; } printf ( "\n" ); printf ( " Weight Sum = %g\n", w_sum ); /* Free memory. */ free ( a ); free ( b ); free ( c ); free ( d ); free ( w ); return; } /******************************************************************************/ void tetrahedron_witherden_rule_test02 ( int p ) /******************************************************************************/ /* Purpose: tetrahedron_witherden_rule_test02() tests a rule of precision P. Licensing: This code is distributed under the GNU GPL license. Modified: 23 May 2023 Author: John Burkardt. */ { double *a; double *b; double *c; double *d; int degree; const int dim_num = 3; double exact; int expon[3]; int h; int i; int ij; double max_error; int more; double quad; double quad_error; int n; int t; double *v; double *w; double *xyz; printf ( "\n" ); printf ( "tetrahedron_witherden_rule_test02():\n" ); printf ( " Test the precision of a quadrature rule for the unit tetrahedron.\n" ); /* Retrieve the quadrature rule. */ n = rule_order ( p ); printf ( "\n" ); printf ( " Number of nodes N = %d\n", n ); w = ( double * ) malloc ( n * sizeof ( double ) ); a = ( double * ) malloc ( n * sizeof ( double ) ); b = ( double * ) malloc ( n * sizeof ( double ) ); c = ( double * ) malloc ( n * sizeof ( double ) ); d = ( double * ) malloc ( n * sizeof ( double ) ); tetrahedron_witherden_rule ( p, n, a, b, c, d, w ); /* Pack the x, y, z vectors as columns of an array. */ xyz = ( double * ) malloc ( n * 3 * sizeof ( double ) ); ij = 0; for ( i = 0; i < n; i++ ) { xyz[ij] = a[i]; ij = ij + 1; xyz[ij] = b[i]; ij = ij + 1; xyz[ij] = c[i]; ij = ij + 1; } printf ( "\n" ); printf ( " Stated precision of rule = %d\n", p ); printf ( " Number of quadrature points = %d\n", n ); printf ( "\n" ); printf ( " Degree Maximum error\n" ); printf ( "\n" ); v = ( double * ) malloc ( n * sizeof ( double ) ); for ( degree = 0; degree <= p + 2; degree++ ) { more = 0; h = 0; t = 0; max_error = 0.0; while ( true ) { comp_next ( degree, dim_num, expon, &more, &h, &t ); v = monomial_value ( dim_num, n, expon, xyz ); quad = 0.0; for ( i = 0; i < n; i++ ) { quad = quad + w[i] * v[i]; } quad = tetrahedron_unit_volume ( ) * quad; exact = tetrahedron_unit_monomial_integral ( expon ); quad_error = fabs ( quad - exact ); max_error = fmax ( max_error, quad_error ); if ( ! more ) { break; } } printf ( " %2d %24.16f\n", degree, max_error ); } free ( a ); free ( b ); free ( c ); free ( d ); free ( v ); free ( w ); free ( xyz ); return; } /******************************************************************************/ void tetrahedron_witherden_rule_test03 ( int p_lo, int p_hi ) /******************************************************************************/ /* Purpose: tetrahedron_witherden_rule_test03() tests absolute and relative error. Licensing: This code is distributed under the GNU GPL license. Modified: 11 July 2023 Author: John Burkardt. Input: int p_lo, p_hi: the lowest and highest precisions to check. */ { double *a; double *b; double *c; double *d; int degree; const int dim_num = 3; double exact; int expon[3]; int h; int i; int ij; double max_abs; double max_rel; int more; int p; double quad; double quad_error; int n; int t; double *v; double *w; double *xyz; printf ( "\n" ); printf ( "tetrahedron_witherden_rule_test03():\n" ); printf ( " Test the precision of quadrature rules for the unit tetrahedron.\n" ); printf ( " Check rules of precision p = %d through %d\n", p_lo, p_hi ); printf ( " for error in approximating integrals of monomials.\n" ); printf ( "\n" ); printf ( " maximum maximum\n" ); printf ( " p absolute relative\n" ); printf ( " error error\n" ); printf ( "\n" ); for ( p = p_lo; p <= p_hi; p++ ) { n = rule_order ( p ); w = ( double * ) malloc ( n * sizeof ( double ) ); a = ( double * ) malloc ( n * sizeof ( double ) ); b = ( double * ) malloc ( n * sizeof ( double ) ); c = ( double * ) malloc ( n * sizeof ( double ) ); d = ( double * ) malloc ( n * sizeof ( double ) ); tetrahedron_witherden_rule ( p, n, a, b, c, d, w ); /* Pack the x, y, z vectors as columns of an array. */ xyz = ( double * ) malloc ( n * 3 * sizeof ( double ) ); ij = 0; for ( i = 0; i < n; i++ ) { xyz[ij] = a[i]; ij = ij + 1; xyz[ij] = b[i]; ij = ij + 1; xyz[ij] = c[i]; ij = ij + 1; } v = ( double * ) malloc ( n * sizeof ( double ) ); max_abs = 0.0; max_rel = 0.0; for ( degree = 0; degree <= p; degree++ ) { more = 0; h = 0; t = 0; while ( true ) { comp_next ( degree, dim_num, expon, &more, &h, &t ); v = monomial_value ( dim_num, n, expon, xyz ); quad = 0.0; for ( i = 0; i < n; i++ ) { quad = quad + w[i] * v[i]; } quad = tetrahedron_unit_volume ( ) * quad; exact = tetrahedron_unit_monomial_integral ( expon ); quad_error = fabs ( quad - exact ); max_abs = fmax ( max_abs, quad_error ); max_rel = fmax ( max_rel, quad_error / fabs ( exact ) ); if ( ! more ) { break; } } } printf ( " %2d %24.16g %24.16g\n", p, max_abs, max_rel ); free ( a ); free ( b ); free ( c ); free ( d ); free ( v ); free ( w ); free ( xyz ); } return; } /******************************************************************************/ void tetrahedron_witherden_rule_test04 ( ) /******************************************************************************/ /* Purpose: tetrahedron_witherden_rule_test04() integrates 1/sqrt(r). Discussion: The integral is taken over the reference tetrahedron with vertices (0,0,0), (1,0,0), (0,1,0), (0,0,1). The test integral is taken from Jaskowiec and Sukumar. Licensing: This code is distributed under the GNU GPL license. Modified: 07 May 2023 Author: John Burkardt. Reference: Jan Jaskowiec, Natarajan Sukumar, High order symmetric cubature rules for tetrahedra and pyramids, International Journal of Numerical Methods in Engineering, Volume 122, Number 1, pages 148-171, 24 August 2020. */ { double *a; double *b; double *c; double *d; double e; double exact; double f; int i; int n; int p; double q; double r; double tetra[4*3] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 }; double volume; double *w; printf ( "\n" ); printf ( "tetrahedron_witherden_rule_test04():\n" ); printf ( " Integrate 1/sqrt(r) over the reference tetrahedron.\n" ); printf ( " Witherden rule #9 fails because a quadrature point is\n" ); printf ( " very near the singularity at the origin.\n" ); exact = 163.0 / 679.0; printf ( " Exact integral value is %g.\n", exact ); volume = tetrahedron_volume ( tetra ); printf ( " Volume of tetrahedron is %g.\n", volume ); printf ( "\n" ); printf ( " P N Q |Q-Exact]\n" ); printf ( "\n" ); for ( p = 0; p <= 10; p++ ) { n = rule_order ( p ); a = ( double * ) malloc ( n * sizeof ( double ) ); b = ( double * ) malloc ( n * sizeof ( double ) ); c = ( double * ) malloc ( n * sizeof ( double ) ); d = ( double * ) malloc ( n * sizeof ( double ) ); w = ( double * ) malloc ( n * sizeof ( double ) ); tetrahedron_witherden_rule ( p, n, a, b, c, d, w ); q = 0.0; for ( i = 0; i < n; i++ ) { r = sqrt ( a[i] * a[i] + b[i] * b[i] + c[i] * c[i] ); f = 1.0 / sqrt ( r ); q = q + w[i] * f; } q = volume * q; e = fabs ( q - exact ); printf ( " %2d %3d %20.16g %20.16g\n", p, n, q, e ); free ( a ); free ( b ); free ( c ); free ( d ); free ( w ); } 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 }