# include # include # include # include "hexahedron_witherden_rule.h" /******************************************************************************/ void comp_next ( int n, int k, int a[], int *more, int *h, int *t ) /******************************************************************************/ /* Purpose: comp_next() computes the compositions of the integer N into K parts. Discussion: A composition of the integer N into K parts is an ordered sequence of K nonnegative integers which sum to N. The compositions (1,2,1) and (1,1,2) are considered to be distinct. The routine computes one composition on each call until there are no more. For instance, one composition of 6 into 3 parts is 3+2+1, another would be 6+0+0. On the first call to this routine, set MORE = FALSE. The routine will compute the first element in the sequence of compositions, and return it, as well as setting MORE = TRUE. If more compositions are desired, call again, and again. Each time, the routine will return with a new composition. However, when the LAST composition in the sequence is computed and returned, the routine will reset MORE to FALSE, signaling that the end of the sequence has been reached. This routine originally used a STATIC statement to maintain the variables H and T. I have decided (based on an wasting an entire morning trying to track down a problem) that it is safer to pass these variables as arguments, even though the user should never alter them. This allows this routine to safely shuffle between several ongoing calculations. There are 28 compositions of 6 into three parts. This routine will produce those compositions in the following order: I A - --------- 1 6 0 0 2 5 1 0 3 4 2 0 4 3 3 0 5 2 4 0 6 1 5 0 7 0 6 0 8 5 0 1 9 4 1 1 10 3 2 1 11 2 3 1 12 1 4 1 13 0 5 1 14 4 0 2 15 3 1 2 16 2 2 2 17 1 3 2 18 0 4 2 19 3 0 3 20 2 1 3 21 1 2 3 22 0 3 3 23 2 0 4 24 1 1 4 25 0 2 4 26 1 0 5 27 0 1 5 28 0 0 6 Licensing: This code is distributed under the MIT license. Modified: 02 July 2008 Author: Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. C version by John Burkardt. Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms for Computers and Calculators, Second Edition, Academic Press, 1978, ISBN: 0-12-519260-6, LC: QA164.N54. Parameters: Input, int N, the integer whose compositions are desired. Input, int K, the number of parts in the composition. Input/output, int A[K], the parts of the composition. Input/output, int *MORE. Set MORE = FALSE on first call. It will be reset to TRUE on return with a new composition. Each new call returns another composition until MORE is set to FALSE when the last composition has been computed and returned. Input/output, int *H, *T, two internal parameters needed for the computation. The user should allocate space for these in the calling program, include them in the calling sequence, but never alter them! */ { int i; if ( !( *more ) ) { *t = n; *h = 0; a[0] = n; for ( i = 1; i < k; i++ ) { a[i] = 0; } } else { if ( 1 < *t ) { *h = 0; } *h = *h + 1; *t = a[*h-1]; a[*h-1] = 0; a[0] = *t - 1; a[*h] = a[*h] + 1; } *more = ( a[k-1] != n ); return; } /******************************************************************************/ double hexahedron_unit_monomial_integral ( int expon[3] ) /******************************************************************************/ /* Purpose: hexahedron_unit_monomial_integral(): monomial integral in a unit hexahedron. Discussion: This function returns the value of the integral of X^ALPHA Y^BETA Z^GAMMA over the unit hexahedron. The unit hexahedron has vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1). Licensing: This code is distributed under the MIT license. Modified: 20 May 2023 Author: John Burkardt Input: int EXPON[3], the exponents. Output: double hexahedron_unit_monomial_integral, the integral of the monomial over the hexahedron. */ { double value; value = 1.0 / ( expon[0] + 1 ) / ( expon[1] + 1 ) / ( expon[2] + 1 ); return value; } /******************************************************************************/ double hexahedron_unit_volume ( ) /******************************************************************************/ /* Purpose: hexahedron_unit_volume() returns the volume of a unit hexahedron. Discussion: The unit hexahedron has vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1). Licensing: This code is distributed under the MIT license. Modified: 02 May 2023 Author: John Burkardt Output: real value: the volume. */ { double value; value = 1.0; return value; } /******************************************************************************/ void hexahedron_witherden_rule ( int p, int n, double x[], double y[], double z[], double w[] ) /******************************************************************************/ /* Purpose: hexahedron_witherden_rule() returns a hexahedron quadrature rule of given precision. Discussion: The unit hexahedron has vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1). Licensing: This code is distributed under the MIT license. Modified: 22 May 2023 Author: John Burkardt Reference: Freddie Witherden, Peter Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers and Mathematics with Applications, Volume 69, pages 1232-1241, 2015. Input: int p: the precision, 0 <= p <= 11. int n: the order of the rule. Output: double x(n), y(n), z(n): the coordinates of quadrature points. double w(n): the quadrature weights. */ { if ( p < 0 ) { printf ( "\n" ); printf ( "hexahedron_witherden_rule(): Fatal error!\n" ); printf ( " Input p < 0.\n" ); exit ( 1 ); } if ( 11 < p ) { printf ( "\n" ); printf ( "hexahedron_witherden_rule(): Fatal error!\n" ); printf ( " Input 11 < p.\n" ); exit ( 1 ); } if ( p <= 1 ) { rule01 ( n, x, y, z, w ); } else if ( p <= 3 ) { rule03 ( n, x, y, z, w ); } else if ( p <= 5 ) { rule05 ( n, x, y, z, w ); } else if ( p <= 7 ) { rule07 ( n, x, y, z, w ); } else if ( p <= 9 ) { rule09 ( n, x, y, z, w ); } else if ( p <= 11 ) { rule11 ( n, x, y, z, w ); } return; } /******************************************************************************/ double *monomial_value ( int m, int n, int e[], double x[] ) /******************************************************************************/ /* Purpose: monomial_value() evaluates a monomial. Discussion: This routine evaluates a monomial of the form product ( 1 <= i <= m ) x(i)^e(i) The combination 0.0^0 is encountered is treated as 1.0. Licensing: This code is distributed under the MIT license. Modified: 17 August 2014 Author: John Burkardt Input: int M, the spatial dimension. int N, the number of evaluation points. int E[M], the exponents. double X[M*N], the point coordinates. Output: double MONOMIAL_VALUE[N], the monomial values. */ { int i; int j; double *v; v = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { v[j] = 1.0; } for ( i = 0; i < m; i++ ) { if ( 0 != e[i] ) { for ( j = 0; j < n; j++ ) { v[j] = v[j] * pow ( x[i+j*m], e[i] ); } } } return v; } /******************************************************************************/ void r8vec_copy ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: r8vec_copy() copies an R8VEC. Licensing: This code is distributed under the MIT license. Modified: 03 July 2005 Author: John Burkardt Input: int N, the number of entries in the vectors. double A1[N], the vector to be copied. Output: double A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ int rule_order ( int p ) /******************************************************************************/ /* Purpose: rule_order() returns the order of a hexahedron quadrature rule of given precision. Discussion: The "order" of a quadrature rule is the number of points involved. Licensing: This code is distributed under the MIT license. Modified: 02 May 2023 Author: John Burkardt Reference: Freddie Witherden, Peter Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers and Mathematics with Applications, Volume 69, pages 1232-1241, 2015. Input: int p: the precision, 0 <= p <= 11. Output: int rule_order: the order of the rule. */ { int order; static int order_save[] = { 1, 1, 6, 6, 14, 14, 34, 34, 58, 58, 90, 90 }; if ( p < 0 ) { printf ( "\n" ); printf ( "rule_order(): Fatal error!\n" ); printf ( " Input p < 0.\n" ); exit ( 1 ); } if ( 11 < p ) { printf ( "\n" ); printf ( "rule_order(): Fatal error!\n" ); printf ( " Input 11 < p.\n" ); exit ( 1 ); } order = order_save[p]; return order; } void rule01 ( int n, double x[], double y[], double z[], double w[] ) /******************************************************************************/ /* Purpose: rule01() returns the hexahedron rule of precision 1. Discussion: We suppose we are given a hexahedron H with vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), We call a rule with n points, returning coordinates x, y, z, and weights w. Then the integral I of f(x,y,z) over H is approximated by Q as follows: Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) Licensing: This code is distributed under the MIT license. Modified: 02 May 2023 Author: John Burkardt Reference: Freddie Witherden, Peter Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers and Mathematics with Applications, Volume 69, pages 1232-1241, 2015. Input: int n: the number of quadrature points for this rule. Output: double x(n), y(n), x(n): the coordinates of quadrature points. double w[n]: the quadrature weights, which sum to 1. */ { static double x_save[] = { 0.5000000000000000 }; static double y_save[] = { 0.5000000000000000 }; static double z_save[] = { 0.5000000000000000 }; double w_save[] = { 1.0000000000000000 }; r8vec_copy ( n, x_save, x ); r8vec_copy ( n, y_save, y ); r8vec_copy ( n, z_save, z ); r8vec_copy ( n, w_save, w ); return; } void rule03 ( int n, double x[], double y[], double z[], double w[] ) /******************************************************************************/ /* Purpose: rule03() returns the hexahedron rule of precision 3. Discussion: We suppose we are given a hexahedron H with vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), We call a rule with n points, returning coordinates x, y, z, and weights w. Then the integral I of f(x,y,z) over H is approximated by Q as follows: Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) Licensing: This code is distributed under the MIT license. Modified: 02 May 2023 Author: John Burkardt Reference: Freddie Witherden, Peter Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers and Mathematics with Applications, Volume 69, pages 1232-1241, 2015. Input: int n: the number of quadrature points for this rule. Output: double x(n), y(n), x(n): the coordinates of quadrature points. double w[n]: the quadrature weights, which sum to 1. */ { static double x_save[] = { 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 1.0000000000000000, 0.5000000000000000 }; static double y_save[] = { 0.5000000000000000, 0.5000000000000000, 1.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.0000000000000000 }; static double z_save[] = { 0.5000000000000000, 1.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000 }; double w_save[] = { 0.1666666666666667, 0.1666666666666667, 0.1666666666666667, 0.1666666666666667, 0.1666666666666667, 0.1666666666666667 }; r8vec_copy ( n, x_save, x ); r8vec_copy ( n, y_save, y ); r8vec_copy ( n, z_save, z ); r8vec_copy ( n, w_save, w ); return; } void rule05 ( int n, double x[], double y[], double z[], double w[] ) /******************************************************************************/ /* Purpose: rule05() returns the hexahedron rule of precision 5. Discussion: We suppose we are given a hexahedron H with vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), We call a rule with n points, returning coordinates x, y, z, and weights w. Then the integral I of f(x,y,z) over H is approximated by Q as follows: Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) Licensing: This code is distributed under the MIT license. Modified: 02 May 2023 Author: John Burkardt Reference: Freddie Witherden, Peter Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers and Mathematics with Applications, Volume 69, pages 1232-1241, 2015. Input: int n: the number of quadrature points for this rule. Output: double x(n), y(n), x(n): the coordinates of quadrature points. double w[n]: the quadrature weights, which sum to 1. */ { static double x_save[] = { 0.1020887871228893, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.8979112128771107, 0.5000000000000000, 0.8793934553196641, 0.1206065446803359, 0.1206065446803359, 0.1206065446803359, 0.1206065446803359, 0.8793934553196641, 0.8793934553196641, 0.8793934553196641 }; static double y_save[] = { 0.5000000000000000, 0.5000000000000000, 0.8979112128771107, 0.5000000000000000, 0.5000000000000000, 0.1020887871228893, 0.1206065446803359, 0.8793934553196641, 0.8793934553196641, 0.1206065446803359, 0.1206065446803359, 0.8793934553196641, 0.8793934553196641, 0.1206065446803359 }; static double z_save[] = { 0.5000000000000000, 0.8979112128771107, 0.5000000000000000, 0.1020887871228893, 0.5000000000000000, 0.5000000000000000, 0.1206065446803359, 0.8793934553196641, 0.1206065446803359, 0.1206065446803359, 0.8793934553196641, 0.1206065446803359, 0.8793934553196641, 0.8793934553196641 }; double w_save[] = { 0.1108033240997230, 0.1108033240997230, 0.1108033240997230, 0.1108033240997230, 0.1108033240997230, 0.1108033240997230, 0.0418975069252078, 0.0418975069252078, 0.0418975069252078, 0.0418975069252078, 0.0418975069252078, 0.0418975069252078, 0.0418975069252078, 0.0418975069252078 }; r8vec_copy ( n, x_save, x ); r8vec_copy ( n, y_save, y ); r8vec_copy ( n, z_save, z ); r8vec_copy ( n, w_save, w ); return; } void rule07 ( int n, double x[], double y[], double z[], double w[] ) /******************************************************************************/ /* Purpose: rule07() returns the hexahedron rule of precision 7. Discussion: We suppose we are given a hexahedron H with vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), We call a rule with n points, returning coordinates x, y, z, and weights w. Then the integral I of f(x,y,z) over H is approximated by Q as follows: Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) Licensing: This code is distributed under the MIT license. Modified: 02 May 2023 Author: John Burkardt Reference: Freddie Witherden, Peter Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers and Mathematics with Applications, Volume 69, pages 1232-1241, 2015. Input: int n: the number of quadrature points for this rule. Output: double x(n), y(n), x(n): the coordinates of quadrature points. double w[n]: the quadrature weights, which sum to 1. */ { static double x_save[] = { 0.0059569194056617, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.9940430805943383, 0.5000000000000000, 0.7039775836791597, 0.2960224163208403, 0.2960224163208403, 0.2960224163208403, 0.2960224163208403, 0.7039775836791597, 0.7039775836791597, 0.7039775836791597, 0.8905514105020593, 0.1094485894979407, 0.1094485894979407, 0.1094485894979407, 0.1094485894979407, 0.8905514105020593, 0.8905514105020593, 0.8905514105020593, 0.0759738621579806, 0.9240261378420194, 0.5000000000000000, 0.9240261378420194, 0.9240261378420194, 0.5000000000000000, 0.5000000000000000, 0.0759738621579806, 0.0759738621579806, 0.9240261378420194, 0.5000000000000000, 0.0759738621579806 }; static double y_save[] = { 0.5000000000000000, 0.5000000000000000, 0.9940430805943383, 0.5000000000000000, 0.5000000000000000, 0.0059569194056617, 0.2960224163208403, 0.7039775836791597, 0.7039775836791597, 0.2960224163208403, 0.2960224163208403, 0.7039775836791597, 0.7039775836791597, 0.2960224163208403, 0.1094485894979407, 0.8905514105020593, 0.8905514105020593, 0.1094485894979407, 0.1094485894979407, 0.8905514105020593, 0.8905514105020593, 0.1094485894979407, 0.0759738621579806, 0.5000000000000000, 0.9240261378420194, 0.9240261378420194, 0.5000000000000000, 0.0759738621579806, 0.0759738621579806, 0.5000000000000000, 0.9240261378420194, 0.0759738621579806, 0.9240261378420194, 0.5000000000000000 }; static double z_save[] = { 0.5000000000000000, 0.9940430805943383, 0.5000000000000000, 0.0059569194056617, 0.5000000000000000, 0.5000000000000000, 0.2960224163208403, 0.7039775836791597, 0.2960224163208403, 0.2960224163208403, 0.7039775836791597, 0.2960224163208403, 0.7039775836791597, 0.7039775836791597, 0.1094485894979407, 0.8905514105020593, 0.1094485894979407, 0.1094485894979407, 0.8905514105020593, 0.1094485894979407, 0.8905514105020593, 0.8905514105020593, 0.5000000000000000, 0.0759738621579806, 0.0759738621579806, 0.5000000000000000, 0.9240261378420194, 0.9240261378420194, 0.0759738621579806, 0.9240261378420194, 0.5000000000000000, 0.5000000000000000, 0.9240261378420194, 0.0759738621579806 }; double w_save[] = { 0.0250162363729459, 0.0250162363729459, 0.0250162363729459, 0.0250162363729459, 0.0250162363729459, 0.0250162363729459, 0.0571442320087316, 0.0571442320087316, 0.0571442320087316, 0.0571442320087316, 0.0571442320087316, 0.0571442320087316, 0.0571442320087316, 0.0571442320087316, 0.0192245175082448, 0.0192245175082448, 0.0192245175082448, 0.0192245175082448, 0.0192245175082448, 0.0192245175082448, 0.0192245175082448, 0.0192245175082448, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761, 0.0199127154688761 }; r8vec_copy ( n, x_save, x ); r8vec_copy ( n, y_save, y ); r8vec_copy ( n, z_save, z ); r8vec_copy ( n, w_save, w ); return; } void rule09 ( int n, double x[], double y[], double z[], double w[] ) /******************************************************************************/ /* Purpose: rule09() returns the hexahedron rule of precision 9. Discussion: We suppose we are given a hexahedron H with vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), We call a rule with n points, returning coordinates x, y, z, and weights w. Then the integral I of f(x,y,z) over H is approximated by Q as follows: Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) Licensing: This code is distributed under the MIT license. Modified: 02 May 2023 Author: John Burkardt Reference: Freddie Witherden, Peter Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers and Mathematics with Applications, Volume 69, pages 1232-1241, 2015. Input: int n: the number of quadrature points for this rule. Output: double x(n), y(n), x(n): the coordinates of quadrature points. double w[n]: the quadrature weights, which sum to 1. */ { static double x_save[] = { 0.1931592652041455, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.8068407347958545, 0.5000000000000000, 0.9350498923309880, 0.0649501076690120, 0.0649501076690120, 0.0649501076690120, 0.0649501076690120, 0.9350498923309880, 0.9350498923309880, 0.9350498923309880, 0.7820554035100150, 0.2179445964899850, 0.2179445964899850, 0.2179445964899850, 0.2179445964899850, 0.7820554035100150, 0.7820554035100150, 0.7820554035100150, 0.0611564383711609, 0.9388435616288391, 0.5000000000000000, 0.9388435616288391, 0.9388435616288391, 0.5000000000000000, 0.5000000000000000, 0.0611564383711609, 0.0611564383711609, 0.9388435616288391, 0.5000000000000000, 0.0611564383711609, 0.7161339513154311, 0.2838660486845689, 0.2838660486845689, 0.0307347890676641, 0.2838660486845689, 0.9692652109323359, 0.7161339513154311, 0.7161339513154311, 0.2838660486845689, 0.7161339513154311, 0.0307347890676641, 0.2838660486845689, 0.2838660486845689, 0.2838660486845689, 0.7161339513154311, 0.9692652109323359, 0.7161339513154311, 0.9692652109323359, 0.9692652109323359, 0.0307347890676641, 0.7161339513154311, 0.2838660486845689, 0.7161339513154311, 0.0307347890676641 }; static double y_save[] = { 0.5000000000000000, 0.5000000000000000, 0.8068407347958545, 0.5000000000000000, 0.5000000000000000, 0.1931592652041455, 0.0649501076690120, 0.9350498923309880, 0.9350498923309880, 0.0649501076690120, 0.0649501076690120, 0.9350498923309880, 0.9350498923309880, 0.0649501076690120, 0.2179445964899850, 0.7820554035100150, 0.7820554035100150, 0.2179445964899850, 0.2179445964899850, 0.7820554035100150, 0.7820554035100150, 0.2179445964899850, 0.0611564383711609, 0.5000000000000000, 0.9388435616288391, 0.9388435616288391, 0.5000000000000000, 0.0611564383711609, 0.0611564383711609, 0.5000000000000000, 0.9388435616288391, 0.0611564383711609, 0.9388435616288391, 0.5000000000000000, 0.9692652109323359, 0.0307347890676641, 0.9692652109323359, 0.2838660486845689, 0.2838660486845689, 0.7161339513154311, 0.0307347890676641, 0.0307347890676641, 0.0307347890676641, 0.7161339513154311, 0.7161339513154311, 0.7161339513154311, 0.7161339513154311, 0.2838660486845689, 0.2838660486845689, 0.2838660486845689, 0.9692652109323359, 0.7161339513154311, 0.2838660486845689, 0.7161339513154311, 0.7161339513154311, 0.9692652109323359, 0.2838660486845689, 0.2838660486845689 }; static double z_save[] = { 0.5000000000000000, 0.8068407347958545, 0.5000000000000000, 0.1931592652041455, 0.5000000000000000, 0.5000000000000000, 0.0649501076690120, 0.9350498923309880, 0.0649501076690120, 0.0649501076690120, 0.9350498923309880, 0.0649501076690120, 0.9350498923309880, 0.9350498923309880, 0.2179445964899850, 0.7820554035100150, 0.2179445964899850, 0.2179445964899850, 0.7820554035100150, 0.2179445964899850, 0.7820554035100150, 0.7820554035100150, 0.5000000000000000, 0.0611564383711609, 0.0611564383711609, 0.5000000000000000, 0.9388435616288391, 0.9388435616288391, 0.0611564383711609, 0.9388435616288391, 0.5000000000000000, 0.5000000000000000, 0.9388435616288391, 0.0611564383711609, 0.2838660486845689, 0.7161339513154311, 0.2838660486845689, 0.7161339513154311, 0.9692652109323359, 0.2838660486845689, 0.2838660486845689, 0.7161339513154311, 0.2838660486845689, 0.0307347890676641, 0.2838660486845689, 0.0307347890676641, 0.9692652109323359, 0.0307347890676641, 0.9692652109323359, 0.2838660486845689, 0.7161339513154311, 0.7161339513154311, 0.7161339513154311, 0.7161339513154311, 0.9692652109323359, 0.7161339513154311, 0.0307347890676641, 0.2838660486845689 }; double w_save[] = { 0.0541593744687068, 0.0541593744687068, 0.0541593744687068, 0.0541593744687068, 0.0541593744687068, 0.0541593744687068, 0.0062685994124186, 0.0062685994124186, 0.0062685994124186, 0.0062685994124186, 0.0062685994124186, 0.0062685994124186, 0.0062685994124186, 0.0062685994124186, 0.0248574797680029, 0.0248574797680029, 0.0248574797680029, 0.0248574797680029, 0.0248574797680029, 0.0248574797680029, 0.0248574797680029, 0.0248574797680029, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0114737257670222, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717, 0.0120146004391717 }; r8vec_copy ( n, x_save, x ); r8vec_copy ( n, y_save, y ); r8vec_copy ( n, z_save, z ); r8vec_copy ( n, w_save, w ); return; } void rule11 ( int n, double x[], double y[], double z[], double w[] ) /******************************************************************************/ /* Purpose: rule11() returns the hexahedron rule of precision 11. Discussion: We suppose we are given a hexahedron H with vertices (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), We call a rule with n points, returning coordinates x, y, z, and weights w. Then the integral I of f(x,y,z) over H is approximated by Q as follows: Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) Licensing: This code is distributed under the MIT license. Modified: 02 May 2023 Author: John Burkardt Reference: Freddie Witherden, Peter Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers and Mathematics with Applications, Volume 69, pages 1232-1241, 2015. Input: int n: the number of quadrature points for this rule. Output: double x(n), y(n), x(n): the coordinates of quadrature points. double w[n]: the quadrature weights, which sum to 1. */ { static double x_save[] = { 0.0936928329501868, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.9063071670498133, 0.5000000000000000, 0.8008376320991313, 0.1991623679008687, 0.1991623679008687, 0.1991623679008687, 0.1991623679008687, 0.8008376320991313, 0.8008376320991313, 0.8008376320991313, 0.9277278805088800, 0.0722721194911200, 0.0722721194911200, 0.0722721194911200, 0.0722721194911200, 0.9277278805088800, 0.9277278805088800, 0.9277278805088800, 0.6566967022580273, 0.3433032977419727, 0.3433032977419727, 0.3433032977419727, 0.3433032977419727, 0.6566967022580273, 0.6566967022580273, 0.6566967022580273, 0.1326658565014960, 0.8673341434985040, 0.5000000000000000, 0.8673341434985040, 0.8673341434985040, 0.5000000000000000, 0.5000000000000000, 0.1326658565014960, 0.1326658565014960, 0.8673341434985040, 0.5000000000000000, 0.1326658565014960, 0.7253999675572547, 0.2746000324427453, 0.2746000324427453, 0.0174501672436448, 0.2746000324427453, 0.9825498327563551, 0.7253999675572547, 0.7253999675572547, 0.2746000324427453, 0.7253999675572547, 0.0174501672436448, 0.2746000324427453, 0.2746000324427453, 0.2746000324427453, 0.7253999675572547, 0.9825498327563551, 0.7253999675572547, 0.9825498327563551, 0.9825498327563551, 0.0174501672436448, 0.7253999675572547, 0.2746000324427453, 0.7253999675572547, 0.0174501672436448, 0.9706224286053016, 0.0293775713946984, 0.0293775713946984, 0.3230485927016850, 0.0293775713946984, 0.6769514072983150, 0.9706224286053016, 0.9706224286053016, 0.0293775713946984, 0.9706224286053016, 0.3230485927016850, 0.0293775713946984, 0.0293775713946984, 0.0293775713946984, 0.9706224286053016, 0.6769514072983150, 0.9706224286053016, 0.6769514072983150, 0.6769514072983150, 0.3230485927016850, 0.9706224286053016, 0.0293775713946984, 0.9706224286053016, 0.3230485927016850 }; static double y_save[] = { 0.5000000000000000, 0.5000000000000000, 0.9063071670498133, 0.5000000000000000, 0.5000000000000000, 0.0936928329501868, 0.1991623679008687, 0.8008376320991313, 0.8008376320991313, 0.1991623679008687, 0.1991623679008687, 0.8008376320991313, 0.8008376320991313, 0.1991623679008687, 0.0722721194911200, 0.9277278805088800, 0.9277278805088800, 0.0722721194911200, 0.0722721194911200, 0.9277278805088800, 0.9277278805088800, 0.0722721194911200, 0.3433032977419727, 0.6566967022580273, 0.6566967022580273, 0.3433032977419727, 0.3433032977419727, 0.6566967022580273, 0.6566967022580273, 0.3433032977419727, 0.1326658565014960, 0.5000000000000000, 0.8673341434985040, 0.8673341434985040, 0.5000000000000000, 0.1326658565014960, 0.1326658565014960, 0.5000000000000000, 0.8673341434985040, 0.1326658565014960, 0.8673341434985040, 0.5000000000000000, 0.9825498327563551, 0.0174501672436448, 0.9825498327563551, 0.2746000324427453, 0.2746000324427453, 0.7253999675572547, 0.0174501672436448, 0.0174501672436448, 0.0174501672436448, 0.7253999675572547, 0.7253999675572547, 0.7253999675572547, 0.7253999675572547, 0.2746000324427453, 0.2746000324427453, 0.2746000324427453, 0.9825498327563551, 0.7253999675572547, 0.2746000324427453, 0.7253999675572547, 0.7253999675572547, 0.9825498327563551, 0.2746000324427453, 0.2746000324427453, 0.6769514072983150, 0.3230485927016850, 0.6769514072983150, 0.0293775713946984, 0.0293775713946984, 0.9706224286053016, 0.3230485927016850, 0.3230485927016850, 0.3230485927016850, 0.9706224286053016, 0.9706224286053016, 0.9706224286053016, 0.9706224286053016, 0.0293775713946984, 0.0293775713946984, 0.0293775713946984, 0.6769514072983150, 0.9706224286053016, 0.0293775713946984, 0.9706224286053016, 0.9706224286053016, 0.6769514072983150, 0.0293775713946984, 0.0293775713946984 }; static double z_save[] = { 0.5000000000000000, 0.9063071670498133, 0.5000000000000000, 0.0936928329501868, 0.5000000000000000, 0.5000000000000000, 0.1991623679008687, 0.8008376320991313, 0.1991623679008687, 0.1991623679008687, 0.8008376320991313, 0.1991623679008687, 0.8008376320991313, 0.8008376320991313, 0.0722721194911200, 0.9277278805088800, 0.0722721194911200, 0.0722721194911200, 0.9277278805088800, 0.0722721194911200, 0.9277278805088800, 0.9277278805088800, 0.3433032977419727, 0.6566967022580273, 0.3433032977419727, 0.3433032977419727, 0.6566967022580273, 0.3433032977419727, 0.6566967022580273, 0.6566967022580273, 0.5000000000000000, 0.1326658565014960, 0.1326658565014960, 0.5000000000000000, 0.8673341434985040, 0.8673341434985040, 0.1326658565014960, 0.8673341434985040, 0.5000000000000000, 0.5000000000000000, 0.8673341434985040, 0.1326658565014960, 0.2746000324427453, 0.7253999675572547, 0.2746000324427453, 0.7253999675572547, 0.9825498327563551, 0.2746000324427453, 0.2746000324427453, 0.7253999675572547, 0.2746000324427453, 0.0174501672436448, 0.2746000324427453, 0.0174501672436448, 0.9825498327563551, 0.0174501672436448, 0.9825498327563551, 0.2746000324427453, 0.7253999675572547, 0.7253999675572547, 0.7253999675572547, 0.7253999675572547, 0.9825498327563551, 0.7253999675572547, 0.0174501672436448, 0.2746000324427453, 0.0293775713946984, 0.9706224286053016, 0.0293775713946984, 0.9706224286053016, 0.6769514072983150, 0.0293775713946984, 0.0293775713946984, 0.9706224286053016, 0.0293775713946984, 0.3230485927016850, 0.0293775713946984, 0.3230485927016850, 0.6769514072983150, 0.3230485927016850, 0.6769514072983150, 0.0293775713946984, 0.9706224286053016, 0.9706224286053016, 0.9706224286053016, 0.9706224286053016, 0.6769514072983150, 0.9706224286053016, 0.3230485927016850, 0.0293775713946984 }; double w_save[] = { 0.0253096342016000, 0.0253096342016000, 0.0253096342016000, 0.0253096342016000, 0.0253096342016000, 0.0253096342016000, 0.0146922934945570, 0.0146922934945570, 0.0146922934945570, 0.0146922934945570, 0.0146922934945570, 0.0146922934945570, 0.0146922934945570, 0.0146922934945570, 0.0055804890098537, 0.0055804890098537, 0.0055804890098537, 0.0055804890098537, 0.0055804890098537, 0.0055804890098537, 0.0055804890098537, 0.0055804890098537, 0.0269990056568711, 0.0269990056568711, 0.0269990056568711, 0.0269990056568711, 0.0269990056568711, 0.0269990056568711, 0.0269990056568711, 0.0269990056568711, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0181499182325145, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0076802492622294, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527, 0.0028267870173527 }; r8vec_copy ( n, x_save, x ); r8vec_copy ( n, y_save, y ); r8vec_copy ( n, z_save, z ); r8vec_copy ( n, w_save, w ); return; }