Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
SpectralFuncs.hpp
Go to the documentation of this file.
00001 #ifndef SPECTRALFUNCS_HPP
00002 #define SPECTRALFUNCS_HPP
00003 
00004 #include <cfloat>
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,
00247                 unsigned mr,
00248                 unsigned nr,
00249                 const real* S,
00250                 unsigned ms,
00251                 unsigned ns,
00252                 const real* u,
00253                 real* v,
00254                 real* work );
00255 void tensor_r2( const real* R,
00256                 unsigned mr,
00257                 unsigned nr,
00258                 const real* S,
00259                 unsigned ms,
00260                 unsigned ns,
00261                 const real* u,
00262                 real* v,
00263                 real* work );
00264 
00265 /* work1 holds mr*ns*nt reals,
00266    work2 holds mr*ms*nt reals */
00267 void tensor_c3( const real* R,
00268                 unsigned mr,
00269                 unsigned nr,
00270                 const real* S,
00271                 unsigned ms,
00272                 unsigned ns,
00273                 const real* T,
00274                 unsigned mt,
00275                 unsigned nt,
00276                 const real* u,
00277                 real* v,
00278                 real* work1,
00279                 real* work2 );
00280 void tensor_r3( const real* R,
00281                 unsigned mr,
00282                 unsigned nr,
00283                 const real* S,
00284                 unsigned ms,
00285                 unsigned ns,
00286                 const real* T,
00287                 unsigned mt,
00288                 unsigned nt,
00289                 const real* u,
00290                 real* v,
00291                 real* work1,
00292                 real* work2 );
00293 
00294 /*--------------------------------------------------------------------------
00295    1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
00296 
00297    the 3d case:
00298    v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
00299    same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
00300      gives v = [ Jr (x) Js (x) Jt ] u
00301      where Jr, Js, Jt are row vectors (interpolation weights)
00302      u is nr x ns x nt in column-major format (inner index is r)
00303      v is a scalar
00304   --------------------------------------------------------------------------*/
00305 
00306 real tensor_i1( const real* Jr, unsigned nr, const real* u );
00307 
00308 /* work holds ns reals */
00309 real tensor_i2( const real* Jr, unsigned nr, const real* Js, unsigned ns, const real* u, real* work );
00310 
00311 /* work holds ns*nt + nt reals */
00312 real tensor_i3( const real* Jr,
00313                 unsigned nr,
00314                 const real* Js,
00315                 unsigned ns,
00316                 const real* Jt,
00317                 unsigned nt,
00318                 const real* u,
00319                 real* work );
00320 
00321 /*--------------------------------------------------------------------------
00322    1-,2-,3-d Tensor Application of Row Vectors
00323              for simultaneous Interpolation and Gradient computation
00324 
00325    the 3d case:
00326    v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
00327      gives v   = [ Jr (x) Js (x) Jt ] u
00328            g_0 = [ Dr (x) Js (x) Jt ] u
00329            g_1 = [ Jr (x) Ds (x) Jt ] u
00330            g_2 = [ Jr (x) Js (x) Dt ] u
00331      where Jr,Dr,Js,Ds,Jt,Dt are row vectors
00332        (interpolation & derivative weights)
00333      u is nr x ns x nt in column-major format (inner index is r)
00334      v is a scalar, g is an array of 3 reals
00335   --------------------------------------------------------------------------*/
00336 
00337 real tensor_ig1( const real* Jr, const real* Dr, unsigned nr, const real* u, real* g );
00338 
00339 /* work holds 2*ns reals */
00340 real tensor_ig2( const real* Jr,
00341                  const real* Dr,
00342                  unsigned nr,
00343                  const real* Js,
00344                  const real* Ds,
00345                  unsigned ns,
00346                  const real* u,
00347                  real* g,
00348                  real* work );
00349 
00350 /* work holds 2*ns*nt + 3*ns reals */
00351 real tensor_ig3( const real* Jr,
00352                  const real* Dr,
00353                  unsigned nr,
00354                  const real* Js,
00355                  const real* Ds,
00356                  unsigned ns,
00357                  const real* Jt,
00358                  const real* Dt,
00359                  unsigned nt,
00360                  const real* u,
00361                  real* g,
00362                  real* work );
00363 
00364 //======================================================
00365 // from findpt.h
00366 //======================================================
00367 
00368 typedef struct
00369 {
00370     const real* xw[2];               /* geometry data */
00371     real* z[2];                      /* lobatto nodes */
00372     lagrange_data ld[2];             /* interpolation, derivative weights & data */
00373     unsigned nptel;                  /* nodes per element */
00374     struct findpt_hash_data_2* hash; /* geometric hashing data */
00375     struct findpt_listel *list, **sorted, **end;
00376     /* pre-allocated list of elements to
00377        check (found by hashing), and
00378        pre-allocated list of pointers into
00379        the first list for sorting */
00380     struct findpt_opt_data_2* od; /* data for the optimization algorithm */
00381     real* od_work;
00382 } findpt_data_2;
00383 
00384 typedef struct
00385 {
00386     const real* xw[3];               /* geometry data */
00387     real* z[3];                      /* lobatto nodes */
00388     lagrange_data ld[3];             /* interpolation, derivative weights & data */
00389     unsigned nptel;                  /* nodes per element */
00390     struct findpt_hash_data_3* hash; /* geometric hashing data */
00391     struct findpt_listel *list, **sorted, **end;
00392     /* pre-allocated list of elements to
00393        check (found by hashing), and
00394        pre-allocated list of pointers into
00395        the first list for sorting */
00396     struct findpt_opt_data_3* od; /* data for the optimization algorithm */
00397     real* od_work;
00398 } findpt_data_3;
00399 
00400 findpt_data_2* findpt_setup_2( const real* const xw[2],
00401                                const unsigned n[2],
00402                                uint nel,
00403                                uint max_hash_size,
00404                                real bbox_tol );
00405 findpt_data_3* findpt_setup_3( const real* const xw[3],
00406                                const unsigned n[3],
00407                                uint nel,
00408                                uint max_hash_size,
00409                                real bbox_tol );
00410 
00411 void findpt_free_2( findpt_data_2* p );
00412 void findpt_free_3( findpt_data_3* p );
00413 
00414 const real* findpt_allbnd_2( const findpt_data_2* p );
00415 const real* findpt_allbnd_3( const findpt_data_3* p );
00416 
00417 typedef int ( *findpt_func )( void*, const real*, int, uint*, real*, real* );
00418 int findpt_2( findpt_data_2* p, const real x[2], int guess, uint* el, real r[2], real* dist );
00419 int findpt_3( findpt_data_3* p, const real x[3], int guess, uint* el, real r[3], real* dist );
00420 
00421 inline void findpt_weights_2( findpt_data_2* p, const real r[2] )
00422 {
00423     lagrange_0( &p->ld[0], r[0] );
00424     lagrange_0( &p->ld[1], r[1] );
00425 }
00426 
00427 inline void findpt_weights_3( findpt_data_3* p, const real r[3] )
00428 {
00429     lagrange_0( &p->ld[0], r[0] );
00430     lagrange_0( &p->ld[1], r[1] );
00431     lagrange_0( &p->ld[2], r[2] );
00432 }
00433 
00434 inline double findpt_eval_2( findpt_data_2* p, const real* u )
00435 {
00436     return tensor_i2( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, u, p->od_work );
00437 }
00438 
00439 inline double findpt_eval_3( findpt_data_3* p, const real* u )
00440 {
00441     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 );
00442 }
00443 
00444 //======================================================
00445 // from extrafindpt.h
00446 //======================================================
00447 
00448 typedef struct
00449 {
00450     unsigned constraints;
00451     unsigned dn, d1, d2;
00452     real *x[3], *fdn[3];
00453 } opt_face_data_3;
00454 
00455 typedef struct
00456 {
00457     unsigned constraints;
00458     unsigned de, d1, d2;
00459     real *x[3], *fd1[3], *fd2[3];
00460 } opt_edge_data_3;
00461 
00462 typedef struct
00463 {
00464     unsigned constraints;
00465     real x[3], jac[9];
00466 } opt_point_data_3;
00467 
00468 typedef struct
00469 {
00470     lagrange_data* ld;
00471     unsigned size[4];
00472     const real* elx[3];
00473     opt_face_data_3 fd;
00474     opt_edge_data_3 ed;
00475     opt_point_data_3 pd;
00476     real* work;
00477     real x[3], jac[9];
00478 } opt_data_3;
00479 
00480 void opt_alloc_3( opt_data_3* p, lagrange_data* ld );
00481 void opt_free_3( opt_data_3* p );
00482 double opt_findpt_3( opt_data_3* p, const real* const elx[3], const real xstar[3], real r[3], unsigned* constr );
00483 void opt_vol_set_intp_3( opt_data_3* p, const real r[3] );
00484 
00485 const unsigned opt_no_constraints_2 = 3 + 1;
00486 const unsigned opt_no_constraints_3 = 9 + 3 + 1;
00487 
00488 /* for 2d spectralQuad */
00489 /*--------------------------------------------------------------------------
00490 
00491    2 - D
00492 
00493   --------------------------------------------------------------------------*/
00494 
00495 typedef struct
00496 {
00497     unsigned constraints;
00498     unsigned de, d1;
00499     real *x[2], *fd1[2];
00500 } opt_edge_data_2;
00501 
00502 typedef struct
00503 {
00504     unsigned constraints;
00505     real x[2], jac[4];
00506 } opt_point_data_2;
00507 
00508 typedef struct
00509 {
00510     lagrange_data* ld;
00511     unsigned size[3];
00512     const real* elx[2];
00513     opt_edge_data_2 ed;
00514     opt_point_data_2 pd;
00515     real* work;
00516     real x[2], jac[4];
00517 } opt_data_2;
00518 void opt_alloc_2( opt_data_2* p, lagrange_data* ld );
00519 void opt_free_2( opt_data_2* p );
00520 double opt_findpt_2( opt_data_2* p, const real* const elx[2], const real xstar[2], real r[2], unsigned* constr );
00521 
00522 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines