MOAB: Mesh Oriented datABase  (version 5.3.0)
poly.c
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines