FindPtFuncs.h

Go to the documentation of this file.

#ifndef FINDPTFUNCS_H
#define FINDPTFUNCS_H

#include "float.h"
#include "math.h"

/*======================================================
/ from types.h
/======================================================
*/

/* integer type to use for everything */
#if   defined(USE_LONG)
#  define INTEGER long
#elif defined(USE_LONG_LONG)
#  define INTEGER long long
#elif defined(USE_SHORT)
#  define INTEGER short
#else
#  define INTEGER int
#endif

/* when defined, use the given type for global indices instead of INTEGER */
#if   defined(USE_GLOBAL_LONG_LONG)
#  define GLOBAL_INT long long
#elif defined(USE_GLOBAL_LONG)
#  define GLOBAL_INT long
#else
#  define GLOBAL_INT long
#endif

/* floating point type to use for everything */
#if   defined(USE_FLOAT)
   typedef float real;
#  define floorr floorf
#  define ceilr  ceilf
#  define sqrtr  sqrtf
#  define fabsr  fabsf
#  define cosr   cosf
#  define sinr   sinf
#  define EPS   (128*FLT_EPSILON)
#  define PI 3.1415926535897932384626433832795028841971693993751058209749445923F
#elif defined(USE_LONG_DOUBLE)
   typedef long double real;
#  define floorr floorl
#  define ceilr  ceill
#  define sqrtr  sqrtl
#  define fabsr  fabsl
#  define cosr   cosl
#  define sinr   sinl
#  define EPS   (128*LDBL_EPSILON)
#  define PI 3.1415926535897932384626433832795028841971693993751058209749445923L
#else
   typedef double real;
#  define floorr floor
#  define ceilr  ceil
#  define sqrtr  sqrt
#  define fabsr  fabs
#  define cosr   cos
#  define sinr   sin
#  define EPS   (128*DBL_EPSILON)
#  define PI 3.1415926535897932384626433832795028841971693993751058209749445923
#endif

/* apparently uint and ulong can be defined already in standard headers */
#define uint uint_
#define ulong ulong_
#define sint sint_
#define slong slong_

typedef   signed INTEGER sint;
typedef unsigned INTEGER uint;
#undef INTEGER

#ifdef GLOBAL_INT
  typedef   signed GLOBAL_INT slong;
  typedef unsigned GLOBAL_INT ulong;
#else
  typedef sint slong;
  typedef uint ulong;
#endif

/*======================================================
/ from poly.h
/======================================================
*/

/* 
  For brevity's sake, some names have been shortened
  Quadrature rules
    Gauss   -> Gauss-Legendre quadrature (open)
    Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends)
  Polynomial bases
    Legendre -> Legendre basis
    Gauss    -> Lagrangian basis using Gauss   quadrature nodes
    Lobatto  -> Lagrangian basis using Lobatto quadrature nodes
*/

/*--------------------------------------------------------------------------
   Legendre Polynomial Matrix Computation
   (compute P_i(x_j) for i = 0, ..., n and a given set of x)
  --------------------------------------------------------------------------*/

/* precondition: n >= 1
   inner index is x index (0 ... m-1);
   outer index is Legendre polynomial number (0 ... n)
 */
void legendre_matrix(const real *x, int m, real *P, int n);

/* precondition: n >= 1
   inner index is Legendre polynomial number (0 ... n)
   outer index is x index (0 ... m-1);
 */
void legendre_matrix_t(const real *x, int m, real *P, int n);

/* precondition: n >= 1
   compute P_i(x) with i = 0 ... n
 */
void legendre_row(real x, real *P, int n);


/*--------------------------------------------------------------------------
   Quadrature Nodes and Weights Calculation
   
   call the _nodes function before calling the _weights function
  --------------------------------------------------------------------------*/

void gauss_nodes(real *z, int n);   /* n nodes (order = 2n-1) */
void lobatto_nodes(real *z, int n); /* n nodes (order = 2n-3) */

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

/*--------------------------------------------------------------------------
   Lagrangian to Legendre Change-of-basis Matrix
  --------------------------------------------------------------------------*/

