# include # include # include # include # include "cube_exactness.h" /******************************************************************************/ void legendre_3d_exactness ( double a[], double b[], int n, double x[], double y[], double z[], double w[], int t ) /******************************************************************************/ /* Purpose: LEGENDRE_3D_EXACTNESS: monomial exactness for the 3D Legendre integral. Licensing: This code is distributed under the MIT license. Modified: 16 August 2014 Author: John Burkardt Parameters: Input, double A[3], the lower limits of integration. Input, double B[3], the upper limits of integration. Input, int N, the number of points in the rule. Input, double X[N], Y[N], Z[N], the quadrature points. Input, double W[N], the quadrature weights. Input, int T, the maximum total degree. 0 <= T. */ { double e; int i; int j; int k; int l; int p[3]; double q; double s; int tt; double *v; v = ( double * ) malloc ( n * sizeof ( double ) ); printf ( "\n" ); printf ( " Quadrature rule for the 3D Legendre integral.\n" ); printf ( " Number of points in rule is %d\n", n ); printf ( "\n" ); printf ( " D I J K Relative Error\n" ); for ( tt = 0; tt <= t; tt++ ) { printf ( " %2d\n", tt ); for ( k = 0; k <= tt; k++ ) { for ( j = 0; j <= tt - k; j++ ) { i = tt - j - k; p[0] = i; p[1] = j; p[2] = k; s = legendre_3d_monomial_integral ( a, b, p ); for ( l = 0; l < n; l++ ) { v[l] = pow ( x[l], p[0] ) * pow ( y[l], p[1] ) * pow ( z[l], p[2] ); } q = r8vec_dot_product ( n, w, v ); if ( s == 0.0 ) { e = fabs ( q ); } else { e = fabs ( q - s ) / fabs ( s ); } printf ( " %6d %6d %6d %24.16f\n", p[0], p[1], p[2], e ); } } } free ( v ); return; } /******************************************************************************/ double legendre_3d_monomial_integral ( double a[], double b[], int p[] ) /******************************************************************************/ /* Purpose: LEGENDRE_3D_MONOMIAL_INTEGRAL the Legendre integral of a monomial. Discussion: The Legendre integral to be evaluated has the form I(f) = integral ( z1 <= z <= z2 ) integral ( y1 <= y <= y2 ) integral ( x1 <= x <= x2 ) x^i y^j z^k dx dy dz Licensing: This code is distributed under the MIT license. Modified: 16 August 2014 Author: John Burkardt Parameters: Input, double A[3], the lower limits of integration. Input, double B[3], the upper limits of integration. Input, int P[3], the exponents of X and Y. Output, double LEGENDRE_3D_MONOMIAL_INTEGRAL, the value of the exact integral. */ { double value; value = ( pow ( b[0], p[0] + 1 ) - pow ( a[0], p[0] + 1 ) ) / ( double ) ( p[0] + 1 ) * ( pow ( b[1], p[1] + 1 ) - pow ( a[1], p[1] + 1 ) ) / ( double ) ( p[1] + 1 ) * ( pow ( b[2], p[2] + 1 ) - pow ( a[2], p[2] + 1 ) ) / ( double ) ( p[2] + 1 ); return value; } /******************************************************************************/ 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 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 }