Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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