/* precondition: n >= 2
   given the Gauss quadrature rule (z,w,n), compute the square matrix J
   for transforming from the Gauss basis to the Legendre basis:
   
      u_legendre(i) = sum_j J(i,j) u_gauss(j)

   computes J   = .5 (2i+1) w  P (z )
             ij              j  i  j
             
   in column major format (inner index is i, the Legendre index)
 */
void gauss_to_legendre(const real *z, const real *w, int n, real *J);

/* precondition: n >= 2
   same as above, but
   in row major format (inner index is j, the Gauss index)
 */
void gauss_to_legendre_t(const real *z, const real *w, int n, real *J);

/* precondition: n >= 3
   given the Lobatto quadrature rule (z,w,n), compute the square matrix J
   for transforming from the Lobatto basis to the Legendre basis:
   
      u_legendre(i) = sum_j J(i,j) u_lobatto(j)

   in column major format (inner index is i, the Legendre index)
 */
void lobatto_to_legendre(const real *z, const real *w, int n, real *J);

/*--------------------------------------------------------------------------
   Lagrangian basis function evaluation
  --------------------------------------------------------------------------*/

/* given the Lagrangian nodes (z,n) and evaluation points (x,m)
   evaluate all Lagrangian basis functions at all points x
   
   inner index of output J is the basis function index (row-major format)
   provide work array with space for 4*n doubles
 */
void lagrange_weights(const real *z, unsigned n,
                      const real *x, unsigned m,
                      real *J, real *work);

/* given the Lagrangian nodes (z,n) and evaluation points (x,m)
   evaluate all Lagrangian basis functions and their derivatives
   
   inner index of outputs J,D is the basis function index (row-major format)
   provide work array with space for 6*n doubles
 */
void lagrange_weights_deriv(const real *z, unsigned n,
                            const real *x, unsigned m,
                            real *J, real *D, real *work);

/*--------------------------------------------------------------------------
   Speedy Lagrangian Interpolation
   
   Usage:
   
     lagrange_data p;
     lagrange_setup(&p,z,n);    *  setup for nodes z[0 ... n-1] *
     
     the weights
       p->J [0 ... n-1]     interpolation weights
       p->D [0 ... n-1]     1st derivative weights
       p->D2[0 ... n-1]     2nd derivative weights
     are computed for a given x with:
       lagrange_0(p,x);  *  compute p->J *
       lagrange_1(p,x);  *  compute p->J, p->D *
       lagrange_2(p,x);  *  compute p->J, p->D, p->D2 *
       lagrange_2u(p);   *  compute p->D2 after call of lagrange_1(p,x); *
     These functions use the z array supplied to setup
       (that pointer should not be freed between calls)
     Weights for x=z[0] and x=z[n-1] are computed during setup; access as:
       p->J_z0, etc. and p->J_zn, etc.

     lagrange_free(&p);  *  deallocate memory allocated by setup *
  --------------------------------------------------------------------------*/

typedef struct {
  unsigned n;                /* number of Lagrange nodes            */
  const real *z;             /* Lagrange nodes (user-supplied)      */
  real *J, *D, *D2;          /* weights for 0th,1st,2nd derivatives */
  real *J_z0, *D_z0, *D2_z0; /* ditto at z[0]   (computed at setup) */
  real *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */
  real *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data           */
} lagrange_data;

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

void lagrange_0(lagrange_data *p, real x) ;
void lagrange_1(lagrange_data *p, real x) ;
void lagrange_2(lagrange_data *p, real x) ;
void lagrange_2u(lagrange_data *p) ;

/*======================================================
/ from tensor.h
/======================================================
*/

/*--------------------------------------------------------------------------
   1-,2-,3-d Tensor Application
   
   the 3d case:
   tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2)
     gives v = [ R (x) S (x) T ] u
     where R is mr x nr, S is ms x ns, T is mt x nt,
       each in row- or column-major format according to f := r | c
     u is nr x ns x nt in column-major format (inner index is r)
     v is mr x ms x mt in column-major format (inner index is r)
  --------------------------------------------------------------------------*/

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

/* work holds mr*ns reals */
void tensor_c2(const real *R, unsigned mr, unsigned nr,
               const real *S, unsigned ms, unsigned ns,
               const real *u, real *v, real *work);
void tensor_r2(const real *R, unsigned mr, unsigned nr,
               const real *S, unsigned ms, unsigned ns,
               const real *u, real *v, real *work);

/* work1 holds mr*ns*nt reals,
   work2 holds mr*ms*nt reals */
void tensor_c3(const real *R, unsigned mr, unsigned nr,
               const real *S, unsigned ms, unsigned ns,
               const real *T, unsigned mt, unsigned nt,
               const real *u, real *v, real *work1, real *work2);
void tensor_r3(const real *R, unsigned mr, unsigned nr,
               const real *S, unsigned ms, unsigned ns,
               const real *T, unsigned mt, unsigned nt,
               const real *u, real *v, real *work1, real *work2);

/*--------------------------------------------------------------------------
   1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
   
   the 3d case:
   v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
   same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
     gives v = [ Jr (x) Js (x) Jt ] u
     where Jr, Js, Jt are row vectors (interpolation weights)
     u is nr x ns x nt in column-major format (inner index is r)
     v is a scalar
  --------------------------------------------------------------------------*/

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

/* work holds ns reals */
real tensor_i2(const real *Jr, unsigned nr,
               const real *Js, unsigned ns,
               const real *u, real *work);

/* work holds ns*nt + nt reals */
real tensor_i3(const real *Jr, unsigned nr,
               const real *Js, unsigned ns,
               const real *Jt, unsigned nt,
               const real *u, real *work);

/*--------------------------------------------------------------------------
   1-,2-,3-d Tensor Application of Row Vectors
             for simultaneous Interpolation and Gradient computation
   
   the 3d case:
   v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
     gives v   = [ Jr (x) Js (x) Jt ] u
           g_0 = [ Dr (x) Js (x) Jt ] u
           g_1 = [ Jr (x) Ds (x) Jt ] u
           g_2 = [ Jr (x) Js (x) Dt ] u
     where Jr,Dr,Js,Ds,Jt,Dt are row vectors
       (interpolation & derivative weights)
     u is nr x ns x nt in column-major format (inner index is r)
     v is a scalar, g is an array of 3 reals
  --------------------------------------------------------------------------*/

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

/* work holds 2*ns reals */
real tensor_ig2(const real *Jr, const real *Dr, unsigned nr,
                const real *Js, const real *Ds, unsigned ns,
                const real *u, real *g, real *work);

/* work holds 2*ns*nt + 3*ns reals */
real tensor_ig3(const real *Jr, const real *Dr, unsigned nr,
                const real *Js, const real *Ds, unsigned ns,
                const real *Jt, const real *Dt, unsigned nt,
                const real *u, real *g, real *work);

/*======================================================
/ from findpt.h
/======================================================
*/

typedef struct {
  real x[2], A[4], axis_bnd[4];
} obbox_2;

typedef struct {
  real x[3], A[9], axis_bnd[6];
} obbox_3;

typedef struct {
  unsigned hash_n;
  real bnd[4]; /* bounds for all elements */
  real fac[2]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
  obbox_2 *obb; /* obb[nel] -- bounding box info for each element */
  uint *offset; /* hash table -- for cell i,j:
                         uint index = j*hash_n+i,
                                  b = offset[index  ],
                                  e = offset[index+1];
                         elements in cell are
                           offset[b], offset[b+1], ..., offset[e-1] */
  unsigned max; /* maximum # of elements in any cell */
} hash_data_2;

typedef struct {
  unsigned hash_n;
  real bnd[6]; /* bounds for all elements */
  real fac[3]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
  obbox_3 *obb; /* obb[nel] -- bounding box info for each element */
  uint *offset; /* hash table -- for cell i,j,k:
                         uint index = (k*hash_n+j)*hash_n+i,
                                  b = offset[index  ],
                                  e = offset[index+1];
                         elements in cell are
                           offset[b], offset[b+1], ..., offset[e-1] */
  unsigned max; /* maximum # of elements in any cell */
} hash_data_3;

typedef struct {
  uint el;
  real r[3];
  real dist;
} findpt_listel;

typedef struct {
  unsigned constraints;
  unsigned de, d1;
  real *x[2], *fd1[2];
} opt_edge_data_2;

typedef struct {
  unsigned constraints;
  real x[2], jac[4];
} opt_point_data_2;

typedef struct {
  lagrange_data *ld;
  unsigned size[3];
  const real *elx[2];
  opt_edge_data_2 ed;
  opt_point_data_2 pd;
  real *work;
  real x[2], jac[4];
} opt_data_2;

typedef struct {
  const real *xw[2];   /* geometry data */
  real *z[2];          /* lobatto nodes */
  lagrange_data ld[2]; /* interpolation, derivative weights & data */
  unsigned nptel;      /* nodes per element */
  hash_data_2 *hash;   /* geometric hashing data */
  findpt_listel *list, **sorted, **end;
                                        /* pre-allocated list of elements to
                                           check (found by hashing), and
                                           pre-allocated list of pointers into
                                           the first list for sorting */
  opt_data_2 *od; /* data for the optimization algorithm */
  real *od_work;
} findpt_data_2;

typedef struct {
  unsigned constraints;
  unsigned dn, d1, d2;
  real *x[3], *fdn[3];
} opt_face_data_3;

typedef struct {
  unsigned constraints;
  unsigned de, d1, d2;
  real *x[3], *fd1[3], *fd2[3];
} opt_edge_data_3;

typedef struct {
  unsigned constraints;
  real x[3], jac[9];
} opt_point_data_3;

typedef struct {
  lagrange_data *ld;
  unsigned size[4];
  const real *elx[3];
  opt_face_data_3 fd;
  opt_edge_data_3 ed;
  opt_point_data_3 pd;
  real *work;
  real x[3], jac[9];
} opt_data_3;

typedef struct {
  const real *xw[3];   /* geometry data */
  real *z[3];          /* lobatto nodes */
  lagrange_data ld[3]; /* interpolation, derivative weights & data */
  unsigned nptel;      /* nodes per element */
  hash_data_3 *hash;   /* geometric hashing data */
  findpt_listel *list, **sorted, **end;
                                        /* pre-allocated list of elements to
                                           check (found by hashing), and
                                           pre-allocated list of pointers into
                                           the first list for sorting */
  opt_data_3 *od; /* data for the optimization algorithm */
  real *od_work;
} findpt_data_3;

findpt_data_2 *findpt_setup_2(
          const real *const xw[2], const unsigned n[2], uint nel,
          uint max_hash_size, real bbox_tol);
findpt_data_3 *findpt_setup_3(
          const real *const xw[3], const unsigned n[3], uint nel,
          uint max_hash_size, real bbox_tol);

void findpt_free_2(findpt_data_2 *p);
void findpt_free_3(findpt_data_3 *p);

const real *findpt_allbnd_2(const findpt_data_2 *p);
const real *findpt_allbnd_3(const findpt_data_3 *p);

typedef int (*findpt_func)(void *, const real *, int, uint *, real *, real *);
int findpt_2(findpt_data_2 *p, const real x[2], int guess,
             uint *el, real r[2], real *dist);
int findpt_3(findpt_data_3 *p, const real x[3], int guess,
             uint *el, real r[3], real *dist);

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

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

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

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

/*======================================================
/ from extrafindpt.h
/======================================================
*/

void opt_alloc_3(opt_data_3 *p, lagrange_data *ld);
void opt_free_3(opt_data_3 *p);
double opt_findpt_3(opt_data_3 *p, const real *const elx[3],
                           const real xstar[3], real r[3], unsigned *constr);
void opt_vol_set_intp_3(opt_data_3 *p, const real r[3]);

/* for 2d spectralQuad */
/*--------------------------------------------------------------------------

   2 - D

  --------------------------------------------------------------------------*/

void opt_alloc_2(opt_data_2 *p, lagrange_data *ld);
void opt_free_2(opt_data_2 *p);
double opt_findpt_2(opt_data_2 *p, const real *const elx[2],
                           const real xstar[2], real r[2], unsigned *constr);

extern const unsigned opt_no_constraints_2;
extern const unsigned opt_no_constraints_3;

/*======================================================
/ from errmem.h
/======================================================
*/

/* requires:
     <stdlib.h> for malloc, calloc, realloc, free
*/

/*--------------------------------------------------------------------------
   Error Reporting
   Memory Allocation Wrappers to Catch Out-of-memory
  --------------------------------------------------------------------------*/

/* #include "malloc.h" */
#include <stdlib.h>

#ifdef __GNUC__
void fail(const char *fmt, ...) __attribute__ ((noreturn));
#define MAYBE_UNUSED __attribute__ ((unused))
#else
void fail(const char *fmt, ...);
#define MAYBE_UNUSED
#endif

