moab
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 /*--------------------------------------------------------------------------
00120    Quadrature Nodes and Weights Calculation
00121    
00122    call the _nodes function before calling the _weights function
00123   --------------------------------------------------------------------------*/
00124 
00125 void gauss_nodes(real *z, int n);   /* n nodes (order = 2n-1) */
00126 void lobatto_nodes(real *z, int n); /* n nodes (order = 2n-3) */
00127 
00128 void gauss_weights(const real *z, real *w, int n);
00129 void lobatto_weights(const real *z, real *w, int n);
00130 
00131 /*--------------------------------------------------------------------------
00132    Lagrangian to Legendre Change-of-basis Matrix
00133   --------------------------------------------------------------------------*/
00134 
00135 /* precondition: n >= 2
00136    given the Gauss quadrature rule (z,w,n), compute the square matrix J
00137    for transforming from the Gauss basis to the Legendre basis:
00138    
00139       u_legendre(i) = sum_j J(i,j) u_gauss(j)
00140 
00141    computes J   = .5 (2i+1) w  P (z )
00142              ij              j  i  j
00143              
00144    in column major format (inner index is i, the Legendre index)
00145  */
00146 void gauss_to_legendre(const real *z, const real *w, int n, real *J);
00147 
00148 /* precondition: n >= 2
00149    same as above, but
00150    in row major format (inner index is j, the Gauss index)
00151  */
00152 void gauss_to_legendre_t(const real *z, const real *w, int n, real *J);
00153 
00154 /* precondition: n >= 3
00155    given the Lobatto quadrature rule (z,w,n), compute the square matrix J
00156    for transforming from the Lobatto basis to the Legendre basis:
00157    
00158       u_legendre(i) = sum_j J(i,j) u_lobatto(j)
00159 
00160    in column major format (inner index is i, the Legendre index)
00161  */
00162 void lobatto_to_legendre(const real *z, const real *w, int n, real *J);
00163 
00164 /*--------------------------------------------------------------------------
00165    Lagrangian basis function evaluation
00166   --------------------------------------------------------------------------*/
00167 
00168 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
00169    evaluate all Lagrangian basis functions at all points x
00170    
00171    inner index of output J is the basis function index (row-major format)
00172    provide work array with space for 4*n doubles
00173  */
00174 void lagrange_weights(const real *z, unsigned n,
00175                       const real *x, unsigned m,
00176                       real *J, real *work);
00177 
00178 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
00179    evaluate all Lagrangian basis functions and their derivatives
00180    
00181    inner index of outputs J,D is the basis function index (row-major format)
00182    provide work array with space for 6*n doubles
00183  */
00184 void lagrange_weights_deriv(const real *z, unsigned n,
00185                             const real *x, unsigned m,
00186                             real *J, real *D, real *work);
00187 
00188 /*--------------------------------------------------------------------------
00189    Speedy Lagrangian Interpolation
00190    
00191    Usage:
00192    
00193      lagrange_data p;
00194      lagrange_setup(&p,z,n);    *  setup for nodes z[0 ... n-1] *
00195      
00196      the weights
00197        p->J [0 ... n-1]     interpolation weights
00198        p->D [0 ... n-1]     1st derivative weights
00199        p->D2[0 ... n-1]     2nd derivative weights
00200      are computed for a given x with:
00201        lagrange_0(p,x);  *  compute p->J *
00202        lagrange_1(p,x);  *  compute p->J, p->D *
00203        lagrange_2(p,x);  *  compute p->J, p->D, p->D2 *
00204        lagrange_2u(p);   *  compute p->D2 after call of lagrange_1(p,x); *
00205      These functions use the z array supplied to setup
00206        (that pointer should not be freed between calls)
00207      Weights for x=z[0] and x=z[n-1] are computed during setup; access as:
00208        p->J_z0, etc. and p->J_zn, etc.
00209 
00210      lagrange_free(&p);  *  deallocate memory allocated by setup *
00211   --------------------------------------------------------------------------*/
00212 
00213 typedef struct {
00214   unsigned n;                /* number of Lagrange nodes            */
00215   const real *z;             /* Lagrange nodes (user-supplied)      */
00216   real *J, *D, *D2;          /* weights for 0th,1st,2nd derivatives */
00217   real *J_z0, *D_z0, *D2_z0; /* ditto at z[0]   (computed at setup) */
00218   real *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */
00219   real *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data           */
00220 } lagrange_data;
00221 
00222 void lagrange_setup(lagrange_data *p, const real *z, unsigned n);
00223 void lagrange_free(lagrange_data *p);
00224 
00225 void lagrange_0(lagrange_data *p, real x) ;
00226 void lagrange_1(lagrange_data *p, real x) ;
00227 void lagrange_2(lagrange_data *p, real x) ;
00228 void lagrange_2u(lagrange_data *p) ;
00229 
00230 //======================================================
00231 // from tensor.h
00232 //======================================================
00233 
00234 /*--------------------------------------------------------------------------
00235    1-,2-,3-d Tensor Application
00236    
00237    the 3d case:
00238    tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2)
00239      gives v = [ R (x) S (x) T ] u
00240      where R is mr x nr, S is ms x ns, T is mt x nt,
00241        each in row- or column-major format according to f := r | c
00242      u is nr x ns x nt in column-major format (inner index is r)
00243      v is mr x ms x mt in column-major format (inner index is r)
00244   --------------------------------------------------------------------------*/
00245 
00246 void tensor_c1(const real *R, unsigned mr, unsigned nr, 
00247                const real *u, real *v);
00248 void tensor_r1(const real *R, unsigned mr, unsigned nr, 
00249                const real *u, real *v);
00250 
00251 /* work holds mr*ns reals */
00252 void tensor_c2(const real *R, unsigned mr, unsigned nr,
00253                const real *S, unsigned ms, unsigned ns,
00254                const real *u, real *v, real *work);
00255 void tensor_r2(const real *R, unsigned mr, unsigned nr,
00256                const real *S, unsigned ms, unsigned ns,
00257                const real *u, real *v, real *work);
00258 
00259 /* work1 holds mr*ns*nt reals,
00260    work2 holds mr*ms*nt reals */
00261 void tensor_c3(const real *R, unsigned mr, unsigned nr,
00262                const real *S, unsigned ms, unsigned ns,
00263                const real *T, unsigned mt, unsigned nt,
00264                const real *u, real *v, real *work1, real *work2);
00265 void tensor_r3(const real *R, unsigned mr, unsigned nr,
00266                const real *S, unsigned ms, unsigned ns,
00267                const real *T, unsigned mt, unsigned nt,
00268                const real *u, real *v, real *work1, real *work2);
00269 
00270 /*--------------------------------------------------------------------------
00271    1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
00272    
00273    the 3d case:
00274    v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
00275    same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
00276      gives v = [ Jr (x) Js (x) Jt ] u
00277      where Jr, Js, Jt are row vectors (interpolation weights)
00278      u is nr x ns x nt in column-major format (inner index is r)
00279      v is a scalar
00280   --------------------------------------------------------------------------*/
00281 
00282 real tensor_i1(const real *Jr, unsigned nr, const real *u);
00283 
00284 /* work holds ns reals */
00285 real tensor_i2(const real *Jr, unsigned nr,
00286                const real *Js, unsigned ns,
00287                const real *u, real *work);
00288 
00289 /* work holds ns*nt + nt reals */
00290 real tensor_i3(const real *Jr, unsigned nr,
00291                const real *Js, unsigned ns,
00292                const real *Jt, unsigned nt,
00293                const real *u, real *work);
00294 
00295 /*--------------------------------------------------------------------------
00296    1-,2-,3-d Tensor Application of Row Vectors
00297              for simultaneous Interpolation and Gradient computation
00298    
00299    the 3d case:
00300    v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
00301      gives v   = [ Jr (x) Js (x) Jt ] u
00302            g_0 = [ Dr (x) Js (x) Jt ] u
00303            g_1 = [ Jr (x) Ds (x) Jt ] u
00304            g_2 = [ Jr (x) Js (x) Dt ] u
00305      where Jr,Dr,Js,Ds,Jt,Dt are row vectors
00306        (interpolation & derivative weights)
00307      u is nr x ns x nt in column-major format (inner index is r)
00308      v is a scalar, g is an array of 3 reals
00309   --------------------------------------------------------------------------*/
00310 
00311 real tensor_ig1(const real *Jr, const real *Dr, unsigned nr,
00312                 const real *u, real *g);
00313 
00314 /* work holds 2*ns reals */
00315 real tensor_ig2(const real *Jr, const real *Dr, unsigned nr,
00316                 const real *Js, const real *Ds, unsigned ns,
00317                 const real *u, real *g, real *work);
00318 
00319 /* work holds 2*ns*nt + 3*ns reals */
00320 real tensor_ig3(const real *Jr, const real *Dr, unsigned nr,
00321                 const real *Js, const real *Ds, unsigned ns,
00322                 const real *Jt, const real *Dt, unsigned nt,
00323                 const real *u, real *g, real *work);
00324 
00325 //======================================================
00326 // from findpt.h
00327 //======================================================
00328 
00329 typedef struct {
00330   const real *xw[2];   /* geometry data */
00331   real *z[2];          /* lobatto nodes */
00332   lagrange_data ld[2]; /* interpolation, derivative weights & data */
00333   unsigned nptel;      /* nodes per element */
00334   struct findpt_hash_data_2 *hash;   /* geometric hashing data */
00335   struct findpt_listel *list, **sorted, **end;
00336                                         /* pre-allocated list of elements to
00337                                            check (found by hashing), and
00338                                            pre-allocated list of pointers into
00339                                            the first list for sorting */
00340   struct findpt_opt_data_2 *od; /* data for the optimization algorithm */
00341   real *od_work;
00342 } findpt_data_2;
00343 
00344 typedef struct {
00345   const real *xw[3];   /* geometry data */
00346   real *z[3];          /* lobatto nodes */
00347   lagrange_data ld[3]; /* interpolation, derivative weights & data */
00348   unsigned nptel;      /* nodes per element */
00349   struct findpt_hash_data_3 *hash;   /* geometric hashing data */
00350   struct findpt_listel *list, **sorted, **end;
00351                                         /* pre-allocated list of elements to
00352                                            check (found by hashing), and
00353                                            pre-allocated list of pointers into
00354                                            the first list for sorting */
00355   struct findpt_opt_data_3 *od; /* data for the optimization algorithm */
00356   real *od_work;
00357 } findpt_data_3;
00358 
00359 findpt_data_2 *findpt_setup_2(
00360           const real *const xw[2], const unsigned n[2], uint nel,
00361           uint max_hash_size, real bbox_tol);
00362 findpt_data_3 *findpt_setup_3(
00363           const real *const xw[3], const unsigned n[3], uint nel,
00364           uint max_hash_size, real bbox_tol);
00365 
00366 void findpt_free_2(findpt_data_2 *p);
00367 void findpt_free_3(findpt_data_3 *p);
00368 
00369 const real *findpt_allbnd_2(const findpt_data_2 *p);
00370 const real *findpt_allbnd_3(const findpt_data_3 *p);
00371 
00372 typedef int (*findpt_func)(void *, const real *, int, uint *, real *, real *);
00373 int findpt_2(findpt_data_2 *p, const real x[2], int guess,
00374              uint *el, real r[2], real *dist);
00375 int findpt_3(findpt_data_3 *p, const real x[3], int guess,
00376              uint *el, real r[3], real *dist);
00377 
00378 inline void findpt_weights_2(findpt_data_2 *p, const real r[2])
00379 {
00380   lagrange_0(&p->ld[0],r[0]);
00381   lagrange_0(&p->ld[1],r[1]);
00382 }
00383 
00384 inline void findpt_weights_3(findpt_data_3 *p, const real r[3])
00385 {
00386   lagrange_0(&p->ld[0],r[0]);
00387   lagrange_0(&p->ld[1],r[1]);
00388   lagrange_0(&p->ld[2],r[2]);
00389 }
00390 
00391 inline double findpt_eval_2(findpt_data_2 *p, const real *u)
00392 {
00393   return tensor_i2(p->ld[0].J,p->ld[0].n,
00394                    p->ld[1].J,p->ld[1].n,
00395                    u, p->od_work);
00396 }
00397 
00398 inline double findpt_eval_3(findpt_data_3 *p, const real *u)
00399 {
00400   return tensor_i3(p->ld[0].J,p->ld[0].n,
00401                    p->ld[1].J,p->ld[1].n,
00402                    p->ld[2].J,p->ld[2].n,
00403                    u, p->od_work);
00404 }
00405 
00406 //======================================================
00407 // from extrafindpt.h
00408 //======================================================
00409 
00410 typedef struct {
00411   unsigned constraints;
00412   unsigned dn, d1, d2;
00413   real *x[3], *fdn[3];
00414 } opt_face_data_3;
00415 
00416 typedef struct {
00417   unsigned constraints;
00418   unsigned de, d1, d2;
00419   real *x[3], *fd1[3], *fd2[3];
00420 } opt_edge_data_3;
00421 
00422 typedef struct {
00423   unsigned constraints;
00424   real x[3], jac[9];
00425 } opt_point_data_3;
00426 
00427 typedef struct {
00428   lagrange_data *ld;
00429   unsigned size[4];
00430   const real *elx[3];
00431   opt_face_data_3 fd;
00432   opt_edge_data_3 ed;
00433   opt_point_data_3 pd;
00434   real *work;
00435   real x[3], jac[9];
00436 } opt_data_3;
00437 
00438 
00439 void opt_alloc_3(opt_data_3 *p, lagrange_data *ld);
00440 void opt_free_3(opt_data_3 *p);
00441 double opt_findpt_3(opt_data_3 *p, const real *const elx[3],
00442                            const real xstar[3], real r[3], unsigned *constr);
00443 void opt_vol_set_intp_3(opt_data_3 *p, const real r[3]);
00444 
00445 const unsigned opt_no_constraints_2 = 3+1;
00446 const unsigned opt_no_constraints_3 = 9+3+1;
00447 
00448 /* for 2d spectralQuad */
00449 /*--------------------------------------------------------------------------
00450 
00451    2 - D
00452 
00453   --------------------------------------------------------------------------*/
00454 
00455 typedef struct {
00456   unsigned constraints;
00457   unsigned de, d1;
00458   real *x[2], *fd1[2];
00459 } opt_edge_data_2;
00460 
00461 typedef struct {
00462   unsigned constraints;
00463   real x[2], jac[4];
00464 } opt_point_data_2;
00465 
00466 typedef struct {
00467   lagrange_data *ld;
00468   unsigned size[3];
00469   const real *elx[2];
00470   opt_edge_data_2 ed;
00471   opt_point_data_2 pd;
00472   real *work;
00473   real x[2], jac[4];
00474 } opt_data_2;
00475 void opt_alloc_2(opt_data_2 *p, lagrange_data *ld);
00476 void opt_free_2(opt_data_2 *p);
00477 double opt_findpt_2(opt_data_2 *p, const real *const elx[2],
00478                            const real xstar[2], real r[2], unsigned *constr);
00479 
00480 #endif
00481 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines