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
src
moab
LocalDiscretization
SpectralFuncs.hpp
Generated on Fri Dec 13 2024 02:06:49 for Mesh Oriented datABase by
1.9.1.