#if 0
{}
#endif

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

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

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

#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__) )
/*
typedef struct { size_t size; void *ptr; } buffer;

static void buffer_init_(buffer *b, size_t size, const char *file) MAYBE_UNUSED;
static void buffer_init_(buffer *b, size_t size, const char *file)
{
  b->size=size, b->ptr=smalloc(size,file);
}
static void buffer_reserve_(buffer *b, size_t min, const char *file)
  MAYBE_UNUSED;
static void buffer_reserve_(buffer *b, size_t min, const char *file)
{
  size_t size = b->size;
  if(size<min) {
    size+=size/2+1;
    if(size<min) size=min;
    b->ptr=srealloc(b->ptr,size,file);
  }
}
static void buffer_free(buffer *b) MAYBE_UNUSED;
static void buffer_free(buffer *b) { free(b->ptr); }

#define buffer_init(b,size) buffer_init_(b,size,__FILE__)
#define buffer_reserve(b,min) buffer_reserve_(b,min,__FILE__)
*/


/*======================================================
/ from minmax.h
/======================================================
*/

/*--------------------------------------------------------------------------
   Min, Max, Norm
  --------------------------------------------------------------------------*/

#ifdef __GNUC__
#define MAYBE_UNUSED __attribute__ ((unused))
#else
#define MAYBE_UNUSED
#endif

#define DECLMINMAX(type, prefix) \
  static type prefix##min_2(type a, type b) MAYBE_UNUSED;               \
  static type prefix##min_2(type a, type b) {                           \
    return b<a?b:a;                                                     \
  }                                                                     \
  static type prefix##max_2(type a, type b) MAYBE_UNUSED;               \
  static type prefix##max_2(type a, type b) {                           \
    return a>b?a:b;                                                     \
  }                                                                     \
  static void prefix##minmax_2(type *min, type *max, type a,            \
                               type b) MAYBE_UNUSED;                    \
  static void prefix##minmax_2(type *min, type *max, type a, type b) {  \
    if(b<a) *min=b, *max=a;                                             \
    else *min=a, *max=b;                                                \
  }                                                                     \
  static type prefix##min_3(type a, type b, type c) MAYBE_UNUSED;       \
  static type prefix##min_3(type a, type b, type c) {                   \
    return b<a?(c<b?c:b):(c<a?c:a);                                     \
  }                                                                     \
  static type prefix##max_3(type a, type b, type c) MAYBE_UNUSED;       \
  static type prefix##max_3(type a, type b, type c) {                   \
    return a>b?(a>c?a:c):(b>c?b:c);                                     \
  }                                                                     \
  static void prefix##minmax_3(type *min, type *max, type a, type b,    \
                               type c) MAYBE_UNUSED;                    \
  static void prefix##minmax_3(type *min, type *max, type a, type b,    \
                               type c) {                                \
    if(b<a) *min=prefix##min_2(b,c), *max=prefix##max_2(a,c);           \
    else    *min=prefix##min_2(a,c), *max=prefix##max_2(b,c);           \
  }

DECLMINMAX(int, i)
DECLMINMAX(unsigned, u)
DECLMINMAX(real, r)
#undef DECLMINMAX

static real r1norm_1(real a) MAYBE_UNUSED;
static real r1norm_1(real a) {
  return fabsr(a);
}
static real r1norm_2(real a, real b) MAYBE_UNUSED;
static real r1norm_2(real a, real b) {
  return fabsr(a)+fabsr(b);
}
static real r1norm_3(real a, real b, real c) MAYBE_UNUSED;
static real r1norm_3(real a, real b, real c) {
  return fabsr(a)+fabsr(b)+fabsr(c);
}
static real r2norm_1(real a) MAYBE_UNUSED;
static real r2norm_1(real a) {
  return sqrtr(a*a);
}
static real r2norm_2(real a, real b) MAYBE_UNUSED;
static real r2norm_2(real a, real b) {
  return sqrtr(a*a+b*b);
}
static real r2norm_3(real a, real b, real c) MAYBE_UNUSED;
static real r2norm_3(real a, real b, real c) {
  return sqrtr(a*a+b*b+c*c);
}

#endif