Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SpectralFuncs.hpp
Go to the documentation of this file.
1 #ifndef SPECTRALFUNCS_HPP 2 #define SPECTRALFUNCS_HPP 3  4 #include <cfloat> 5  6 //====================================================== 7 // from types.h 8 //====================================================== 9  10 /* integer type to use for everything */ 11 #if defined( USE_LONG ) 12 #define INTEGER long 13 #elif defined( USE_LONG_LONG ) 14 #define INTEGER long long 15 #elif defined( USE_SHORT ) 16 #define INTEGER short 17 #else 18 #define INTEGER int 19 #endif 20  21 /* when defined, use the given type for global indices instead of INTEGER */ 22 #if defined( USE_GLOBAL_LONG_LONG ) 23 #define GLOBAL_INT long long 24 #elif defined( USE_GLOBAL_LONG ) 25 #define GLOBAL_INT long 26 #else 27 #define GLOBAL_INT long 28 #endif 29  30 /* floating point type to use for everything */ 31 #if defined( USE_FLOAT ) 32 typedef float real; 33 #define floorr floorf 34 #define ceilr ceilf 35 #define sqrtr sqrtf 36 #define fabsr fabsf 37 #define cosr cosf 38 #define sinr sinf 39 #define EPS ( 128 * FLT_EPSILON ) 40 #define PI 3.1415926535897932384626433832795028841971693993751058209749445923F 41 #elif defined( USE_LONG_DOUBLE ) 42 typedef long double real; 43 #define floorr floorl 44 #define ceilr ceill 45 #define sqrtr sqrtl 46 #define fabsr fabsl 47 #define cosr cosl 48 #define sinr sinl 49 #define EPS ( 128 * LDBL_EPSILON ) 50 #define PI 3.1415926535897932384626433832795028841971693993751058209749445923L 51 #else 52 typedef double real; 53 #define floorr floor 54 #define ceilr ceil 55 #define sqrtr sqrt 56 #define fabsr fabs 57 #define cosr cos 58 #define sinr sin 59 #define EPS ( 128 * DBL_EPSILON ) 60 #define PI 3.1415926535897932384626433832795028841971693993751058209749445923 61 #endif 62  63 /* apparently uint and ulong can be defined already in standard headers */ 64 #define uint uint_ 65 #define ulong ulong_ 66 #define sint sint_ 67 #define slong slong_ 68  69 typedef signed INTEGER sint; 70 typedef unsigned INTEGER uint; 71 #undef INTEGER 72  73 #ifdef GLOBAL_INT 74 typedef signed GLOBAL_INT slong; 75 typedef unsigned GLOBAL_INT ulong; 76 #else 77 typedef sint slong; 78 typedef uint ulong; 79 #endif 80  81 //====================================================== 82 // from poly.h 83 //====================================================== 84  85 /* 86  For brevity's sake, some names have been shortened 87  Quadrature rules 88  Gauss -> Gauss-Legendre quadrature (open) 89  Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends) 90  Polynomial bases 91  Legendre -> Legendre basis 92  Gauss -> Lagrangian basis using Gauss quadrature nodes 93  Lobatto -> Lagrangian basis using Lobatto quadrature nodes 94 */ 95  96 /*-------------------------------------------------------------------------- 97  Legendre Polynomial Matrix Computation 98  (compute P_i(x_j) for i = 0, ..., n and a given set of x) 99  --------------------------------------------------------------------------*/ 100  101 /* precondition: n >= 1 102  inner index is x index (0 ... m-1); 103  outer index is Legendre polynomial number (0 ... n) 104  */ 105 void legendre_matrix( const real* x, int m, real* P, int n ); 106  107 /* precondition: n >= 1 108  inner index is Legendre polynomial number (0 ... n) 109  outer index is x index (0 ... m-1); 110  */ 111 void legendre_matrix_t( const real* x, int m, real* P, int n ); 112  113 /* precondition: n >= 1 114  compute P_i(x) with i = 0 ... n 115  */ 116 void legendre_row( real x, real* P, int n ); 117  118 /*-------------------------------------------------------------------------- 119  Quadrature Nodes and Weights Calculation 120  121  call the _nodes function before calling the _weights function 122  --------------------------------------------------------------------------*/ 123  124 void gauss_nodes( real* z, int n ); /* n nodes (order = 2n-1) */ 125 void lobatto_nodes( real* z, int n ); /* n nodes (order = 2n-3) */ 126  127 void gauss_weights( const real* z, real* w, int n ); 128 void lobatto_weights( const real* z, real* w, int n ); 129  130 /*-------------------------------------------------------------------------- 131  Lagrangian to Legendre Change-of-basis Matrix 132  --------------------------------------------------------------------------*/ 133  134 /* precondition: n >= 2 135  given the Gauss quadrature rule (z,w,n), compute the square matrix J 136  for transforming from the Gauss basis to the Legendre basis: 137  138  u_legendre(i) = sum_j J(i,j) u_gauss(j) 139  140  computes J = .5 (2i+1) w P (z ) 141  ij j i j 142  143  in column major format (inner index is i, the Legendre index) 144  */ 145 void gauss_to_legendre( const real* z, const real* w, int n, real* J ); 146  147 /* precondition: n >= 2 148  same as above, but 149  in row major format (inner index is j, the Gauss index) 150  */ 151 void gauss_to_legendre_t( const real* z, const real* w, int n, real* J ); 152  153 /* precondition: n >= 3 154  given the Lobatto quadrature rule (z,w,n), compute the square matrix J 155  for transforming from the Lobatto basis to the Legendre basis: 156  157  u_legendre(i) = sum_j J(i,j) u_lobatto(j) 158  159  in column major format (inner index is i, the Legendre index) 160  */ 161 void lobatto_to_legendre( const real* z, const real* w, int n, real* J ); 162  163 /*-------------------------------------------------------------------------- 164  Lagrangian basis function evaluation 165  --------------------------------------------------------------------------*/ 166  167 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 168  evaluate all Lagrangian basis functions at all points x 169  170  inner index of output J is the basis function index (row-major format) 171  provide work array with space for 4*n doubles 172  */ 173 void lagrange_weights( const real* z, unsigned n, const real* x, unsigned m, real* J, real* work ); 174  175 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 176  evaluate all Lagrangian basis functions and their derivatives 177  178  inner index of outputs J,D is the basis function index (row-major format) 179  provide work array with space for 6*n doubles 180  */ 181 void lagrange_weights_deriv( const real* z, unsigned n, const real* x, unsigned m, real* J, real* D, real* work ); 182  183 /*-------------------------------------------------------------------------- 184  Speedy Lagrangian Interpolation 185  186  Usage: 187  188  lagrange_data p; 189  lagrange_setup(&p,z,n); * setup for nodes z[0 ... n-1] * 190  191  the weights 192  p->J [0 ... n-1] interpolation weights 193  p->D [0 ... n-1] 1st derivative weights 194  p->D2[0 ... n-1] 2nd derivative weights 195  are computed for a given x with: 196  lagrange_0(p,x); * compute p->J * 197  lagrange_1(p,x); * compute p->J, p->D * 198  lagrange_2(p,x); * compute p->J, p->D, p->D2 * 199  lagrange_2u(p); * compute p->D2 after call of lagrange_1(p,x); * 200  These functions use the z array supplied to setup 201  (that pointer should not be freed between calls) 202  Weights for x=z[0] and x=z[n-1] are computed during setup; access as: 203  p->J_z0, etc. and p->J_zn, etc. 204  205  lagrange_free(&p); * deallocate memory allocated by setup * 206  --------------------------------------------------------------------------*/ 207  208 typedef struct 209 { 210  unsigned n; /* number of Lagrange nodes */ 211  const real* z; /* Lagrange nodes (user-supplied) */ 212  real *J, *D, *D2; /* weights for 0th,1st,2nd derivatives */ 213  real *J_z0, *D_z0, *D2_z0; /* ditto at z[0] (computed at setup) */ 214  real *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */ 215  real *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data */ 216 } lagrange_data; 217  218 void lagrange_setup( lagrange_data* p, const real* z, unsigned n ); 219 void lagrange_free( lagrange_data* p ); 220  221 void lagrange_0( lagrange_data* p, real x ); 222 void lagrange_1( lagrange_data* p, real x ); 223 void lagrange_2( lagrange_data* p, real x ); 224 void lagrange_2u( lagrange_data* p ); 225  226 //====================================================== 227 // from tensor.h 228 //====================================================== 229  230 /*-------------------------------------------------------------------------- 231  1-,2-,3-d Tensor Application 232  233  the 3d case: 234  tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2) 235  gives v = [ R (x) S (x) T ] u 236  where R is mr x nr, S is ms x ns, T is mt x nt, 237  each in row- or column-major format according to f := r | c 238  u is nr x ns x nt in column-major format (inner index is r) 239  v is mr x ms x mt in column-major format (inner index is r) 240  --------------------------------------------------------------------------*/ 241  242 void tensor_c1( const real* R, unsigned mr, unsigned nr, const real* u, real* v ); 243 void tensor_r1( const real* R, unsigned mr, unsigned nr, const real* u, real* v ); 244  245 /* work holds mr*ns reals */ 246 void tensor_c2( const real* R, 247  unsigned mr, 248  unsigned nr, 249  const real* S, 250  unsigned ms, 251  unsigned ns, 252  const real* u, 253  real* v, 254  real* work ); 255 void tensor_r2( const real* R, 256  unsigned mr, 257  unsigned nr, 258  const real* S, 259  unsigned ms, 260  unsigned ns, 261  const real* u, 262  real* v, 263  real* work ); 264  265 /* work1 holds mr*ns*nt reals, 266  work2 holds mr*ms*nt reals */ 267 void tensor_c3( const real* R, 268  unsigned mr, 269  unsigned nr, 270  const real* S, 271  unsigned ms, 272  unsigned ns, 273  const real* T, 274  unsigned mt, 275  unsigned nt, 276  const real* u, 277  real* v, 278  real* work1, 279  real* work2 ); 280 void tensor_r3( const real* R, 281  unsigned mr, 282  unsigned nr, 283  const real* S, 284  unsigned ms, 285  unsigned ns, 286  const real* T, 287  unsigned mt, 288  unsigned nt, 289  const real* u, 290  real* v, 291  real* work1, 292  real* work2 ); 293  294 /*-------------------------------------------------------------------------- 295  1-,2-,3-d Tensor Application of Row Vectors (for Interpolation) 296  297  the 3d case: 298  v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work) 299  same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2): 300  gives v = [ Jr (x) Js (x) Jt ] u 301  where Jr, Js, Jt are row vectors (interpolation weights) 302  u is nr x ns x nt in column-major format (inner index is r) 303  v is a scalar 304  --------------------------------------------------------------------------*/ 305  306 real tensor_i1( const real* Jr, unsigned nr, const real* u ); 307  308 /* work holds ns reals */ 309 real tensor_i2( const real* Jr, unsigned nr, const real* Js, unsigned ns, const real* u, real* work ); 310  311 /* work holds ns*nt + nt reals */ 312 real tensor_i3( const real* Jr, 313  unsigned nr, 314  const real* Js, 315  unsigned ns, 316  const real* Jt, 317  unsigned nt, 318  const real* u, 319  real* work ); 320  321 /*-------------------------------------------------------------------------- 322  1-,2-,3-d Tensor Application of Row Vectors 323  for simultaneous Interpolation and Gradient computation 324  325  the 3d case: 326  v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work) 327  gives v = [ Jr (x) Js (x) Jt ] u 328  g_0 = [ Dr (x) Js (x) Jt ] u 329  g_1 = [ Jr (x) Ds (x) Jt ] u 330  g_2 = [ Jr (x) Js (x) Dt ] u 331  where Jr,Dr,Js,Ds,Jt,Dt are row vectors 332  (interpolation & derivative weights) 333  u is nr x ns x nt in column-major format (inner index is r) 334  v is a scalar, g is an array of 3 reals 335  --------------------------------------------------------------------------*/ 336  337 real tensor_ig1( const real* Jr, const real* Dr, unsigned nr, const real* u, real* g ); 338  339 /* work holds 2*ns reals */ 340 real tensor_ig2( const real* Jr, 341  const real* Dr, 342  unsigned nr, 343  const real* Js, 344  const real* Ds, 345  unsigned ns, 346  const real* u, 347  real* g, 348  real* work ); 349  350 /* work holds 2*ns*nt + 3*ns reals */ 351 real tensor_ig3( const real* Jr, 352  const real* Dr, 353  unsigned nr, 354  const real* Js, 355  const real* Ds, 356  unsigned ns, 357  const real* Jt, 358  const real* Dt, 359  unsigned nt, 360  const real* u, 361  real* g, 362  real* work ); 363  364 //====================================================== 365 // from findpt.h 366 //====================================================== 367  368 typedef struct 369 { 370  const real* xw[2]; /* geometry data */ 371  real* z[2]; /* lobatto nodes */ 372  lagrange_data ld[2]; /* interpolation, derivative weights & data */ 373  unsigned nptel; /* nodes per element */ 374  struct findpt_hash_data_2* hash; /* geometric hashing data */ 375  struct findpt_listel *list, **sorted, **end; 376  /* pre-allocated list of elements to 377  check (found by hashing), and 378  pre-allocated list of pointers into 379  the first list for sorting */ 380  struct findpt_opt_data_2* od; /* data for the optimization algorithm */ 381  real* od_work; 382 } findpt_data_2; 383  384 typedef struct 385 { 386  const real* xw[3]; /* geometry data */ 387  real* z[3]; /* lobatto nodes */ 388  lagrange_data ld[3]; /* interpolation, derivative weights & data */ 389  unsigned nptel; /* nodes per element */ 390  struct findpt_hash_data_3* hash; /* geometric hashing data */ 391  struct findpt_listel *list, **sorted, **end; 392  /* pre-allocated list of elements to 393  check (found by hashing), and 394  pre-allocated list of pointers into 395  the first list for sorting */ 396  struct findpt_opt_data_3* od; /* data for the optimization algorithm */ 397  real* od_work; 398 } findpt_data_3; 399  400 findpt_data_2* findpt_setup_2( const real* const xw[2], 401  const unsigned n[2], 402  uint nel, 403  uint max_hash_size, 404  real bbox_tol ); 405 findpt_data_3* findpt_setup_3( const real* const xw[3], 406  const unsigned n[3], 407  uint nel, 408  uint max_hash_size, 409  real bbox_tol ); 410  411 void findpt_free_2( findpt_data_2* p ); 412 void findpt_free_3( findpt_data_3* p ); 413  414 const real* findpt_allbnd_2( const findpt_data_2* p ); 415 const real* findpt_allbnd_3( const findpt_data_3* p ); 416  417 typedef int ( *findpt_func )( void*, const real*, int, uint*, real*, real* ); 418 int findpt_2( findpt_data_2* p, const real x[2], int guess, uint* el, real r[2], real* dist ); 419 int findpt_3( findpt_data_3* p, const real x[3], int guess, uint* el, real r[3], real* dist ); 420  421 inline void findpt_weights_2( findpt_data_2* p, const real r[2] ) 422 { 423  lagrange_0( &p->ld[0], r[0] ); 424  lagrange_0( &p->ld[1], r[1] ); 425 } 426  427 inline void findpt_weights_3( findpt_data_3* p, const real r[3] ) 428 { 429  lagrange_0( &p->ld[0], r[0] ); 430  lagrange_0( &p->ld[1], r[1] ); 431  lagrange_0( &p->ld[2], r[2] ); 432 } 433  434 inline double findpt_eval_2( findpt_data_2* p, const real* u ) 435 { 436  return tensor_i2( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, u, p->od_work ); 437 } 438  439 inline double findpt_eval_3( findpt_data_3* p, const real* u ) 440 { 441  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 ); 442 } 443  444 //====================================================== 445 // from extrafindpt.h 446 //====================================================== 447  448 typedef struct 449 { 450  unsigned constraints; 451  unsigned dn, d1, d2; 452  real *x[3], *fdn[3]; 453 } opt_face_data_3; 454  455 typedef struct 456 { 457  unsigned constraints; 458  unsigned de, d1, d2; 459  real *x[3], *fd1[3], *fd2[3]; 460 } opt_edge_data_3; 461  462 typedef struct 463 { 464  unsigned constraints; 465  real x[3], jac[9]; 466 } opt_point_data_3; 467  468 typedef struct 469 { 470  lagrange_data* ld; 471  unsigned size[4]; 472  const real* elx[3]; 473  opt_face_data_3 fd; 474  opt_edge_data_3 ed; 475  opt_point_data_3 pd; 476  real* work; 477  real x[3], jac[9]; 478 } opt_data_3; 479  480 void opt_alloc_3( opt_data_3* p, lagrange_data* ld ); 481 void opt_free_3( opt_data_3* p ); 482 double opt_findpt_3( opt_data_3* p, const real* const elx[3], const real xstar[3], real r[3], unsigned* constr ); 483 void opt_vol_set_intp_3( opt_data_3* p, const real r[3] ); 484  485 const unsigned opt_no_constraints_2 = 3 + 1; 486 const unsigned opt_no_constraints_3 = 9 + 3 + 1; 487  488 /* for 2d spectralQuad */ 489 /*-------------------------------------------------------------------------- 490  491  2 - D 492  493  --------------------------------------------------------------------------*/ 494  495 typedef struct 496 { 497  unsigned constraints; 498  unsigned de, d1; 499  real *x[2], *fd1[2]; 500 } opt_edge_data_2; 501  502 typedef struct 503 { 504  unsigned constraints; 505  real x[2], jac[4]; 506 } opt_point_data_2; 507  508 typedef struct 509 { 510  lagrange_data* ld; 511  unsigned size[3]; 512  const real* elx[2]; 513  opt_edge_data_2 ed; 514  opt_point_data_2 pd; 515  real* work; 516  real x[2], jac[4]; 517 } opt_data_2; 518 void opt_alloc_2( opt_data_2* p, lagrange_data* ld ); 519 void opt_free_2( opt_data_2* p ); 520 double opt_findpt_2( opt_data_2* p, const real* const elx[2], const real xstar[2], real r[2], unsigned* constr ); 521  522 #endif