MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include <stdio.h> 00002 #include <stdlib.h> 00003 #include <stdarg.h> 00004 #include <math.h> /* for cos, fabs */ 00005 #include <string.h> /* for memcpy */ 00006 #include <float.h> 00007 00008 #include "moab/FindPtFuncs.h" 00009 00010 /* 00011 For brevity's sake, some names have been shortened 00012 Quadrature rules 00013 Gauss -> Gauss-Legendre quadrature (open) 00014 Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends) 00015 Polynomial bases 00016 Legendre -> Legendre basis 00017 Gauss -> Lagrangian basis using Gauss quadrature nodes 00018 Lobatto -> Lagrangian basis using Lobatto quadrature nodes 00019 */ 00020 00021 /*-------------------------------------------------------------------------- 00022 Legendre Polynomial Matrix Computation 00023 (compute P_i(x_j) for i = 0, ..., n and a given set of x) 00024 --------------------------------------------------------------------------*/ 00025 00026 /* precondition: n >= 1 00027 inner index is x index (0 ... m-1); 00028 outer index is Legendre polynomial number (0 ... n) 00029 */ 00030 void legendre_matrix( const realType* x, int m, realType* P, int n ) 00031 { 00032 int i, j; 00033 realType *Pjm1 = P, *Pj = Pjm1 + m, *Pjp1 = Pj + m; 00034 for( i = 0; i < m; ++i ) 00035 Pjm1[ i ] = 1; 00036 for( i = 0; i < m; ++i ) 00037 Pj[ i ] = x[ i ]; 00038 for( j = 1; j < n; ++j ) 00039 { 00040 realType c = 1 / ( realType )( j + 1 ), a = c * ( 2 * j + 1 ), b = c * j; 00041 for( i = 0; i < m; ++i ) 00042 Pjp1[ i ] = a * x[ i ] * Pj[ i ] - b * Pjm1[ i ]; 00043 Pjp1 += m, Pj += m, Pjm1 += m; 00044 } 00045 } 00046 00047 /* precondition: n >= 1, n even */ 00048 static void legendre_row_even( realType x, realType* P, int n ) 00049 { 00050 int i; 00051 P[ 0 ] = 1, P[ 1 ] = x; 00052 for( i = 1; i <= n - 2; i += 2 ) 00053 { 00054 P[ i + 1 ] = ( ( 2 * i + 1 ) * x * P[ i ] - i * P[ i - 1 ] ) / ( i + 1 ); 00055 P[ i + 2 ] = ( ( 2 * i + 3 ) * x * P[ i - 1 ] - ( i + 1 ) * P[ i ] ) / ( i + 2 ); 00056 } 00057 P[ n ] = ( ( 2 * n - 1 ) * x * P[ n - 1 ] - ( n - 1 ) * P[ n - 2 ] ) / n; 00058 } 00059 00060 /* precondition: n >= 1, n odd */ 00061 static void legendre_row_odd( realType x, realType* P, int n ) 00062 { 00063 int i; 00064 P[ 0 ] = 1, P[ 1 ] = x; 00065 for( i = 1; i <= n - 2; i += 2 ) 00066 { 00067 P[ i + 1 ] = ( ( 2 * i + 1 ) * x * P[ i ] - i * P[ i - 1 ] ) / ( i + 1 ); 00068 P[ i + 2 ] = ( ( 2 * i + 3 ) * x * P[ i - 1 ] - ( i + 1 ) * P[ i ] ) / ( i + 2 ); 00069 } 00070 } 00071 00072 /* precondition: n >= 1 00073 compute P_i(x) with i = 0 ... n 00074 */ 00075 void legendre_row( realType x, realType* P, int n ) 00076 { 00077 if( n & 1 ) 00078 legendre_row_odd( x, P, n ); 00079 else 00080 legendre_row_even( x, P, n ); 00081 } 00082 00083 /* precondition: n >= 1 00084 inner index is Legendre polynomial number (0 ... n) 00085 outer index is x index (0 ... m-1); 00086 */ 00087 void legendre_matrix_t( const realType* x, int m, realType* P, int n ) 00088 { 00089 int i; 00090 if( n & 1 ) 00091 for( i = 0; i < m; ++i, P += n + 1 ) 00092 legendre_row_odd( x[ i ], P, n ); 00093 else 00094 for( i = 0; i < m; ++i, P += n + 1 ) 00095 legendre_row_even( x[ i ], P, n ); 00096 } 00097 00098 /*-------------------------------------------------------------------------- 00099 Legendre Polynomial Computation 00100 compute P_n(x) or P_n'(x) or P_n''(x) 00101 --------------------------------------------------------------------------*/ 00102 00103 /* precondition: n >= 0 */ 00104 static realType legendre( int n, realType x ) 00105 { 00106 realType p[ 2 ]; 00107 int i; 00108 p[ 0 ] = 1; 00109 p[ 1 ] = x; 00110 for( i = 1; i < n; i += 2 ) 00111 { 00112 p[ 0 ] = ( ( 2 * i + 1 ) * x * p[ 1 ] - i * p[ 0 ] ) / ( i + 1 ); 00113 p[ 1 ] = ( ( 2 * i + 3 ) * x * p[ 0 ] - ( i + 1 ) * p[ 1 ] ) / ( i + 2 ); 00114 } 00115 return p[ n & 1 ]; 00116 } 00117 00118 /* precondition: n > 0 */ 00119 static realType legendre_d1( int n, realType x ) 00120 { 00121 realType p[ 2 ]; 00122 int i; 00123 p[ 0 ] = 3 * x; 00124 p[ 1 ] = 1; 00125 for( i = 2; i < n; i += 2 ) 00126 { 00127 p[ 1 ] = ( ( 2 * i + 1 ) * x * p[ 0 ] - ( i + 1 ) * p[ 1 ] ) / i; 00128 p[ 0 ] = ( ( 2 * i + 3 ) * x * p[ 1 ] - ( i + 2 ) * p[ 0 ] ) / ( i + 1 ); 00129 } 00130 return p[ n & 1 ]; 00131 } 00132 00133 /* precondition: n > 1 */ 00134 static realType legendre_d2( int n, realType x ) 00135 { 00136 realType p[ 2 ]; 00137 int i; 00138 p[ 0 ] = 3; 00139 p[ 1 ] = 15 * x; 00140 for( i = 3; i < n; i += 2 ) 00141 { 00142 p[ 0 ] = ( ( 2 * i + 1 ) * x * p[ 1 ] - ( i + 2 ) * p[ 0 ] ) / ( i - 1 ); 00143 p[ 1 ] = ( ( 2 * i + 3 ) * x * p[ 0 ] - ( i + 3 ) * p[ 1 ] ) / i; 00144 } 00145 return p[ n & 1 ]; 00146 } 00147 00148 /*-------------------------------------------------------------------------- 00149 Quadrature Nodes and Weights Calculation 00150 compute the n Gauss-Legendre nodes and weights or 00151 the n Gauss-Lobatto-Legendre nodes and weights 00152 --------------------------------------------------------------------------*/ 00153 00154 /* n nodes */ 00155 void gauss_nodes( realType* z, int n ) 00156 { 00157 int i, j; 00158 for( i = 0; i <= n / 2 - 1; ++i ) 00159 { 00160 realType ox, x = mbcos( ( 2 * n - 2 * i - 1 ) * ( MOAB_POLY_PI / 2 ) / n ); 00161 do 00162 { 00163 ox = x; 00164 x -= legendre( n, x ) / legendre_d1( n, x ); 00165 } while( mbabs( x - ox ) > -x * MOAB_POLY_EPS ); 00166 z[ i ] = x - legendre( n, x ) / legendre_d1( n, x ); 00167 } 00168 if( n & 1 ) z[ n / 2 ] = 0; 00169 for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i ) 00170 z[ j ] = -z[ i ]; 00171 } 00172 00173 /* n inner lobatto nodes (excluding -1,1) */ 00174 static void lobatto_nodes_aux( realType* z, int n ) 00175 { 00176 int i, j, np = n + 1; 00177 for( i = 0; i <= n / 2 - 1; ++i ) 00178 { 00179 realType ox, x = mbcos( ( n - i ) * MOAB_POLY_PI / np ); 00180 do 00181 { 00182 ox = x; 00183 x -= legendre_d1( np, x ) / legendre_d2( np, x ); 00184 } while( mbabs( x - ox ) > -x * MOAB_POLY_EPS ); 00185 z[ i ] = x - legendre_d1( np, x ) / legendre_d2( np, x ); 00186 } 00187 if( n & 1 ) z[ n / 2 ] = 0; 00188 for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i ) 00189 z[ j ] = -z[ i ]; 00190 } 00191 00192 /* n lobatto nodes */ 00193 void lobatto_nodes( realType* z, int n ) 00194 { 00195 z[ 0 ] = -1, z[ n - 1 ] = 1; 00196 lobatto_nodes_aux( &z[ 1 ], n - 2 ); 00197 } 00198 00199 void gauss_weights( const realType* z, realType* w, int n ) 00200 { 00201 int i, j; 00202 for( i = 0; i <= ( n - 1 ) / 2; ++i ) 00203 { 00204 realType d = ( n + 1 ) * legendre( n + 1, z[ i ] ); 00205 w[ i ] = 2 * ( 1 - z[ i ] * z[ i ] ) / ( d * d ); 00206 } 00207 for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i ) 00208 w[ j ] = w[ i ]; 00209 } 00210 00211 void lobatto_weights( const realType* z, realType* w, int n ) 00212 { 00213 int i, j; 00214 for( i = 0; i <= ( n - 1 ) / 2; ++i ) 00215 { 00216 realType d = legendre( n - 1, z[ i ] ); 00217 w[ i ] = 2 / ( ( n - 1 ) * n * d * d ); 00218 } 00219 for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i ) 00220 w[ j ] = w[ i ]; 00221 } 00222 00223 /*-------------------------------------------------------------------------- 00224 Lagrangian to Legendre Change-of-basis Matrix 00225 where the nodes of the Lagrangian basis are the GL or GLL quadrature nodes 00226 --------------------------------------------------------------------------*/ 00227 00228 /* precondition: n >= 2 00229 given the Gauss quadrature rule (z,w,n), compute the square matrix J 00230 for transforming from the Gauss basis to the Legendre basis: 00231 00232 u_legendre(i) = sum_j J(i,j) u_gauss(j) 00233 00234 computes J = .5 (2i+1) w P (z ) 00235 ij j i j 00236 00237 in column major format (inner index is i, the Legendre index) 00238 */ 00239 void gauss_to_legendre( const realType* z, const realType* w, int n, realType* J ) 00240 { 00241 int i, j; 00242 legendre_matrix_t( z, n, J, n - 1 ); 00243 for( j = 0; j < n; ++j ) 00244 { 00245 realType ww = w[ j ]; 00246 for( i = 0; i < n; ++i ) 00247 *J++ *= ( 2 * i + 1 ) * ww / 2; 00248 } 00249 } 00250 00251 /* precondition: n >= 2 00252 same as above, but 00253 in row major format (inner index is j, the Gauss index) 00254 */ 00255 void gauss_to_legendre_t( const realType* z, const realType* w, int n, realType* J ) 00256 { 00257 int i, j; 00258 legendre_matrix( z, n, J, n - 1 ); 00259 for( i = 0; i < n; ++i ) 00260 { 00261 realType ii = ( realType )( 2 * i + 1 ) / 2; 00262 for( j = 0; j < n; ++j ) 00263 *J++ *= ii * w[ j ]; 00264 } 00265 } 00266 00267 /* precondition: n >= 3 00268 given the Lobatto quadrature rule (z,w,n), compute the square matrix J 00269 for transforming from the Gauss basis to the Legendre basis: 00270 00271 u_legendre(i) = sum_j J(i,j) u_lobatto(j) 00272 00273 in column major format (inner index is i, the Legendre index) 00274 */ 00275 void lobatto_to_legendre( const realType* z, const realType* w, int n, realType* J ) 00276 { 00277 int i, j, m = ( n + 1 ) / 2; 00278 realType *p = J, *q; 00279 realType ww, sum; 00280 if( n & 1 ) 00281 for( j = 0; j < m; ++j, p += n ) 00282 legendre_row_odd( z[ j ], p, n - 2 ); 00283 else 00284 for( j = 0; j < m; ++j, p += n ) 00285 legendre_row_even( z[ j ], p, n - 2 ); 00286 p = J; 00287 for( j = 0; j < m; ++j ) 00288 { 00289 ww = w[ j ], sum = 0; 00290 for( i = 0; i < n - 1; ++i ) 00291 *p *= ( 2 * i + 1 ) * ww / 2, sum += *p++; 00292 *p++ = -sum; 00293 } 00294 q = J + ( n / 2 - 1 ) * n; 00295 if( n & 1 ) 00296 for( ; j < n; ++j, p += n, q -= n ) 00297 { 00298 for( i = 0; i < n - 1; i += 2 ) 00299 p[ i ] = q[ i ], p[ i + 1 ] = -q[ i + 1 ]; 00300 p[ i ] = q[ i ]; 00301 } 00302 else 00303 for( ; j < n; ++j, p += n, q -= n ) 00304 { 00305 for( i = 0; i < n - 1; i += 2 ) 00306 p[ i ] = q[ i ], p[ i + 1 ] = -q[ i + 1 ]; 00307 } 00308 } 00309 00310 /*-------------------------------------------------------------------------- 00311 Lagrangian to Lagrangian change-of-basis matrix, and derivative matrix 00312 --------------------------------------------------------------------------*/ 00313 00314 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 00315 evaluate all Lagrangian basis functions at all points x 00316 00317 inner index of output J is the basis function index (row-major format) 00318 provide work array with space for 4*n reals 00319 */ 00320 void lagrange_weights( const realType* z, unsigned n, const realType* x, unsigned m, realType* J, realType* work ) 00321 { 00322 unsigned i, j; 00323 realType *w = work, *d = w + n, *u = d + n, *v = u + n; 00324 for( i = 0; i < n; ++i ) 00325 { 00326 realType ww = 1, zi = z[ i ]; 00327 for( j = 0; j < i; ++j ) 00328 ww *= zi - z[ j ]; 00329 for( ++j; j < n; ++j ) 00330 ww *= zi - z[ j ]; 00331 w[ i ] = 1 / ww; 00332 } 00333 u[ 0 ] = v[ n - 1 ] = 1; 00334 for( i = 0; i < m; ++i ) 00335 { 00336 realType xi = x[ i ]; 00337 for( j = 0; j < n; ++j ) 00338 d[ j ] = xi - z[ j ]; 00339 for( j = 0; j < n - 1; ++j ) 00340 u[ j + 1 ] = d[ j ] * u[ j ]; 00341 for( j = n - 1; j; --j ) 00342 v[ j - 1 ] = d[ j ] * v[ j ]; 00343 for( j = 0; j < n; ++j ) 00344 *J++ = w[ j ] * u[ j ] * v[ j ]; 00345 } 00346 } 00347 00348 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 00349 evaluate all Lagrangian basis functions and their derivatives 00350 00351 inner index of outputs J,D is the basis function index (row-major format) 00352 provide work array with space for 6*n reals 00353 */ 00354 void lagrange_weights_deriv( const realType* z, unsigned n, const realType* x, unsigned m, realType* J, realType* D, 00355 realType* work ) 00356 { 00357 unsigned i, j; 00358 realType *w = work, *d = w + n, *u = d + n, *v = u + n, *up = v + n, *vp = up + n; 00359 for( i = 0; i < n; ++i ) 00360 { 00361 realType ww = 1, zi = z[ i ]; 00362 for( j = 0; j < i; ++j ) 00363 ww *= zi - z[ j ]; 00364 for( ++j; j < n; ++j ) 00365 ww *= zi - z[ j ]; 00366 w[ i ] = 1 / ww; 00367 } 00368 u[ 0 ] = v[ n - 1 ] = 1; 00369 up[ 0 ] = vp[ n - 1 ] = 0; 00370 for( i = 0; i < m; ++i ) 00371 { 00372 realType xi = x[ i ]; 00373 for( j = 0; j < n; ++j ) 00374 d[ j ] = xi - z[ j ]; 00375 for( j = 0; j < n - 1; ++j ) 00376 u[ j + 1 ] = d[ j ] * u[ j ], up[ j + 1 ] = d[ j ] * up[ j ] + u[ j ]; 00377 for( j = n - 1; j; --j ) 00378 v[ j - 1 ] = d[ j ] * v[ j ], vp[ j - 1 ] = d[ j ] * vp[ j ] + v[ j ]; 00379 for( j = 0; j < n; ++j ) 00380 *J++ = w[ j ] * u[ j ] * v[ j ], *D++ = w[ j ] * ( up[ j ] * v[ j ] + u[ j ] * vp[ j ] ); 00381 } 00382 } 00383 00384 /*-------------------------------------------------------------------------- 00385 Speedy Lagrangian Interpolation 00386 00387 Usage: 00388 00389 lagrange_data p; 00390 lagrange_setup(&p,z,n); // setup for nodes z[0 ... n-1] 00391 00392 the weights 00393 p->J [0 ... n-1] interpolation weights 00394 p->D [0 ... n-1] 1st derivative weights 00395 p->D2[0 ... n-1] 2nd derivative weights 00396 are computed for a given x with: 00397 lagrange_0(p,x); // compute p->J 00398 lagrange_1(p,x); // compute p->J, p->D 00399 lagrange_2(p,x); // compute p->J, p->D, p->D2 00400 lagrange_2u(p); // compute p->D2 after call of lagrange_1(p,x); 00401 These functions use the z array supplied to setup 00402 (that pointer should not be freed between calls) 00403 Weights for x=z[0] and x=z[n-1] are computed during setup; access as: 00404 p->J_z0, etc. and p->J_zn, etc. 00405 00406 lagrange_free(&p); // deallocate memory allocated by setup 00407 --------------------------------------------------------------------------*/ 00408 00409 void lagrange_0( lagrange_data* p, realType x ) 00410 { 00411 unsigned i, n = p->n; 00412 for( i = 0; i < n; ++i ) 00413 p->d[ i ] = x - p->z[ i ]; 00414 for( i = 0; i < n - 1; ++i ) 00415 p->u0[ i + 1 ] = p->d[ i ] * p->u0[ i ]; 00416 for( i = n - 1; i; --i ) 00417 p->v0[ i - 1 ] = p->d[ i ] * p->v0[ i ]; 00418 for( i = 0; i < n; ++i ) 00419 p->J[ i ] = p->w[ i ] * p->u0[ i ] * p->v0[ i ]; 00420 } 00421 00422 void lagrange_1( lagrange_data* p, realType x ) 00423 { 00424 unsigned i, n = p->n; 00425 for( i = 0; i < n; ++i ) 00426 p->d[ i ] = x - p->z[ i ]; 00427 for( i = 0; i < n - 1; ++i ) 00428 p->u0[ i + 1 ] = p->d[ i ] * p->u0[ i ], p->u1[ i + 1 ] = p->d[ i ] * p->u1[ i ] + p->u0[ i ]; 00429 for( i = n - 1; i; --i ) 00430 p->v0[ i - 1 ] = p->d[ i ] * p->v0[ i ], p->v1[ i - 1 ] = p->d[ i ] * p->v1[ i ] + p->v0[ i ]; 00431 for( i = 0; i < n; ++i ) 00432 p->J[ i ] = p->w[ i ] * p->u0[ i ] * p->v0[ i ], 00433 p->D[ i ] = p->w[ i ] * ( p->u1[ i ] * p->v0[ i ] + p->u0[ i ] * p->v1[ i ] ); 00434 } 00435 00436 void lagrange_2( lagrange_data* p, realType x ) 00437 { 00438 unsigned i, n = p->n; 00439 for( i = 0; i < n; ++i ) 00440 p->d[ i ] = x - p->z[ i ]; 00441 for( i = 0; i < n - 1; ++i ) 00442 p->u0[ i + 1 ] = p->d[ i ] * p->u0[ i ], p->u1[ i + 1 ] = p->d[ i ] * p->u1[ i ] + p->u0[ i ], 00443 p->u2[ i + 1 ] = p->d[ i ] * p->u2[ i ] + 2 * p->u1[ i ]; 00444 for( i = n - 1; i; --i ) 00445 p->v0[ i - 1 ] = p->d[ i ] * p->v0[ i ], p->v1[ i - 1 ] = p->d[ i ] * p->v1[ i ] + p->v0[ i ], 00446 p->v2[ i - 1 ] = p->d[ i ] * p->v2[ i ] + 2 * p->v1[ i ]; 00447 for( i = 0; i < n; ++i ) 00448 p->J[ i ] = p->w[ i ] * p->u0[ i ] * p->v0[ i ], 00449 p->D[ i ] = p->w[ i ] * ( p->u1[ i ] * p->v0[ i ] + p->u0[ i ] * p->v1[ i ] ), 00450 p->D2[ i ] = 00451 p->w[ i ] * ( p->u2[ i ] * p->v0[ i ] + 2 * p->u1[ i ] * p->v1[ i ] + p->u0[ i ] * p->v2[ i ] ); 00452 } 00453 00454 void lagrange_2u( lagrange_data* p ) 00455 { 00456 unsigned i, n = p->n; 00457 for( i = 0; i < n - 1; ++i ) 00458 p->u2[ i + 1 ] = p->d[ i ] * p->u2[ i ] + 2 * p->u1[ i ]; 00459 for( i = n - 1; i; --i ) 00460 p->v2[ i - 1 ] = p->d[ i ] * p->v2[ i ] + 2 * p->v1[ i ]; 00461 for( i = 0; i < n; ++i ) 00462 p->D2[ i ] = p->w[ i ] * ( p->u2[ i ] * p->v0[ i ] + 2 * p->u1[ i ] * p->v1[ i ] + p->u0[ i ] * p->v2[ i ] ); 00463 } 00464 00465 void lagrange_setup( lagrange_data* p, const realType* z, unsigned n ) 00466 { 00467 unsigned i, j; 00468 p->n = n, p->z = z; 00469 p->w = tmalloc( realType, 17 * n ); 00470 p->d = p->w + n; 00471 p->J = p->d + n, p->D = p->J + n, p->D2 = p->D + n; 00472 p->u0 = p->D2 + n, p->v0 = p->u0 + n; 00473 p->u1 = p->v0 + n, p->v1 = p->u1 + n; 00474 p->u2 = p->v1 + n, p->v2 = p->u2 + n; 00475 p->J_z0 = p->v2 + n, p->D_z0 = p->J_z0 + n, p->D2_z0 = p->D_z0 + n; 00476 p->J_zn = p->D2_z0 + n, p->D_zn = p->J_zn + n, p->D2_zn = p->D_zn + n; 00477 for( i = 0; i < n; ++i ) 00478 { 00479 realType ww = 1, zi = z[ i ]; 00480 for( j = 0; j < i; ++j ) 00481 ww *= zi - z[ j ]; 00482 for( ++j; j < n; ++j ) 00483 ww *= zi - z[ j ]; 00484 p->w[ i ] = 1 / ww; 00485 } 00486 p->u0[ 0 ] = p->v0[ n - 1 ] = 1; 00487 p->u1[ 0 ] = p->v1[ n - 1 ] = 0; 00488 p->u2[ 0 ] = p->v2[ n - 1 ] = 0; 00489 lagrange_2( p, z[ 0 ] ); 00490 memcpy( p->J_z0, p->J, 3 * n * sizeof( realType ) ); 00491 lagrange_2( p, z[ n - 1 ] ); 00492 memcpy( p->J_zn, p->J, 3 * n * sizeof( realType ) ); 00493 } 00494 00495 void lagrange_free( lagrange_data* p ) 00496 { 00497 free( p->w ); 00498 }