Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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

Defines

#define MOAB_POLY_EPS   ( 128 * DBL_EPSILON )
#define MOAB_POLY_PI   3.1415926535897932384626433832795028841971693993751058209749445923
#define mbsqrt   sqrt
#define mbabs   fabs
#define mbcos   cos
#define mbsin   sin
#define mbfloor   floor
#define mbceil   ceil
#define INTEGER   int
#define GLOBAL_INT   long
#define MAYBE_UNUSED
#define tmalloc(type, count)   ( (type*)smalloc( ( count ) * sizeof( type ), __FILE__ ) )
#define tcalloc(type, count)   ( (type*)scalloc( ( count ), sizeof( type ), __FILE__ ) )
#define trealloc(type, ptr, count)   ( (type*)srealloc( ( ptr ), ( count ) * sizeof( type ), __FILE__ ) )
#define MAYBE_UNUSED
#define DECLMINMAX(type, prefix)

Typedefs

typedef double realType
typedef signed INTEGER sint
typedef unsigned INTEGER uint
typedef signed GLOBAL_INT slong
typedef unsigned GLOBAL_INT ulong
typedef int(* findpt_func )(void *, const realType *, int, uint *, realType *, realType *)

Functions

void legendre_matrix (const realType *x, int m, realType *P, int n)
void legendre_matrix_t (const realType *x, int m, realType *P, int n)
void legendre_row (realType x, realType *P, int n)
void gauss_nodes (realType *z, int n)
void lobatto_nodes (realType *z, int n)
void gauss_weights (const realType *z, realType *w, int n)
void lobatto_weights (const realType *z, realType *w, int n)
void gauss_to_legendre (const realType *z, const realType *w, int n, realType *J)
void gauss_to_legendre_t (const realType *z, const realType *w, int n, realType *J)
void lobatto_to_legendre (const realType *z, const realType *w, int n, realType *J)
void lagrange_weights (const realType *z, unsigned n, const realType *x, unsigned m, realType *J, realType *work)
void lagrange_weights_deriv (const realType *z, unsigned n, const realType *x, unsigned m, realType *J, realType *D, realType *work)
void lagrange_setup (lagrange_data *p, const realType *z, unsigned n)
void lagrange_free (lagrange_data *p)
void lagrange_0 (lagrange_data *p, realType x)
void lagrange_1 (lagrange_data *p, realType x)
void lagrange_2 (lagrange_data *p, realType x)
void lagrange_2u (lagrange_data *p)
void tensor_c1 (const realType *R, unsigned mr, unsigned nr, const realType *u, realType *v)
void tensor_r1 (const realType *R, unsigned mr, unsigned nr, const realType *u, realType *v)
void tensor_c2 (const realType *R, unsigned mr, unsigned nr, const realType *S, unsigned ms, unsigned ns, const realType *u, realType *v, realType *work)
void tensor_r2 (const realType *R, unsigned mr, unsigned nr, const realType *S, unsigned ms, unsigned ns, const realType *u, realType *v, realType *work)
void tensor_c3 (const realType *R, unsigned mr, unsigned nr, const realType *S, unsigned ms, unsigned ns, const realType *T, unsigned mt, unsigned nt, const realType *u, realType *v, realType *work1, realType *work2)
void tensor_r3 (const realType *R, unsigned mr, unsigned nr, const realType *S, unsigned ms, unsigned ns, const realType *T, unsigned mt, unsigned nt, const realType *u, realType *v, realType *work1, realType *work2)
realType tensor_i1 (const realType *Jr, unsigned nr, const realType *u)
realType tensor_i2 (const realType *Jr, unsigned nr, const realType *Js, unsigned ns, const realType *u, realType *work)
realType tensor_i3 (const realType *Jr, unsigned nr, const realType *Js, unsigned ns, const realType *Jt, unsigned nt, const realType *u, realType *work)
realType tensor_ig1 (const realType *Jr, const realType *Dr, unsigned nr, const realType *u, realType *g)
realType tensor_ig2 (const realType *Jr, const realType *Dr, unsigned nr, const realType *Js, const realType *Ds, unsigned ns, const realType *u, realType *g, realType *work)
realType tensor_ig3 (const realType *Jr, const realType *Dr, unsigned nr, const realType *Js, const realType *Ds, unsigned ns, const realType *Jt, const realType *Dt, unsigned nt, const realType *u, realType *g, realType *work)
findpt_data_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

Define Documentation

#define DECLMINMAX (   type,
  prefix 
)

Definition at line 679 of file FindPtFuncs.h.

#define GLOBAL_INT   long

Definition at line 47 of file FindPtFuncs.h.

#define INTEGER   int

Definition at line 38 of file FindPtFuncs.h.

#define MAYBE_UNUSED

Definition at line 676 of file FindPtFuncs.h.

#define MAYBE_UNUSED

Definition at line 676 of file FindPtFuncs.h.

#define mbabs   fabs
#define mbceil   ceil

Definition at line 20 of file FindPtFuncs.h.

Referenced by iceil().

#define mbcos   cos

Definition at line 17 of file FindPtFuncs.h.

Referenced by gauss_nodes(), lob_bnd_base_setup(), and lobatto_nodes_aux().

#define mbfloor   floor

Definition at line 19 of file FindPtFuncs.h.

Referenced by ifloor().

#define mbsin   sin

