MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#include <cfloat>
Go to the source code of this file.
Classes | |
struct | lagrange_data |
struct | findpt_data_2 |
struct | findpt_data_3 |
struct | opt_face_data_3 |
struct | opt_edge_data_3 |
struct | opt_point_data_3 |
struct | opt_data_3 |
struct | opt_edge_data_2 |
struct | opt_point_data_2 |
struct | opt_data_2 |
Defines | |
#define | INTEGER int |
#define | GLOBAL_INT long |
#define | floorr floor |
#define | ceilr ceil |
#define | sqrtr sqrt |
#define | fabsr fabs |
#define | cosr cos |
#define | sinr sin |
#define | EPS ( 128 * DBL_EPSILON ) |
#define | PI 3.1415926535897932384626433832795028841971693993751058209749445923 |
#define | uint uint_ |
#define | ulong ulong_ |
#define | sint sint_ |
#define | slong slong_ |
Typedefs | |
typedef double | real |
typedef signed INTEGER | sint |
typedef unsigned INTEGER | uint |
typedef signed GLOBAL_INT | slong |
typedef unsigned GLOBAL_INT | ulong |
typedef int(* | findpt_func )(void *, const real *, int, uint *, real *, real *) |
Functions | |
void | legendre_matrix (const real *x, int m, real *P, int n) |
void | legendre_matrix_t (const real *x, int m, real *P, int n) |
void | legendre_row (real x, real *P, int n) |
void | gauss_nodes (real *z, int n) |
void | lobatto_nodes (real *z, int n) |
void | gauss_weights (const real *z, real *w, int n) |
void | lobatto_weights (const real *z, real *w, int n) |
void | gauss_to_legendre (const real *z, const real *w, int n, real *J) |
void | gauss_to_legendre_t (const real *z, const real *w, int n, real *J) |
void | lobatto_to_legendre (const real *z, const real *w, int n, real *J) |
void | lagrange_weights (const real *z, unsigned n, const real *x, unsigned m, real *J, real *work) |
void | lagrange_weights_deriv (const real *z, unsigned n, const real *x, unsigned m, real *J, real *D, real *work) |
void | lagrange_setup (lagrange_data *p, const real *z, unsigned n) |
void | lagrange_free (lagrange_data *p) |
void | lagrange_0 (lagrange_data *p, real x) |
void | lagrange_1 (lagrange_data *p, real x) |
void | lagrange_2 (lagrange_data *p, real x) |
void | lagrange_2u (lagrange_data *p) |
void | tensor_c1 (const real *R, unsigned mr, unsigned nr, const real *u, real *v) |
void | tensor_r1 (const real *R, unsigned mr, unsigned nr, const real *u, real *v) |
void | tensor_c2 (const real *R, unsigned mr, unsigned nr, const real *S, unsigned ms, unsigned ns, const real *u, real *v, real *work) |
void | tensor_r2 (const real *R, unsigned mr, unsigned nr, const real *S, unsigned ms, unsigned ns, const real *u, real *v, real *work) |
void | tensor_c3 (const real *R, unsigned mr, unsigned nr, const real *S, unsigned ms, unsigned ns, const real *T, unsigned mt, unsigned nt, const real *u, real *v, real *work1, real *work2) |
void | tensor_r3 (const real *R, unsigned mr, unsigned nr, const real *S, unsigned ms, unsigned ns, const real *T, unsigned mt, unsigned nt, const real *u, real *v, real *work1, real *work2) |
real | tensor_i1 (const real *Jr, unsigned nr, const real *u) |
real | tensor_i2 (const real *Jr, unsigned nr, const real *Js, unsigned ns, const real *u, real *work) |
real | tensor_i3 (const real *Jr, unsigned nr, const real *Js, unsigned ns, const real *Jt, unsigned nt, const real *u, real *work) |
real | tensor_ig1 (const real *Jr, const real *Dr, unsigned nr, const real *u, real *g) |
real | tensor_ig2 (const real *Jr, const real *Dr, unsigned nr, const real *Js, const real *Ds, unsigned ns, const real *u, real *g, real *work) |
real | tensor_ig3 (const real *Jr, const real *Dr, unsigned nr, const real *Js, const real *Ds, unsigned ns, const real *Jt, const real *Dt, unsigned nt, const real *u, real *g, real *work) |
findpt_data_2 * | findpt_setup_2 (const real *const xw[2], const unsigned n[2], uint nel, uint max_hash_size, real bbox_tol) |
findpt_data_3 * | findpt_setup_3 (const real *const xw[3], const unsigned n[3], uint nel, uint max_hash_size, real bbox_tol) |
void | findpt_free_2 (findpt_data_2 *p) |
void | findpt_free_3 (findpt_data_3 *p) |
const real * | findpt_allbnd_2 (const findpt_data_2 *p) |
const real * | findpt_allbnd_3 (const findpt_data_3 *p) |
int | findpt_2 (findpt_data_2 *p, const real x[2], int guess, uint *el, real r[2], real *dist) |
int | findpt_3 (findpt_data_3 *p, const real x[3], int guess, uint *el, real r[3], real *dist) |
void | findpt_weights_2 (findpt_data_2 *p, const real r[2]) |
void | findpt_weights_3 (findpt_data_3 *p, const real r[3]) |
double | findpt_eval_2 (findpt_data_2 *p, const real *u) |
double | findpt_eval_3 (findpt_data_3 *p, const real *u) |
void | opt_alloc_3 (opt_data_3 *p, lagrange_data *ld) |
void | opt_free_3 (opt_data_3 *p) |
double | opt_findpt_3 (opt_data_3 *p, const real *const elx[3], const real xstar[3], real r[3], unsigned *constr) |
void | opt_vol_set_intp_3 (opt_data_3 *p, const real r[3]) |
void | opt_alloc_2 (opt_data_2 *p, lagrange_data *ld) |
void | opt_free_2 (opt_data_2 *p) |
double | opt_findpt_2 (opt_data_2 *p, const real *const elx[2], const real xstar[2], real r[2], unsigned *constr) |
Variables | |
const unsigned | opt_no_constraints_2 = 3 + 1 |
const unsigned | opt_no_constraints_3 = 9 + 3 + 1 |
#define ceilr ceil |
Definition at line 54 of file SpectralFuncs.hpp.
#define cosr cos |
Definition at line 57 of file SpectralFuncs.hpp.
#define EPS ( 128 * DBL_EPSILON ) |
Definition at line 59 of file SpectralFuncs.hpp.
#define fabsr fabs |
Definition at line 56 of file SpectralFuncs.hpp.
#define floorr floor |
Definition at line 53 of file SpectralFuncs.hpp.
#define GLOBAL_INT long |
Definition at line 27 of file SpectralFuncs.hpp.
#define INTEGER int |
Definition at line 18 of file SpectralFuncs.hpp.
#define PI 3.1415926535897932384626433832795028841971693993751058209749445923 |
Definition at line 60 of file SpectralFuncs.hpp.
Referenced by main().
#define sinr sin |
Definition at line 58 of file SpectralFuncs.hpp.
#define sint sint_ |
Definition at line 66 of file SpectralFuncs.hpp.
Referenced by moab::Coupler::consolidate_tuples(), moab::TupleList::initialize(), moab::pack_tuples(), moab::TupleList::permute(), moab::TupleList::resize(), moab::TupleList::sort(), and moab::unpack_tuples().
#define slong slong_ |
Definition at line 67 of file SpectralFuncs.hpp.
Referenced by moab::TupleList::initialize(), moab::pack_tuples(), moab::TupleList::resize(), moab::TupleList::sort(), and moab::unpack_tuples().
#define sqrtr sqrt |
Definition at line 55 of file SpectralFuncs.hpp.
#define uint uint_ |
Definition at line 64 of file SpectralFuncs.hpp.
Referenced by hash_build_2(), hash_build_3(), hash_opt_size_2(), and hash_opt_size_3().
#define ulong ulong_ |
Definition at line 65 of file SpectralFuncs.hpp.
Referenced by moab::pack_tuples(), and moab::unpack_tuples().
typedef int( * findpt_func)(void *, const real *, int, uint *, real *, real *) |
Definition at line 417 of file SpectralFuncs.hpp.
typedef double real |
Definition at line 52 of file SpectralFuncs.hpp.
Definition at line 69 of file SpectralFuncs.hpp.
typedef signed GLOBAL_INT slong |
Definition at line 74 of file SpectralFuncs.hpp.
Definition at line 70 of file SpectralFuncs.hpp.
typedef unsigned GLOBAL_INT ulong |
Definition at line 75 of file SpectralFuncs.hpp.
Definition at line 2251 of file findpt.c.
References findpt_data_2::end, findpt_guess_2(), findpt_hash_2(), findpt_pass_2(), and findpt_data_2::sorted.
{ if( guess && findpt_guess_2( p, x, *el, r, dist ) ) return 0; findpt_hash_2( p, x ); if( p->sorted == p->end ) return -1; return findpt_pass_2( p, x, el, r, dist ); }
Definition at line 2259 of file findpt.c.
References findpt_data_3::end, findpt_guess_3(), findpt_hash_3(), findpt_pass_3(), and findpt_data_3::sorted.
{ if( guess && findpt_guess_3( p, x, *el, r, dist ) ) return 0; findpt_hash_3( p, x ); #if DIAGNOSTICS printf( "hashing leaves %d elements to consider\n", p->end - p->sorted ); #endif if( p->sorted == p->end ) return -1; return findpt_pass_3( p, x, el, r, dist ); }
const real* findpt_allbnd_2 | ( | const findpt_data_2 * | p | ) |
Definition at line 2107 of file findpt.c.
References hash_data_2::bnd, and findpt_data_2::hash.
const real* findpt_allbnd_3 | ( | const findpt_data_3 * | p | ) |
Definition at line 2112 of file findpt.c.
References hash_data_3::bnd, and findpt_data_3::hash.
double findpt_eval_2 | ( | findpt_data_2 * | p, |
const real * | u | ||
) | [inline] |
Definition at line 434 of file SpectralFuncs.hpp.
References lagrange_data::J, findpt_data_2::ld, lagrange_data::n, findpt_data_2::od_work, and tensor_i2().
double findpt_eval_3 | ( | findpt_data_3 * | p, |
const real * | u | ||
) | [inline] |
Definition at line 439 of file SpectralFuncs.hpp.
References lagrange_data::J, findpt_data_3::ld, lagrange_data::n, findpt_data_3::od_work, and tensor_i3().
void findpt_free_2 | ( | findpt_data_2 * | p | ) |
Definition at line 2079 of file findpt.c.
References findpt_data_2::hash, hash_free_2(), findpt_data_2::list, findpt_data_2::od, opt_free_2(), findpt_data_2::sorted, and findpt_data_2::z.
{ unsigned d; opt_free_2( p->od ); free( p->od ); hash_free_2( p->hash ); free( p->hash ); free( p->list ); free( p->sorted ); for( d = 0; d < 2; ++d ) free( p->z[d] ); free( p ); }
void findpt_free_3 | ( | findpt_data_3 * | p | ) |
Definition at line 2093 of file findpt.c.
References findpt_data_3::hash, hash_free_3(), findpt_data_3::list, findpt_data_3::od, opt_free_3(), findpt_data_3::sorted, and findpt_data_3::z.
{ unsigned d; opt_free_3( p->od ); free( p->od ); hash_free_3( p->hash ); free( p->hash ); free( p->list ); free( p->sorted ); for( d = 0; d < 3; ++d ) free( p->z[d] ); free( p ); }
findpt_data_2* findpt_setup_2 | ( | const real *const | xw[2], |
const unsigned | n[2], | ||
uint | nel, | ||
uint | max_hash_size, | ||
real | bbox_tol | ||
) |
Definition at line 2011 of file findpt.c.
References findpt_data_2::hash, hash_build_2(), lagrange_setup(), findpt_data_2::ld, findpt_data_2::list, lobatto_nodes(), hash_data_2::max, findpt_data_2::nptel, findpt_data_2::od, findpt_data_2::od_work, opt_alloc_2(), findpt_data_2::sorted, tmalloc, opt_data_2::work, findpt_data_2::xw, and findpt_data_2::z.
{ unsigned d; findpt_data_2* p = tmalloc( findpt_data_2, 1 ); p->hash = tmalloc( hash_data_2, 1 ); p->od = tmalloc( opt_data_2, 1 ); for( d = 0; d < 2; ++d ) p->xw[d] = xw[d]; p->nptel = n[0] * n[1]; hash_build_2( p->hash, xw, n, nel, max_hash_size, bbox_tol ); for( d = 0; d < 2; ++d ) { p->z[d] = tmalloc( realType, n[d] ); lobatto_nodes( p->z[d], n[d] ); lagrange_setup( &p->ld[d], p->z[d], n[d] ); } p->list = tmalloc( findpt_listel, p->hash->max ); p->sorted = tmalloc( findpt_listel*, p->hash->max ); opt_alloc_2( p->od, p->ld ); p->od_work = p->od->work; return p; }
findpt_data_3* findpt_setup_3 | ( | const real *const | xw[3], |
const unsigned | n[3], | ||
uint | nel, | ||
uint | max_hash_size, | ||
real | bbox_tol | ||
) |
Definition at line 2045 of file findpt.c.
References findpt_data_3::hash, hash_build_3(), lagrange_setup(), findpt_data_3::ld, findpt_data_3::list, lobatto_nodes(), hash_data_3::max, findpt_data_3::nptel, findpt_data_3::od, findpt_data_3::od_work, opt_alloc_3(), findpt_data_3::sorted, tmalloc, opt_data_3::work, findpt_data_3::xw, and findpt_data_3::z.
{ unsigned d; findpt_data_3* p = tmalloc( findpt_data_3, 1 ); p->hash = tmalloc( hash_data_3, 1 ); p->od = tmalloc( opt_data_3, 1 ); for( d = 0; d < 3; ++d ) p->xw[d] = xw[d]; p->nptel = n[0] * n[1] * n[2]; hash_build_3( p->hash, xw, n, nel, max_hash_size, bbox_tol ); for( d = 0; d < 3; ++d ) { p->z[d] = tmalloc( realType, n[d] ); lobatto_nodes( p->z[d], n[d] ); lagrange_setup( &p->ld[d], p->z[d], n[d] ); } p->list = tmalloc( findpt_listel, p->hash->max ); p->sorted = tmalloc( findpt_listel*, p->hash->max ); opt_alloc_3( p->od, p->ld ); p->od_work = p->od->work; return p; }
void findpt_weights_2 | ( | findpt_data_2 * | p, |
const real | r[2] | ||
) | [inline] |
Definition at line 421 of file SpectralFuncs.hpp.
References lagrange_0(), and findpt_data_2::ld.
{ lagrange_0( &p->ld[0], r[0] ); lagrange_0( &p->ld[1], r[1] ); }
void findpt_weights_3 | ( | findpt_data_3 * | p, |
const real | r[3] | ||
) | [inline] |
Definition at line 427 of file SpectralFuncs.hpp.
References lagrange_0(), and findpt_data_3::ld.
{ lagrange_0( &p->ld[0], r[0] ); lagrange_0( &p->ld[1], r[1] ); lagrange_0( &p->ld[2], r[2] ); }
void gauss_nodes | ( | real * | 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 real * | z, |
const real * | w, | ||
int | n, | ||
real * | 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 real * | z, |
const real * | w, | ||
int | n, | ||
real * | 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 real * | z, |
real * | w, | ||
int | n | ||
) |
Definition at line 199 of file poly.c.
References legendre().
void lagrange_0 | ( | lagrange_data * | p, |
real | 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, |
real | 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, |
real | 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 real * | 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 real * | z, |
unsigned | n, | ||
const real * | x, | ||
unsigned | m, | ||
real * | J, | ||
real * | 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 real * | z, |
unsigned | n, | ||
const real * | x, | ||
unsigned | m, | ||
real * | J, | ||
real * | D, | ||
real * | 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] ); } }
void legendre_matrix | ( | const real * | x, |
int | m, | ||
real * | 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 real * | x, |
int | m, | ||
real * | 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 | ( | real | x, |
real * | 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 ); }
void lobatto_nodes | ( | real * | 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 ); }
void lobatto_to_legendre | ( | const real * | z, |
const real * | w, | ||
int | n, | ||
real * | 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 real * | z, |
real * | w, | ||
int | n | ||
) |
Definition at line 211 of file poly.c.
References legendre().
Referenced by hash_getbb_2(), and hash_getbb_3().
void opt_alloc_2 | ( | opt_data_2 * | p, |
lagrange_data * | ld | ||
) |
Definition at line 1662 of file findpt.c.
References opt_data_2::ed, opt_edge_data_2::fd1, opt_data_2::ld, lagrange_data::n, nr, opt_data_2::size, tmalloc, umax_2, opt_data_2::work, and opt_edge_data_2::x.
Referenced by findpt_setup_2(), and moab::Element::SpectralQuad::Init().
{ const unsigned nr = ld[0].n, ns = ld[1].n, ne = umax_2( nr, ns ), nw = 2 * ns; p->size[0] = 1; p->size[1] = nr; p->size[2] = nr * ns; p->ld = ld; p->work = tmalloc( realType, 4 * ne + nw ); p->ed.x[0] = p->work + nw; p->ed.x[1] = p->ed.x[0] + ne; p->ed.fd1[0] = p->ed.x[1] + ne; p->ed.fd1[1] = p->ed.fd1[0] + ne; }
void opt_alloc_3 | ( | opt_data_3 * | p, |
lagrange_data * | ld | ||
) |
Definition at line 1251 of file findpt.c.
References opt_data_3::ed, opt_data_3::fd, opt_edge_data_3::fd1, opt_edge_data_3::fd2, opt_face_data_3::fdn, opt_data_3::ld, lagrange_data::n, nr, opt_data_3::size, tmalloc, opt_data_3::work, opt_face_data_3::x, and opt_edge_data_3::x.
Referenced by findpt_setup_3(), moab::ElemUtil::hex_findpt(), moab::Element::SpectralHex::Init(), moab::SpectralHex::initFcn(), and moab::element_utility::Spectral_hex_map< moab::Matrix3 >::initialize_spectral_hex().
{ const unsigned nr = ld[0].n, ns = ld[1].n, nt = ld[2].n, nf = umax_3( nr * ns, nr * nt, ns * nt ), ne = umax_3( nr, ns, nt ), nw = 2 * ns * nt + 3 * ns; p->size[0] = 1; p->size[1] = nr; p->size[2] = nr * ns; p->size[3] = p->size[2] * nt; p->ld = ld; p->work = tmalloc( realType, 6 * nf + 9 * ne + nw ); p->fd.x[0] = p->work + nw; p->fd.x[1] = p->fd.x[0] + nf; p->fd.x[2] = p->fd.x[1] + nf; p->fd.fdn[0] = p->fd.x[2] + nf; p->fd.fdn[1] = p->fd.fdn[0] + nf; p->fd.fdn[2] = p->fd.fdn[1] + nf; p->ed.x[0] = p->fd.fdn[2] + nf; p->ed.x[1] = p->ed.x[0] + ne; p->ed.x[2] = p->ed.x[1] + ne; p->ed.fd1[0] = p->ed.x[2] + ne; p->ed.fd1[1] = p->ed.fd1[0] + ne; p->ed.fd1[2] = p->ed.fd1[1] + ne; p->ed.fd2[0] = p->ed.fd1[2] + ne; p->ed.fd2[1] = p->ed.fd2[0] + ne; p->ed.fd2[2] = p->ed.fd2[1] + ne; }
double opt_findpt_2 | ( | opt_data_2 * | p, |
const real *const | elx[2], | ||
const real | xstar[2], | ||
real | r[2], | ||
unsigned * | constr | ||
) |
Definition at line 1818 of file findpt.c.
References opt_edge_data_2::constraints, opt_point_data_2::constraints, opt_edge_data_2::d1, opt_edge_data_2::de, DIAGNOSTICS, opt_data_2::ed, opt_data_2::elx, opt_data_2::jac, mat_app_2c(), MOAB_POLY_EPS, opt_area_set_intp_2(), opt_constr_num_2, opt_constr_pack_2(), opt_constr_unpack_2(), opt_edge_hess_2(), opt_edge_set_2(), opt_edge_set_intp_2(), opt_no_constraints_2, opt_point_set_2(), opt_point_set_intp_2(), opt_data_2::pd, r1norm_2(), r2norm_2(), tinyla_solve_2(), and opt_data_2::x.
Referenced by findpt_guess_2(), findpt_pass_2(), moab::Element::SpectralQuad::ievaluate(), and moab::SpectralQuad::reverseEvalFcn().
{ realType dr[2], resid[2], steep[2]; unsigned c = *constr, ac, d, cc[2], step = 0; p->elx[0] = elx[0], p->elx[1] = elx[1]; p->ed.constraints = opt_no_constraints_2; p->pd.constraints = opt_no_constraints_2; #if DIAGNOSTICS printf( "opt_findpt: xstar = %g, %g\n", xstar[0], xstar[1] ); #endif do { ++step; if( step == 50 ) /*fail("%s: opt_findpt_2 did not converge\n",__FILE__);*/ return 1.e+30; #if DIAGNOSTICS printf( " iteration %u\n", step ); printf( " %d constraint(s) active\n", (int)opt_constr_num_2[c] ); #endif /* update face/edge/point data if necessary, and evaluate x(r) as well as the jacobian */ switch( opt_constr_num_2[c] ) { case 0: opt_area_set_intp_2( p, r ); break; case 1: opt_edge_set_intp_2( p, r, c ); break; case 2: opt_point_set_intp_2( p, c ); break; } #if DIAGNOSTICS printf( " r = %g, %g\n", r[0], r[1] ); printf( " x = %g, %g\n", p->x[0], p->x[1] ); #endif /* compute residual */ for( d = 0; d < 2; ++d ) resid[d] = xstar[d] - p->x[d]; #if DIAGNOSTICS printf( " resid = %g, %g\n", resid[0], resid[1] ); printf( " 2-norm = %g\n", r2norm_2( resid[0], resid[1] ) ); #endif /* check constraints against steepest descent direction */ ac = c; if( opt_constr_num_2[c] ) { opt_constr_unpack_2( c, cc ); mat_app_2c( steep, p->jac, resid ); /* steepest descent = J^T r */ #if DIAGNOSTICS printf( " steepest descent = %g, %g\n", steep[0], steep[1] ); #endif for( d = 0; d < 2; ++d ) if( ( cc[d] == 0 && steep[d] > 0 ) || ( cc[d] == 2 && steep[d] < 0 ) ) cc[d] = 1; ac = opt_constr_pack_2( cc ); } /* update face/edge/point data if necessary */ if( ac != c ) { c = ac; #if DIAGNOSTICS printf( " relaxed to %d constraints\n", (int)opt_constr_num_2[c] ); #endif switch( opt_constr_num_2[c] ) { case 1: opt_edge_set_2( p, r, c ); break; case 2: opt_point_set_2( p, c ); break; } } /* compute Newton step */ switch( opt_constr_num_2[c] ) { case 0: tinyla_solve_2( dr, p->jac, resid ); break; case 1: { const unsigned de = p->ed.de, d1 = p->ed.d1; realType fac, H[2]; const realType* J = p->jac + de; opt_edge_hess_2( p, H ); fac = J[0] * J[0] + J[2] * J[2] - ( resid[0] * H[0] + resid[1] * H[1] ); dr[de] = steep[de] / fac; dr[d1] = 0; } break; case 2: dr[0] = dr[1] = 0; break; } #if DIAGNOSTICS printf( " dr = %g, %g\n", dr[0], dr[1] ); #endif /* project new iteration onto [-1,1]^2 */ opt_constr_unpack_2( c, cc ); for( d = 0; d < 2; ++d ) { if( cc[d] != 1 ) continue; r[d] += dr[d]; if( r[d] <= -1 ) dr[d] -= r[d] + 1, r[d] = -1, cc[d] = 0; else if( r[d] >= 1 ) dr[d] -= r[d] - 1, r[d] = 1, cc[d] = 2; } c = opt_constr_pack_2( cc ); } while( r1norm_2( dr[0], dr[1] ) > 2 * MOAB_POLY_EPS ); *constr = c; return r2norm_2( resid[0], resid[1] ); }
double opt_findpt_3 | ( | opt_data_3 * | p, |
const real *const | elx[3], | ||
const real | xstar[3], | ||
real | r[3], | ||
unsigned * | constr | ||
) |
Definition at line 1512 of file findpt.c.
References opt_face_data_3::constraints, opt_edge_data_3::constraints, opt_point_data_3::constraints, opt_face_data_3::d1, opt_edge_data_3::d1, opt_face_data_3::d2, opt_edge_data_3::d2, opt_edge_data_3::de, DIAGNOSTICS, opt_face_data_3::dn, opt_data_3::ed, opt_data_3::elx, opt_data_3::fd, opt_data_3::jac, mat_app_3c(), MOAB_POLY_EPS, opt_constr_num_3, opt_constr_pack_3(), opt_constr_unpack_3(), opt_edge_hess_3(), opt_edge_set_3(), opt_edge_set_intp_3(), opt_face_hess_3(), opt_face_set_3(), opt_face_set_intp_3(), opt_no_constraints_3, opt_point_set_3(), opt_point_set_intp_3(), opt_vol_set_intp_3(), opt_data_3::pd, r1norm_3(), r2norm_3(), tinyla_solve_3(), tinyla_solve_sym_2(), and opt_data_3::x.
Referenced by findpt_guess_3(), findpt_pass_3(), moab::ElemUtil::hex_findpt(), and moab::Element::SpectralHex::ievaluate().
{ realType dr[3], resid[3], steep[3]; unsigned c = *constr, ac, d, cc[3], step = 0; p->elx[0] = elx[0], p->elx[1] = elx[1], p->elx[2] = elx[2]; p->fd.constraints = opt_no_constraints_3; p->ed.constraints = opt_no_constraints_3; p->pd.constraints = opt_no_constraints_3; #if DIAGNOSTICS printf( "opt_findpt: xstar = %g, %g, %g\n", xstar[0], xstar[1], xstar[2] ); #endif do { ++step; if( step == 50 ) /*fail("%s: opt_findpt_3 did not converge\n",__FILE__);*/ return 1.e+30; #if DIAGNOSTICS printf( " iteration %u\n", step ); printf( " %d constraint(s) active\n", (int)opt_constr_num_3[c] ); #endif /* update face/edge/point data if necessary, and evaluate x(r) as well as the jacobian */ switch( opt_constr_num_3[c] ) { case 0: opt_vol_set_intp_3( p, r ); break; case 1: opt_face_set_intp_3( p, r, c ); break; case 2: opt_edge_set_intp_3( p, r, c ); break; case 3: opt_point_set_intp_3( p, c ); break; } #if DIAGNOSTICS printf( " r = %g, %g, %g\n", r[0], r[1], r[2] ); printf( " x = %g, %g, %g\n", p->x[0], p->x[1], p->x[2] ); #endif /* compute residual */ for( d = 0; d < 3; ++d ) resid[d] = xstar[d] - p->x[d]; #if DIAGNOSTICS printf( " resid = %g, %g, %g\n", resid[0], resid[1], resid[2] ); printf( " 2-norm = %g\n", r2norm_3( resid[0], resid[1], resid[2] ) ); #endif /* check constraints against steepest descent direction */ ac = c; if( opt_constr_num_3[c] ) { opt_constr_unpack_3( c, cc ); mat_app_3c( steep, p->jac, resid ); /* steepest descent = J^T r */ #if DIAGNOSTICS printf( " steepest descent = %g, %g, %g\n", steep[0], steep[1], steep[2] ); #endif for( d = 0; d < 3; ++d ) if( ( cc[d] == 0 && steep[d] > 0 ) || ( cc[d] == 2 && steep[d] < 0 ) ) cc[d] = 1; ac = opt_constr_pack_3( cc ); } /* update face/edge/point data if necessary */ if( ac != c ) { c = ac; #if DIAGNOSTICS printf( " relaxed to %d constraints\n", (int)opt_constr_num_3[c] ); #endif switch( opt_constr_num_3[c] ) { case 1: opt_face_set_3( p, r, c ); break; case 2: opt_edge_set_3( p, r, c ); break; case 3: opt_point_set_3( p, c ); break; } } /* compute Newton step */ switch( opt_constr_num_3[c] ) { case 0: tinyla_solve_3( dr, p->jac, resid ); break; case 1: { const unsigned dn = p->fd.dn, d1 = p->fd.d1, d2 = p->fd.d2; realType A[4], H[9]; const realType* J = p->jac; opt_face_hess_3( p, H ); A[0] = J[d1] * J[d1] + J[3 + d1] * J[3 + d1] + J[6 + d1] * J[6 + d1]; A[1] = J[d2] * J[d2] + J[3 + d2] * J[3 + d2] + J[6 + d2] * J[6 + d2]; A[2] = J[d1] * J[d2] + J[3 + d1] * J[3 + d2] + J[6 + d1] * J[6 + d2]; A[0] -= resid[0] * H[0] + resid[1] * H[3] + resid[2] * H[6]; A[1] -= resid[0] * H[1] + resid[1] * H[4] + resid[2] * H[7]; A[2] -= resid[0] * H[2] + resid[1] * H[5] + resid[2] * H[8]; tinyla_solve_sym_2( &dr[d1], &dr[d2], A, steep[d1], steep[d2] ); dr[dn] = 0; } break; case 2: { const unsigned de = p->ed.de, d1 = p->ed.d1, d2 = p->ed.d2; realType fac, H[3]; const realType* J = p->jac + de; opt_edge_hess_3( p, H ); fac = J[0] * J[0] + J[3] * J[3] + J[6] * J[6] - ( resid[0] * H[0] + resid[1] * H[1] + resid[2] * H[2] ); dr[de] = steep[de] / fac; dr[d1] = 0, dr[d2] = 0; } break; case 3: dr[0] = dr[1] = dr[2] = 0; break; } #if DIAGNOSTICS printf( " dr = %g, %g, %g\n", dr[0], dr[1], dr[2] ); #endif /* project new iteration onto [-1,1]^3 */ opt_constr_unpack_3( c, cc ); for( d = 0; d < 3; ++d ) { if( cc[d] != 1 ) continue; r[d] += dr[d]; if( r[d] <= -1 ) dr[d] -= r[d] + 1, r[d] = -1, cc[d] = 0; else if( r[d] >= 1 ) dr[d] -= r[d] - 1, r[d] = 1, cc[d] = 2; } c = opt_constr_pack_3( cc ); } while( r1norm_3( dr[0], dr[1], dr[2] ) > 3 * MOAB_POLY_EPS ); *constr = c; #if 0 printf("opt_findpt_3 converged in %u iterations\n", step); #endif return r2norm_3( resid[0], resid[1], resid[2] ); }
void opt_free_2 | ( | opt_data_2 * | p | ) |
Definition at line 1676 of file findpt.c.
References opt_data_2::work.
Referenced by findpt_free_2(), and moab::Element::SpectralQuad::freedata().
{ free( p->work ); }
void opt_free_3 | ( | opt_data_3 * | p | ) |
Definition at line 1278 of file findpt.c.
References opt_data_3::work.
Referenced by findpt_free_3(), moab::element_utility::Spectral_hex_map< moab::Matrix3 >::free_data(), moab::Element::SpectralHex::freedata(), and moab::ElemUtil::hex_findpt().
{ free( p->work ); }
void opt_vol_set_intp_3 | ( | opt_data_3 * | p, |
const real | r[3] | ||
) |
Definition at line 1301 of file findpt.c.
References opt_vol_intp_3(), and opt_vol_set_3().
Referenced by moab::element_utility::Spectral_hex_map< moab::Matrix3 >::integrate_scalar_field(), moab::Element::SpectralHex::integrate_scalar_field(), moab::element_utility::Spectral_hex_map< moab::Matrix3 >::jacobian(), moab::Element::SpectralHex::jacobian(), and opt_findpt_3().
{ opt_vol_set_3( p, r ); opt_vol_intp_3( p ); }
Definition at line 255 of file tensor.c.
References inner().
Referenced by opt_edge_hess_2(), opt_edge_hess_3(), opt_edge_intp_2(), and opt_edge_intp_3().
real tensor_i2 | ( | const real * | Jr, |
unsigned | nr, | ||
const real * | Js, | ||
unsigned | ns, | ||
const real * | u, | ||
real * | work | ||
) |
Definition at line 261 of file tensor.c.
References inner(), and mxv_r().
Referenced by moab::SpectralQuad::evalFcn(), moab::Element::SpectralQuad::evaluate(), moab::Element::SpectralQuad::evaluate_scalar_field(), findpt_eval_2(), opt_face_hess_3(), and opt_face_intp_3().
real tensor_i3 | ( | const real * | Jr, |
unsigned | nr, | ||
const real * | Js, | ||
unsigned | ns, | ||
const real * | Jt, | ||
unsigned | nt, | ||
const real * | u, | ||
real * | work | ||
) |
Definition at line 273 of file tensor.c.
References inner(), and mxv_r().
Referenced by moab::element_utility::Spectral_hex_map< moab::Matrix3 >::evaluate(), moab::Element::SpectralHex::evaluate(), moab::element_utility::Spectral_hex_map< moab::Matrix3 >::evaluate_scalar_field(), moab::Element::SpectralHex::evaluate_scalar_field(), findpt_eval_3(), and moab::ElemUtil::hex_eval().
Definition at line 304 of file tensor.c.
References inner().
Referenced by opt_edge_intp_2(), and opt_edge_intp_3().
real tensor_ig2 | ( | const real * | Jr, |
const real * | Dr, | ||
unsigned | nr, | ||
const real * | Js, | ||
const real * | Ds, | ||
unsigned | ns, | ||
const real * | u, | ||
real * | g, | ||
real * | work | ||
) |
Definition at line 311 of file tensor.c.
References inner(), and mxv_r().
Referenced by obbox_calc_2(), opt_area_intp_2(), opt_face_hess_3(), and opt_face_intp_3().
real tensor_ig3 | ( | const real * | Jr, |
const real * | Dr, | ||
unsigned | nr, | ||
const real * | Js, | ||
const real * | Ds, | ||
unsigned | ns, | ||
const real * | Jt, | ||
const real * | Dt, | ||
unsigned | nt, | ||
const real * | u, | ||
real * | g, | ||
real * | work | ||
) |
Definition at line 330 of file tensor.c.
References inner(), and mxv_r().
Referenced by obbox_calc_3(), and opt_vol_intp_3().
{ unsigned nsnt = ns * nt; realType *a = work, *ar = a + nsnt, *b = ar + nsnt, *br = b + ns, *bs = br + ns; mxv_r( a, nsnt, u, Jr, nr ); mxv_r( ar, nsnt, u, Dr, nr ); mxv_r( b, nt, a, Js, ns ); mxv_r( br, nt, ar, Js, ns ); mxv_r( bs, nt, a, Ds, ns ); g[0] = inner( Jt, br, nt ); g[1] = inner( Jt, bs, nt ); g[2] = inner( Dt, b, nt ); return inner( Jt, b, nt ); }
const unsigned opt_no_constraints_2 = 3 + 1 |
Definition at line 485 of file SpectralFuncs.hpp.
const unsigned opt_no_constraints_3 = 9 + 3 + 1 |
Definition at line 486 of file SpectralFuncs.hpp.