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
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 }