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