Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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,
00355                              unsigned n,
00356                              const realType* x,
00357                              unsigned m,
00358                              realType* J,
00359                              realType* D,
00360                              realType* work )
00361 {
00362     unsigned i, j;
00363     realType *w = work, *d = w + n, *u = d + n, *v = u + n, *up = v + n, *vp = up + n;
00364     for( i = 0; i < n; ++i )
00365     {
00366         realType ww = 1, zi = z[i];
00367         for( j = 0; j < i; ++j )
00368             ww *= zi - z[j];
00369         for( ++j; j < n; ++j )
00370             ww *= zi - z[j];
00371         w[i] = 1 / ww;
00372     }
00373     u[0] = v[n - 1] = 1;
00374     up[0] = vp[n - 1] = 0;
00375     for( i = 0; i < m; ++i )
00376     {
00377         realType xi = x[i];
00378         for( j = 0; j < n; ++j )
00379             d[j] = xi - z[j];
00380         for( j = 0; j < n - 1; ++j )
00381             u[j + 1] = d[j] * u[j], up[j + 1] = d[j] * up[j] + u[j];
00382         for( j = n - 1; j; --j )
00383             v[j - 1] = d[j] * v[j], vp[j - 1] = d[j] * vp[j] + v[j];
00384         for( j = 0; j < n; ++j )
00385             *J++ = w[j] * u[j] * v[j], *D++ = w[j] * ( up[j] * v[j] + u[j] * vp[j] );
00386     }
00387 }
00388 
00389 /*--------------------------------------------------------------------------
00390    Speedy Lagrangian Interpolation
00391 
00392    Usage:
00393 
00394      lagrange_data p;
00395      lagrange_setup(&p,z,n);    // setup for nodes z[0 ... n-1]
00396 
00397      the weights
00398        p->J [0 ... n-1]     interpolation weights
00399        p->D [0 ... n-1]     1st derivative weights
00400        p->D2[0 ... n-1]     2nd derivative weights
00401      are computed for a given x with:
00402        lagrange_0(p,x);  // compute p->J
00403        lagrange_1(p,x);  // compute p->J, p->D
00404        lagrange_2(p,x);  // compute p->J, p->D, p->D2
00405        lagrange_2u(p);   // compute p->D2 after call of lagrange_1(p,x);
00406      These functions use the z array supplied to setup
00407        (that pointer should not be freed between calls)
00408      Weights for x=z[0] and x=z[n-1] are computed during setup; access as:
00409        p->J_z0, etc. and p->J_zn, etc.
00410 
00411      lagrange_free(&p);  // deallocate memory allocated by setup
00412   --------------------------------------------------------------------------*/
00413 
00414 void lagrange_0( lagrange_data* p, realType x )
00415 {
00416     unsigned i, n = p->n;
00417     for( i = 0; i < n; ++i )
00418         p->d[i] = x - p->z[i];
00419     for( i = 0; i < n - 1; ++i )
00420         p->u0[i + 1] = p->d[i] * p->u0[i];
00421     for( i = n - 1; i; --i )
00422         p->v0[i - 1] = p->d[i] * p->v0[i];
00423     for( i = 0; i < n; ++i )
00424         p->J[i] = p->w[i] * p->u0[i] * p->v0[i];
00425 }
00426 
00427 void lagrange_1( lagrange_data* p, realType x )
00428 {
00429     unsigned i, n = p->n;
00430     for( i = 0; i < n; ++i )
00431         p->d[i] = x - p->z[i];
00432     for( i = 0; i < n - 1; ++i )
00433         p->u0[i + 1] = p->d[i] * p->u0[i], p->u1[i + 1] = p->d[i] * p->u1[i] + p->u0[i];
00434     for( i = n - 1; i; --i )
00435         p->v0[i - 1] = p->d[i] * p->v0[i], p->v1[i - 1] = p->d[i] * p->v1[i] + p->v0[i];
00436     for( i = 0; i < n; ++i )
00437         p->J[i] = p->w[i] * p->u0[i] * p->v0[i], p->D[i] = p->w[i] * ( p->u1[i] * p->v0[i] + p->u0[i] * p->v1[i] );
00438 }
00439 
00440 void lagrange_2( lagrange_data* p, realType x )
00441 {
00442     unsigned i, n = p->n;
00443     for( i = 0; i < n; ++i )
00444         p->d[i] = x - p->z[i];
00445     for( i = 0; i < n - 1; ++i )
00446         p->u0[i + 1] = p->d[i] * p->u0[i], p->u1[i + 1] = p->d[i] * p->u1[i] + p->u0[i],
00447                   p->u2[i + 1] = p->d[i] * p->u2[i] + 2 * p->u1[i];
00448     for( i = n - 1; i; --i )
00449         p->v0[i - 1] = p->d[i] * p->v0[i], p->v1[i - 1] = p->d[i] * p->v1[i] + p->v0[i],
00450                   p->v2[i - 1] = p->d[i] * p->v2[i] + 2 * p->v1[i];
00451     for( i = 0; i < n; ++i )
00452         p->J[i] = p->w[i] * p->u0[i] * p->v0[i], p->D[i] = p->w[i] * ( p->u1[i] * p->v0[i] + p->u0[i] * p->v1[i] ),
00453         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] );
00454 }
00455 
00456 void lagrange_2u( lagrange_data* p )
00457 {
00458     unsigned i, n = p->n;
00459     for( i = 0; i < n - 1; ++i )
00460         p->u2[i + 1] = p->d[i] * p->u2[i] + 2 * p->u1[i];
00461     for( i = n - 1; i; --i )
00462         p->v2[i - 1] = p->d[i] * p->v2[i] + 2 * p->v1[i];
00463     for( i = 0; i < n; ++i )
00464         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] );
00465 }
00466 
00467 void lagrange_setup( lagrange_data* p, const realType* z, unsigned n )
00468 {
00469     unsigned i, j;
00470     p->n = n, p->z = z;
00471     p->w = tmalloc( realType, 17 * n );
00472     p->d = p->w + n;
00473     p->J = p->d + n, p->D = p->J + n, p->D2 = p->D + n;
00474     p->u0 = p->D2 + n, p->v0 = p->u0 + n;
00475     p->u1 = p->v0 + n, p->v1 = p->u1 + n;
00476     p->u2 = p->v1 + n, p->v2 = p->u2 + n;
00477     p->J_z0 = p->v2 + n, p->D_z0 = p->J_z0 + n, p->D2_z0 = p->D_z0 + n;
00478     p->J_zn = p->D2_z0 + n, p->D_zn = p->J_zn + n, p->D2_zn = p->D_zn + n;
00479     for( i = 0; i < n; ++i )
00480     {
00481         realType ww = 1, zi = z[i];
00482         for( j = 0; j < i; ++j )
00483             ww *= zi - z[j];
00484         for( ++j; j < n; ++j )
00485             ww *= zi - z[j];
00486         p->w[i] = 1 / ww;
00487     }
00488     p->u0[0] = p->v0[n - 1] = 1;
00489     p->u1[0] = p->v1[n - 1] = 0;
00490     p->u2[0] = p->v2[n - 1] = 0;
00491     lagrange_2( p, z[0] );
00492     memcpy( p->J_z0, p->J, 3 * n * sizeof( realType ) );
00493     lagrange_2( p, z[n - 1] );
00494     memcpy( p->J_zn, p->J, 3 * n * sizeof( realType ) );
00495 }
00496 
00497 void lagrange_free( lagrange_data* p )
00498 {
00499     free( p->w );
00500 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines