MOAB: Mesh Oriented datABase
(version 5.2.1)
|
#include "float.h"
#include "math.h"
#include <stdlib.h>
Go to the source code of this file.
Classes | |
struct | lagrange_data |
struct | obbox_2 |
struct | obbox_3 |
struct | hash_data_2 |
struct | hash_data_3 |
struct | findpt_listel |
struct | opt_edge_data_2 |
struct | opt_point_data_2 |
struct | opt_data_2 |
struct | findpt_data_2 |
struct | opt_face_data_3 |
struct | opt_edge_data_3 |
struct | opt_point_data_3 |
struct | opt_data_3 |
struct | findpt_data_3 |
Defines | |
#define | MOAB_POLY_EPS (128*DBL_EPSILON) |
#define | MOAB_POLY_PI 3.1415926535897932384626433832795028841971693993751058209749445923 |
#define | mbsqrt sqrt |
#define | mbabs fabs |
#define | mbcos cos |
#define | mbsin sin |
#define | mbfloor floor |
#define | mbceil ceil |
#define | INTEGER int |
#define | GLOBAL_INT long |
#define | MAYBE_UNUSED |
#define | tmalloc(type, count) ((type*) smalloc((count)*sizeof(type),__FILE__) ) |
#define | tcalloc(type, count) ((type*) scalloc((count),sizeof(type),__FILE__) ) |
#define | trealloc(type, ptr, count) ((type*) srealloc((ptr),(count)*sizeof(type),__FILE__) ) |
#define | MAYBE_UNUSED |
#define | DECLMINMAX(type, prefix) |
Typedefs | |
typedef double | realType |
typedef signed INTEGER | sint |
typedef unsigned INTEGER | uint |
typedef signed GLOBAL_INT | slong |
typedef unsigned GLOBAL_INT | ulong |
typedef int(* | findpt_func )(void *, const realType *, int, uint *, realType *, realType *) |
Functions | |
void | legendre_matrix (const realType *x, int m, realType *P, int n) |
void | legendre_matrix_t (const realType *x, int m, realType *P, int n) |
void | legendre_row (realType x, realType *P, int n) |
void | gauss_nodes (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_setup (lagrange_data *p, const realType *z, unsigned n) |
void | lagrange_free (lagrange_data *p) |
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 | tensor_c1 (const realType *R, unsigned mr, unsigned nr, const realType *u, realType *v) |
void | tensor_r1 (const realType *R, unsigned mr, unsigned nr, const realType *u, realType *v) |
void | tensor_c2 (const realType *R, unsigned mr, unsigned nr, const realType *S, unsigned ms, unsigned ns, const realType *u, realType *v, realType *work) |
void | tensor_r2 (const realType *R, unsigned mr, unsigned nr, const realType *S, unsigned ms, unsigned ns, const realType *u, realType *v, realType *work) |
void | tensor_c3 (const realType *R, unsigned mr, unsigned nr, const realType *S, unsigned ms, unsigned ns, const realType *T, unsigned mt, unsigned nt, const realType *u, realType *v, realType *work1, realType *work2) |
void | tensor_r3 (const realType *R, unsigned mr, unsigned nr, const realType *S, unsigned ms, unsigned ns, const realType *T, unsigned mt, unsigned nt, const realType *u, realType *v, realType *work1, realType *work2) |
realType | tensor_i1 (const realType *Jr, unsigned nr, const realType *u) |
realType | tensor_i2 (const realType *Jr, unsigned nr, const realType *Js, unsigned ns, const realType *u, realType *work) |
realType | tensor_i3 (const realType *Jr, unsigned nr, const realType *Js, unsigned ns, const realType *Jt, unsigned nt, const realType *u, realType *work) |
realType | tensor_ig1 (const realType *Jr, const realType *Dr, unsigned nr, const realType *u, realType *g) |
realType | tensor_ig2 (const realType *Jr, const realType *Dr, unsigned nr, const realType *Js, const realType *Ds, unsigned ns, const realType *u, realType *g, realType *work) |
realType | tensor_ig3 (const realType *Jr, const realType *Dr, unsigned nr, const realType *Js, const realType *Ds, unsigned ns, const realType *Jt, const realType *Dt, unsigned nt, const realType *u, realType *g, realType *work) |
findpt_data_2 * | findpt_setup_2 (const realType *const xw[2], const unsigned n[2], uint nel, uint max_hash_size, realType bbox_tol) |
findpt_data_3 * | findpt_setup_3 (const realType *const xw[3], const unsigned n[3], uint nel, uint max_hash_size, realType bbox_tol) |
void | findpt_free_2 (findpt_data_2 *p) |
void | findpt_free_3 (findpt_data_3 *p) |
const realType * | findpt_allbnd_2 (const findpt_data_2 *p) |
const realType * | findpt_allbnd_3 (const findpt_data_3 *p) |
int | findpt_2 (findpt_data_2 *p, const realType x[2], int guess, uint *el, realType r[2], realType *dist) |
int | findpt_3 (findpt_data_3 *p, const realType x[3], int guess, uint *el, realType r[3], realType *dist) |
void | findpt_weights_2 (findpt_data_2 *p, const realType r[2]) |
void | findpt_weights_3 (findpt_data_3 *p, const realType r[3]) |
double | findpt_eval_2 (findpt_data_2 *p, const realType *u) |
double | findpt_eval_3 (findpt_data_3 *p, const realType *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 realType *const elx[3], const realType xstar[3], realType r[3], unsigned *constr) |
void | opt_vol_set_intp_3 (opt_data_3 *p, const realType 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 realType *const elx[2], const realType xstar[2], realType r[2], unsigned *constr) |
void | fail (const char *fmt,...) |
static void * | smalloc (size_t size, const char *file) MAYBE_UNUSED |
static void * | scalloc (size_t nmemb, size_t size, const char *file) MAYBE_UNUSED |
static void * | srealloc (void *ptr, size_t size, const char *file) MAYBE_UNUSED |
static realType | r1norm_1 (realType a) MAYBE_UNUSED |
static realType | r1norm_2 (realType a, realType b) MAYBE_UNUSED |
static realType | r1norm_3 (realType a, realType b, realType c) MAYBE_UNUSED |
static realType | r2norm_1 (realType a) MAYBE_UNUSED |
static realType | r2norm_2 (realType a, realType b) MAYBE_UNUSED |
static realType | r2norm_3 (realType a, realType b, realType c) MAYBE_UNUSED |
Variables | |
const unsigned | opt_no_constraints_2 |
const unsigned | opt_no_constraints_3 |
#define DECLMINMAX | ( | type, | |
prefix | |||
) |
static type prefix##min_2(type a, type b) MAYBE_UNUSED; \ static type prefix##min_2(type a, type b) { \ return b<a?b:a; \ } \ static type prefix##max_2(type a, type b) MAYBE_UNUSED; \ static type prefix##max_2(type a, type b) { \ return a>b?a:b; \ } \ static void prefix##minmax_2(type *min, type *max, type a, \ type b) MAYBE_UNUSED; \ static void prefix##minmax_2(type *min, type *max, type a, type b) { \ if(b<a) *min=b, *max=a; \ else *min=a, *max=b; \ } \ static type prefix##min_3(type a, type b, type c) MAYBE_UNUSED; \ static type prefix##min_3(type a, type b, type c) { \ return b<a?(c<b?c:b):(c<a?c:a); \ } \ static type prefix##max_3(type a, type b, type c) MAYBE_UNUSED; \ static type prefix##max_3(type a, type b, type c) { \ return a>b?(a>c?a:c):(b>c?b:c); \ } \ static void prefix##minmax_3(type *min, type *max, type a, type b, \ type c) MAYBE_UNUSED; \ static void prefix##minmax_3(type *min, type *max, type a, type b, \ type c) { \ if(b<a) *min=prefix##min_2(b,c), *max=prefix##max_2(a,c); \ else *min=prefix##min_2(a,c), *max=prefix##max_2(b,c); \ }
Definition at line 612 of file FindPtFuncs.h.
#define GLOBAL_INT long |
Definition at line 47 of file FindPtFuncs.h.
#define INTEGER int |
Definition at line 38 of file FindPtFuncs.h.
#define MAYBE_UNUSED |
Definition at line 609 of file FindPtFuncs.h.
#define MAYBE_UNUSED |
Definition at line 609 of file FindPtFuncs.h.
#define mbabs fabs |
Definition at line 16 of file FindPtFuncs.h.
Referenced by gauss_nodes(), lobatto_nodes_aux(), obbox_test_2(), obbox_test_3(), r1norm_1(), r1norm_2(), and r1norm_3().
#define mbceil ceil |
Definition at line 20 of file FindPtFuncs.h.
Referenced by iceil().
#define mbcos cos |
Definition at line 17 of file FindPtFuncs.h.
Referenced by gauss_nodes(), lob_bnd_base_setup(), and lobatto_nodes_aux().
#define mbfloor floor |
Definition at line 19 of file FindPtFuncs.h.
Referenced by ifloor().
#define mbsin sin |
Definition at line 18 of file FindPtFuncs.h.
#define mbsqrt sqrt |
Definition at line 15 of file FindPtFuncs.h.
Referenced by r2norm_1(), r2norm_2(), and r2norm_3().
#define MOAB_POLY_EPS (128*DBL_EPSILON) |
Definition at line 12 of file FindPtFuncs.h.
Referenced by gauss_nodes(), lobatto_nodes_aux(), moab::ElemUtil::nat_coords_trilinear_hex2(), opt_findpt_2(), and opt_findpt_3().
#define MOAB_POLY_PI 3.1415926535897932384626433832795028841971693993751058209749445923 |
Definition at line 13 of file FindPtFuncs.h.
Referenced by gauss_nodes(), lob_bnd_base_setup(), and lobatto_nodes_aux().
Definition at line 566 of file FindPtFuncs.h.
Referenced by hash_build_2(), and hash_build_3().
Definition at line 564 of file FindPtFuncs.h.
Referenced by findpt_setup_2(), findpt_setup_3(), hash_build_2(), hash_build_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(), moab::element_utility::Spectral_hex_map< moab::Matrix3 >::initialize_spectral_hex(), lagrange_setup(), lob_bnd_base_alloc(), lob_bnd_base_setup(), lob_bnd_ext_alloc(), obbox_data_alloc_2(), obbox_data_alloc_3(), obbox_setup_2(), obbox_setup_3(), opt_alloc_2(), and opt_alloc_3().
Definition at line 568 of file FindPtFuncs.h.
typedef int(* findpt_func)(void *, const realType *, int, uint *, realType *, realType *) |
Definition at line 470 of file FindPtFuncs.h.
typedef double realType |
Definition at line 56 of file FindPtFuncs.h.
Definition at line 61 of file FindPtFuncs.h.
typedef signed GLOBAL_INT slong |
Definition at line 71 of file FindPtFuncs.h.
Definition at line 65 of file FindPtFuncs.h.
typedef unsigned GLOBAL_INT ulong |
Definition at line 79 of file FindPtFuncs.h.
void fail | ( | const char * | fmt, |
... | |||
) |
Definition at line 6 of file errmem.c.
Referenced by moab::TupleList::buffer::buffer_init_(), moab::TupleList::buffer::buffer_reserve_(), NumericalTestOF::evaluate(), moab::gs_data::gs_data_op_many(), moab::gs_data::gs_data_op_vec(), moab::TupleList::initialize(), moab::ReadNCDF::load_file(), main(), ZoltanPartitioner::mbGlobalSuccess(), moab::WriteSLAC::open_file(), moab::WriteNCDF::open_file(), moab::ReadNCDF::read_exodus_header(), moab::ReadNCDF::read_nodes(), moab::ReadNCDF::read_polyhedra_faces(), moab::ReadNCDF::read_tag_values(), moab::TupleList::resize(), runner_run_tests(), scalloc(), smalloc(), srealloc(), test_parallel_partition(), moab::WriteNCDF::write_BCs(), moab::WriteNCDF::write_elementblocks(), moab::WriteNCDF::write_exodus_integer_variable(), moab::WriteSLAC::write_file(), moab::WriteNCDF::write_file(), moab::WriteSLAC::write_matsets(), moab::WriteSLAC::write_nodes(), moab::WriteNCDF::write_nodes(), moab::WriteNCDF::write_poly_faces(), and moab::WriteNCDF::write_qa_string().
{ va_list ap; va_start( ap, fmt ); vfprintf( stderr, fmt, ap ); va_end( ap ); exit( 1 ); }
int findpt_2 | ( | findpt_data_2 * | p, |
const realType | x[2], | ||
int | guess, | ||
uint * | el, | ||
realType | r[2], | ||
realType * | dist | ||
) |
Definition at line 2186 of file findpt.c.
{ 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 ); }
int findpt_3 | ( | findpt_data_3 * | p, |
const realType | x[3], | ||
int | guess, | ||
uint * | el, | ||
realType | r[3], | ||
realType * | dist | ||
) |
Definition at line 2194 of file findpt.c.
{ 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 realType* findpt_allbnd_2 | ( | const findpt_data_2 * | p | ) |
const realType* findpt_allbnd_3 | ( | const findpt_data_3 * | p | ) |
double findpt_eval_2 | ( | findpt_data_2 * | p, |
const realType * | u | ||
) | [inline] |
Definition at line 2218 of file findpt.c.
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 realType * | u | ||
) | [inline] |
Definition at line 2223 of file findpt.c.
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 | ) |
void findpt_free_3 | ( | findpt_data_3 * | p | ) |
findpt_data_2* findpt_setup_2 | ( | const realType *const | xw[2], |
const unsigned | n[2], | ||
uint | nel, | ||
uint | max_hash_size, | ||
realType | bbox_tol | ||
) |
Definition at line 1952 of file findpt.c.
{ 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 realType *const | xw[3], |
const unsigned | n[3], | ||
uint | nel, | ||
uint | max_hash_size, | ||
realType | bbox_tol | ||
) |
Definition at line 1983 of file findpt.c.
{ 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 realType | r[2] | ||
) | [inline] |
Definition at line 2205 of file findpt.c.
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 realType | r[3] | ||
) | [inline] |
Definition at line 2211 of file findpt.c.
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 | ( | realType * | z, |
int | n | ||
) |
Definition at line 155 of file poly.c.
{ 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 | ||
) |
void gauss_to_legendre_t | ( | const realType * | z, |
const realType * | w, | ||
int | n, | ||
realType * | J | ||
) |
void gauss_weights | ( | const realType * | z, |
realType * | w, | ||
int | n | ||
) |
void lagrange_0 | ( | lagrange_data * | p, |
realType | x | ||
) |
void lagrange_1 | ( | lagrange_data * | p, |
realType | x | ||
) |
Definition at line 422 of file poly.c.
{ 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.
{ 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.
{ 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_free | ( | lagrange_data * | p | ) |
void lagrange_setup | ( | lagrange_data * | p, |
const realType * | z, | ||
unsigned | n | ||
) |
Definition at line 465 of file poly.c.
{ 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.
{ 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 realType * | x, |
int | m, | ||
realType * | P, | ||
int | n | ||
) |
Definition at line 30 of file poly.c.
{ 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.
{ 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.
{ if( n & 1 ) legendre_row_odd( x, P, n ); else legendre_row_even( x, P, n ); }
void lobatto_nodes | ( | realType * | z, |
int | n | ||
) |
void lobatto_to_legendre | ( | const realType * | z, |
const realType * | w, | ||
int | n, | ||
realType * | J | ||
) |
Definition at line 275 of file poly.c.
{ 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 | ||
) |
void opt_alloc_2 | ( | opt_data_2 * | p, |
lagrange_data * | ld | ||
) |
Definition at line 1605 of file findpt.c.
{ 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 1196 of file findpt.c.
{ 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 realType *const | elx[2], | ||
const realType | xstar[2], | ||
realType | r[2], | ||
unsigned * | constr | ||
) |
Definition at line 1762 of file findpt.c.
{ 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 realType *const | elx[3], | ||
const realType | xstar[3], | ||
realType | r[3], | ||
unsigned * | constr | ||
) |
Definition at line 1457 of file findpt.c.
{ 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 | ) |
void opt_free_3 | ( | opt_data_3 * | p | ) |
void opt_vol_set_intp_3 | ( | opt_data_3 * | p, |
const realType | r[3] | ||
) |
Definition at line 1246 of file findpt.c.
{ opt_vol_set_3( p, r ); opt_vol_intp_3( p ); }
Definition at line 653 of file FindPtFuncs.h.
References mbabs.
Referenced by findpt_hash_2(), and opt_findpt_2().
Definition at line 657 of file FindPtFuncs.h.
References mbabs.
Referenced by findpt_hash_3(), and opt_findpt_3().
Definition at line 665 of file FindPtFuncs.h.
References mbsqrt.
Referenced by findpt_pass_2(), and opt_findpt_2().
Definition at line 669 of file FindPtFuncs.h.
References mbsqrt.
Referenced by findpt_pass_3(), and opt_findpt_3().
static void * scalloc | ( | size_t | nmemb, |
size_t | size, | ||
const char * | file | ||
) | [static] |
static void * smalloc | ( | size_t | size, |
const char * | file | ||
) | [static] |
static void * srealloc | ( | void * | ptr, |
size_t | size, | ||
const char * | file | ||
) | [static] |
realType tensor_ig3 | ( | const realType * | Jr, |
const realType * | Dr, | ||
unsigned | nr, | ||
const realType * | Js, | ||
const realType * | Ds, | ||
unsigned | ns, | ||
const realType * | Jt, | ||
const realType * | Dt, | ||
unsigned | nt, | ||
const realType * | u, | ||
realType * | g, | ||
realType * | work | ||
) |
Definition at line 277 of file tensor.c.
{ 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 |
Definition at line 10 of file findpt.c.
Referenced by findpt_guess_2(), findpt_pass_2(), and opt_findpt_2().
const unsigned opt_no_constraints_3 |
Definition at line 11 of file findpt.c.
Referenced by findpt_guess_3(), findpt_pass_3(), moab::ElemUtil::hex_findpt(), moab::Element::SpectralHex::ievaluate(), moab::Element::SpectralQuad::ievaluate(), opt_findpt_3(), and moab::SpectralQuad::reverseEvalFcn().