MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include "moab/FindPtFuncs.h"
Go to the source code of this file.
void gauss_nodes | ( | realType * | z, |
int | n | ||
) |
Definition at line 155 of file poly.c.
References legendre(), legendre_d1(), mbabs, mbcos, MOAB_POLY_EPS, and MOAB_POLY_PI.
{ 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().
{ 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().
{ 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().
void lagrange_0 | ( | lagrange_data * | p, |
realType | x | ||
) |
Definition at line 414 of file poly.c.
References lagrange_data::d, lagrange_data::J, lagrange_data::n, lagrange_data::u0, lagrange_data::v0, lagrange_data::w, and lagrange_data::z.
Referenced by moab::SpectralQuad::evalFcn(), moab::element_utility::Spectral_hex_map< moab::Matrix3 >::evaluate(), moab::Element::SpectralHex::evaluate(), moab::Element::SpectralQuad::evaluate(), moab::element_utility::Spectral_hex_map< moab::Matrix3 >::evaluate_scalar_field(), moab::Element::SpectralHex::evaluate_scalar_field(), moab::Element::SpectralQuad::evaluate_scalar_field(), findpt_weights_2(), findpt_weights_3(), and moab::ElemUtil::hex_eval().
void lagrange_1 | ( | lagrange_data * | p, |
realType | x | ||
) |
Definition at line 427 of file poly.c.
References lagrange_data::D, lagrange_data::d, lagrange_data::J, 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 440 of file poly.c.
References lagrange_data::D, lagrange_data::d, lagrange_data::D2, lagrange_data::J, 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 456 of file poly.c.
References lagrange_data::d, lagrange_data::D2, 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().
void lagrange_free | ( | lagrange_data * | p | ) |
Definition at line 497 of file poly.c.
References lagrange_data::w.
Referenced by moab::element_utility::Spectral_hex_map< moab::Matrix3 >::free_data(), moab::Element::SpectralHex::freedata(), moab::Element::SpectralQuad::freedata(), moab::ElemUtil::hex_eval(), and moab::ElemUtil::hex_findpt().
{ free( p->w ); }
void lagrange_setup | ( | lagrange_data * | p, |
const realType * | z, | ||
unsigned | n | ||
) |
Definition at line 467 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(), 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, 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.
{ 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.
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] ); } }
Definition at line 104 of file poly.c.
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.
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.
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.
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.
Referenced by legendre_matrix_t(), legendre_row(), and lobatto_to_legendre().
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().
void lobatto_nodes | ( | realType * | z, |
int | n | ||
) |
Definition at line 193 of file poly.c.
References lobatto_nodes_aux().
Referenced by findpt_setup_2(), findpt_setup_3(), hash_getbb_2(), hash_getbb_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().
{ z[0] = -1, z[n - 1] = 1; lobatto_nodes_aux( &z[1], n - 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, 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(), 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().
Referenced by hash_getbb_2(), and hash_getbb_3().