Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
Classes | Macros | Typedefs | Functions | Variables
FindPtFuncs.h File Reference
#include "float.h"
#include "math.h"
#include <stdlib.h>
+ Include dependency graph for FindPtFuncs.h:
+ This graph shows which files directly or indirectly include this file:

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
 

Macros

#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
 

Macro Definition Documentation

◆ DECLMINMAX

#define DECLMINMAX (   type,
  prefix 
)

Definition at line 679 of file FindPtFuncs.h.

◆ GLOBAL_INT

#define GLOBAL_INT   long

Definition at line 47 of file FindPtFuncs.h.

◆ INTEGER

#define INTEGER   int

Definition at line 38 of file FindPtFuncs.h.

◆ MAYBE_UNUSED [1/2]

#define MAYBE_UNUSED

Definition at line 676 of file FindPtFuncs.h.

◆ MAYBE_UNUSED [2/2]

#define MAYBE_UNUSED

Definition at line 676 of file FindPtFuncs.h.

◆ mbabs

#define mbabs   fabs

Definition at line 16 of file FindPtFuncs.h.

◆ mbceil

#define mbceil   ceil

Definition at line 20 of file FindPtFuncs.h.

◆ mbcos

#define mbcos   cos

Definition at line 17 of file FindPtFuncs.h.

◆ mbfloor

#define mbfloor   floor

Definition at line 19 of file FindPtFuncs.h.

◆ mbsin

#define mbsin   sin

Definition at line 18 of file FindPtFuncs.h.

◆ mbsqrt

#define mbsqrt   sqrt

Definition at line 15 of file FindPtFuncs.h.

◆ MOAB_POLY_EPS

#define MOAB_POLY_EPS   ( 128 * DBL_EPSILON )

Definition at line 12 of file FindPtFuncs.h.

◆ MOAB_POLY_PI

#define MOAB_POLY_PI   3.1415926535897932384626433832795028841971693993751058209749445923

Definition at line 13 of file FindPtFuncs.h.

◆ tcalloc

#define tcalloc (   type,
  count 
)    ( (type*)scalloc( ( count ), sizeof( type ), __FILE__ ) )

Definition at line 636 of file FindPtFuncs.h.

◆ tmalloc

#define tmalloc (   type,
  count 
)    ( (type*)smalloc( ( count ) * sizeof( type ), __FILE__ ) )

Definition at line 635 of file FindPtFuncs.h.

◆ trealloc

#define trealloc (   type,
  ptr,
  count 
)    ( (type*)srealloc( ( ptr ), ( count ) * sizeof( type ), __FILE__ ) )

Definition at line 637 of file FindPtFuncs.h.

Typedef Documentation

◆ findpt_func

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

Definition at line 538 of file FindPtFuncs.h.

◆ realType

typedef double realType

Definition at line 56 of file FindPtFuncs.h.

◆ sint

typedef signed INTEGER sint

Definition at line 61 of file FindPtFuncs.h.

◆ slong

typedef signed GLOBAL_INT slong

Definition at line 71 of file FindPtFuncs.h.

◆ uint

typedef unsigned INTEGER uint

Definition at line 65 of file FindPtFuncs.h.

◆ ulong

typedef unsigned GLOBAL_INT ulong

Definition at line 79 of file FindPtFuncs.h.

Function Documentation

◆ fail()

void fail ( const char *  fmt,
  ... 
)

Definition at line 6 of file errmem.c.

7 {
8  va_list ap;
9  va_start( ap, fmt );
10  vfprintf( stderr, fmt, ap );
11  va_end( ap );
12  exit( 1 );
13 }

Referenced by main(), ZoltanPartitioner::mbGlobalSuccess(), scalloc(), smalloc(), and srealloc().

◆ findpt_2()

int findpt_2 ( findpt_data_2 *  p,
const realType  x[2],
int  guess,
uint *  el,
realType  r[2],
realType *  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 }

◆ findpt_3()

int findpt_3 ( findpt_data_3 *  p,
const realType  x[3],
int  guess,
uint *  el,
realType  r[3],
realType *  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 }

◆ findpt_allbnd_2()

const realType* findpt_allbnd_2 ( const findpt_data_2 *  p)

Definition at line 2107 of file findpt.c.

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

◆ findpt_allbnd_3()

const realType* findpt_allbnd_3 ( const findpt_data_3 *  p)

Definition at line 2112 of file findpt.c.

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

◆ findpt_eval_2()

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

Definition at line 2283 of file findpt.c.

2284 {
2285  return tensor_i2( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, u, p->od_work );
2286 }

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 realType *  u 
)
inline

Definition at line 2288 of file findpt.c.

2289 {
2290  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 );
2291 }

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 }

◆ 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 }

◆ findpt_setup_2()

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 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 }

◆ findpt_setup_3()

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 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 }

◆ findpt_weights_2()

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

Definition at line 2270 of file findpt.c.

2271 {
2272  lagrange_0( &p->ld[0], r[0] );
2273  lagrange_0( &p->ld[1], r[1] );
2274 }

References lagrange_0(), and findpt_data_2::ld.

◆ findpt_weights_3()

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

Definition at line 2276 of file findpt.c.

2277 {
2278  lagrange_0( &p->ld[0], r[0] );
2279  lagrange_0( &p->ld[1], r[1] );
2280  lagrange_0( &p->ld[2], r[2] );
2281 }

References lagrange_0(), and findpt_data_3::ld.

◆ gauss_nodes()

void gauss_nodes ( realType *  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 }

◆ gauss_to_legendre()

void gauss_to_legendre ( const realType *  z,
const realType *  w,
int  n,
realType *  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 }

◆ gauss_to_legendre_t()

void gauss_to_legendre_t ( const realType *  z,
const realType *  w,
int  n,
realType *  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 }

◆ gauss_weights()

void gauss_weights ( const realType *  z,
realType *  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 }

◆ lagrange_0()

void lagrange_0 ( lagrange_data *  p,
realType  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 }

◆ lagrange_1()

void lagrange_1 ( lagrange_data *  p,
realType  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 }

◆ lagrange_2()

void lagrange_2 ( lagrange_data *  p,
realType  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 }

◆ 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 }

◆ lagrange_free()

void lagrange_free ( lagrange_data *  p)

Definition at line 497 of file poly.c.

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

◆ lagrange_setup()

void lagrange_setup ( lagrange_data *  p,
const realType *  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 }

◆ lagrange_weights()

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.

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 realType *  z,
unsigned  n,
const realType *  x,
unsigned  m,
realType *  J,
realType *  D,
realType *  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 }

◆ legendre_matrix()

void legendre_matrix ( const realType *  x,
int  m,
realType *  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 }

◆ legendre_matrix_t()

void legendre_matrix_t ( const realType *  x,
int  m,
realType *  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 }

◆ legendre_row()

void legendre_row ( realType  x,
realType *  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 }

◆ lobatto_nodes()

void lobatto_nodes ( realType *  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 }

◆ lobatto_to_legendre()

void lobatto_to_legendre ( const realType *  z,
const realType *  w,
int  n,
realType *  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 }

◆ lobatto_weights()

void lobatto_weights ( const realType *  z,
realType *  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 }

◆ 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 }

◆ 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 }

◆ opt_findpt_2()

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 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 }

◆ opt_findpt_3()

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 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 }

◆ 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 }

◆ 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 }

◆ opt_vol_set_intp_3()

void opt_vol_set_intp_3 ( opt_data_3 *  p,
const realType  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 }

◆ r1norm_1()

static realType r1norm_1 ( realType  a)
static

Definition at line 723 of file FindPtFuncs.h.

724 {
725  return mbabs( a );
726 }

References mbabs.

◆ r1norm_2()

static realType r1norm_2 ( realType  a,
realType  b 
)
static

Definition at line 728 of file FindPtFuncs.h.

729 {
730  return mbabs( a ) + mbabs( b );
731 }

References mbabs.

Referenced by findpt_hash_2(), and opt_findpt_2().

◆ r1norm_3()

static realType r1norm_3 ( realType  a,
realType  b,
realType  c 
)
static

Definition at line 733 of file FindPtFuncs.h.

734 {
735  return mbabs( a ) + mbabs( b ) + mbabs( c );
736 }

References mbabs.

Referenced by findpt_hash_3(), and opt_findpt_3().

◆ r2norm_1()

static realType r2norm_1 ( realType  a)
static

Definition at line 738 of file FindPtFuncs.h.

739 {
740  return mbsqrt( a * a );
741 }

References mbsqrt.

◆ r2norm_2()

static realType r2norm_2 ( realType  a,
realType  b 
)
static

Definition at line 743 of file FindPtFuncs.h.

744 {
745  return mbsqrt( a * a + b * b );
746 }

References mbsqrt.

Referenced by findpt_pass_2(), and opt_findpt_2().

◆ r2norm_3()

static realType r2norm_3 ( realType  a,
realType  b,
realType  c 
)
static

Definition at line 748 of file FindPtFuncs.h.

749 {
750  return mbsqrt( a * a + b * b + c * c );
751 }

References mbsqrt.

Referenced by findpt_pass_3(), and opt_findpt_3().

◆ scalloc()

static void * scalloc ( size_t  nmemb,
size_t  size,
const char *  file 
)
static

Definition at line 620 of file FindPtFuncs.h.

621 {
622  void* res = calloc( nmemb, size );
623  if( !res && nmemb ) fail( "%s: allocation of %d bytes failed\n", file, (int)size * nmemb );
624  return res;
625 }

References fail(), and size.

◆ smalloc()

static void * smalloc ( size_t  size,
const char *  file 
)
static

Definition at line 612 of file FindPtFuncs.h.

613 {
614  void* res = malloc( size );
615  if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
616  return res;
617 }

References fail(), and size.

◆ srealloc()

static void * srealloc ( void *  ptr,
size_t  size,
const char *  file 
)
static

Definition at line 628 of file FindPtFuncs.h.

629 {
630  void* res = realloc( ptr, size );
631  if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
632  return res;
633 }

References fail(), and size.

◆ tensor_c1()

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

Definition at line 155 of file tensor.c.

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

◆ tensor_c2()

void tensor_c2 ( const realType *  R,
unsigned  mr,
unsigned  nr,
const realType *  S,
unsigned  ms,
unsigned  ns,
const realType *  u,
realType *  v,
realType *  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 }

◆ tensor_c3()

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 
)

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 }

◆ tensor_i1()

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

Definition at line 255 of file tensor.c.

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

◆ tensor_i2()

realType tensor_i2 ( const realType *  Jr,
unsigned  nr,
const realType *  Js,
unsigned  ns,
const realType *  u,
realType *  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 }

◆ tensor_i3()

realType tensor_i3 ( const realType *  Jr,
unsigned  nr,
const realType *  Js,
unsigned  ns,
const realType *  Jt,
unsigned  nt,
const realType *  u,
realType *  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 }

◆ tensor_ig1()

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

Definition at line 304 of file tensor.c.

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

◆ tensor_ig2()

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 
)

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 }

◆ tensor_ig3()

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 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 }

◆ tensor_r1()

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

Definition at line 160 of file tensor.c.

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

◆ tensor_r2()

void tensor_r2 ( const realType *  R,
unsigned  mr,
unsigned  nr,
const realType *  S,
unsigned  ms,
unsigned  ns,
const realType *  u,
realType *  v,
realType *  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 }

◆ tensor_r3()

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 
)

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 }

Variable Documentation

◆ opt_no_constraints_2

const unsigned opt_no_constraints_2
extern

Definition at line 10 of file findpt.c.

Referenced by findpt_guess_2(), findpt_pass_2(), and opt_findpt_2().

◆ opt_no_constraints_3

const unsigned opt_no_constraints_3
extern

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