MOAB: Mesh Oriented datABase  (version 5.2.1)
SpectralFuncs.hpp
Go to the documentation of this file.
00001 #ifndef SPECTRALFUNCS_HPP
00002 #define SPECTRALFUNCS_HPP
00003 
00004 #include "float.h"
00005 
00006 //======================================================
00007 // from types.h
00008 //======================================================
00009 
00010 /* integer type to use for everything */
00011 #if defined( USE_LONG )
00012 #define INTEGER long
00013 #elif defined( USE_LONG_LONG )
00014 #define INTEGER long long
00015 #elif defined( USE_SHORT )
00016 #define INTEGER short
00017 #else
00018 #define INTEGER int
00019 #endif
00020 
00021 /* when defined, use the given type for global indices instead of INTEGER */
00022 #if defined( USE_GLOBAL_LONG_LONG )
00023 #define GLOBAL_INT long long
00024 #elif defined( USE_GLOBAL_LONG )
00025 #define GLOBAL_INT long
00026 #else
00027 #define GLOBAL_INT long
00028 #endif
00029 
00030 /* floating point type to use for everything */
00031 #if defined( USE_FLOAT )
00032 typedef float real;
00033 #define floorr floorf
00034 #define ceilr  ceilf
00035 #define sqrtr  sqrtf
00036 #define fabsr  fabsf
00037 #define cosr   cosf
00038 #define sinr   sinf
00039 #define EPS    ( 128 * FLT_EPSILON )
00040 #define PI     3.1415926535897932384626433832795028841971693993751058209749445923F
00041 #elif defined( USE_LONG_DOUBLE )
00042 typedef long double real;
00043 #define floorr floorl
00044 #define ceilr  ceill
00045 #define sqrtr  sqrtl
00046 #define fabsr  fabsl
00047 #define cosr   cosl
00048 #define sinr   sinl
00049 #define EPS    ( 128 * LDBL_EPSILON )
00050 #define PI     3.1415926535897932384626433832795028841971693993751058209749445923L
00051 #else
00052 typedef double real;
00053 #define floorr floor
00054 #define ceilr  ceil
00055 #define sqrtr  sqrt
00056 #define fabsr  fabs
00057 #define cosr   cos
00058 #define sinr   sin
00059 #define EPS    ( 128 * DBL_EPSILON )
00060 #define PI     3.1415926535897932384626433832795028841971693993751058209749445923
00061 #endif
00062 
00063 /* apparently uint and ulong can be defined already in standard headers */
00064 #define uint  uint_
00065 #define ulong ulong_
00066 #define sint  sint_
00067 #define slong slong_
00068 
00069 typedef signed INTEGER sint;
00070 typedef unsigned INTEGER uint;
00071 #undef INTEGER
00072 
00073 #ifdef GLOBAL_INT
00074 typedef signed GLOBAL_INT slong;
00075 typedef unsigned GLOBAL_INT ulong;
00076 #else
00077 typedef sint slong;
00078 typedef uint ulong;
00079 #endif
00080 
00081 //======================================================
00082 // from poly.h
00083 //======================================================
00084 
00085 /*
00086   For brevity's sake, some names have been shortened
00087   Quadrature rules
00088     Gauss   -> Gauss-Legendre quadrature (open)
00089     Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends)
00090   Polynomial bases
00091     Legendre -> Legendre basis
00092     Gauss    -> Lagrangian basis using Gauss   quadrature nodes
00093     Lobatto  -> Lagrangian basis using Lobatto quadrature nodes
00094 */
00095 
00096 /*--------------------------------------------------------------------------
00097    Legendre Polynomial Matrix Computation
00098    (compute P_i(x_j) for i = 0, ..., n and a given set of x)
00099   --------------------------------------------------------------------------*/
00100 
00101 /* precondition: n >= 1
00102    inner index is x index (0 ... m-1);
00103    outer index is Legendre polynomial number (0 ... n)
00104  */
00105 void legendre_matrix( const real* x, int m, real* P, int n );
00106 
00107 /* precondition: n >= 1
00108    inner index is Legendre polynomial number (0 ... n)
00109    outer index is x index (0 ... m-1);
00110  */
00111 void legendre_matrix_t( const real* x, int m, real* P, int n );
00112 
00113 /* precondition: n >= 1
00114    compute P_i(x) with i = 0 ... n
00115  */
00116 void legendre_row( real x, real* P, int n );
00117 
00118 /*--------------------------------------------------------------------------
00119    Quadrature Nodes and Weights Calculation
00120 
00121    call the _nodes function before calling the _weights function
00122   --------------------------------------------------------------------------*/
00123 
00124 void gauss_nodes( real* z, int n );   /* n nodes (order = 2n-1) */
00125 void lobatto_nodes( real* z, int n ); /* n nodes (order = 2n-3) */
00126 
00127 void gauss_weights( const real* z, real* w, int n );
00128 void lobatto_weights( const real* z, real* w, int n );
00129 
00130 /*--------------------------------------------------------------------------
00131    Lagrangian to Legendre Change-of-basis Matrix
00132   --------------------------------------------------------------------------*/
00133 
00134 /* precondition: n >= 2
00135    given the Gauss quadrature rule (z,w,n), compute the square matrix J
00136    for transforming from the Gauss basis to the Legendre basis:
00137 
00138       u_legendre(i) = sum_j J(i,j) u_gauss(j)
00139 
00140    computes J   = .5 (2i+1) w  P (z )
00141              ij              j  i  j
00142 
00143    in column major format (inner index is i, the Legendre index)
00144  */
00145 void gauss_to_legendre( const real* z, const real* w, int n, real* J );
00146 
00147 /* precondition: n >= 2
00148    same as above, but
00149    in row major format (inner index is j, the Gauss index)
00150  */
00151 void gauss_to_legendre_t( const real* z, const real* w, int n, real* J );
00152 
00153 /* precondition: n >= 3
00154    given the Lobatto quadrature rule (z,w,n), compute the square matrix J
00155    for transforming from the Lobatto basis to the Legendre basis:
00156 
00157       u_legendre(i) = sum_j J(i,j) u_lobatto(j)
00158 
00159    in column major format (inner index is i, the Legendre index)
00160  */
00161 void lobatto_to_legendre( const real* z, const real* w, int n, real* J );
00162 
00163 /*--------------------------------------------------------------------------
00164    Lagrangian basis function evaluation
00165   --------------------------------------------------------------------------*/
00166 
00167 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
00168    evaluate all Lagrangian basis functions at all points x
00169 
00170    inner index of output J is the basis function index (row-major format)
00171    provide work array with space for 4*n doubles
00172  */
00173 void lagrange_weights( const real* z, unsigned n, const real* x, unsigned m, real* J, real* work );
00174 
00175 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
00176    evaluate all Lagrangian basis functions and their derivatives
00177 
00178    inner index of outputs J,D is the basis function index (row-major format)
00179    provide work array with space for 6*n doubles
00180  */
00181 void lagrange_weights_deriv( const real* z, unsigned n, const real* x, unsigned m, real* J, real* D, real* work );
00182 
00183 /*--------------------------------------------------------------------------
00184    Speedy Lagrangian Interpolation
00185 
00186    Usage:
00187 
00188      lagrange_data p;
00189      lagrange_setup(&p,z,n);    *  setup for nodes z[0 ... n-1] *
00190 
00191      the weights
00192        p->J [0 ... n-1]     interpolation weights
00193        p->D [0 ... n-1]     1st derivative weights
00194        p->D2[0 ... n-1]     2nd derivative weights
00195      are computed for a given x with:
00196        lagrange_0(p,x);  *  compute p->J *
00197        lagrange_1(p,x);  *  compute p->J, p->D *
00198        lagrange_2(p,x);  *  compute p->J, p->D, p->D2 *
00199        lagrange_2u(p);   *  compute p->D2 after call of lagrange_1(p,x); *
00200      These functions use the z array supplied to setup
00201        (that pointer should not be freed between calls)
00202      Weights for x=z[0] and x=z[n-1] are computed during setup; access as:
00203        p->J_z0, etc. and p->J_zn, etc.
00204 
00205      lagrange_free(&p);  *  deallocate memory allocated by setup *
00206   --------------------------------------------------------------------------*/
00207 
00208 typedef struct
00209 {
00210     unsigned n;                                /* number of Lagrange nodes            */
00211     const real* z;                             /* Lagrange nodes (user-supplied)      */
00212     real *J, *D, *D2;                          /* weights for 0th,1st,2nd derivatives */
00213     real *J_z0, *D_z0, *D2_z0;                 /* ditto at z[0]   (computed at setup) */
00214     real *J_zn, *D_zn, *D2_zn;                 /* ditto at z[n-1] (computed at setup) */
00215     real *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data           */
00216 } lagrange_data;
00217 
00218 void lagrange_setup( lagrange_data* p, const real* z, unsigned n );
00219 void lagrange_free( lagrange_data* p );
00220 
00221 void lagrange_0( lagrange_data* p, real x );
00222 void lagrange_1( lagrange_data* p, real x );
00223 void lagrange_2( lagrange_data* p, real x );
00224 void lagrange_2u( lagrange_data* p );
00225 
00226 //======================================================
00227 // from tensor.h
00228 //======================================================
00229 
00230 /*--------------------------------------------------------------------------
00231    1-,2-,3-d Tensor Application
00232 
00233    the 3d case:
00234    tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2)
00235      gives v = [ R (x) S (x) T ] u
00236      where R is mr x nr, S is ms x ns, T is mt x nt,
00237        each in row- or column-major format according to f := r | c
00238      u is nr x ns x nt in column-major format (inner index is r)
00239      v is mr x ms x mt in column-major format (inner index is r)
00240   --------------------------------------------------------------------------*/
00241 
00242 void tensor_c1( const real* R, unsigned mr, unsigned nr, const real* u, real* v );
00243 void tensor_r1( const real* R, unsigned mr, unsigned nr, const real* u, real* v );
00244 
00245 /* work holds mr*ns reals */
00246 void tensor_c2( const real* R, unsigned mr, unsigned nr, const real* S, unsigned ms, unsigned ns, const real* u,
00247                 real* v, real* work );
00248 void tensor_r2( const real* R, unsigned mr, unsigned nr, const real* S, unsigned ms, unsigned ns, const real* u,
00249                 real* v, real* work );
00250 
00251 /* work1 holds mr*ns*nt reals,
00252    work2 holds mr*ms*nt reals */
00253 void tensor_c3( const real* R, unsigned mr, unsigned nr, const real* S, unsigned ms, unsigned ns, const real* T,
00254                 unsigned mt, unsigned nt, const real* u, real* v, real* work1, real* work2 );
00255 void tensor_r3( const real* R, unsigned mr, unsigned nr, const real* S, unsigned ms, unsigned ns, const real* T,
00256                 unsigned mt, unsigned nt, const real* u, real* v, real* work1, real* work2 );
00257 
00258 /*--------------------------------------------------------------------------
00259    1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
00260 
00261    the 3d case:
00262    v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
00263    same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
00264      gives v = [ Jr (x) Js (x) Jt ] u
00265      where Jr, Js, Jt are row vectors (interpolation weights)
00266      u is nr x ns x nt in column-major format (inner index is r)
00267      v is a scalar
00268   --------------------------------------------------------------------------*/
00269 
00270 real tensor_i1( const real* Jr, unsigned nr, const real* u );
00271 
00272 /* work holds ns reals */
00273 real tensor_i2( const real* Jr, unsigned nr, const real* Js, unsigned ns, const real* u, real* work );
00274 
00275 /* work holds ns*nt + nt reals */
00276 real tensor_i3( const real* Jr, unsigned nr, const real* Js, unsigned ns, const real* Jt, unsigned nt, const real* u,
00277                 real* work );
00278 
00279 /*--------------------------------------------------------------------------
00280    1-,2-,3-d Tensor Application of Row Vectors
00281              for simultaneous Interpolation and Gradient computation
00282 
00283    the 3d case:
00284    v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
00285      gives v   = [ Jr (x) Js (x) Jt ] u
00286            g_0 = [ Dr (x) Js (x) Jt ] u
00287            g_1 = [ Jr (x) Ds (x) Jt ] u
00288            g_2 = [ Jr (x) Js (x) Dt ] u
00289      where Jr,Dr,Js,Ds,Jt,Dt are row vectors
00290        (interpolation & derivative weights)
00291      u is nr x ns x nt in column-major format (inner index is r)
00292      v is a scalar, g is an array of 3 reals
00293   --------------------------------------------------------------------------*/
00294 
00295 real tensor_ig1( const real* Jr, const real* Dr, unsigned nr, const real* u, real* g );
00296 
00297 /* work holds 2*ns reals */
00298 real tensor_ig2( const real* Jr, const real* Dr, unsigned nr, const real* Js, const real* Ds, unsigned ns,
00299                  const real* u, real* g, real* work );
00300 
00301 /* work holds 2*ns*nt + 3*ns reals */
00302 real tensor_ig3( const real* Jr, const real* Dr, unsigned nr, const real* Js, const real* Ds, unsigned ns,
00303                  const real* Jt, const real* Dt, unsigned nt, const real* u, real* g, real* work );
00304 
00305 //======================================================
00306 // from findpt.h
00307 //======================================================
00308 
00309 typedef struct
00310 {
00311     const real* xw[2];               /* geometry data */
00312     real* z[2];                      /* lobatto nodes */
00313     lagrange_data ld[2];             /* interpolation, derivative weights & data */
00314     unsigned nptel;                  /* nodes per element */
00315     struct findpt_hash_data_2* hash; /* geometric hashing data */
00316     struct findpt_listel *list, **sorted, **end;
00317     /* pre-allocated list of elements to
00318        check (found by hashing), and
00319        pre-allocated list of pointers into
00320        the first list for sorting */
00321     struct findpt_opt_data_2* od; /* data for the optimization algorithm */
00322     real* od_work;
00323 } findpt_data_2;
00324 
00325 typedef struct
00326 {
00327     const real* xw[3];               /* geometry data */
00328     real* z[3];                      /* lobatto nodes */
00329     lagrange_data ld[3];             /* interpolation, derivative weights & data */
00330     unsigned nptel;                  /* nodes per element */
00331     struct findpt_hash_data_3* hash; /* geometric hashing data */
00332     struct findpt_listel *list, **sorted, **end;
00333     /* pre-allocated list of elements to
00334        check (found by hashing), and
00335        pre-allocated list of pointers into
00336        the first list for sorting */
00337     struct findpt_opt_data_3* od; /* data for the optimization algorithm */
00338     real* od_work;
00339 } findpt_data_3;
00340 
00341 findpt_data_2* findpt_setup_2( const real* const xw[2], const unsigned n[2], uint nel, uint max_hash_size,
00342                                real bbox_tol );
00343 findpt_data_3* findpt_setup_3( const real* const xw[3], const unsigned n[3], uint nel, uint max_hash_size,
00344                                real bbox_tol );
00345 
00346 void findpt_free_2( findpt_data_2* p );
00347 void findpt_free_3( findpt_data_3* p );
00348 
00349 const real* findpt_allbnd_2( const findpt_data_2* p );
00350 const real* findpt_allbnd_3( const findpt_data_3* p );
00351 
00352 typedef int ( *findpt_func )( void*, const real*, int, uint*, real*, real* );
00353 int findpt_2( findpt_data_2* p, const real x[2], int guess, uint* el, real r[2], real* dist );
00354 int findpt_3( findpt_data_3* p, const real x[3], int guess, uint* el, real r[3], real* dist );
00355 
00356 inline void findpt_weights_2( findpt_data_2* p, const real r[2] )
00357 {
00358     lagrange_0( &p->ld[0], r[0] );
00359     lagrange_0( &p->ld[1], r[1] );
00360 }
00361 
00362 inline void findpt_weights_3( findpt_data_3* p, const real r[3] )
00363 {
00364     lagrange_0( &p->ld[0], r[0] );
00365     lagrange_0( &p->ld[1], r[1] );
00366     lagrange_0( &p->ld[2], r[2] );
00367 }
00368 
00369 inline double findpt_eval_2( findpt_data_2* p, const real* u )
00370 {
00371     return tensor_i2( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, u, p->od_work );
00372 }
00373 
00374 inline double findpt_eval_3( findpt_data_3* p, const real* u )
00375 {
00376     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 );
00377 }
00378 
00379 //======================================================
00380 // from extrafindpt.h
00381 //======================================================
00382 
00383 typedef struct
00384 {
00385     unsigned constraints;
00386     unsigned dn, d1, d2;
00387     real *x[3], *fdn[3];
00388 } opt_face_data_3;
00389 
00390 typedef struct
00391 {
00392     unsigned constraints;
00393     unsigned de, d1, d2;
00394     real *x[3], *fd1[3], *fd2[3];
00395 } opt_edge_data_3;
00396 
00397 typedef struct
00398 {
00399     unsigned constraints;
00400     real x[3], jac[9];
00401 } opt_point_data_3;
00402 
00403 typedef struct
00404 {
00405     lagrange_data* ld;
00406     unsigned size[4];
00407     const real* elx[3];
00408     opt_face_data_3 fd;
00409     opt_edge_data_3 ed;
00410     opt_point_data_3 pd;
00411     real* work;
00412     real x[3], jac[9];
00413 } opt_data_3;
00414 
00415 void opt_alloc_3( opt_data_3* p, lagrange_data* ld );
00416 void opt_free_3( opt_data_3* p );
00417 double opt_findpt_3( opt_data_3* p, const real* const elx[3], const real xstar[3], real r[3], unsigned* constr );
00418 void opt_vol_set_intp_3( opt_data_3* p, const real r[3] );
00419 
00420 const unsigned opt_no_constraints_2 = 3 + 1;
00421 const unsigned opt_no_constraints_3 = 9 + 3 + 1;
00422 
00423 /* for 2d spectralQuad */
00424 /*--------------------------------------------------------------------------
00425 
00426    2 - D
00427 
00428   --------------------------------------------------------------------------*/
00429 
00430 typedef struct
00431 {
00432     unsigned constraints;
00433     unsigned de, d1;
00434     real *x[2], *fd1[2];
00435 } opt_edge_data_2;
00436 
00437 typedef struct
00438 {
00439     unsigned constraints;
00440     real x[2], jac[4];
00441 } opt_point_data_2;
00442 
00443 typedef struct
00444 {
00445     lagrange_data* ld;
00446     unsigned size[3];
00447     const real* elx[2];
00448     opt_edge_data_2 ed;
00449     opt_point_data_2 pd;
00450     real* work;
00451     real x[2], jac[4];
00452 } opt_data_2;
00453 void opt_alloc_2( opt_data_2* p, lagrange_data* ld );
00454 void opt_free_2( opt_data_2* p );
00455 double opt_findpt_2( opt_data_2* p, const real* const elx[2], const real xstar[2], real r[2], unsigned* constr );
00456 
00457 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines