Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Macros | Typedefs | Functions | Variables
SpectralFuncs.hpp File Reference
#include <cfloat>
+ Include dependency graph for SpectralFuncs.hpp:
+ This graph shows which files directly or indirectly include this file:

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
 

Macros

#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
 

Macro Definition Documentation

◆ ceilr

#define ceilr   ceil

Definition at line 54 of file SpectralFuncs.hpp.

◆ cosr

#define cosr   cos

Definition at line 57 of file SpectralFuncs.hpp.

◆ EPS

#define EPS   ( 128 * DBL_EPSILON )

Definition at line 59 of file SpectralFuncs.hpp.

◆ fabsr

#define fabsr   fabs

Definition at line 56 of file SpectralFuncs.hpp.

◆ floorr

#define floorr   floor

Definition at line 53 of file SpectralFuncs.hpp.

◆ GLOBAL_INT

#define GLOBAL_INT   long

Definition at line 27 of file SpectralFuncs.hpp.

◆ INTEGER

#define INTEGER   int

Definition at line 18 of file SpectralFuncs.hpp.

◆ PI

#define PI   3.1415926535897932384626433832795028841971693993751058209749445923

Definition at line 60 of file SpectralFuncs.hpp.

◆ sinr

#define sinr   sin

Definition at line 58 of file SpectralFuncs.hpp.

◆ sint

#define sint   sint_

Definition at line 66 of file SpectralFuncs.hpp.

◆ slong

#define slong   slong_

Definition at line 67 of file SpectralFuncs.hpp.

◆ sqrtr

#define sqrtr   sqrt

Definition at line 55 of file SpectralFuncs.hpp.

◆ uint

#define uint   uint_

Definition at line 64 of file SpectralFuncs.hpp.

◆ ulong

#define ulong   ulong_

Definition at line 65 of file SpectralFuncs.hpp.

Typedef Documentation

◆ findpt_func

typedef int( * findpt_func) (void *, const real *, int, uint *, real *, real *)

Definition at line 417 of file SpectralFuncs.hpp.

◆ real

typedef double real

Definition at line 52 of file SpectralFuncs.hpp.

◆ sint

typedef signed INTEGER sint

Definition at line 69 of file SpectralFuncs.hpp.

◆ slong

typedef signed GLOBAL_INT slong

Definition at line 74 of file SpectralFuncs.hpp.

◆ uint

typedef unsigned INTEGER uint

Definition at line 70 of file SpectralFuncs.hpp.

◆ ulong

typedef unsigned GLOBAL_INT ulong

Definition at line 75 of file SpectralFuncs.hpp.

Function Documentation

◆ findpt_2()

int findpt_2 ( findpt_data_2 *  p,
const real  x[2],
int  guess,
uint *  el,
real  r[2],
real *  dist 
)

Definition at line 2251 of file findpt.c.

 2252 {
 2253     if( guess && findpt_guess_2( p, x, *el, r, dist ) ) return 0;
 2254     findpt_hash_2( p, x );
 2255     if( p->sorted == p->end ) return -1;
 2256     return findpt_pass_2( p, x, el, r, dist );
 2257 }

References findpt_data_2::end, findpt_guess_2(), findpt_hash_2(), findpt_pass_2(), and findpt_data_2::sorted.

◆ findpt_3()

int findpt_3 ( findpt_data_3 *  p,
const real  x[3],
int  guess,
uint *  el,
real  r[3],
real *  dist 
)

Definition at line 2259 of file findpt.c.

 2260 {
 2261     if( guess && findpt_guess_3( p, x, *el, r, dist ) ) return 0;
 2262     findpt_hash_3( p, x );
 2263 #if DIAGNOSTICS
 2264     printf( "hashing leaves %d elements to consider\n", p->end - p->sorted );
 2265 #endif
 2266     if( p->sorted == p->end ) return -1;
 2267     return findpt_pass_3( p, x, el, r, dist );
 2268 }

References findpt_data_3::end, findpt_guess_3(), findpt_hash_3(), findpt_pass_3(), and findpt_data_3::sorted.

◆ findpt_allbnd_2()

const real* findpt_allbnd_2 ( const findpt_data_2 *  p)

Definition at line 2107 of file findpt.c.

 2108 {
 2109     return p->hash->bnd;
 2110 }

References hash_data_2::bnd, and findpt_data_2::hash.

◆ findpt_allbnd_3()

const real* findpt_allbnd_3 ( const findpt_data_3 *  p)

Definition at line 2112 of file findpt.c.

 2113 {
 2114     return p->hash->bnd;
 2115 }

References hash_data_3::bnd, and findpt_data_3::hash.

◆ findpt_eval_2()

double findpt_eval_2 ( findpt_data_2 *  p,
const real *  u 
)
inline

Definition at line 434 of file SpectralFuncs.hpp.

  435 {
  436     return tensor_i2( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, u, p->od_work );
  437 }

References lagrange_data::J, findpt_data_2::ld, lagrange_data::n, findpt_data_2::od_work, and tensor_i2().

◆ findpt_eval_3()

double findpt_eval_3 ( findpt_data_3 *  p,
const real *  u 
)
inline

Definition at line 439 of file SpectralFuncs.hpp.

  440 {
  441     return tensor_i3( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, p->ld[2].J, p->ld[2].n, u, p->od_work );
  442 }

References lagrange_data::J, findpt_data_3::ld, lagrange_data::n, findpt_data_3::od_work, and tensor_i3().

◆ findpt_free_2()

void findpt_free_2 ( findpt_data_2 *  p)

Definition at line 2079 of file findpt.c.

 2080 {
 2081     unsigned d;
 2082     opt_free_2( p->od );
 2083     free( p->od );
 2084     hash_free_2( p->hash );
 2085     free( p->hash );
 2086     free( p->list );
 2087     free( p->sorted );
 2088     for( d = 0; d < 2; ++d )
 2089         free( p->z[d] );
 2090     free( p );
 2091 }

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.

◆ findpt_free_3()

void findpt_free_3 ( findpt_data_3 *  p)

Definition at line 2093 of file findpt.c.

 2094 {
 2095     unsigned d;
 2096     opt_free_3( p->od );
 2097     free( p->od );
 2098     hash_free_3( p->hash );
 2099     free( p->hash );
 2100     free( p->list );
 2101     free( p->sorted );
 2102     for( d = 0; d < 3; ++d )
 2103         free( p->z[d] );
 2104     free( p );
 2105 }

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.

◆ findpt_setup_2()

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.

 2016 {
 2017     unsigned d;
 2018     findpt_data_2* p = tmalloc( findpt_data_2, 1 );
 2019  
 2020     p->hash = tmalloc( hash_data_2, 1 );
 2021     p->od   = tmalloc( opt_data_2, 1 );
 2022  
 2023     for( d = 0; d < 2; ++d )
 2024         p->xw[d] = xw[d];
 2025     p->nptel = n[0] * n[1];
 2026  
 2027     hash_build_2( p->hash, xw, n, nel, max_hash_size, bbox_tol );
 2028  
 2029     for( d = 0; d < 2; ++d )
 2030     {
 2031         p->z[d] = tmalloc( realType, n[d] );
 2032         lobatto_nodes( p->z[d], n[d] );
 2033         lagrange_setup( &p->ld[d], p->z[d], n[d] );
 2034     }
 2035  
 2036     p->list   = tmalloc( findpt_listel, p->hash->max );
 2037     p->sorted = tmalloc( findpt_listel*, p->hash->max );
 2038  
 2039     opt_alloc_2( p->od, p->ld );
 2040     p->od_work = p->od->work;
 2041  
 2042     return p;
 2043 }

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.

◆ findpt_setup_3()

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.

 2050 {
 2051     unsigned d;
 2052     findpt_data_3* p = tmalloc( findpt_data_3, 1 );
 2053  
 2054     p->hash = tmalloc( hash_data_3, 1 );
 2055     p->od   = tmalloc( opt_data_3, 1 );
 2056  
 2057     for( d = 0; d < 3; ++d )
 2058         p->xw[d] = xw[d];
 2059     p->nptel = n[0] * n[1] * n[2];
 2060  
 2061     hash_build_3( p->hash, xw, n, nel, max_hash_size, bbox_tol );
 2062  
 2063     for( d = 0; d < 3; ++d )
 2064     {
 2065         p->z[d] = tmalloc( realType, n[d] );
 2066         lobatto_nodes( p->z[d], n[d] );
 2067         lagrange_setup( &p->ld[d], p->z[d], n[d] );
 2068     }
 2069  
 2070     p->list   = tmalloc( findpt_listel, p->hash->max );
 2071     p->sorted = tmalloc( findpt_listel*, p->hash->max );
 2072  
 2073     opt_alloc_3( p->od, p->ld );
 2074     p->od_work = p->od->work;
 2075  
 2076     return p;
 2077 }

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.

◆ findpt_weights_2()

void findpt_weights_2 ( findpt_data_2 *  p,
const real  r[2] 
)
inline

Definition at line 421 of file SpectralFuncs.hpp.

  422 {
  423     lagrange_0( &p->ld[0], r[0] );
  424     lagrange_0( &p->ld[1], r[1] );
  425 }

References lagrange_0(), and findpt_data_2::ld.

◆ findpt_weights_3()

void findpt_weights_3 ( findpt_data_3 *  p,
const real  r[3] 
)
inline

Definition at line 427 of file SpectralFuncs.hpp.

  428 {
  429     lagrange_0( &p->ld[0], r[0] );
  430     lagrange_0( &p->ld[1], r[1] );
  431     lagrange_0( &p->ld[2], r[2] );
  432 }

References lagrange_0(), and findpt_data_3::ld.

◆ gauss_nodes()

void gauss_nodes ( real *  z,
int  n 
)

Definition at line 155 of file poly.c.

  156 {
  157     int i, j;
  158     for( i = 0; i <= n / 2 - 1; ++i )
  159     {
  160         realType ox, x = mbcos( ( 2 * n - 2 * i - 1 ) * ( MOAB_POLY_PI / 2 ) / n );
  161         do
  162         {
  163             ox = x;
  164             x -= legendre( n, x ) / legendre_d1( n, x );
  165         } while( mbabs( x - ox ) > -x * MOAB_POLY_EPS );
  166         z[i] = x - legendre( n, x ) / legendre_d1( n, x );
  167     }
  168     if( n & 1 ) z[n / 2] = 0;
  169     for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
  170         z[j] = -z[i];
  171 }

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

◆ gauss_to_legendre()

void gauss_to_legendre ( const real *  z,
const real *  w,
int  n,
real *  J 
)

Definition at line 239 of file poly.c.

  240 {
  241     int i, j;
  242     legendre_matrix_t( z, n, J, n - 1 );
  243     for( j = 0; j < n; ++j )
  244     {
  245         realType ww = w[j];
  246         for( i = 0; i < n; ++i )
  247             *J++ *= ( 2 * i + 1 ) * ww / 2;
  248     }
  249 }

References legendre_matrix_t().

◆ gauss_to_legendre_t()

void gauss_to_legendre_t ( const real *  z,
const real *  w,
int  n,
real *  J 
)

Definition at line 255 of file poly.c.

  256 {
  257     int i, j;
  258     legendre_matrix( z, n, J, n - 1 );
  259     for( i = 0; i < n; ++i )
  260     {
  261         realType ii = (realType)( 2 * i + 1 ) / 2;
  262         for( j = 0; j < n; ++j )
  263             *J++ *= ii * w[j];
  264     }
  265 }

References legendre_matrix().

◆ gauss_weights()

void gauss_weights ( const real *  z,
real *  w,
int  n 
)

Definition at line 199 of file poly.c.

  200 {
  201     int i, j;
  202     for( i = 0; i <= ( n - 1 ) / 2; ++i )
  203     {
  204         realType d = ( n + 1 ) * legendre( n + 1, z[i] );
  205         w[i]       = 2 * ( 1 - z[i] * z[i] ) / ( d * d );
  206     }
  207     for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
  208         w[j] = w[i];
  209 }

References legendre().

◆ lagrange_0()

void lagrange_0 ( lagrange_data *  p,
real  x 
)

Definition at line 414 of file poly.c.

  415 {
  416     unsigned i, n = p->n;
  417     for( i = 0; i < n; ++i )
  418         p->d[i] = x - p->z[i];
  419     for( i = 0; i < n - 1; ++i )
  420         p->u0[i + 1] = p->d[i] * p->u0[i];
  421     for( i = n - 1; i; --i )
  422         p->v0[i - 1] = p->d[i] * p->v0[i];
  423     for( i = 0; i < n; ++i )
  424         p->J[i] = p->w[i] * p->u0[i] * p->v0[i];
  425 }

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::SpectralHex::evaluate(), moab::Element::SpectralQuad::evaluate(), moab::element_utility::Spectral_hex_map< _Matrix >::evaluate(), moab::Element::SpectralHex::evaluate_scalar_field(), moab::Element::SpectralQuad::evaluate_scalar_field(), moab::element_utility::Spectral_hex_map< _Matrix >::evaluate_scalar_field(), findpt_weights_2(), findpt_weights_3(), and moab::ElemUtil::hex_eval().

◆ lagrange_1()

void lagrange_1 ( lagrange_data *  p,
real  x 
)

Definition at line 427 of file poly.c.

  428 {
  429     unsigned i, n = p->n;
  430     for( i = 0; i < n; ++i )
  431         p->d[i] = x - p->z[i];
  432     for( i = 0; i < n - 1; ++i )
  433         p->u0[i + 1] = p->d[i] * p->u0[i], p->u1[i + 1] = p->d[i] * p->u1[i] + p->u0[i];
  434     for( i = n - 1; i; --i )
  435         p->v0[i - 1] = p->d[i] * p->v0[i], p->v1[i - 1] = p->d[i] * p->v1[i] + p->v0[i];
  436     for( i = 0; i < n; ++i )
  437         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] );
  438 }

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().

◆ lagrange_2()

void lagrange_2 ( lagrange_data *  p,
real  x 
)

Definition at line 440 of file poly.c.

  441 {
  442     unsigned i, n = p->n;
  443     for( i = 0; i < n; ++i )
  444         p->d[i] = x - p->z[i];
  445     for( i = 0; i < n - 1; ++i )
  446         p->u0[i + 1] = p->d[i] * p->u0[i], p->u1[i + 1] = p->d[i] * p->u1[i] + p->u0[i],
  447                   p->u2[i + 1] = p->d[i] * p->u2[i] + 2 * p->u1[i];
  448     for( i = n - 1; i; --i )
  449         p->v0[i - 1] = p->d[i] * p->v0[i], p->v1[i - 1] = p->d[i] * p->v1[i] + p->v0[i],
  450                   p->v2[i - 1] = p->d[i] * p->v2[i] + 2 * p->v1[i];
  451     for( i = 0; i < n; ++i )
  452         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] ),
  453         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] );
  454 }

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().

◆ lagrange_2u()

void lagrange_2u ( lagrange_data *  p)

Definition at line 456 of file poly.c.

  457 {
  458     unsigned i, n = p->n;
  459     for( i = 0; i < n - 1; ++i )
  460         p->u2[i + 1] = p->d[i] * p->u2[i] + 2 * p->u1[i];
  461     for( i = n - 1; i; --i )
  462         p->v2[i - 1] = p->d[i] * p->v2[i] + 2 * p->v1[i];
  463     for( i = 0; i < n; ++i )
  464         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] );
  465 }

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().

◆ lagrange_free()

void lagrange_free ( lagrange_data *  p)

Definition at line 497 of file poly.c.

  498 {
  499     free( p->w );
  500 }

References lagrange_data::w.

Referenced by moab::element_utility::Spectral_hex_map< _Matrix >::free_data(), moab::Element::SpectralHex::freedata(), moab::Element::SpectralQuad::freedata(), moab::ElemUtil::hex_eval(), and moab::ElemUtil::hex_findpt().

◆ lagrange_setup()

void lagrange_setup ( lagrange_data *  p,
const real *  z,
unsigned  n 
)

Definition at line 467 of file poly.c.

  468 {
  469     unsigned i, j;
  470     p->n = n, p->z = z;
  471     p->w = tmalloc( realType, 17 * n );
  472     p->d = p->w + n;
  473     p->J = p->d + n, p->D = p->J + n, p->D2 = p->D + n;
  474     p->u0 = p->D2 + n, p->v0 = p->u0 + n;
  475     p->u1 = p->v0 + n, p->v1 = p->u1 + n;
  476     p->u2 = p->v1 + n, p->v2 = p->u2 + n;
  477     p->J_z0 = p->v2 + n, p->D_z0 = p->J_z0 + n, p->D2_z0 = p->D_z0 + n;
  478     p->J_zn = p->D2_z0 + n, p->D_zn = p->J_zn + n, p->D2_zn = p->D_zn + n;
  479     for( i = 0; i < n; ++i )
  480     {
  481         realType ww = 1, zi = z[i];
  482         for( j = 0; j < i; ++j )
  483             ww *= zi - z[j];
  484         for( ++j; j < n; ++j )
  485             ww *= zi - z[j];
  486         p->w[i] = 1 / ww;
  487     }
  488     p->u0[0] = p->v0[n - 1] = 1;
  489     p->u1[0] = p->v1[n - 1] = 0;
  490     p->u2[0] = p->v2[n - 1] = 0;
  491     lagrange_2( p, z[0] );
  492     memcpy( p->J_z0, p->J, 3 * n * sizeof( realType ) );
  493     lagrange_2( p, z[n - 1] );
  494     memcpy( p->J_zn, p->J, 3 * n * sizeof( realType ) );
  495 }

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< _Matrix >::initialize_spectral_hex().

◆ lagrange_weights()

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.

  321 {
  322     unsigned i, j;
  323     realType *w = work, *d = w + n, *u = d + n, *v = u + n;
  324     for( i = 0; i < n; ++i )
  325     {
  326         realType ww = 1, zi = z[i];
  327         for( j = 0; j < i; ++j )
  328             ww *= zi - z[j];
  329         for( ++j; j < n; ++j )
  330             ww *= zi - z[j];
  331         w[i] = 1 / ww;
  332     }
  333     u[0] = v[n - 1] = 1;
  334     for( i = 0; i < m; ++i )
  335     {
  336         realType xi = x[i];
  337         for( j = 0; j < n; ++j )
  338             d[j] = xi - z[j];
  339         for( j = 0; j < n - 1; ++j )
  340             u[j + 1] = d[j] * u[j];
  341         for( j = n - 1; j; --j )
  342             v[j - 1] = d[j] * v[j];
  343         for( j = 0; j < n; ++j )
  344             *J++ = w[j] * u[j] * v[j];
  345     }
  346 }

◆ lagrange_weights_deriv()

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.

  361 {
  362     unsigned i, j;
  363     realType *w = work, *d = w + n, *u = d + n, *v = u + n, *up = v + n, *vp = up + n;
  364     for( i = 0; i < n; ++i )
  365     {
  366         realType ww = 1, zi = z[i];
  367         for( j = 0; j < i; ++j )
  368             ww *= zi - z[j];
  369         for( ++j; j < n; ++j )
  370             ww *= zi - z[j];
  371         w[i] = 1 / ww;
  372     }
  373     u[0] = v[n - 1] = 1;
  374     up[0] = vp[n - 1] = 0;
  375     for( i = 0; i < m; ++i )
  376     {
  377         realType xi = x[i];
  378         for( j = 0; j < n; ++j )
  379             d[j] = xi - z[j];
  380         for( j = 0; j < n - 1; ++j )
  381             u[j + 1] = d[j] * u[j], up[j + 1] = d[j] * up[j] + u[j];
  382         for( j = n - 1; j; --j )
  383             v[j - 1] = d[j] * v[j], vp[j - 1] = d[j] * vp[j] + v[j];
  384         for( j = 0; j < n; ++j )
  385             *J++ = w[j] * u[j] * v[j], *D++ = w[j] * ( up[j] * v[j] + u[j] * vp[j] );
  386     }
  387 }

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

◆ legendre_matrix()

void legendre_matrix ( const real *  x,
int  m,
real *  P,
int  n 
)

Definition at line 30 of file poly.c.

   31 {
   32     int i, j;
   33     realType *Pjm1 = P, *Pj = Pjm1 + m, *Pjp1 = Pj + m;
   34     for( i = 0; i < m; ++i )
   35         Pjm1[i] = 1;
   36     for( i = 0; i < m; ++i )
   37         Pj[i] = x[i];
   38     for( j = 1; j < n; ++j )
   39     {
   40         realType c = 1 / (realType)( j + 1 ), a = c * ( 2 * j + 1 ), b = c * j;
   41         for( i = 0; i < m; ++i )
   42             Pjp1[i] = a * x[i] * Pj[i] - b * Pjm1[i];
   43         Pjp1 += m, Pj += m, Pjm1 += m;
   44     }
   45 }

Referenced by gauss_to_legendre_t().

◆ legendre_matrix_t()

void legendre_matrix_t ( const real *  x,
int  m,
real *  P,
int  n 
)

Definition at line 87 of file poly.c.

   88 {
   89     int i;
   90     if( n & 1 )
   91         for( i = 0; i < m; ++i, P += n + 1 )
   92             legendre_row_odd( x[i], P, n );
   93     else
   94         for( i = 0; i < m; ++i, P += n + 1 )
   95             legendre_row_even( x[i], P, n );
   96 }

References legendre_row_even(), and legendre_row_odd().

Referenced by gauss_to_legendre().

◆ legendre_row()

void legendre_row ( real  x,
real *  P,
int  n 
)

Definition at line 75 of file poly.c.

   76 {
   77     if( n & 1 )
   78         legendre_row_odd( x, P, n );
   79     else
   80         legendre_row_even( x, P, n );
   81 }

References legendre_row_even(), and legendre_row_odd().

◆ lobatto_nodes()

void lobatto_nodes ( real *  z,
int  n 
)

Definition at line 193 of file poly.c.

  194 {
  195     z[0] = -1, z[n - 1] = 1;
  196     lobatto_nodes_aux( &z[1], n - 2 );
  197 }

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< _Matrix >::initialize_spectral_hex().

◆ lobatto_to_legendre()

void lobatto_to_legendre ( const real *  z,
const real *  w,
int  n,
real *  J 
)

Definition at line 275 of file poly.c.

  276 {
  277     int i, j, m = ( n + 1 ) / 2;
  278     realType *p = J, *q;
  279     realType ww, sum;
  280     if( n & 1 )
  281         for( j = 0; j < m; ++j, p += n )
  282             legendre_row_odd( z[j], p, n - 2 );
  283     else
  284         for( j = 0; j < m; ++j, p += n )
  285             legendre_row_even( z[j], p, n - 2 );
  286     p = J;
  287     for( j = 0; j < m; ++j )
  288     {
  289         ww = w[j], sum = 0;
  290         for( i = 0; i < n - 1; ++i )
  291             *p *= ( 2 * i + 1 ) * ww / 2, sum += *p++;
  292         *p++ = -sum;
  293     }
  294     q = J + ( n / 2 - 1 ) * n;
  295     if( n & 1 )
  296         for( ; j < n; ++j, p += n, q -= n )
  297         {
  298             for( i = 0; i < n - 1; i += 2 )
  299                 p[i] = q[i], p[i + 1] = -q[i + 1];
  300             p[i] = q[i];
  301         }
  302     else
  303         for( ; j < n; ++j, p += n, q -= n )
  304         {
  305             for( i = 0; i < n - 1; i += 2 )
  306                 p[i] = q[i], p[i + 1] = -q[i + 1];
  307         }
  308 }

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

◆ lobatto_weights()

void lobatto_weights ( const real *  z,
real *  w,
int  n 
)

Definition at line 211 of file poly.c.

  212 {
  213     int i, j;
  214     for( i = 0; i <= ( n - 1 ) / 2; ++i )
  215     {
  216         realType d = legendre( n - 1, z[i] );
  217         w[i]       = 2 / ( ( n - 1 ) * n * d * d );
  218     }
  219     for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
  220         w[j] = w[i];
  221 }

References legendre().

Referenced by hash_getbb_2(), and hash_getbb_3().

◆ opt_alloc_2()

void opt_alloc_2 ( opt_data_2 *  p,
lagrange_data *  ld 
)

Definition at line 1662 of file findpt.c.

 1663 {
 1664     const unsigned nr = ld[0].n, ns = ld[1].n, ne = umax_2( nr, ns ), nw = 2 * ns;
 1665     p->size[0]   = 1;
 1666     p->size[1]   = nr;
 1667     p->size[2]   = nr * ns;
 1668     p->ld        = ld;
 1669     p->work      = tmalloc( realType, 4 * ne + nw );
 1670     p->ed.x[0]   = p->work + nw;
 1671     p->ed.x[1]   = p->ed.x[0] + ne;
 1672     p->ed.fd1[0] = p->ed.x[1] + ne;
 1673     p->ed.fd1[1] = p->ed.fd1[0] + ne;
 1674 }

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().

◆ opt_alloc_3()

void opt_alloc_3 ( opt_data_3 *  p,
lagrange_data *  ld 
)

Definition at line 1251 of file findpt.c.

 1252 {
 1253     const unsigned nr = ld[0].n, ns = ld[1].n, nt = ld[2].n, nf = umax_3( nr * ns, nr * nt, ns * nt ),
 1254                    ne = umax_3( nr, ns, nt ), nw = 2 * ns * nt + 3 * ns;
 1255     p->size[0]   = 1;
 1256     p->size[1]   = nr;
 1257     p->size[2]   = nr * ns;
 1258     p->size[3]   = p->size[2] * nt;
 1259     p->ld        = ld;
 1260     p->work      = tmalloc( realType, 6 * nf + 9 * ne + nw );
 1261     p->fd.x[0]   = p->work + nw;
 1262     p->fd.x[1]   = p->fd.x[0] + nf;
 1263     p->fd.x[2]   = p->fd.x[1] + nf;
 1264     p->fd.fdn[0] = p->fd.x[2] + nf;
 1265     p->fd.fdn[1] = p->fd.fdn[0] + nf;
 1266     p->fd.fdn[2] = p->fd.fdn[1] + nf;
 1267     p->ed.x[0]   = p->fd.fdn[2] + nf;
 1268     p->ed.x[1]   = p->ed.x[0] + ne;
 1269     p->ed.x[2]   = p->ed.x[1] + ne;
 1270     p->ed.fd1[0] = p->ed.x[2] + ne;
 1271     p->ed.fd1[1] = p->ed.fd1[0] + ne;
 1272     p->ed.fd1[2] = p->ed.fd1[1] + ne;
 1273     p->ed.fd2[0] = p->ed.fd1[2] + ne;
 1274     p->ed.fd2[1] = p->ed.fd2[0] + ne;
 1275     p->ed.fd2[2] = p->ed.fd2[1] + ne;
 1276 }

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< _Matrix >::initialize_spectral_hex().

◆ opt_findpt_2()

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.

 1823 {
 1824     realType dr[2], resid[2], steep[2];
 1825  
 1826     unsigned c = *constr, ac, d, cc[2], step = 0;
 1827  
 1828     p->elx[0] = elx[0], p->elx[1] = elx[1];
 1829  
 1830     p->ed.constraints = opt_no_constraints_2;
 1831     p->pd.constraints = opt_no_constraints_2;
 1832  
 1833 #if DIAGNOSTICS
 1834     printf( "opt_findpt: xstar = %g, %g\n", xstar[0], xstar[1] );
 1835 #endif
 1836  
 1837     do
 1838     {
 1839         ++step;
 1840         if( step == 50 ) /*fail("%s: opt_findpt_2 did not converge\n",__FILE__);*/
 1841             return 1.e+30;
 1842 #if DIAGNOSTICS
 1843         printf( "  iteration %u\n", step );
 1844         printf( "    %d constraint(s) active\n", (int)opt_constr_num_2[c] );
 1845 #endif
 1846         /* update face/edge/point data if necessary,
 1847            and evaluate x(r) as well as the jacobian */
 1848         switch( opt_constr_num_2[c] )
 1849         {
 1850             case 0:
 1851                 opt_area_set_intp_2( p, r );
 1852                 break;
 1853             case 1:
 1854                 opt_edge_set_intp_2( p, r, c );
 1855                 break;
 1856             case 2:
 1857                 opt_point_set_intp_2( p, c );
 1858                 break;
 1859         }
 1860 #if DIAGNOSTICS
 1861         printf( "    r = %g, %g\n", r[0], r[1] );
 1862         printf( "    x = %g, %g\n", p->x[0], p->x[1] );
 1863 #endif
 1864         /* compute residual */
 1865         for( d = 0; d < 2; ++d )
 1866             resid[d] = xstar[d] - p->x[d];
 1867 #if DIAGNOSTICS
 1868         printf( "    resid = %g, %g\n", resid[0], resid[1] );
 1869         printf( "    2-norm = %g\n", r2norm_2( resid[0], resid[1] ) );
 1870 #endif
 1871         /* check constraints against steepest descent direction */
 1872         ac = c;
 1873         if( opt_constr_num_2[c] )
 1874         {
 1875             opt_constr_unpack_2( c, cc );
 1876             mat_app_2c( steep, p->jac, resid ); /* steepest descent = J^T r */
 1877 #if DIAGNOSTICS
 1878             printf( "    steepest descent = %g, %g\n", steep[0], steep[1] );
 1879 #endif
 1880             for( d = 0; d < 2; ++d )
 1881                 if( ( cc[d] == 0 && steep[d] > 0 ) || ( cc[d] == 2 && steep[d] < 0 ) ) cc[d] = 1;
 1882             ac = opt_constr_pack_2( cc );
 1883         }
 1884         /* update face/edge/point data if necessary */
 1885         if( ac != c )
 1886         {
 1887             c = ac;
 1888 #if DIAGNOSTICS
 1889             printf( "    relaxed to %d constraints\n", (int)opt_constr_num_2[c] );
 1890 #endif
 1891             switch( opt_constr_num_2[c] )
 1892             {
 1893                 case 1:
 1894                     opt_edge_set_2( p, r, c );
 1895                     break;
 1896                 case 2:
 1897                     opt_point_set_2( p, c );
 1898                     break;
 1899             }
 1900         }
 1901         /* compute Newton step */
 1902         switch( opt_constr_num_2[c] )
 1903         {
 1904             case 0:
 1905                 tinyla_solve_2( dr, p->jac, resid );
 1906                 break;
 1907             case 1: {
 1908                 const unsigned de = p->ed.de, d1 = p->ed.d1;
 1909                 realType fac, H[2];
 1910                 const realType* J = p->jac + de;
 1911                 opt_edge_hess_2( p, H );
 1912                 fac    = J[0] * J[0] + J[2] * J[2] - ( resid[0] * H[0] + resid[1] * H[1] );
 1913                 dr[de] = steep[de] / fac;
 1914                 dr[d1] = 0;
 1915             }
 1916             break;
 1917             case 2:
 1918                 dr[0] = dr[1] = 0;
 1919                 break;
 1920         }
 1921 #if DIAGNOSTICS
 1922         printf( "    dr = %g, %g\n", dr[0], dr[1] );
 1923 #endif
 1924         /* project new iteration onto [-1,1]^2 */
 1925         opt_constr_unpack_2( c, cc );
 1926         for( d = 0; d < 2; ++d )
 1927         {
 1928             if( cc[d] != 1 ) continue;
 1929             r[d] += dr[d];
 1930             if( r[d] <= -1 )
 1931                 dr[d] -= r[d] + 1, r[d] = -1, cc[d] = 0;
 1932             else if( r[d] >= 1 )
 1933                 dr[d] -= r[d] - 1, r[d] = 1, cc[d] = 2;
 1934         }
 1935         c = opt_constr_pack_2( cc );
 1936     } while( r1norm_2( dr[0], dr[1] ) > 2 * MOAB_POLY_EPS );
 1937     *constr = c;
 1938     return r2norm_2( resid[0], resid[1] );
 1939 }

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().

◆ opt_findpt_3()

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.

 1517 {
 1518     realType dr[3], resid[3], steep[3];
 1519  
 1520     unsigned c = *constr, ac, d, cc[3], step = 0;
 1521  
 1522     p->elx[0] = elx[0], p->elx[1] = elx[1], p->elx[2] = elx[2];
 1523  
 1524     p->fd.constraints = opt_no_constraints_3;
 1525     p->ed.constraints = opt_no_constraints_3;
 1526     p->pd.constraints = opt_no_constraints_3;
 1527  
 1528 #if DIAGNOSTICS
 1529     printf( "opt_findpt: xstar = %g, %g, %g\n", xstar[0], xstar[1], xstar[2] );
 1530 #endif
 1531  
 1532     do
 1533     {
 1534         ++step;
 1535         if( step == 50 ) /*fail("%s: opt_findpt_3 did not converge\n",__FILE__);*/
 1536             return 1.e+30;
 1537 #if DIAGNOSTICS
 1538         printf( "  iteration %u\n", step );
 1539         printf( "    %d constraint(s) active\n", (int)opt_constr_num_3[c] );
 1540 #endif
 1541         /* update face/edge/point data if necessary,
 1542            and evaluate x(r) as well as the jacobian */
 1543         switch( opt_constr_num_3[c] )
 1544         {
 1545             case 0:
 1546                 opt_vol_set_intp_3( p, r );
 1547                 break;
 1548             case 1:
 1549                 opt_face_set_intp_3( p, r, c );
 1550                 break;
 1551             case 2:
 1552                 opt_edge_set_intp_3( p, r, c );
 1553                 break;
 1554             case 3:
 1555                 opt_point_set_intp_3( p, c );
 1556                 break;
 1557         }
 1558 #if DIAGNOSTICS
 1559         printf( "    r = %g, %g, %g\n", r[0], r[1], r[2] );
 1560         printf( "    x = %g, %g, %g\n", p->x[0], p->x[1], p->x[2] );
 1561 #endif
 1562         /* compute residual */
 1563         for( d = 0; d < 3; ++d )
 1564             resid[d] = xstar[d] - p->x[d];
 1565 #if DIAGNOSTICS
 1566         printf( "    resid = %g, %g, %g\n", resid[0], resid[1], resid[2] );
 1567         printf( "    2-norm = %g\n", r2norm_3( resid[0], resid[1], resid[2] ) );
 1568 #endif
 1569         /* check constraints against steepest descent direction */
 1570         ac = c;
 1571         if( opt_constr_num_3[c] )
 1572         {
 1573             opt_constr_unpack_3( c, cc );
 1574             mat_app_3c( steep, p->jac, resid ); /* steepest descent = J^T r */
 1575 #if DIAGNOSTICS
 1576             printf( "    steepest descent = %g, %g, %g\n", steep[0], steep[1], steep[2] );
 1577 #endif
 1578             for( d = 0; d < 3; ++d )
 1579                 if( ( cc[d] == 0 && steep[d] > 0 ) || ( cc[d] == 2 && steep[d] < 0 ) ) cc[d] = 1;
 1580             ac = opt_constr_pack_3( cc );
 1581         }
 1582         /* update face/edge/point data if necessary */
 1583         if( ac != c )
 1584         {
 1585             c = ac;
 1586 #if DIAGNOSTICS
 1587             printf( "    relaxed to %d constraints\n", (int)opt_constr_num_3[c] );
 1588 #endif
 1589             switch( opt_constr_num_3[c] )
 1590             {
 1591                 case 1:
 1592                     opt_face_set_3( p, r, c );
 1593                     break;
 1594                 case 2:
 1595                     opt_edge_set_3( p, r, c );
 1596                     break;
 1597                 case 3:
 1598                     opt_point_set_3( p, c );
 1599                     break;
 1600             }
 1601         }
 1602         /* compute Newton step */
 1603         switch( opt_constr_num_3[c] )
 1604         {
 1605             case 0:
 1606                 tinyla_solve_3( dr, p->jac, resid );
 1607                 break;
 1608             case 1: {
 1609                 const unsigned dn = p->fd.dn, d1 = p->fd.d1, d2 = p->fd.d2;
 1610                 realType A[4], H[9];
 1611                 const realType* J = p->jac;
 1612                 opt_face_hess_3( p, H );
 1613                 A[0] = J[d1] * J[d1] + J[3 + d1] * J[3 + d1] + J[6 + d1] * J[6 + d1];
 1614                 A[1] = J[d2] * J[d2] + J[3 + d2] * J[3 + d2] + J[6 + d2] * J[6 + d2];
 1615                 A[2] = J[d1] * J[d2] + J[3 + d1] * J[3 + d2] + J[6 + d1] * J[6 + d2];
 1616                 A[0] -= resid[0] * H[0] + resid[1] * H[3] + resid[2] * H[6];
 1617                 A[1] -= resid[0] * H[1] + resid[1] * H[4] + resid[2] * H[7];
 1618                 A[2] -= resid[0] * H[2] + resid[1] * H[5] + resid[2] * H[8];
 1619                 tinyla_solve_sym_2( &dr[d1], &dr[d2], A, steep[d1], steep[d2] );
 1620                 dr[dn] = 0;
 1621             }
 1622             break;
 1623             case 2: {
 1624                 const unsigned de = p->ed.de, d1 = p->ed.d1, d2 = p->ed.d2;
 1625                 realType fac, H[3];
 1626                 const realType* J = p->jac + de;
 1627                 opt_edge_hess_3( p, H );
 1628                 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] );
 1629                 dr[de] = steep[de] / fac;
 1630                 dr[d1] = 0, dr[d2] = 0;
 1631             }
 1632             break;
 1633             case 3:
 1634                 dr[0] = dr[1] = dr[2] = 0;
 1635                 break;
 1636         }
 1637 #if DIAGNOSTICS
 1638         printf( "    dr = %g, %g, %g\n", dr[0], dr[1], dr[2] );
 1639 #endif
 1640         /* project new iteration onto [-1,1]^3 */
 1641         opt_constr_unpack_3( c, cc );
 1642         for( d = 0; d < 3; ++d )
 1643         {
 1644             if( cc[d] != 1 ) continue;
 1645             r[d] += dr[d];
 1646             if( r[d] <= -1 )
 1647                 dr[d] -= r[d] + 1, r[d] = -1, cc[d] = 0;
 1648             else if( r[d] >= 1 )
 1649                 dr[d] -= r[d] - 1, r[d] = 1, cc[d] = 2;
 1650         }
 1651         c = opt_constr_pack_3( cc );
 1652     } while( r1norm_3( dr[0], dr[1], dr[2] ) > 3 * MOAB_POLY_EPS );
 1653     *constr = c;
 1654 #if 0
 1655   printf("opt_findpt_3 converged in %u iterations\n", step);
 1656 #endif
 1657     return r2norm_3( resid[0], resid[1], resid[2] );
 1658 }

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().

◆ opt_free_2()

void opt_free_2 ( opt_data_2 *  p)

Definition at line 1676 of file findpt.c.

 1677 {
 1678     free( p->work );
 1679 }

References opt_data_2::work.

Referenced by findpt_free_2(), and moab::Element::SpectralQuad::freedata().

◆ opt_free_3()

void opt_free_3 ( opt_data_3 *  p)

Definition at line 1278 of file findpt.c.

 1279 {
 1280     free( p->work );
 1281 }

References opt_data_3::work.

Referenced by findpt_free_3(), moab::element_utility::Spectral_hex_map< _Matrix >::free_data(), moab::Element::SpectralHex::freedata(), and moab::ElemUtil::hex_findpt().

◆ opt_vol_set_intp_3()

void opt_vol_set_intp_3 ( opt_data_3 *  p,
const real  r[3] 
)

Definition at line 1301 of file findpt.c.

 1302 {
 1303     opt_vol_set_3( p, r );
 1304     opt_vol_intp_3( p );
 1305 }

References opt_vol_intp_3(), and opt_vol_set_3().

Referenced by moab::Element::SpectralHex::integrate_scalar_field(), moab::element_utility::Spectral_hex_map< _Matrix >::integrate_scalar_field(), moab::Element::SpectralHex::jacobian(), moab::element_utility::Spectral_hex_map< _Matrix >::jacobian(), and opt_findpt_3().

◆ tensor_c1()

void tensor_c1 ( const real *  R,
unsigned  mr,
unsigned  nr,
const real *  u,
real *  v 
)

Definition at line 155 of file tensor.c.

  156 {
  157     mxv_c( v, mr, R, u, nr );
  158 }

References mxv_c(), nr, and moab::R.

◆ tensor_c2()

void tensor_c2 ( const real *  R,
unsigned  mr,
unsigned  nr,
const real *  S,
unsigned  ms,
unsigned  ns,
const real *  u,
real *  v,
real *  work 
)

Definition at line 166 of file tensor.c.

  175 {
  176     mxm_cc( R, mr, u, nr, W, ns );
  177     mxm_cr( W, mr, S, ns, v, ms );
  178 }

References mxm_cc(), mxm_cr(), nr, and moab::R.

◆ tensor_c3()

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 
)

Definition at line 197 of file tensor.c.

  210 {
  211     unsigned n, mrns = mr * ns, mrms = mr * ms;
  212     realType* Zp = Z;
  213     mxm_cc( R, mr, u, nr, W, ns * nt );
  214     for( n = 0; n < nt; ++n, W += mrns, Zp += mrms )
  215         mxm_cr( W, mr, S, ns, Zp, ms );
  216     mxm_cr( Z, mrms, T, nt, v, mt );
  217 }

References mxm_cc(), mxm_cr(), nr, and moab::R.

◆ tensor_i1()

real tensor_i1 ( const real *  Jr,
unsigned  nr,
const real *  u 
)

Definition at line 255 of file tensor.c.

  256 {
  257     return inner( Jr, u, nr );
  258 }

References inner(), and nr.

Referenced by opt_edge_hess_2(), opt_edge_hess_3(), opt_edge_intp_2(), and opt_edge_intp_3().

◆ tensor_i2()

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.

  267 {
  268     mxv_r( work, ns, u, Jr, nr );
  269     return inner( Js, work, ns );
  270 }

References inner(), mxv_r(), and nr.

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().

◆ tensor_i3()

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.

  281 {
  282     realType* work2 = work + nt;
  283     mxv_r( work2, ns * nt, u, Jr, nr );
  284     mxv_r( work, nt, work2, Js, ns );
  285     return inner( Jt, work, nt );
  286 }

References inner(), mxv_r(), and nr.

Referenced by moab::Element::SpectralHex::evaluate(), moab::element_utility::Spectral_hex_map< _Matrix >::evaluate(), moab::Element::SpectralHex::evaluate_scalar_field(), moab::element_utility::Spectral_hex_map< _Matrix >::evaluate_scalar_field(), findpt_eval_3(), and moab::ElemUtil::hex_eval().

◆ tensor_ig1()

real tensor_ig1 ( const real *  Jr,
const real *  Dr,
unsigned  nr,
const real *  u,
real *  g 
)

Definition at line 304 of file tensor.c.

  305 {
  306     *g = inner( Dr, u, nr );
  307     return inner( Jr, u, nr );
  308 }

References inner(), and nr.

Referenced by opt_edge_intp_2(), and opt_edge_intp_3().

◆ tensor_ig2()

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.

  320 {
  321     realType *a = work, *ar = a + ns;
  322     mxv_r( a, ns, u, Jr, nr );
  323     mxv_r( ar, ns, u, Dr, nr );
  324     g[0] = inner( Js, ar, ns );
  325     g[1] = inner( Ds, a, ns );
  326     return inner( Js, a, ns );
  327 }

References inner(), mxv_r(), and nr.

Referenced by obbox_calc_2(), opt_area_intp_2(), opt_face_hess_3(), and opt_face_intp_3().

◆ tensor_ig3()

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.

  342 {
  343     unsigned nsnt = ns * nt;
  344     realType *a = work, *ar = a + nsnt, *b = ar + nsnt, *br = b + ns, *bs = br + ns;
  345     mxv_r( a, nsnt, u, Jr, nr );
  346     mxv_r( ar, nsnt, u, Dr, nr );
  347     mxv_r( b, nt, a, Js, ns );
  348     mxv_r( br, nt, ar, Js, ns );
  349     mxv_r( bs, nt, a, Ds, ns );
  350     g[0] = inner( Jt, br, nt );
  351     g[1] = inner( Jt, bs, nt );
  352     g[2] = inner( Dt, b, nt );
  353     return inner( Jt, b, nt );
  354 }

References inner(), mxv_r(), and nr.

Referenced by obbox_calc_3(), and opt_vol_intp_3().

◆ tensor_r1()

void tensor_r1 ( const real *  R,
unsigned  mr,
unsigned  nr,
const real *  u,
real *  v 
)

Definition at line 160 of file tensor.c.

  161 {
  162     mxv_r( v, mr, R, u, nr );
  163 }

References mxv_r(), nr, and moab::R.

◆ tensor_r2()

void tensor_r2 ( const real *  R,
unsigned  mr,
unsigned  nr,
const real *  S,
unsigned  ms,
unsigned  ns,
const real *  u,
real *  v,
real *  work 
)

Definition at line 181 of file tensor.c.

  190 {
  191     mxm_rc( R, mr, u, nr, W, ns );
  192     mxm_cc( W, mr, S, ns, v, ms );
  193 }

References mxm_cc(), mxm_rc(), nr, and moab::R.

◆ tensor_r3()

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 
)

Definition at line 221 of file tensor.c.

  234 {
  235     unsigned n, mrns = mr * ns, mrms = mr * ms;
  236     realType* Zp = Z;
  237     mxm_rc( R, mr, u, nr, W, ns * nt );
  238     for( n = 0; n < nt; ++n, W += mrns, Zp += mrms )
  239         mxm_cc( W, mr, S, ns, Zp, ms );
  240     mxm_cc( Z, mrms, T, nt, v, mt );
  241 }

References mxm_cc(), mxm_rc(), nr, and moab::R.

Variable Documentation

◆ opt_no_constraints_2

const unsigned opt_no_constraints_2 = 3 + 1

Definition at line 485 of file SpectralFuncs.hpp.

◆ opt_no_constraints_3

const unsigned opt_no_constraints_3 = 9 + 3 + 1

Definition at line 486 of file SpectralFuncs.hpp.