C 代码 返回用于近似积分的费利帕正交规则 在 3D 立方体的内部.rar

# include <stdlib.h> # include <stdio.h> # include <time.h> # include <string.h> # include <math.h> # include "cube_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 GNU LGPL 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 cube_monomial ( double a[], double b[], int expon[3] ) /******************************************************************************/ /* Purpose: CUBE_MONOMIAL integrates a monomial over a cube in 3D. Discussion: This routine integrates a monomial of the form product ( 1 <= dim <= 3 ) x(dim)^expon(dim) The combination 0^0 should be treated as 1. The integration region is: A(1) <= X <= B(1) A(2) <= Y <= B(2) A(3) <= Z <= B(3) Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2014 Author: John Burkardt Parameters: Input, double A[3], B[3], the lower and upper limits. Input, int EXPON[3], the exponents. Output, double CUBE_MONOMIAL, the integral of the monomial. */ { int i; double value; for ( i = 0; i < 3; i++ ) { if ( expon[i] == -1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CUBE_MONOMIAL - Fatal error!\n" ); fprintf ( stderr, " Exponent of -1 encountered.\n" ); exit ( 1 ); } } value = 1.0; for ( i = 0; i < 3; i++ ) { value = value * ( pow ( b[i], expon[i] + 1 ) - pow ( a[i], expon[i] + 1 ) ) / ( double ) ( expon[i] + 1 ); } return value; } /******************************************************************************/ void cube_monomial_test ( int degree_max ) /******************************************************************************/ /* Purpose: CUBE_MONOMIAL_TEST tests CUBE_MONOMIAL. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2014 Author: John Burkardt Parameters: Input, int DEGREE_MAX, the maximum total degree of the monomials to check. */ { double a[3] = { -1.0, -1.0, -1.0 }; int alpha; double b[3] = { +1.0, +1.0, +1.0 }; int beta; int expon[3]; int gamma; double value; printf ( "\n" ); printf ( "CUBE_MONOMIAL_TEST\n" ); printf ( " For the unit hexahedron,\n" ); printf ( " CUBE_MONOMIAL returns the exact value of the\n" ); printf ( " integral of X^ALPHA Y^BETA Z^GAMMA\n" ); printf ( "\n" ); printf ( " Volume = %g\n", cube_volume ( a, b ) ); printf ( "\n" ); printf ( " ALPHA BETA GAMMA INTEGRAL\n" ); printf ( "\n" ); for ( alpha = 0; alpha <= degree_max; alpha++ ) { expon[0] = alpha; for ( beta = 0; beta <= degree_max - alpha; beta++ ) { expon[1] = beta; for ( gamma = 0; gamma <= degree_max - alpha - beta; gamma++ ) { expon[2] = gamma; value = cube_monomial ( a, b, expon ); printf ( " %8d %8d %8d %14.g\n", expon[0], expon[1], expon[2], value ); } } } return; } /******************************************************************************/ void cube_quad_test ( int degree_max ) /******************************************************************************/ /* Purpose: CUBE_QUAD_TEST tests the rules for a cube in 3D. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2014 Author: John Burkardt Parameters: Input, int DEGREE_MAX, the maximum total degree of the monomials to check. */ { # define DIM_NUM 3 double a[3] = { -1.0, -1.0, -1.0 }; double b[3] = { +1.0, +1.0, +1.0 }; int dim; int dim_num = DIM_NUM; int expon[DIM_NUM]; int h; int k; int more; int order; int order_1d[DIM_NUM]; double quad; int t; double *v; double *w; double *xyz; printf ( "\n" ); printf ( "CUBE_QUAD_TEST\n" ); printf ( " For the unit hexahedron,\n" ); printf ( " we approximate monomial integrals with:\n" ); printf ( " CUBE_RULE, which returns N1 by N2 by N3 point rules..\n" ); more = 0; for ( ; ; ) { subcomp_next ( degree_max, dim_num, expon, &more, &h, &t ); printf ( "\n" ); printf ( " Monomial exponents: " ); for ( dim = 0; dim < dim_num; dim++ ) { printf (