# include # include # include # include # include # include "pyramid_felippa_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 STATICE 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 *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 Parameters: Input, int M, the spatial dimension. Input, int N, the number of evaluation points. Input, int E[M], the exponents. Input, 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; } /******************************************************************************/ double pyramid_unit_monomial ( int expon[3] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_MONOMIAL: monomial integral in a unit pyramid. Discussion: This function returns the value of the integral of X^ALPHA Y^BETA Z^GAMMA over the unit pyramid. The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. Licensing: This code is distributed under the MIT license. Modified: 24 March 2008 Author: John Burkardt Reference: Arthur Stroud, Approximate Calculation of Multiple Integrals, Prentice Hall, 1971, ISBN: 0130438936, LC: QA311.S85. Parameters: Input, int EXPON[3], the exponents. Output, double PYRAMID_UNIT_MONOMIAL, the volume of the pyramid. */ { int i; int i_hi; double value; value = 0.0; if ( ( expon[0] % 2 ) == 0 && ( expon[1] % 2 ) == 0 ) { i_hi = 2 + expon[0] + expon[1]; for ( i = 0; i <= i_hi; i++ ) { value = value + r8_mop ( i ) * r8_choose ( i_hi, i ) / ( double ) ( i + expon[2] + 1 ); } value = value * 2.0 / ( double ) ( expon[0] + 1 ) * 2.0 / ( double ) ( expon[1] + 1 ); } return value; } /******************************************************************************/ void pyramid_unit_o01 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O01 returns a 1 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Arthur Stroud, Approximate Calculation of Multiple Integrals, Prentice Hall, 1971, ISBN: 0130438936, LC: QA311.S85. Parameters: Output, double W[1], the weights. Output, double XYZ[3*1], the abscissas. */ { int order = 1; double w_save[1] = { 1.0 }; double xyz_save[3*1] = { 0.00, 0.00, 0.25 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o05 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O05 returns a 5 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Carlos Felippa, A compendium of FEM integration formulas for symbolic work, Engineering Computation, Volume 21, Number 8, 2004, pages 867-890. Parameters: Output, double W[5], the weights. Output, double XYZ[3*5], the abscissas. */ { int order = 5; double w_save[5] = { 0.21093750000000000000, 0.21093750000000000000, 0.21093750000000000000, 0.21093750000000000000, 0.15625000000000000000 }; double xyz_save[3*5] = { -0.48686449556014765641, -0.48686449556014765641, 0.16666666666666666667, 0.48686449556014765641, -0.48686449556014765641, 0.16666666666666666667, 0.48686449556014765641, 0.48686449556014765641, 0.16666666666666666667, -0.48686449556014765641, 0.48686449556014765641, 0.16666666666666666667, 0.00000000000000000000, 0.00000000000000000000, 0.70000000000000000000 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o06 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O06 returns a 6 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Carlos Felippa, A compendium of FEM integration formulas for symbolic work, Engineering Computation, Volume 21, Number 8, 2004, pages 867-890. Parameters: Output, double W[6], the weights. Output, double XYZ[3*6], the abscissas. */ { int order = 6; double w_save[6] = { 0.21000000000000000000, 0.21000000000000000000, 0.21000000000000000000, 0.21000000000000000000, 0.06000000000000000000, 0.10000000000000000000 }; double xyz_save[3*6] = { -0.48795003647426658968, -0.48795003647426658968, 0.16666666666666666667, 0.48795003647426658968, -0.48795003647426658968, 0.16666666666666666667, 0.48795003647426658968, 0.48795003647426658968, 0.16666666666666666667, -0.48795003647426658968, 0.48795003647426658968, 0.16666666666666666667, 0.00000000000000000000, 0.00000000000000000000, 0.58333333333333333333, 0.00000000000000000000, 0.00000000000000000000, 0.75000000000000000000 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o08 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O08 returns an 8 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Carlos Felippa, A compendium of FEM integration formulas for symbolic work, Engineering Computation, Volume 21, Number 8, 2004, pages 867-890. Parameters: Output, double W[8], the weights. Output, double XYZ[3*8], the abscissas. */ { int order = 8; double w_save[8] = { 0.075589411559869072938, 0.075589411559869072938, 0.075589411559869072938, 0.075589411559869072938, 0.17441058844013092706, 0.17441058844013092706, 0.17441058844013092706, 0.17441058844013092706 }; double xyz_save[3*8] = { -0.26318405556971359557, -0.26318405556971359557, 0.54415184401122528880, 0.26318405556971359557, -0.26318405556971359557, 0.54415184401122528880, 0.26318405556971359557, 0.26318405556971359557, 0.54415184401122528880, -0.26318405556971359557, 0.26318405556971359557, 0.54415184401122528880, -0.50661630334978742377, -0.50661630334978742377, 0.12251482265544137787, 0.50661630334978742377, -0.50661630334978742377, 0.12251482265544137787, 0.50661630334978742377, 0.50661630334978742377, 0.12251482265544137787, -0.50661630334978742377, 0.50661630334978742377, 0.12251482265544137787 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o08b ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O08B returns an 8 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Carlos Felippa, A compendium of FEM integration formulas for symbolic work, Engineering Computation, Volume 21, Number 8, 2004, pages 867-890. Parameters: Output, double W[8], the weights. Output, double XYZ[3*8], the abscissas. */ { int order = 8; double w_save[8] = { 0.16438287736328777572, 0.16438287736328777572, 0.16438287736328777572, 0.16438287736328777572, 0.085617122636712224276, 0.085617122636712224276, 0.085617122636712224276, 0.085617122636712224276 }; double xyz_save[3*8] = { -0.51197009372656270107, -0.51197009372656270107, 0.11024490204163285720, 0.51197009372656270107, -0.51197009372656270107, 0.11024490204163285720, 0.51197009372656270107, 0.51197009372656270107, 0.11024490204163285720, -0.51197009372656270107, 0.51197009372656270107, 0.11024490204163285720, -0.28415447557052037456, -0.28415447557052037456, 0.518326526529795714229, 0.28415447557052037456, -0.28415447557052037456, 0.518326526529795714229, 0.28415447557052037456, 0.28415447557052037456, 0.518326526529795714229, -0.28415447557052037456, 0.28415447557052037456, 0.518326526529795714229 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o09 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O09 returns a 9 point quadrature rule for the unit pyramid. / Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Carlos Felippa, A compendium of FEM integration formulas for symbolic work, Engineering Computation, Volume 21, Number 8, 2004, pages 867-890. Parameters: Output, double W[9], the weights. Output, double XYZ[3*9], the abscissas. */ { int order = 9; double w_save[9] = { 0.13073389672275944791, 0.13073389672275944791, 0.13073389672275944791, 0.13073389672275944791, 0.10989110327724055209, 0.10989110327724055209, 0.10989110327724055209, 0.10989110327724055209, 0.03750000000000000000 }; double xyz_save[3*9] = { -0.52966422253852215131, -0.52966422253852215131, 0.08176876558246862335, 0.52966422253852215131, -0.52966422253852215131, 0.08176876558246862335, 0.52966422253852215131, 0.52966422253852215131, 0.08176876558246862335, -0.52966422253852215131, 0.52966422253852215131, 0.08176876558246862335, -0.34819753825720418039, -0.34819753825720418039, 0.400374091560388519511, 0.34819753825720418039, -0.34819753825720418039, 0.400374091560388519511, 0.34819753825720418039, 0.34819753825720418039, 0.400374091560388519511, -0.34819753825720418039, 0.34819753825720418039, 0.400374091560388519511, 0.00000000000000000000, 0.00000000000000000000, 0.83333333333333333333 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o13 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O13 returns a 13 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Carlos Felippa, A compendium of FEM integration formulas for symbolic work, Engineering Computation, Volume 21, Number 8, 2004, pages 867-890. Parameters: Output, double W[13], the weights. Output, double XYZ[3*13], the abscissas. */ { int order = 13; double w_save[13] = { 0.063061594202898550725, 0.063061594202898550725, 0.063061594202898550725, 0.063061594202898550725, 0.042101946815575556199, 0.042101946815575556199, 0.042101946815575556199, 0.042101946815575556199, 0.13172030707666776585, 0.13172030707666776585, 0.13172030707666776585, 0.13172030707666776585, 0.05246460761943250889 }; double xyz_save[3*13] = { -0.38510399211870384331, -0.38510399211870384331, 0.428571428571428571429, 0.38510399211870384331, -0.38510399211870384331, 0.428571428571428571429, 0.38510399211870384331, 0.38510399211870384331, 0.428571428571428571429, -0.38510399211870384331, 0.38510399211870384331, 0.428571428571428571429, -0.40345831960728204766, 0.00000000000000000000, 0.33928571428571428571, 0.40345831960728204766, 0.00000000000000000000, 0.33928571428571428571, 0.00000000000000000000, -0.40345831960728204766, 0.33928571428571428571, 0.00000000000000000000, 0.40345831960728204766, 0.33928571428571428571, -0.53157877436961973359, -0.53157877436961973359, 0.08496732026143790850, 0.53157877436961973359, -0.53157877436961973359, 0.08496732026143790850, 0.53157877436961973359, 0.53157877436961973359, 0.08496732026143790850, -0.53157877436961973359, 0.53157877436961973359, 0.08496732026143790850, 0.00000000000000000000, 0.00000000000000000000, 0.76219701803768503595 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o18 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O18 returns an 18 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Carlos Felippa, A compendium of FEM integration formulas for symbolic work, Engineering Computation, Volume 21, Number 8, 2004, pages 867-890. Parameters: Output, double W[18], the weights. Output, double XYZ[3*18], the abscissas. */ { int order = 18; double w_save[18] = { 0.023330065296255886709, 0.037328104474009418735, 0.023330065296255886709, 0.037328104474009418735, 0.059724967158415069975, 0.037328104474009418735, 0.023330065296255886709, 0.037328104474009418735, 0.023330065296255886709, 0.05383042853090460712, 0.08612868564944737139, 0.05383042853090460712, 0.08612868564944737139, 0.13780589703911579422, 0.08612868564944737139, 0.05383042853090460712, 0.08612868564944737139, 0.05383042853090460712 }; double xyz_save[3*18] = { -0.35309846330877704481, -0.35309846330877704481, 0.544151844011225288800, 0.00000000000000000000, -0.35309846330877704481, 0.544151844011225288800, 0.35309846330877704481, -0.35309846330877704481, 0.544151844011225288800, -0.35309846330877704481, 0.00000000000000000000, 0.544151844011225288800, 0.00000000000000000000, 0.00000000000000000000, 0.544151844011225288800, 0.35309846330877704481, 0.00000000000000000000, 0.544151844011225288800, -0.35309846330877704481, 0.35309846330877704481, 0.544151844011225288800, 0.00000000000000000000, 0.35309846330877704481, 0.544151844011225288800, 0.35309846330877704481, 0.35309846330877704481, 0.544151844011225288800, -0.67969709567986745790, -0.67969709567986745790, 0.12251482265544137787, 0.00000000000000000000, -0.67969709567986745790, 0.12251482265544137787, 0.67969709567986745790, -0.67969709567986745790, 0.12251482265544137787, -0.67969709567986745790, 0.00000000000000000000, 0.12251482265544137787, 0.00000000000000000000, 0.00000000000000000000, 0.12251482265544137787, 0.67969709567986745790, 0.00000000000000000000, 0.12251482265544137787, -0.67969709567986745790, 0.67969709567986745790, 0.12251482265544137787, 0.00000000000000000000, 0.67969709567986745790, 0.12251482265544137787, 0.67969709567986745790, 0.67969709567986745790, 0.12251482265544137787 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o27 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O27 returns a 27 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Carlos Felippa, A compendium of FEM integration formulas for symbolic work, Engineering Computation, Volume 21, Number 8, 2004, pages 867-890. Parameters: Output, double W[27], the weights. Output, double XYZ[3*27], the abscissas. */ { int order = 27; double w_save[27] = { 0.036374157653908938268, 0.05819865224625430123, 0.036374157653908938268, 0.05819865224625430123, 0.09311784359400688197, 0.05819865224625430123, 0.036374157653908938268, 0.05819865224625430123, 0.036374157653908938268, 0.033853303069413431019, 0.054165284911061489631, 0.033853303069413431019, 0.054165284911061489631, 0.08666445585769838341, 0.054165284911061489631, 0.033853303069413431019, 0.054165284911061489631, 0.033853303069413431019, 0.006933033103838124540, 0.011092852966140999264, 0.006933033103838124540, 0.011092852966140999264, 0.017748564745825598822, 0.011092852966140999264, 0.006933033103838124540, 0.011092852966140999264, 0.006933033103838124540 }; double xyz_save[3*27] = { -0.7180557413198889387, -0.7180557413198889387, 0.07299402407314973216, 0.00000000000000000000, -0.7180557413198889387, 0.07299402407314973216, 0.7180557413198889387, -0.7180557413198889387, 0.07299402407314973216, -0.7180557413198889387, 0.00000000000000000000, 0.07299402407314973216, 0.00000000000000000000, 0.00000000000000000000, 0.07299402407314973216, 0.7180557413198889387, 0.00000000000000000000, 0.07299402407314973216, -0.7180557413198889387, 0.7180557413198889387, 0.07299402407314973216, 0.00000000000000000000, 0.7180557413198889387, 0.07299402407314973216, 0.7180557413198889387, 0.7180557413198889387, 0.07299402407314973216, -0.50580870785392503961, -0.50580870785392503961, 0.34700376603835188472, 0.00000000000000000000, -0.50580870785392503961, 0.34700376603835188472, 0.50580870785392503961, -0.50580870785392503961, 0.34700376603835188472, -0.50580870785392503961, 0.00000000000000000000, 0.34700376603835188472, 0.00000000000000000000, 0.00000000000000000000, 0.34700376603835188472, 0.50580870785392503961, 0.00000000000000000000, 0.34700376603835188472, -0.50580870785392503961, 0.50580870785392503961, 0.34700376603835188472, 0.00000000000000000000, 0.50580870785392503961, 0.34700376603835188472, 0.50580870785392503961, 0.50580870785392503961, 0.34700376603835188472, -0.22850430565396735360, -0.22850430565396735360, 0.70500220988849838312, 0.00000000000000000000, -0.22850430565396735360, 0.70500220988849838312, 0.22850430565396735360, -0.22850430565396735360, 0.70500220988849838312, -0.22850430565396735360, 0.00000000000000000000, 0.70500220988849838312, 0.00000000000000000000, 0.00000000000000000000, 0.70500220988849838312, 0.22850430565396735360, 0.00000000000000000000, 0.70500220988849838312, -0.22850430565396735360, 0.22850430565396735360, 0.70500220988849838312, 0.00000000000000000000, 0.22850430565396735360, 0.70500220988849838312, 0.22850430565396735360, 0.22850430565396735360, 0.70500220988849838312 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ void pyramid_unit_o48 ( double w[], double xyz[] ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_O48 returns a 48 point quadrature rule for the unit pyramid. Discussion: The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. When Z is zero, the integration region is a square lying in the (X,Y) plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the radius of the square diminishes, and when Z reaches 1, the square has contracted to the single point (0,0,1). Licensing: This code is distributed under the MIT license. Modified: 02 April 2009 Author: John Burkardt Reference: Arthur Stroud, Approximate Calculation of Multiple Integrals, Prentice Hall, 1971, ISBN: 0130438936, LC: QA311.S85. Parameters: Output, double W[48], the weights. Output, double XYZ[3*48], the abscissas. */ { int order = 48; double w_save[48] = { 2.01241939442682455E-002, 2.01241939442682455E-002, 2.01241939442682455E-002, 2.01241939442682455E-002, 2.60351137043010779E-002, 2.60351137043010779E-002, 2.60351137043010779E-002, 2.60351137043010779E-002, 1.24557795239745531E-002, 1.24557795239745531E-002, 1.24557795239745531E-002, 1.24557795239745531E-002, 1.87873998794808156E-003, 1.87873998794808156E-003, 1.87873998794808156E-003, 1.87873998794808156E-003, 4.32957927807745280E-002, 4.32957927807745280E-002, 4.32957927807745280E-002, 4.32957927807745280E-002, 1.97463249834127288E-002, 1.97463249834127288E-002, 1.97463249834127288E-002, 1.97463249834127288E-002, 5.60127223523590526E-002, 5.60127223523590526E-002, 5.60127223523590526E-002, 5.60127223523590526E-002, 2.55462562927473852E-002, 2.55462562927473852E-002, 2.55462562927473852E-002, 2.55462562927473852E-002, 2.67977366291788643E-002, 2.67977366291788643E-002, 2.67977366291788643E-002, 2.67977366291788643E-002, 1.22218992265373354E-002, 1.22218992265373354E-002, 1.22218992265373354E-002, 1.22218992265373354E-002, 4.04197740453215038E-003, 4.04197740453215038E-003, 4.04197740453215038E-003, 4.04197740453215038E-003, 1.84346316995826843E-003, 1.84346316995826843E-003, 1.84346316995826843E-003, 1.84346316995826843E-003 }; double xyz_save[3*48] = { 0.88091731624450909E+00, 0.00000000000000000E+00, 4.85005494469969989E-02, -0.88091731624450909E+00, 0.00000000000000000E+00, 4.85005494469969989E-02, 0.00000000000000000E+00, 0.88091731624450909E+00, 4.85005494469969989E-02, 0.00000000000000000E+00, -0.88091731624450909E+00, 4.85005494469969989E-02, 0.70491874112648223E+00, 0.00000000000000000E+00, 0.23860073755186201E+00, -0.70491874112648223E+00, 0.00000000000000000E+00, 0.23860073755186201E+00, 0.00000000000000000E+00, 0.70491874112648223E+00, 0.23860073755186201E+00, 0.00000000000000000E+00, -0.70491874112648223E+00, 0.23860073755186201E+00, 0.44712732143189760E+00, 0.00000000000000000E+00, 0.51704729510436798E+00, -0.44712732143189760E+00, 0.00000000000000000E+00, 0.51704729510436798E+00, 0.00000000000000000E+00, 0.44712732143189760E+00, 0.51704729510436798E+00, 0.00000000000000000E+00, -0.44712732143189760E+00, 0.51704729510436798E+00, 0.18900486065123448E+00, 0.00000000000000000E+00, 0.79585141789677305E+00, -0.18900486065123448E+00, 0.00000000000000000E+00, 0.79585141789677305E+00, 0.00000000000000000E+00, 0.18900486065123448E+00, 0.79585141789677305E+00, 0.00000000000000000E+00, -0.18900486065123448E+00, 0.79585141789677305E+00, 0.36209733410322176E+00, 0.36209733410322176E+00, 4.85005494469969989E-02, -0.36209733410322176E+00, 0.36209733410322176E+00, 4.85005494469969989E-02, -0.36209733410322176E+00, -0.36209733410322176E+00, 4.85005494469969989E-02, 0.36209733410322176E+00, -0.36209733410322176E+00, 4.85005494469969989E-02, 0.76688932060387538E+00, 0.76688932060387538E+00, 4.85005494469969989E-02, -0.76688932060387538E+00, 0.76688932060387538E+00, 4.85005494469969989E-02, -0.76688932060387538E+00, -0.76688932060387538E+00, 4.85005494469969989E-02, 0.76688932060387538E+00, -0.76688932060387538E+00, 4.85005494469969989E-02, 0.28975386476618070E+00, 0.28975386476618070E+00, 0.23860073755186201E+00, -0.28975386476618070E+00, 0.28975386476618070E+00, 0.23860073755186201E+00, -0.28975386476618070E+00, -0.28975386476618070E+00, 0.23860073755186201E+00, 0.28975386476618070E+00, -0.28975386476618070E+00, 0.23860073755186201E+00, 0.61367241226233160E+00, 0.61367241226233160E+00, 0.23860073755186201E+00, -0.61367241226233160E+00, 0.61367241226233160E+00, 0.23860073755186201E+00, -0.61367241226233160E+00, -0.61367241226233160E+00, 0.23860073755186201E+00, 0.61367241226233160E+00, -0.61367241226233160E+00, 0.23860073755186201E+00, 0.18378979287798017E+00, 0.18378979287798017E+00, 0.51704729510436798E+00, -0.18378979287798017E+00, 0.18378979287798017E+00, 0.51704729510436798E+00, -0.18378979287798017E+00, -0.18378979287798017E+00, 0.51704729510436798E+00, 0.18378979287798017E+00, -0.18378979287798017E+00, 0.51704729510436798E+00, 0.38925011625173161E+00, 0.38925011625173161E+00, 0.51704729510436798E+00, -0.38925011625173161E+00, 0.38925011625173161E+00, 0.51704729510436798E+00, -0.38925011625173161E+00, -0.38925011625173161E+00, 0.51704729510436798E+00, 0.38925011625173161E+00, -0.38925011625173161E+00, 0.51704729510436798E+00, 7.76896479525748113E-02, 7.76896479525748113E-02, 0.79585141789677305E+00, -7.76896479525748113E-02, 7.76896479525748113E-02, 0.79585141789677305E+00, -7.76896479525748113E-02, -7.76896479525748113E-02, 0.79585141789677305E+00, 7.76896479525748113E-02, -7.76896479525748113E-02, 0.79585141789677305E+00, 0.16453962988669860E+00, 0.16453962988669860E+00, 0.79585141789677305E+00, -0.16453962988669860E+00, 0.16453962988669860E+00, 0.79585141789677305E+00, -0.16453962988669860E+00, -0.16453962988669860E+00, 0.79585141789677305E+00, 0.16453962988669860E+00, -0.16453962988669860E+00, 0.79585141789677305E+00 }; r8vec_copy ( order, w_save, w ); r8vec_copy ( 3 * order, xyz_save, xyz ); return; } /******************************************************************************/ double pyramid_unit_volume ( ) /******************************************************************************/ /* Purpose: PYRAMID_UNIT_VOLUME: volume of a unit pyramid with square base. Discussion: The volume of this unit pyramid is 4/3. The integration region is: - ( 1 - Z ) <= X <= 1 - Z - ( 1 - Z ) <= Y <= 1 - Z 0 <= Z <= 1. Licensing: This code is distributed under the MIT license. Modified: 22 March 2008 Author: John Burkardt Parameters: Output, double PYRAMID_UNIT_VOLUME, the volume of the pyramid. */ { double volume; volume = 4.0 / 3.0; return volume; } /******************************************************************************/ double r8_choose ( int n, int k ) /******************************************************************************/ /* Purpose: R8_CHOOSE computes the binomial coefficient C(N,K) as an R8. Discussion: The value is calculated in such a way as to avoid overflow and roundoff. The calculation is done in R8 arithmetic. The formula used is: C(N,K) = N! / ( K! * (N-K)! ) Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Reference: ML Wolfson, HV Wright, Algorithm 160: Combinatorial of M Things Taken N at a Time, Communications of the ACM, Volume 6, Number 4, April 1963, page 161. Parameters: Input, int N, K, the values of N and K. Output, double R8_CHOOSE, the number of combinations of N things taken K at a time. */ { int i; int mn; int mx; double value; if ( k < n - k ) { mn = k; mx = n - k; } else { mn = n - k; mx = k; } if ( mn < 0 ) { value = 0.0; } else if ( mn == 0 ) { value = 1.0; } else { value = ( double ) ( mx + 1 ); for ( i = 2; i <= mn; i++ ) { value = ( value * ( double ) ( mx + i ) ) / ( double ) i; } } return value; } /******************************************************************************/ double r8_mop ( int i ) /******************************************************************************/ /* Purpose: R8_MOP returns the I-th power of -1 as an R8 value. Discussion: An R8 is an double value. Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int I, the power of -1. Output, double R8_MOP, the I-th power of -1. */ { double value; if ( ( i % 2 ) == 0 ) { value = + 1.0; } else { value = - 1.0; } return value; } /******************************************************************************/ void r8vec_copy ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_COPY copies an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 03 July 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], the vector to be copied. Input, double A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ double r8vec_dot_product ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's. Licensing: This code is distributed under the MIT license. Modified: 26 July 2007 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], A2[N], the two vectors to be considered. Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a1[i] * a2[i]; } return value; } /******************************************************************************/ void subcomp_next ( int n, int k, int a[], int *more, int *h, int *t ) /******************************************************************************/ /* Purpose: SUBCOMP_NEXT computes the next subcomposition of 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 a value of N. A subcomposition of the integer N into K parts is a composition of M into K parts, where 0 <= M <= N. Licensing: This code is distributed under the MIT license. Modified: 27 March 2009 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 subcompositions are desired. Input, int K, the number of parts in the subcomposition. Input/output, int A[K], the parts of the subcomposition. Input/output, int *MORE, set by the user to start the computation, and by the routine to terminate it. 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; static int more2 = 0; static int n2 = 0; /* The first computation. */ if ( !( *more ) ) { n2 = 0; for ( i = 0; i < k; i++ ) { a[i] = 0; } more2 = 0; *h = 0; *t = 0; *more = 1; } /* Do the next element at the current value of N. */ else if ( more2 ) { comp_next ( n2, k, a, &more2, h, t ); } else { more2 = 0; n2 = n2 + 1; comp_next ( n2, k, a, &more2, h, t ); } /* Termination occurs if MORE2 = FALSE and N2 = N. */ if ( !more2 && n2 == n ) { *more = 0; } return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* 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: 24 September 2003 Author: John Burkardt Parameters: None */ { # 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 ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }