MOAB: Mesh Oriented datABase  (version 5.2.1)
poly.c File Reference
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include "moab/FindPtFuncs.h"
+ Include dependency graph for poly.c:

Go to the source code of this file.

Functions

void legendre_matrix (const realType *x, int m, realType *P, int n)
static void legendre_row_even (realType x, realType *P, int n)
static void legendre_row_odd (realType x, realType *P, int n)
void legendre_row (realType x, realType *P, int n)
void legendre_matrix_t (const realType *x, int m, realType *P, int n)
static realType legendre (int n, realType x)
static realType legendre_d1 (int n, realType x)
static realType legendre_d2 (int n, realType x)
void gauss_nodes (realType *z, int n)
static void lobatto_nodes_aux (realType *z, int n)
void lobatto_nodes (realType *z, int n)
void gauss_weights (const realType *z, realType *w, int n)
void lobatto_weights (const realType *z, realType *w, int n)
void gauss_to_legendre (const realType *z, const realType *w, int n, realType *J)
void gauss_to_legendre_t (const realType *z, const realType *w, int n, realType *J)
void lobatto_to_legendre (const realType *z, const realType *w, int n, realType *J)
void lagrange_weights (const realType *z, unsigned n, const realType *x, unsigned m, realType *J, realType *work)
void lagrange_weights_deriv (const realType *z, unsigned n, const realType *x, unsigned m, realType *J, realType *D, realType *work)
void lagrange_0 (lagrange_data *p, realType x)
void lagrange_1 (lagrange_data *p, realType x)
void lagrange_2 (lagrange_data *p, realType x)
void lagrange_2u (lagrange_data *p)
void lagrange_setup (lagrange_data *p, const realType *z, unsigned n)
void lagrange_free (lagrange_data *p)

Function Documentation

void gauss_nodes ( realType z,
int  n 
)

Definition at line 155 of file poly.c.

References legendre(), legendre_d1(), mbabs, mbcos, MOAB_POLY_EPS, MOAB_POLY_PI, and n.

{
    int i, j;
    for( i = 0; i <= n / 2 - 1; ++i )
    {
        realType ox, x = mbcos( ( 2 * n - 2 * i - 1 ) * ( MOAB_POLY_PI / 2 ) / n );
        do
        {
            ox = x;
            x -= legendre( n, x ) / legendre_d1( n, x );
        } while( mbabs( x - ox ) > -x * MOAB_POLY_EPS );
        z[ i ] = x - legendre( n, x ) / legendre_d1( n, x );
    }
    if( n & 1 ) z[ n / 2 ] = 0;
    for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
        z[ j ] = -z[ i ];
}
void gauss_to_legendre ( const realType z,
const realType w,
int  n,
realType J 
)

Definition at line 239 of file poly.c.

References legendre_matrix_t(), and n.

{
    int i, j;
    legendre_matrix_t( z, n, J, n - 1 );
    for( j = 0; j < n; ++j )
    {
        realType ww = w[ j ];
        for( i = 0; i < n; ++i )
            *J++ *= ( 2 * i + 1 ) * ww / 2;
    }
}
void gauss_to_legendre_t ( const realType z,
const realType w,
int  n,
realType J 
)

Definition at line 255 of file poly.c.

References legendre_matrix(), and n.

{
    int i, j;
    legendre_matrix( z, n, J, n - 1 );
    for( i = 0; i < n; ++i )
    {
        realType ii = ( realType )( 2 * i + 1 ) / 2;
        for( j = 0; j < n; ++j )
            *J++ *= ii * w[ j ];
    }
}
void gauss_weights ( const realType z,
realType w,
int  n 
)

Definition at line 199 of file poly.c.

References legendre(), and n.

{
    int i, j;
    for( i = 0; i <= ( n - 1 ) / 2; ++i )
    {
        realType d = ( n + 1 ) * legendre( n + 1, z[ i ] );
        w[ i ] = 2 * ( 1 - z[ i ] * z[ i ] ) / ( d * d );
    }
    for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
        w[ j ] = w[ i ];
}
void lagrange_1 ( lagrange_data p,
realType  x 
)

Definition at line 422 of file poly.c.

References lagrange_data::D, lagrange_data::d, lagrange_data::J, n, lagrange_data::n, lagrange_data::u0, lagrange_data::u1, lagrange_data::v0, lagrange_data::v1, lagrange_data::w, and lagrange_data::z.

Referenced by opt_area_set_2(), opt_edge_set_2(), opt_edge_set_3(), opt_face_set_3(), and opt_vol_set_3().

{
    unsigned i, n = p->n;
    for( i = 0; i < n; ++i )
        p->d[ i ] = x - p->z[ i ];
    for( i = 0; i < n - 1; ++i )
        p->u0[ i + 1 ] = p->d[ i ] * p->u0[ i ], p->u1[ i + 1 ] = p->d[ i ] * p->u1[ i ] + p->u0[ i ];
    for( i = n - 1; i; --i )
        p->v0[ i - 1 ] = p->d[ i ] * p->v0[ i ], p->v1[ i - 1 ] = p->d[ i ] * p->v1[ i ] + p->v0[ i ];
    for( i = 0; i < n; ++i )
        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 ] );
}
void lagrange_2 ( lagrange_data p,
realType  x 
)

Definition at line 436 of file poly.c.

References lagrange_data::D, lagrange_data::d, lagrange_data::D2, lagrange_data::J, n, lagrange_data::n, lagrange_data::u0, lagrange_data::u1, lagrange_data::u2, lagrange_data::v0, lagrange_data::v1, lagrange_data::v2, lagrange_data::w, and lagrange_data::z.

Referenced by lagrange_setup().

{
    unsigned i, n = p->n;
    for( i = 0; i < n; ++i )
        p->d[ i ] = x - p->z[ i ];
    for( i = 0; i < n - 1; ++i )
        p->u0[ i + 1 ] = p->d[ i ] * p->u0[ i ], p->u1[ i + 1 ] = p->d[ i ] * p->u1[ i ] + p->u0[ i ],
                     p->u2[ i + 1 ] = p->d[ i ] * p->u2[ i ] + 2 * p->u1[ i ];
    for( i = n - 1; i; --i )
        p->v0[ i - 1 ] = p->d[ i ] * p->v0[ i ], p->v1[ i - 1 ] = p->d[ i ] * p->v1[ i ] + p->v0[ i ],
                     p->v2[ i - 1 ] = p->d[ i ] * p->v2[ i ] + 2 * p->v1[ i ];
    for( i = 0; i < n; ++i )
        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 ] ),
                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 ] );
}
void lagrange_2u ( lagrange_data p)

Definition at line 454 of file poly.c.

References lagrange_data::d, lagrange_data::D2, n, lagrange_data::n, lagrange_data::u0, lagrange_data::u1, lagrange_data::u2, lagrange_data::v0, lagrange_data::v1, lagrange_data::v2, and lagrange_data::w.

Referenced by opt_edge_hess_2(), opt_edge_hess_3(), and opt_face_hess_3().

{
    unsigned i, n = p->n;
    for( i = 0; i < n - 1; ++i )
        p->u2[ i + 1 ] = p->d[ i ] * p->u2[ i ] + 2 * p->u1[ i ];
    for( i = n - 1; i; --i )
        p->v2[ i - 1 ] = p->d[ i ] * p->v2[ i ] + 2 * p->v1[ i ];
    for( i = 0; i < n; ++i )
        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 ] );
}
void lagrange_setup ( lagrange_data p,
const realType z,
unsigned  n 
)

Definition at line 465 of file poly.c.

References lagrange_data::D, lagrange_data::d, lagrange_data::D2, lagrange_data::D2_z0, lagrange_data::D2_zn, lagrange_data::D_z0, lagrange_data::D_zn, lagrange_data::J, lagrange_data::J_z0, lagrange_data::J_zn, lagrange_2(), n, lagrange_data::n, tmalloc, lagrange_data::u0, lagrange_data::u1, lagrange_data::u2, lagrange_data::v0, lagrange_data::v1, lagrange_data::v2, lagrange_data::w, z, and lagrange_data::z.

Referenced by findpt_setup_2(), findpt_setup_3(), moab::ElemUtil::hex_eval(), moab::ElemUtil::hex_findpt(), moab::Element::SpectralHex::Init(), moab::Element::SpectralQuad::Init(), moab::SpectralHex::initFcn(), and moab::element_utility::Spectral_hex_map< moab::Matrix3 >::initialize_spectral_hex().

{
    unsigned i, j;
    p->n = n, p->z = z;
    p->w = tmalloc( realType, 17 * n );
    p->d = p->w + n;
    p->J = p->d + n, p->D = p->J + n, p->D2 = p->D + n;
    p->u0 = p->D2 + n, p->v0 = p->u0 + n;
    p->u1 = p->v0 + n, p->v1 = p->u1 + n;
    p->u2 = p->v1 + n, p->v2 = p->u2 + n;
    p->J_z0 = p->v2 + n, p->D_z0 = p->J_z0 + n, p->D2_z0 = p->D_z0 + n;
    p->J_zn = p->D2_z0 + n, p->D_zn = p->J_zn + n, p->D2_zn = p->D_zn + n;
    for( i = 0; i < n; ++i )
    {
        realType ww = 1, zi = z[ i ];
        for( j = 0; j < i; ++j )
            ww *= zi - z[ j ];
        for( ++j; j < n; ++j )
            ww *= zi - z[ j ];
        p->w[ i ] = 1 / ww;
    }
    p->u0[ 0 ] = p->v0[ n - 1 ] = 1;
    p->u1[ 0 ] = p->v1[ n - 1 ] = 0;
    p->u2[ 0 ] = p->v2[ n - 1 ] = 0;
    lagrange_2( p, z[ 0 ] );
    memcpy( p->J_z0, p->J, 3 * n * sizeof( realType ) );
    lagrange_2( p, z[ n - 1 ] );
    memcpy( p->J_zn, p->J, 3 * n * sizeof( realType ) );
}
void lagrange_weights ( const realType z,
unsigned  n,
const realType x,
unsigned  m,
realType J,
realType work 
)

Definition at line 320 of file poly.c.

References n, u, and MBMesquite::xi.

{
    unsigned  i, j;
    realType *w = work, *d = w + n, *u = d + n, *v = u + n;
    for( i = 0; i < n; ++i )
    {
        realType ww = 1, zi = z[ i ];
        for( j = 0; j < i; ++j )
            ww *= zi - z[ j ];
        for( ++j; j < n; ++j )
            ww *= zi - z[ j ];
        w[ i ] = 1 / ww;
    }
    u[ 0 ] = v[ n - 1 ] = 1;
    for( i = 0; i < m; ++i )
    {
        realType xi = x[ i ];
        for( j = 0; j < n; ++j )
            d[ j ] = xi - z[ j ];
        for( j = 0; j < n - 1; ++j )
            u[ j + 1 ] = d[ j ] * u[ j ];
        for( j = n - 1; j; --j )
            v[ j - 1 ] = d[ j ] * v[ j ];
        for( j = 0; j < n; ++j )
            *J++ = w[ j ] * u[ j ] * v[ j ];
    }
}
void lagrange_weights_deriv ( const realType z,
unsigned  n,
const realType x,
unsigned  m,
realType J,
realType D,
realType work 
)

Definition at line 354 of file poly.c.

References n, u, and MBMesquite::xi.

Referenced by lob_bnd_base_setup(), obbox_setup_2(), and obbox_setup_3().

{
    unsigned  i, j;
    realType *w = work, *d = w + n, *u = d + n, *v = u + n, *up = v + n, *vp = up + n;
    for( i = 0; i < n; ++i )
    {
        realType ww = 1, zi = z[ i ];
        for( j = 0; j < i; ++j )
            ww *= zi - z[ j ];
        for( ++j; j < n; ++j )
            ww *= zi - z[ j ];
        w[ i ] = 1 / ww;
    }
    u[ 0 ] = v[ n - 1 ] = 1;
    up[ 0 ] = vp[ n - 1 ] = 0;
    for( i = 0; i < m; ++i )
    {
        realType xi = x[ i ];
        for( j = 0; j < n; ++j )
            d[ j ] = xi - z[ j ];
        for( j = 0; j < n - 1; ++j )
            u[ j + 1 ] = d[ j ] * u[ j ], up[ j + 1 ] = d[ j ] * up[ j ] + u[ j ];
        for( j = n - 1; j; --j )
            v[ j - 1 ] = d[ j ] * v[ j ], vp[ j - 1 ] = d[ j ] * vp[ j ] + v[ j ];
        for( j = 0; j < n; ++j )
            *J++ = w[ j ] * u[ j ] * v[ j ], *D++ = w[ j ] * ( up[ j ] * v[ j ] + u[ j ] * vp[ j ] );
    }
}
static realType legendre ( int  n,
realType  x 
) [static]

Definition at line 104 of file poly.c.

References n.

Referenced by gauss_nodes(), gauss_weights(), and lobatto_weights().

{
    realType p[ 2 ];
    int      i;
    p[ 0 ] = 1;
    p[ 1 ] = x;
    for( i = 1; i < n; i += 2 )
    {
        p[ 0 ] = ( ( 2 * i + 1 ) * x * p[ 1 ] - i * p[ 0 ] ) / ( i + 1 );
        p[ 1 ] = ( ( 2 * i + 3 ) * x * p[ 0 ] - ( i + 1 ) * p[ 1 ] ) / ( i + 2 );
    }
    return p[ n & 1 ];
}
static realType legendre_d1 ( int  n,
realType  x 
) [static]

Definition at line 119 of file poly.c.

References n.

Referenced by gauss_nodes(), and lobatto_nodes_aux().

{
    realType p[ 2 ];
    int      i;
    p[ 0 ] = 3 * x;
    p[ 1 ] = 1;
    for( i = 2; i < n; i += 2 )
    {
        p[ 1 ] = ( ( 2 * i + 1 ) * x * p[ 0 ] - ( i + 1 ) * p[ 1 ] ) / i;
        p[ 0 ] = ( ( 2 * i + 3 ) * x * p[ 1 ] - ( i + 2 ) * p[ 0 ] ) / ( i + 1 );
    }
    return p[ n & 1 ];
}
static realType legendre_d2 ( int  n,
realType  x 
) [static]

Definition at line 134 of file poly.c.

References n.

Referenced by lobatto_nodes_aux().

{
    realType p[ 2 ];
    int      i;
    p[ 0 ] = 3;
    p[ 1 ] = 15 * x;
    for( i = 3; i < n; i += 2 )
    {
        p[ 0 ] = ( ( 2 * i + 1 ) * x * p[ 1 ] - ( i + 2 ) * p[ 0 ] ) / ( i - 1 );
        p[ 1 ] = ( ( 2 * i + 3 ) * x * p[ 0 ] - ( i + 3 ) * p[ 1 ] ) / i;
    }
    return p[ n & 1 ];
}
void legendre_matrix ( const realType x,
int  m,
realType P,
int  n 
)

Definition at line 30 of file poly.c.

References b, n, and P.

Referenced by gauss_to_legendre_t().

{
    int       i, j;
    realType *Pjm1 = P, *Pj = Pjm1 + m, *Pjp1 = Pj + m;
    for( i = 0; i < m; ++i )
        Pjm1[ i ] = 1;
    for( i = 0; i < m; ++i )
        Pj[ i ] = x[ i ];
    for( j = 1; j < n; ++j )
    {
        realType c = 1 / ( realType )( j + 1 ), a = c * ( 2 * j + 1 ), b = c * j;
        for( i = 0; i < m; ++i )
            Pjp1[ i ] = a * x[ i ] * Pj[ i ] - b * Pjm1[ i ];
        Pjp1 += m, Pj += m, Pjm1 += m;
    }
}
void legendre_matrix_t ( const realType x,
int  m,
realType P,
int  n 
)

Definition at line 87 of file poly.c.

References legendre_row_even(), and legendre_row_odd().

Referenced by gauss_to_legendre().

{
    int i;
    if( n & 1 )
        for( i = 0; i < m; ++i, P += n + 1 )
            legendre_row_odd( x[ i ], P, n );
    else
        for( i = 0; i < m; ++i, P += n + 1 )
            legendre_row_even( x[ i ], P, n );
}
void legendre_row ( realType  x,
realType P,
int  n 
)

Definition at line 75 of file poly.c.

References legendre_row_even(), and legendre_row_odd().

{
    if( n & 1 )
        legendre_row_odd( x, P, n );
    else
        legendre_row_even( x, P, n );
}
static void legendre_row_even ( realType  x,
realType P,
int  n 
) [static]

Definition at line 48 of file poly.c.

References n.

Referenced by legendre_matrix_t(), legendre_row(), and lobatto_to_legendre().

{
    int i;
    P[ 0 ] = 1, P[ 1 ] = x;
    for( i = 1; i <= n - 2; i += 2 )
    {
        P[ i + 1 ] = ( ( 2 * i + 1 ) * x * P[ i ] - i * P[ i - 1 ] ) / ( i + 1 );
        P[ i + 2 ] = ( ( 2 * i + 3 ) * x * P[ i - 1 ] - ( i + 1 ) * P[ i ] ) / ( i + 2 );
    }
    P[ n ] = ( ( 2 * n - 1 ) * x * P[ n - 1 ] - ( n - 1 ) * P[ n - 2 ] ) / n;
}
static void legendre_row_odd ( realType  x,
realType P,
int  n 
) [static]

Definition at line 61 of file poly.c.

Referenced by legendre_matrix_t(), legendre_row(), and lobatto_to_legendre().

{
    int i;
    P[ 0 ] = 1, P[ 1 ] = x;
    for( i = 1; i <= n - 2; i += 2 )
    {
        P[ i + 1 ] = ( ( 2 * i + 1 ) * x * P[ i ] - i * P[ i - 1 ] ) / ( i + 1 );
        P[ i + 2 ] = ( ( 2 * i + 3 ) * x * P[ i - 1 ] - ( i + 1 ) * P[ i ] ) / ( i + 2 );
    }
}
static void lobatto_nodes_aux ( realType z,
int  n 
) [static]

Definition at line 174 of file poly.c.

References legendre_d1(), legendre_d2(), mbabs, mbcos, MOAB_POLY_EPS, MOAB_POLY_PI, n, and np.

Referenced by lobatto_nodes().

{
    int i, j, np = n + 1;
    for( i = 0; i <= n / 2 - 1; ++i )
    {
        realType ox, x = mbcos( ( n - i ) * MOAB_POLY_PI / np );
        do
        {
            ox = x;
            x -= legendre_d1( np, x ) / legendre_d2( np, x );
        } while( mbabs( x - ox ) > -x * MOAB_POLY_EPS );
        z[ i ] = x - legendre_d1( np, x ) / legendre_d2( np, x );
    }
    if( n & 1 ) z[ n / 2 ] = 0;
    for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
        z[ j ] = -z[ i ];
}
void lobatto_to_legendre ( const realType z,
const realType w,
int  n,
realType J 
)

Definition at line 275 of file poly.c.

References legendre_row_even(), legendre_row_odd(), n, and moab::sum().

{
    int       i, j, m = ( n + 1 ) / 2;
    realType *p = J, *q;
    realType  ww, sum;
    if( n & 1 )
        for( j = 0; j < m; ++j, p += n )
            legendre_row_odd( z[ j ], p, n - 2 );
    else
        for( j = 0; j < m; ++j, p += n )
            legendre_row_even( z[ j ], p, n - 2 );
    p = J;
    for( j = 0; j < m; ++j )
    {
        ww = w[ j ], sum = 0;
        for( i = 0; i < n - 1; ++i )
            *p *= ( 2 * i + 1 ) * ww / 2, sum += *p++;
        *p++ = -sum;
    }
    q = J + ( n / 2 - 1 ) * n;
    if( n & 1 )
        for( ; j < n; ++j, p += n, q -= n )
        {
            for( i = 0; i < n - 1; i += 2 )
                p[ i ] = q[ i ], p[ i + 1 ] = -q[ i + 1 ];
            p[ i ] = q[ i ];
        }
    else
        for( ; j < n; ++j, p += n, q -= n )
        {
            for( i = 0; i < n - 1; i += 2 )
                p[ i ] = q[ i ], p[ i + 1 ] = -q[ i + 1 ];
        }
}
void lobatto_weights ( const realType z,
realType w,
int  n 
)

Definition at line 211 of file poly.c.

References legendre(), and n.

Referenced by hash_getbb_2(), and hash_getbb_3().

{
    int i, j;
    for( i = 0; i <= ( n - 1 ) / 2; ++i )
    {
        realType d = legendre( n - 1, z[ i ] );
        w[ i ] = 2 / ( ( n - 1 ) * n * d * d );
    }
    for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
        w[ j ] = w[ i ];
}
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines