Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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_2findpt_setup_2 (const realType *const xw[2], const unsigned n[2], uint nel, uint max_hash_size, realType bbox_tol)
 
findpt_data_3findpt_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 realTypefindpt_allbnd_2 (const findpt_data_2 *p)
 
const realTypefindpt_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 
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 
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