# include # include # include # include using namespace std; # include "tetrahedron_jaskowiec_rule.hpp" int main ( ); void tetrahedron_jaskowiec_rule_test01 ( int p ); void tetrahedron_jaskowiec_rule_test02 ( int p ); void tetrahedron_jaskowiec_rule_test03 ( int p_lo, int p_hi ); void tetrahedron_jaskowiec_rule_test04 ( ); void timestamp ( ); //****************************************************************************80 int main ( ) //****************************************************************************80 // // Purpose: // // tetrahedron_jaskowiec_rule_test() tests tetrahedron_jaskowiec_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 ( ); cout << "\n"; cout << "tetrahedron_jaskowiec_rule_test():\n"; cout << " C++ version\n"; cout << " Test tetrahedron_jaskowiec_rule().\n"; p = 5; tetrahedron_jaskowiec_rule_test01 ( p ); p = 5; tetrahedron_jaskowiec_rule_test02 ( p ); p_lo = 0; p_hi = 20; tetrahedron_jaskowiec_rule_test03 ( p_lo, p_hi ); tetrahedron_jaskowiec_rule_test04 ( ); // // Terminate. // cout << "\n"; cout << "tetrahedron_jaskowiec_rule_test()\n"; cout << " Normal end of execution.\n"; cout << "\n"; timestamp ( ); return 0; } //****************************************************************************80 void tetrahedron_jaskowiec_rule_test01 ( int p ) //****************************************************************************80 // // Purpose: // // tetrahedron_jaskowiec_rule_test01() prints the rule of precision P. // // Licensing: // // This code is distributed under the GNU GPL license. // // Modified: // // 26 May 2023 // // Author: // // John Burkardt. // { double *a; double *b; double *c; double *d; int i; int n; double *w; double w_sum; cout << "\n"; cout << "tetrahedron_jaskowiec_rule_test01():\n"; cout << " Quadrature rule for the tetrahedron,\n"; cout << " given in barycentric coordinates.\n"; cout << " Precision p = " << p << "\n"; // // Retrieve the rule. // n = rule_order ( p ); cout << "\n"; cout << " Number of nodes N = " << n << "\n"; a = new double[n]; b = new double[n]; c = new double[n]; d = new double[n]; w = new double[n]; tetrahedron_jaskowiec_rule ( p, n, a, b, c, d, w ); // // Print the rule. // cout << "\n"; cout << " I W A B C D\n"; cout << "\n"; for ( i = 0; i < n; i++ ) { cout << " " << setw(4) << i << " " << setw(10) << w[i] << " " << setw(10) << a[i] << " " << setw(10) << b[i] << " " << setw(10) << c[i] << " " << setw(10) << d[i] << "\n"; } // // Verify that weights sum to 1. // w_sum = 0.0; for ( i = 0; i < n; i++ ) { w_sum = w_sum + w[i]; } cout << "\n"; cout << " Weight Sum = " << w_sum << "\n"; // // Free memory. // delete [] a; delete [] b; delete [] c; delete [] d; delete [] w; return; } //****************************************************************************80 void tetrahedron_jaskowiec_rule_test02 ( int p ) //****************************************************************************80 // // Purpose: // // tetrahedron_jaskowiec_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; bool more; double quad; double quad_error; int n; int t; double *v; double *w; double *xyz; cout << "\n"; cout << "tetrahedron_jaskowiec_rule_test02():\n"; cout << " Test the precision of a quadrature rule for the unit tetrahedron.\n"; // // Retrieve the quadrature rule. // n = rule_order ( p ); cout << "\n"; cout << " Number of nodes N = " << n << "\n"; a = new double[n]; b = new double[n]; c = new double[n]; d = new double[n]; w = new double[n]; tetrahedron_jaskowiec_rule ( p, n, a, b, c, d, w ); // // Pack the x, y, z vectors as columns of an array. // xyz = new double[n*3]; 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; } cout << "\n"; cout << " Stated precision of rule = " << p << "\n"; cout << " Number of quadrature points = " << n << "\n"; cout << "\n"; cout << " Degree Maximum error\n"; cout << "\n"; v = new double[n]; 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; } } cout << " " << setw(2) << degree << " " << setw(24) << max_error << "\n"; } // // Free memory. // delete [] a; delete [] b; delete [] c; delete [] d; delete [] v; delete [] w; delete [] xyz; return; } //****************************************************************************80 void tetrahedron_jaskowiec_rule_test03 ( int p_lo, int p_hi ) //****************************************************************************80 // // Purpose: // // tetrahedron_jaskowiec_rule_test03() tests absolute and relative precision. // // Licensing: // // This code is distributed under the GNU GPL license. // // Modified: // // 11 July 2023 // // Author: // // John Burkardt. // // Input: // // integer p_lo, p_hi: the lowest and highest rules 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; bool more; int p; double quad; double quad_error; int n; int t; double *v; double *w; double *xyz; cout << "\n"; cout << "tetrahedron_jaskowiec_rule_test03():\n"; cout << " Test the precision of quadrature rules for the unit tetrahedron.\n"; cout << " Check rules of precision p = " << p_lo << " through " << p_hi << "\n"; cout << " for error in approximating integrals of monomials.\n"; cout << "\n"; cout << " maximum maximum\n"; cout << " p absolute relative\n"; cout << " error error\n"; cout << "\n"; for ( p = p_lo; p <= p_hi; p++ ) { n = rule_order ( p ); a = new double[n]; b = new double[n]; c = new double[n]; d = new double[n]; w = new double[n]; tetrahedron_jaskowiec_rule ( p, n, a, b, c, d, w ); // // Pack the x, y, z vectors as columns of an array. // xyz = new double[n*3]; 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 = new double[n]; 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; } } } cout << " " << setw(2) << p << " " << setw(24) << max_abs << " " << setw(24) << max_rel << "\n"; delete [] a; delete [] b; delete [] c; delete [] d; delete [] v; delete [] w; delete [] xyz; } return; } //****************************************************************************80 void tetrahedron_jaskowiec_rule_test04 ( ) //****************************************************************************80 // // Purpose: // // tetrahedron_jaskowiec_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: // // 08 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; cout << "\n"; cout << "tetrahedron_jaskowiec_rule_test04():\n"; cout << " Integrate 1/sqrt(r) over the reference tetrahedron.\n"; exact = 163.0 / 679.0; cout << " Exact integral value is " << exact << "\n"; volume = tetrahedron_volume ( tetra ); cout << " Volume of tetrahedron is " << volume << "\n"; cout << "\n"; cout << " P N Q |Q-Exact]\n"; cout << "\n"; for ( p = 0; p <= 20; p++ ) { n = rule_order ( p ); a = new double[n]; b = new double[n]; c = new double[n]; d = new double[n]; w = new double[n]; tetrahedron_jaskowiec_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 ); cout << " " << setw(2) << p << " " << setw(3) << n << " " << setw(20) << q << " " << setw(20) << e << "\n"; delete [] a; delete [] b; delete [] c; delete [] d; delete [] w; } return; } //****************************************************************************80 void timestamp ( ) //****************************************************************************80 // // Purpose: // // timestamp() prints the current YMDHMS date as a time stamp. // // Example: // // 31 May 2001 09:45:54 AM // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 19 March 2018 // // Author: // // John Burkardt // { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct std::tm *tm_ptr; std::time_t now; now = std::time ( NULL ); tm_ptr = std::localtime ( &now ); std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr ); std::cout << time_buffer << "\n"; return; # undef TIME_SIZE }