Definition at line 18 of file FindPtFuncs.h.

#define mbsqrt   sqrt

Definition at line 15 of file FindPtFuncs.h.

Referenced by r2norm_1(), r2norm_2(), and r2norm_3().

#define MOAB_POLY_EPS   ( 128 * DBL_EPSILON )
#define MOAB_POLY_PI   3.1415926535897932384626433832795028841971693993751058209749445923

Definition at line 13 of file FindPtFuncs.h.

Referenced by gauss_nodes(), lob_bnd_base_setup(), and lobatto_nodes_aux().

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

Definition at line 636 of file FindPtFuncs.h.

Referenced by hash_build_2(), and hash_build_3().

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

Definition at line 637 of file FindPtFuncs.h.


Typedef Documentation

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

Definition at line 538 of file FindPtFuncs.h.

typedef double realType

Definition at line 56 of file FindPtFuncs.h.

typedef signed INTEGER sint

Definition at line 61 of file FindPtFuncs.h.

typedef signed GLOBAL_INT slong

Definition at line 71 of file FindPtFuncs.h.

typedef unsigned INTEGER uint

Definition at line 65 of file FindPtFuncs.h.

typedef unsigned GLOBAL_INT ulong

Definition at line 79 of file FindPtFuncs.h.


Function Documentation

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.

{
    if( guess && findpt_guess_2( p, x, *el, r, dist ) ) return 0;
    findpt_hash_2( p, x );
    if( p->sorted == p->end ) return -1;
    return findpt_pass_2( p, x, el, r, dist );
}
int findpt_3 ( findpt_data_3 p,
const realType  x[3],
int  guess,
uint el,
realType  r[3],
realType dist 
)

Definition at line 2259 of file findpt.c.

{
    if( guess && findpt_guess_3( p, x, *el, r, dist ) ) return 0;
    findpt_hash_3( p, x );
#if DIAGNOSTICS
    printf( "hashing leaves %d elements to consider\n", p->end - p->sorted );
#endif
    if( p->sorted == p->end ) return -1;
    return findpt_pass_3( p, x, el, r, dist );
}
const realType* findpt_allbnd_2 ( const findpt_data_2 p)

Definition at line 2107 of file findpt.c.

{
    return p->hash->bnd;
}
const realType* findpt_allbnd_3 ( const findpt_data_3 p)

Definition at line 2112 of file findpt.c.

{
    return p->hash->bnd;
}
double findpt_eval_2 ( findpt_data_2 p,
const realType u 
) [inline]

Definition at line 2283 of file findpt.c.

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

{
    return tensor_i2( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, u, p->od_work );
}
double findpt_eval_3 ( findpt_data_3 p,
const realType u 
) [inline]

Definition at line 2288 of file findpt.c.

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

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

Definition at line 2079 of file findpt.c.

{
    unsigned d;
    opt_free_2( p->od );
    free( p->od );
    hash_free_2( p->hash );
    free( p->hash );
    free( p->list );
    free( p->sorted );
    for( d = 0; d < 2; ++d )
        free( p->z[d] );
    free( p );
}

Definition at line 2093 of file findpt.c.

{
    unsigned d;
    opt_free_3( p->od );
    free( p->od );
    hash_free_3( p->hash );
    free( p->hash );
    free( p->list );
    free( p->sorted );
    for( d = 0; d < 3; ++d )
        free( p->z[d] );
    free( p );
}
findpt_data_2* findpt_setup_2 ( const 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.

{
    unsigned d;
    findpt_data_2* p = tmalloc( findpt_data_2, 1 );

    p->hash = tmalloc( hash_data_2, 1 );
    p->od   = tmalloc( opt_data_2, 1 );

    for( d = 0; d < 2; ++d )
        p->xw[d] = xw[d];
    p->nptel = n[0] * n[1];

    hash_build_2( p->hash, xw, n, nel, max_hash_size, bbox_tol );

    for( d = 0; d < 2; ++d )
    {
        p->z[d] = tmalloc( realType, n[d] );
        lobatto_nodes( p->z[d], n[d] );
        lagrange_setup( &p->ld[d], p->z[d], n[d] );
    }

    p->list   = tmalloc( findpt_listel, p->hash->max );
    p->sorted = tmalloc( findpt_listel*, p->hash->max );

    opt_alloc_2( p->od, p->ld );
    p->od_work = p->od->work;

    return p;
}
findpt_data_3* findpt_setup_3 ( const realType *const  xw[3],
const unsigned  n[3],
uint  nel,
uint  max_hash_size,
realType  bbox_tol 
)

Definition at line 2045 of file findpt.c.

{
    unsigned d;
    findpt_data_3* p = tmalloc( findpt_data_3, 1 );

    p->hash = tmalloc( hash_data_3, 1 );
    p->od   = tmalloc( opt_data_3, 1 );

    for( d = 0; d < 3; ++d )
        p->xw[d] = xw[d];
    p->nptel = n[0] * n[1] * n[2];

    hash_build_3( p->hash, xw, n, nel, max_hash_size, bbox_tol );

    for( d = 0; d < 3; ++d )
    {
        p->z[d] = tmalloc( realType, n[d] );
        lobatto_nodes( p->z[d], n[d] );
        lagrange_setup( &p->ld[d], p->z[d], n[d] );
    }

    p->list   = tmalloc( findpt_listel, p->hash->max );
    p->sorted = tmalloc( findpt_listel*, p->hash->max );

    opt_alloc_3( p->od, p->ld );
    p->od_work = p->od->work;

    return p;
}
void findpt_weights_2 ( findpt_data_2 p,
const realType  r[2] 
) [inline]

Definition at line 2270 of file findpt.c.

References lagrange_0(), and findpt_data_2::ld.

{
    lagrange_0( &p->ld[0], r[0] );
    lagrange_0( &p->ld[1], r[1] );
}
void findpt_weights_3 ( findpt_data_3 p,
const realType  r[3] 
) [inline]

Definition at line 2276 of file findpt.c.

References lagrange_0(), and findpt_data_3::ld.

{
    lagrange_0( &p->ld[0], r[0] );
    lagrange_0( &p->ld[1], r[1] );
    lagrange_0( &p->ld[2], r[2] );
}
void gauss_nodes ( realType z,
int  n 
)

Definition at line 155 of file poly.c.

{
    int i, j;
    for( i = 0; i <= n / 2 - 1; ++i )
    {
        realType ox, x = mbcos( ( 2 * n - 2 * i - 1 ) * ( MOAB_POLY_PI / 2 ) / n );
        do
        {
            ox = x;
            x -= legendre( n, x ) / legendre_d1( n, x );
        } while( mbabs( x - ox ) > -x * MOAB_POLY_EPS );
        z[i] = x - legendre( n, x ) / legendre_d1( n, x );
    }
    if( n & 1 ) z[n / 2] = 0;
    for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
        z[j] = -z[i];
}
void gauss_to_legendre ( const realType z,
const realType w,
int  n,
realType J 
)

Definition at line 239 of file poly.c.

{
    int i, j;
    legendre_matrix_t( z, n, J, n - 1 );
    for( j = 0; j < n; ++j )
    {
        realType ww = w[j];
        for( i = 0; i < n; ++i )
            *J++ *= ( 2 * i + 1 ) * ww / 2;
    }
}
void gauss_to_legendre_t ( const realType z,
const realType w,
int  n,
realType J 
)

Definition at line 255 of file poly.c.

{
    int i, j;
    legendre_matrix( z, n, J, n - 1 );
    for( i = 0; i < n; ++i )
    {
        realType ii = (realType)( 2 * i + 1 ) / 2;
        for( j = 0; j < n; ++j )
            *J++ *= ii * w[j];
    }
}
void gauss_weights ( const realType z,
realType w,
int  n 
)

Definition at line 199 of file poly.c.

{
    int i, j;
    for( i = 0; i <= ( n - 1 ) / 2; ++i )
    {
        realType d = ( n + 1 ) * legendre( n + 1, z[i] );
        w[i]       = 2 * ( 1 - z[i] * z[i] ) / ( d * d );
    }
    for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
        w[j] = w[i];
}
void lagrange_0 ( lagrange_data p,
realType  x 
)

Definition at line 414 of file poly.c.

{
    unsigned i, n = p->n;
    for( i = 0; i < n; ++i )
        p->d[i] = x - p->z[i];
    for( i = 0; i < n - 1; ++i )
        p->u0[i + 1] = p->d[i] * p->u0[i];
    for( i = n - 1; i; --i )
        p->v0[i - 1] = p->d[i] * p->v0[i];
    for( i = 0; i < n; ++i )
        p->J[i] = p->w[i] * p->u0[i] * p->v0[i];
}
void lagrange_1 ( lagrange_data p,
realType  x 
)

Definition at line 427 of file poly.c.

{
    unsigned i, n = p->n;
    for( i = 0; i < n; ++i )
        p->d[i] = x - p->z[i];
    for( i = 0; i < n - 1; ++i )
        p->u0[i + 1] = p->d[i] * p->u0[i], p->u1[i + 1] = p->d[i] * p->u1[i] + p->u0[i];
    for( i = n - 1; i; --i )
        p->v0[i - 1] = p->d[i] * p->v0[i], p->v1[i - 1] = p->d[i] * p->v1[i] + p->v0[i];
    for( i = 0; i < n; ++i )
        p->J[i] = p->w[i] * p->u0[i] * p->v0[i], p->D[i] = p->w[i] * ( p->u1[i] * p->v0[i] + p->u0[i] * p->v1[i] );
}
void lagrange_2 ( lagrange_data p,
realType  x 
)

Definition at line 440 of file poly.c.

{
    unsigned i, n = p->n;
    for( i = 0; i < n; ++i )
        p->d[i] = x - p->z[i];
    for( i = 0; i < n - 1; ++i )
        p->u0[i + 1] = p->d[i] * p->u0[i], p->u1[i + 1] = p->d[i] * p->u1[i] + p->u0[i],
                  p->u2[i + 1] = p->d[i] * p->u2[i] + 2 * p->u1[i];
    for( i = n - 1; i; --i )
        p->v0[i - 1] = p->d[i] * p->v0[i], p->v1[i - 1] = p->d[i] * p->v1[i] + p->v0[i],
                  p->v2[i - 1] = p->d[i] * p->v2[i] + 2 * p->v1[i];
    for( i = 0; i < n; ++i )
        p->J[i] = p->w[i] * p->u0[i] * p->v0[i], p->D[i] = p->w[i] * ( p->u1[i] * p->v0[i] + p->u0[i] * p->v1[i] ),
        p->D2[i] = p->w[i] * ( p->u2[i] * p->v0[i] + 2 * p->u1[i] * p->v1[i] + p->u0[i] * p->v2[i] );
}
void lagrange_2u ( lagrange_data p)

Definition at line 456 of file poly.c.

{
    unsigned i, n = p->n;
    for( i = 0; i < n - 1; ++i )
        p->u2[i + 1] = p->d[i] * p->u2[i] + 2 * p->u1[i];
    for( i = n - 1; i; --i )
        p->v2[i - 1] = p->d[i] * p->v2[i] + 2 * p->v1[i];
    for( i = 0; i < n; ++i )
        p->D2[i] = p->w[i] * ( p->u2[i] * p->v0[i] + 2 * p->u1[i] * p->v1[i] + p->u0[i] * p->v2[i] );
}

Definition at line 497 of file poly.c.

{
    free( p->w );
}
void lagrange_setup ( lagrange_data p,
const realType z,
unsigned  n 
)

Definition at line 467 of file poly.c.

{
    unsigned i, j;
    p->n = n, p->z = z;
    p->w = tmalloc( realType, 17 * n );
    p->d = p->w + n;
    p->J = p->d + n, p->D = p->J + n, p->D2 = p->D + n;
    p->u0 = p->D2 + n, p->v0 = p->u0 + n;
    p->u1 = p->v0 + n, p->v1 = p->u1 + n;
    p->u2 = p->v1 + n, p->v2 = p->u2 + n;
    p->J_z0 = p->v2 + n, p->D_z0 = p->J_z0 + n, p->D2_z0 = p->D_z0 + n;
    p->J_zn = p->D2_z0 + n, p->D_zn = p->J_zn + n, p->D2_zn = p->D_zn + n;
    for( i = 0; i < n; ++i )
    {
        realType ww = 1, zi = z[i];
        for( j = 0; j < i; ++j )
            ww *= zi - z[j];
        for( ++j; j < n; ++j )
            ww *= zi - z[j];
        p->w[i] = 1 / ww;
    }
    p->u0[0] = p->v0[n - 1] = 1;
    p->u1[0] = p->v1[n - 1] = 0;
    p->u2[0] = p->v2[n - 1] = 0;
    lagrange_2( p, z[0] );
    memcpy( p->J_z0, p->J, 3 * n * sizeof( realType ) );
    lagrange_2( p, z[n - 1] );
    memcpy( p->J_zn, p->J, 3 * n * sizeof( realType ) );
}
void lagrange_weights ( const realType z,
unsigned  n,
const realType x,
unsigned  m,
realType J,
realType work 
)

Definition at line 320 of file poly.c.

{
    unsigned i, j;
    realType *w = work, *d = w + n, *u = d + n, *v = u + n;
    for( i = 0; i < n; ++i )
    {
        realType ww = 1, zi = z[i];
        for( j = 0; j < i; ++j )
            ww *= zi - z[j];
        for( ++j; j < n; ++j )
            ww *= zi - z[j];
        w[i] = 1 / ww;
    }
    u[0] = v[n - 1] = 1;
    for( i = 0; i < m; ++i )
    {
        realType xi = x[i];
        for( j = 0; j < n; ++j )
            d[j] = xi - z[j];
        for( j = 0; j < n - 1; ++j )
            u[j + 1] = d[j] * u[j];
        for( j = n - 1; j; --j )
            v[j - 1] = d[j] * v[j];
        for( j = 0; j < n; ++j )
            *J++ = w[j] * u[j] * v[j];
    }
}
void lagrange_weights_deriv ( const realType z,
unsigned  n,
const realType x,
unsigned  m,
realType J,
realType D,
realType work 
)

Definition at line 354 of file poly.c.

{
    unsigned i, j;
    realType *w = work, *d = w + n, *u = d + n, *v = u + n, *up = v + n, *vp = up + n;
    for( i = 0; i < n; ++i )
    {
        realType ww = 1, zi = z[i];
        for( j = 0; j < i; ++j )
            ww *= zi - z[j];
        for( ++j; j < n; ++j )
            ww *= zi - z[j];
        w[i] = 1 / ww;
    }
    u[0] = v[n - 1] = 1;
    up[0] = vp[n - 1] = 0;
    for( i = 0; i < m; ++i )
    {
        realType xi = x[i];
        for( j = 0; j < n; ++j )
            d[j] = xi - z[j];
        for( j = 0; j < n - 1; ++j )
            u[j + 1] = d[j] * u[j], up[j + 1] = d[j] * up[j] + u[j];
        for( j = n - 1; j; --j )
            v[j - 1] = d[j] * v[j], vp[j - 1] = d[j] * vp[j] + v[j];
        for( j = 0; j < n; ++j )
            *J++ = w[j] * u[j] * v[j], *D++ = w[j] * ( up[j] * v[j] + u[j] * vp[j] );
    }
}
void legendre_matrix ( const realType x,
int  m,
realType P,
int  n 
)

Definition at line 30 of file poly.c.

{
    int i, j;
    realType *Pjm1 = P, *Pj = Pjm1 + m, *Pjp1 = Pj + m;
    for( i = 0; i < m; ++i )
        Pjm1[i] = 1;
    for( i = 0; i < m; ++i )
        Pj[i] = x[i];
    for( j = 1; j < n; ++j )
    {
        realType c = 1 / (realType)( j + 1 ), a = c * ( 2 * j + 1 ), b = c * j;
        for( i = 0; i < m; ++i )
            Pjp1[i] = a * x[i] * Pj[i] - b * Pjm1[i];
        Pjp1 += m, Pj += m, Pjm1 += m;
    }
}
void legendre_matrix_t ( const realType x,
int  m,
realType P,
int  n 
)

Definition at line 87 of file poly.c.

{
    int i;
    if( n & 1 )
        for( i = 0; i < m; ++i, P += n + 1 )
            legendre_row_odd( x[i], P, n );
    else
        for( i = 0; i < m; ++i, P += n + 1 )
            legendre_row_even( x[i], P, n );
}
void legendre_row ( realType  x,
realType P,
int  n 
)

Definition at line 75 of file poly.c.

{
    if( n & 1 )
        legendre_row_odd( x, P, n );
    else
        legendre_row_even( x, P, n );
}
void lobatto_nodes ( realType z,
int  n 
)

Definition at line 193 of file poly.c.

{
    z[0] = -1, z[n - 1] = 1;
    lobatto_nodes_aux( &z[1], n - 2 );
}
void lobatto_to_legendre ( const realType z,
const realType w,
int  n,
realType J 
)

Definition at line 275 of file poly.c.

{
    int i, j, m = ( n + 1 ) / 2;
    realType *p = J, *q;
    realType ww, sum;
    if( n & 1 )
        for( j = 0; j < m; ++j, p += n )
            legendre_row_odd( z[j], p, n - 2 );
    else
        for( j = 0; j < m; ++j, p += n )
            legendre_row_even( z[j], p, n - 2 );
    p = J;
    for( j = 0; j < m; ++j )
    {
        ww = w[j], sum = 0;
        for( i = 0; i < n - 1; ++i )
            *p *= ( 2 * i + 1 ) * ww / 2, sum += *p++;
        *p++ = -sum;
    }
    q = J + ( n / 2 - 1 ) * n;
    if( n & 1 )
        for( ; j < n; ++j, p += n, q -= n )
        {
            for( i = 0; i < n - 1; i += 2 )
                p[i] = q[i], p[i + 1] = -q[i + 1];
            p[i] = q[i];
        }
    else
        for( ; j < n; ++j, p += n, q -= n )
        {
            for( i = 0; i < n - 1; i += 2 )
                p[i] = q[i], p[i + 1] = -q[i + 1];
        }
}
void lobatto_weights ( const realType z,
realType w,
int  n 
)

Definition at line 211 of file poly.c.

{
    int i, j;
    for( i = 0; i <= ( n - 1 ) / 2; ++i )
    {
        realType d = legendre( n - 1, z[i] );
        w[i]       = 2 / ( ( n - 1 ) * n * d * d );
    }
    for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
        w[j] = w[i];
}
void opt_alloc_2 ( opt_data_2 p,
lagrange_data ld 
)

Definition at line 1662 of file findpt.c.

{
    const unsigned nr = ld[0].n, ns = ld[1].n, ne = umax_2( nr, ns ), nw = 2 * ns;
    p->size[0]   = 1;
    p->size[1]   = nr;
    p->size[2]   = nr * ns;
    p->ld        = ld;
    p->work      = tmalloc( realType, 4 * ne + nw );
    p->ed.x[0]   = p->work + nw;
    p->ed.x[1]   = p->ed.x[0] + ne;
    p->ed.fd1[0] = p->ed.x[1] + ne;
    p->ed.fd1[1] = p->ed.fd1[0] + ne;
}
void opt_alloc_3 ( opt_data_3 p,
lagrange_data ld 
)

Definition at line 1251 of file findpt.c.

{
    const unsigned nr = ld[0].n, ns = ld[1].n, nt = ld[2].n, nf = umax_3( nr * ns, nr * nt, ns * nt ),
                   ne = umax_3( nr, ns, nt ), nw = 2 * ns * nt + 3 * ns;
    p->size[0]   = 1;
    p->size[1]   = nr;
    p->size[2]   = nr * ns;
    p->size[3]   = p->size[2] * nt;
    p->ld        = ld;
    p->work      = tmalloc( realType, 6 * nf + 9 * ne + nw );
    p->fd.x[0]   = p->work + nw;
    p->fd.x[1]   = p->fd.x[0] + nf;
    p->fd.x[2]   = p->fd.x[1] + nf;
    p->fd.fdn[0] = p->fd.x[2] + nf;
    p->fd.fdn[1] = p->fd.fdn[0] + nf;
    p->fd.fdn[2] = p->fd.fdn[1] + nf;
    p->ed.x[0]   = p->fd.fdn[2] + nf;
    p->ed.x[1]   = p->ed.x[0] + ne;
    p->ed.x[2]   = p->ed.x[1] + ne;
    p->ed.fd1[0] = p->ed.x[2] + ne;
    p->ed.fd1[1] = p->ed.fd1[0] + ne;
    p->ed.fd1[2] = p->ed.fd1[1] + ne;
    p->ed.fd2[0] = p->ed.fd1[2] + ne;
    p->ed.fd2[1] = p->ed.fd2[0] + ne;
    p->ed.fd2[2] = p->ed.fd2[1] + ne;
}
double opt_findpt_2 ( opt_data_2 p,
const realType *const  elx[2],
const realType  xstar[2],
realType  r[2],
unsigned *  constr 
)

Definition at line 1818 of file findpt.c.

{
    realType dr[2], resid[2], steep[2];

    unsigned c = *constr, ac, d, cc[2], step = 0;

    p->elx[0] = elx[0], p->elx[1] = elx[1];

    p->ed.constraints = opt_no_constraints_2;
    p->pd.constraints = opt_no_constraints_2;

#if DIAGNOSTICS
    printf( "opt_findpt: xstar = %g, %g\n", xstar[0], xstar[1] );
#endif

    do
    {
        ++step;
        if( step == 50 ) /*fail("%s: opt_findpt_2 did not converge\n",__FILE__);*/
            return 1.e+30;
#if DIAGNOSTICS
        printf( "  iteration %u\n", step );
        printf( "    %d constraint(s) active\n", (int)opt_constr_num_2[c] );
#endif
        /* update face/edge/point data if necessary,
           and evaluate x(r) as well as the jacobian */
        switch( opt_constr_num_2[c] )
        {
            case 0:
                opt_area_set_intp_2( p, r );
                break;
            case 1:
                opt_edge_set_intp_2( p, r, c );
                break;
            case 2:
                opt_point_set_intp_2( p, c );
                break;
        }
#if DIAGNOSTICS
        printf( "    r = %g, %g\n", r[0], r[1] );
        printf( "    x = %g, %g\n", p->x[0], p->x[1] );
#endif
        /* compute residual */
        for( d = 0; d < 2; ++d )
            resid[d] = xstar[d] - p->x[d];
#if DIAGNOSTICS
        printf( "    resid = %g, %g\n", resid[0], resid[1] );
        printf( "    2-norm = %g\n", r2norm_2( resid[0], resid[1] ) );
#endif
        /* check constraints against steepest descent direction */
        ac = c;
        if( opt_constr_num_2[c] )
        {
            opt_constr_unpack_2( c, cc );
            mat_app_2c( steep, p->jac, resid ); /* steepest descent = J^T r */
#if DIAGNOSTICS
            printf( "    steepest descent = %g, %g\n", steep[0], steep[1] );
#endif
            for( d = 0; d < 2; ++d )
                if( ( cc[d] == 0 && steep[d] > 0 ) || ( cc[d] == 2 && steep[d] < 0 ) ) cc[d] = 1;
            ac = opt_constr_pack_2( cc );
        }
        /* update face/edge/point data if necessary */
        if( ac != c )
        {
            c = ac;
#if DIAGNOSTICS
            printf( "    relaxed to %d constraints\n", (int)opt_constr_num_2[c] );
#endif
            switch( opt_constr_num_2[c] )
            {
                case 1:
                    opt_edge_set_2( p, r, c );
                    break;
                case 2:
                    opt_point_set_2( p, c );
                    break;
            }
        }
        /* compute Newton step */
        switch( opt_constr_num_2[c] )
        {
            case 0:
                tinyla_solve_2( dr, p->jac, resid );
                break;
            case 1: {
                const unsigned de = p->ed.de, d1 = p->ed.d1;
                realType fac, H[2];
                const realType* J = p->jac + de;
                opt_edge_hess_2( p, H );
                fac    = J[0] * J[0] + J[2] * J[2] - ( resid[0] * H[0] + resid[1] * H[1] );
                dr[de] = steep[de] / fac;
                dr[d1] = 0;
            }
            break;
            case 2:
                dr[0] = dr[1] = 0;
                break;
        }
#if DIAGNOSTICS
        printf( "    dr = %g, %g\n", dr[0], dr[1] );
#endif
        /* project new iteration onto [-1,1]^2 */
        opt_constr_unpack_2( c, cc );
        for( d = 0; d < 2; ++d )
        {
            if( cc[d] != 1 ) continue;
            r[d] += dr[d];
            if( r[d] <= -1 )
                dr[d] -= r[d] + 1, r[d] = -1, cc[d] = 0;
            else if( r[d] >= 1 )
                dr[d] -= r[d] - 1, r[d] = 1, cc[d] = 2;
        }
        c = opt_constr_pack_2( cc );
    } while( r1norm_2( dr[0], dr[1] ) > 2 * MOAB_POLY_EPS );
    *constr = c;
    return r2norm_2( resid[0], resid[1] );
}
double opt_findpt_3 ( opt_data_3 p,
const realType *const  elx[3],
const realType  xstar[3],
realType  r[3],
unsigned *  constr 
)

Definition at line 1512 of file findpt.c.

{
    realType dr[3], resid[3], steep[3];

    unsigned c = *constr, ac, d, cc[3], step = 0;

    p->elx[0] = elx[0], p->elx[1] = elx[1], p->elx[2] = elx[2];

    p->fd.constraints = opt_no_constraints_3;
    p->ed.constraints = opt_no_constraints_3;
    p->pd.constraints = opt_no_constraints_3;

#if DIAGNOSTICS
    printf( "opt_findpt: xstar = %g, %g, %g\n", xstar[0], xstar[1], xstar[2] );
#endif

    do
    {
        ++step;
        if( step == 50 ) /*fail("%s: opt_findpt_3 did not converge\n",__FILE__);*/
            return 1.e+30;
#if DIAGNOSTICS
        printf( "  iteration %u\n", step );
        printf( "    %d constraint(s) active\n", (int)opt_constr_num_3[c] );
#endif
        /* update face/edge/point data if necessary,
           and evaluate x(r) as well as the jacobian */
        switch( opt_constr_num_3[c] )
        {
            case 0:
                opt_vol_set_intp_3( p, r );
                break;
            case 1:
                opt_face_set_intp_3( p, r, c );
                break;
            case 2:
                opt_edge_set_intp_3( p, r, c );
                break;
            case 3:
                opt_point_set_intp_3( p, c );
                break;
        }
#if DIAGNOSTICS
        printf( "    r = %g, %g, %g\n", r[0], r[1], r[2] );
        printf( "    x = %g, %g, %g\n", p->x[0], p->x[1], p->x[2] );
#endif
        /* compute residual */
        for( d = 0; d < 3; ++d )
            resid[d] = xstar[d] - p->x[d];
#if DIAGNOSTICS
        printf( "    resid = %g, %g, %g\n", resid[0], resid[1], resid[2] );
        printf( "    2-norm = %g\n", r2norm_3( resid[0], resid[1], resid[2] ) );
#endif
        /* check constraints against steepest descent direction */
        ac = c;
        if( opt_constr_num_3[c] )
        {
            opt_constr_unpack_3( c, cc );
            mat_app_3c( steep, p->jac, resid ); /* steepest descent = J^T r */
#if DIAGNOSTICS
            printf( "    steepest descent = %g, %g, %g\n", steep[0], steep[1], steep[2] );
#endif
            for( d = 0; d < 3; ++d )
                if( ( cc[d] == 0 && steep[d] > 0 ) || ( cc[d] == 2 && steep[d] < 0 ) ) cc[d] = 1;
            ac = opt_constr_pack_3( cc );
        }
        /* update face/edge/point data if necessary */
        if( ac != c )
        {
            c = ac;
#if DIAGNOSTICS
            printf( "    relaxed to %d constraints\n", (int)opt_constr_num_3[c] );
#endif
            switch( opt_constr_num_3[c] )
            {
                case 1:
                    opt_face_set_3( p, r, c );
                    break;
                case 2:
                    opt_edge_set_3( p, r, c );
                    break;
                case 3:
                    opt_point_set_3( p, c );
                    break;
            }
        }
        /* compute Newton step */
        switch( opt_constr_num_3[c] )
        {
            case 0:
                tinyla_solve_3( dr, p->jac, resid );
                break;
            case 1: {
                const unsigned dn = p->fd.dn, d1 = p->fd.d1, d2 = p->fd.d2;
                realType A[4], H[9];
                const realType* J = p->jac;
                opt_face_hess_3( p, H );
                A[0] = J[d1] * J[d1] + J[3 + d1] * J[3 + d1] + J[6 + d1] * J[6 + d1];
                A[1] = J[d2] * J[d2] + J[3 + d2] * J[3 + d2] + J[6 + d2] * J[6 + d2];
                A[2] = J[d1] * J[d2] + J[3 + d1] * J[3 + d2] + J[6 + d1] * J[6 + d2];
                A[0] -= resid[0] * H[0] + resid[1] * H[3] + resid[2] * H[6];
                A[1] -= resid[0] * H[1] + resid[1] * H[4] + resid[2] * H[7];
                A[2] -= resid[0] * H[2] + resid[1] * H[5] + resid[2] * H[8];
                tinyla_solve_sym_2( &dr[d1], &dr[d2], A, steep[d1], steep[d2] );
                dr[dn] = 0;
            }
            break;
            case 2: {
                const unsigned de = p->ed.de, d1 = p->ed.d1, d2 = p->ed.d2;
                realType fac, H[3];
                const realType* J = p->jac + de;
                opt_edge_hess_3( p, H );
                fac = J[0] * J[0] + J[3] * J[3] + J[6] * J[6] - ( resid[0] * H[0] + resid[1] * H[1] + resid[2] * H[2] );
                dr[de] = steep[de] / fac;
                dr[d1] = 0, dr[d2] = 0;
            }
            break;
            case 3:
                dr[0] = dr[1] = dr[2] = 0;
                break;
        }
#if DIAGNOSTICS
        printf( "    dr = %g, %g, %g\n", dr[0], dr[1], dr[2] );
#endif
        /* project new iteration onto [-1,1]^3 */
        opt_constr_unpack_3( c, cc );
        for( d = 0; d < 3; ++d )
        {
            if( cc[d] != 1 ) continue;
            r[d] += dr[d];
            if( r[d] <= -1 )
                dr[d] -= r[d] + 1, r[d] = -1, cc[d] = 0;
            else if( r[d] >= 1 )
                dr[d] -= r[d] - 1, r[d] = 1, cc[d] = 2;
        }
        c = opt_constr_pack_3( cc );
    } while( r1norm_3( dr[0], dr[1], dr[2] ) > 3 * MOAB_POLY_EPS );
    *constr = c;
#if 0
  printf("opt_findpt_3 converged in %u iterations\n", step);
#endif
    return r2norm_3( resid[0], resid[1], resid[2] );
}
void opt_free_2 ( opt_data_2 p)

Definition at line 1676 of file findpt.c.

{
    free( p->work );
}
void opt_free_3 ( opt_data_3 p)

Definition at line 1278 of file findpt.c.

{
    free( p->work );
}
void opt_vol_set_intp_3 ( opt_data_3 p,
const realType  r[3] 
)

Definition at line 1301 of file findpt.c.

{
    opt_vol_set_3( p, r );
    opt_vol_intp_3( p );
}
static realType r1norm_1 ( realType  a) [static]

Definition at line 723 of file FindPtFuncs.h.

References mbabs.

{
    return mbabs( a );
}
static realType r1norm_2 ( realType  a,
realType  b 
) [static]

Definition at line 728 of file FindPtFuncs.h.

References mbabs.

Referenced by findpt_hash_2(), and opt_findpt_2().

{
    return mbabs( a ) + mbabs( b );
}
static realType r1norm_3 ( realType  a,
realType  b,
realType  c 
) [static]

Definition at line 733 of file FindPtFuncs.h.

References mbabs.

Referenced by findpt_hash_3(), and opt_findpt_3().

{
    return mbabs( a ) + mbabs( b ) + mbabs( c );
}
static realType r2norm_1 ( realType  a) [static]

Definition at line 738 of file FindPtFuncs.h.

References mbsqrt.

{
    return mbsqrt( a * a );
}
static realType r2norm_2 ( realType  a,
realType  b 
) [static]

Definition at line 743 of file FindPtFuncs.h.

References mbsqrt.

Referenced by findpt_pass_2(), and opt_findpt_2().

{
    return mbsqrt( a * a + b * b );
}
static realType r2norm_3 ( realType  a,
realType  b,
realType  c 
) [static]

Definition at line 748 of file FindPtFuncs.h.

References mbsqrt.

Referenced by findpt_pass_3(), and opt_findpt_3().

{
    return mbsqrt( a * a + b * b + c * c );
}
static void * scalloc ( size_t  nmemb,
size_t  size,
const char *  file 
) [static]

Definition at line 620 of file FindPtFuncs.h.

References fail().

{
    void* res = calloc( nmemb, size );
    if( !res && nmemb ) fail( "%s: allocation of %d bytes failed\n", file, (int)size * nmemb );
    return res;
}
static void * smalloc ( size_t  size,
const char *  file 
) [static]

Definition at line 612 of file FindPtFuncs.h.

References fail().

{
    void* res = malloc( size );
    if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
    return res;
}
static void * srealloc ( void *  ptr,
size_t  size,
const char *  file 
) [static]

Definition at line 628 of file FindPtFuncs.h.

References fail().

{
    void* res = realloc( ptr, size );
    if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
    return res;
}
void tensor_c1 ( const realType R,
unsigned  mr,
unsigned  nr,
const realType u,
realType v 
)

Definition at line 155 of file tensor.c.

{
    mxv_c( v, mr, R, u, nr );
}
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.

{
    mxm_cc( R, mr, u, nr, W, ns );
    mxm_cr( W, mr, S, ns, v, ms );
}
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.

{
    unsigned n, mrns = mr * ns, mrms = mr * ms;
    realType* Zp = Z;
    mxm_cc( R, mr, u, nr, W, ns * nt );
    for( n = 0; n < nt; ++n, W += mrns, Zp += mrms )
        mxm_cr( W, mr, S, ns, Zp, ms );
    mxm_cr( Z, mrms, T, nt, v, mt );
}
realType tensor_i1 ( const realType Jr,
unsigned  nr,
const realType u 
)

Definition at line 255 of file tensor.c.

{
    return inner( Jr, u, nr );
}
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.

{
    mxv_r( work, ns, u, Jr, nr );
    return inner( Js, work, ns );
}
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.

{
    realType* work2 = work + nt;
    mxv_r( work2, ns * nt, u, Jr, nr );
    mxv_r( work, nt, work2, Js, ns );
    return inner( Jt, work, nt );
}
realType tensor_ig1 ( const realType Jr,
const realType Dr,
unsigned  nr,
const realType u,
realType g 
)

Definition at line 304 of file tensor.c.

{
    *g = inner( Dr, u, nr );
    return inner( Jr, u, nr );
}
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.

{
    realType *a = work, *ar = a + ns;
    mxv_r( a, ns, u, Jr, nr );
    mxv_r( ar, ns, u, Dr, nr );
    g[0] = inner( Js, ar, ns );
    g[1] = inner( Ds, a, ns );
    return inner( Js, a, ns );
}
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.

{
    unsigned nsnt = ns * nt;
    realType *a = work, *ar = a + nsnt, *b = ar + nsnt, *br = b + ns, *bs = br + ns;
    mxv_r( a, nsnt, u, Jr, nr );
    mxv_r( ar, nsnt, u, Dr, nr );
    mxv_r( b, nt, a, Js, ns );
    mxv_r( br, nt, ar, Js, ns );
    mxv_r( bs, nt, a, Ds, ns );
    g[0] = inner( Jt, br, nt );
    g[1] = inner( Jt, bs, nt );
    g[2] = inner( Dt, b, nt );
    return inner( Jt, b, nt );
}
void tensor_r1 ( const realType R,
unsigned  mr,
unsigned  nr,
const realType u,
realType v 
)

Definition at line 160 of file tensor.c.

{
    mxv_r( v, mr, R, u, nr );
}
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.

{
    mxm_rc( R, mr, u, nr, W, ns );
    mxm_cc( W, mr, S, ns, v, ms );
}
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.

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

Variable Documentation

const unsigned opt_no_constraints_2

Definition at line 10 of file findpt.c.

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

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines