DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
QuadratureHub.hpp
Go to the documentation of this file.
1#pragma once
2
3// =============================================================================
4// QuadratureHub.hpp - Master include for all quadrature rules
5// =============================================================================
6// This file includes all individual quadrature data files and provides:
7// 1. GetQuadratureScheme() - Select appropriate scheme for given order
8// 2. GetQuadraturePoint() - Retrieve quadrature points and weights
9//
10// The design follows the modular pattern used for Elements:
11// - QuadratureBase.hpp: Shared constants and type aliases
12// - Line.hpp, Tri.hpp, Tet.hpp, etc.: Individual parametric space data
13// - GaussJacobi.hpp: Special quadrature for pyramid elements
14// - QuadratureHub.hpp: Unified interface (this file)
15// =============================================================================
16
17#include "Elements.hpp"
18#include "Quadratures/QuadratureBase.hpp"
19#include "Quadratures/Line.hpp"
20#include "Quadratures/Tri.hpp"
21#include "Quadratures/Tet.hpp"
22#include "Quadratures/GaussJacobi.hpp"
23
24namespace DNDS::Geom::Elem
25{
26 // =========================================================================
27 // Scheme Selection
28 // =========================================================================
29 // Returns the number of quadrature points for the given parametric space
30 // and integration order. Returns 0 if the order is not supported.
31 //
32 // The returned value is both:
33 // - The scheme identifier (INT_SCHEME_*)
34 // - The number of quadrature points
35 // =========================================================================
36
37 inline constexpr t_index GetQuadratureScheme(ParamSpace ps, int int_order)
38 {
39 if (ps == LineSpace)
40 switch (int_order)
41 {
42 case 0:
43 case 1:
44 return INT_SCHEME_Line_1;
45 case 2:
46 case 3:
47 return INT_SCHEME_Line_2;
48 case 4:
49 case 5:
50 return INT_SCHEME_Line_3;
51 case 6:
52 return INT_SCHEME_Line_4;
53 default:
54 return 0;
55 }
56
57 if (ps == TriSpace)
58 switch (int_order)
59 {
60 case 0:
61 case 1:
62 return INT_SCHEME_Tri_1;
63 case 2:
64 return INT_SCHEME_Tri_3;
65 case 3:
66 case 4:
67 return INT_SCHEME_Tri_6;
68 case 5:
69 return INT_SCHEME_Tri_7;
70 case 6:
71 return INT_SCHEME_Tri_12;
72 default:
73 return 0;
74 }
75
76 if (ps == QuadSpace)
77 switch (int_order)
78 {
79 case 0:
80 case 1:
81 return INT_SCHEME_Quad_1;
82 case 2:
83 case 3:
84 return INT_SCHEME_Quad_4;
85 case 4:
86 case 5:
87 return INT_SCHEME_Quad_9;
88 case 6:
89 return INT_SCHEME_Quad_16;
90 default:
91 return 0;
92 }
93
94 if (ps == TetSpace)
95 switch (int_order)
96 {
97 case 0:
98 case 1:
99 return INT_SCHEME_Tet_1;
100 case 2:
101 return INT_SCHEME_Tet_4;
102 case 3:
103 return INT_SCHEME_Tet_8;
104 case 4:
105 case 5:
106 return INT_SCHEME_Tet_14;
107 case 6:
108 return INT_SCHEME_Tet_24;
109 default:
110 return 0;
111 }
112
113 if (ps == HexSpace)
114 switch (int_order)
115 {
116 case 0:
117 case 1:
118 return INT_SCHEME_Hex_1;
119 case 2:
120 case 3:
121 return INT_SCHEME_Hex_8;
122 case 4:
123 case 5:
124 return INT_SCHEME_Hex_27;
125 case 6:
126 return INT_SCHEME_Hex_64;
127 default:
128 return 0;
129 }
130
131 if (ps == PyramidSpace)
132 switch (int_order)
133 {
134 case 0:
135 case 1:
136 return INT_SCHEME_Pyramid_1;
137 case 2:
138 case 3:
139 return INT_SCHEME_Pyramid_8;
140 case 4:
141 case 5:
142 return INT_SCHEME_Pyramid_27;
143 case 6:
144 return INT_SCHEME_Pyramid_64;
145 default:
146 return 0;
147 }
148
149 if (ps == PrismSpace)
150 switch (int_order)
151 {
152 case 0:
153 case 1:
154 return INT_SCHEME_Prism_1;
155 case 2:
156 return INT_SCHEME_Prism_6;
157 case 3:
158 case 4:
159 return INT_SCHEME_Prism_18;
160 case 5:
161 return INT_SCHEME_Prism_21;
162 case 6:
163 return INT_SCHEME_Prism_48;
164 default:
165 return 0;
166 }
167 return 0;
168 }
169
170 // =========================================================================
171 // Quadrature Point Access
172 // =========================================================================
173 // Retrieves the iG-th quadrature point (pParam) and weight (w) for the
174 // given parametric space and scheme.
175 //
176 // @param ps Parametric space (LineSpace, TriSpace, etc.)
177 // @param scheme Integration scheme (from GetQuadratureScheme)
178 // @param iG Quadrature point index (0 <= iG < scheme)
179 // @param pParam Output: parametric coordinates [xi, eta, zeta]
180 // @param w Output: quadrature weight
181 //
182 // @warning pParam should be initialized (with 0) before calling
183 // =========================================================================
184
185 template <class TPoint>
186 inline void GetQuadraturePoint(ParamSpace ps, t_index scheme, int iG, TPoint &pParam, t_real &w)
187 {
188 int scheme_size = scheme;
189 DNDS_assert(iG < scheme_size);
190
191 // -----------------------------------------------------------------
192 // Line Space - Gauss-Legendre
193 // -----------------------------------------------------------------
194 if (ps == LineSpace)
195 {
196 switch (scheme)
197 {
198 case INT_SCHEME_Line_1:
199 pParam[0] = __GaussLegendre_1[0][iG];
200 w = __GaussLegendre_1[1][iG];
201 return;
202 case INT_SCHEME_Line_2:
203 pParam[0] = __GaussLegendre_2[0][iG];
204 w = __GaussLegendre_2[1][iG];
205 return;
206 case INT_SCHEME_Line_3:
207 pParam[0] = __GaussLegendre_3[0][iG];
208 w = __GaussLegendre_3[1][iG];
209 return;
210 case INT_SCHEME_Line_4:
211 pParam[0] = __GaussLegendre_4[0][iG];
212 w = __GaussLegendre_4[1][iG];
213 return;
214 default:
215 w = 1e100;
216 }
217 }
218
219 // -----------------------------------------------------------------
220 // Triangle Space - Hammer Rules
221 // -----------------------------------------------------------------
222 if (ps == TriSpace)
223 {
224 switch (scheme)
225 {
226 case INT_SCHEME_Tri_1:
227 pParam[0] = __HammerTri_1[0][iG];
228 pParam[1] = __HammerTri_1[1][iG];
229 w = __HammerTri_1[2][iG];
230 return;
231 case INT_SCHEME_Tri_3:
232 pParam[0] = __HammerTri_3[0][iG];
233 pParam[1] = __HammerTri_3[1][iG];
234 w = __HammerTri_3[2][iG];
235 return;
236 case INT_SCHEME_Tri_6:
237 pParam[0] = __HammerTri_6[0][iG];
238 pParam[1] = __HammerTri_6[1][iG];
239 w = __HammerTri_6[2][iG];
240 return;
241 case INT_SCHEME_Tri_7:
242 pParam[0] = __HammerTri_7[0][iG];
243 pParam[1] = __HammerTri_7[1][iG];
244 w = __HammerTri_7[2][iG];
245 return;
246 case INT_SCHEME_Tri_12:
247 pParam[0] = __HammerTri_12[0][iG];
248 pParam[1] = __HammerTri_12[1][iG];
249 w = __HammerTri_12[2][iG];
250 return;
251 default:
252 w = 1e100;
253 }
254 }
255
256 // -----------------------------------------------------------------
257 // Tetrahedron Space - Hammer Rules
258 // -----------------------------------------------------------------
259 if (ps == TetSpace)
260 {
261 switch (scheme)
262 {
263 case INT_SCHEME_Tet_1:
264 pParam[0] = __HammerTet_1[0][iG];
265 pParam[1] = __HammerTet_1[1][iG];
266 pParam[2] = __HammerTet_1[2][iG];
267 w = __HammerTet_1[3][iG];
268 return;
269 case INT_SCHEME_Tet_4:
270 pParam[0] = __HammerTet_4[0][iG];
271 pParam[1] = __HammerTet_4[1][iG];
272 pParam[2] = __HammerTet_4[2][iG];
273 w = __HammerTet_4[3][iG];
274 return;
275 case INT_SCHEME_Tet_8:
276 pParam[0] = __HammerTet_8[0][iG];
277 pParam[1] = __HammerTet_8[1][iG];
278 pParam[2] = __HammerTet_8[2][iG];
279 w = __HammerTet_8[3][iG];
280 return;
281 case INT_SCHEME_Tet_14:
282 pParam[0] = __HammerTet_14[0][iG];
283 pParam[1] = __HammerTet_14[1][iG];
284 pParam[2] = __HammerTet_14[2][iG];
285 w = __HammerTet_14[3][iG];
286 return;
287 case INT_SCHEME_Tet_24:
288 pParam[0] = __HammerTet_24[0][iG];
289 pParam[1] = __HammerTet_24[1][iG];
290 pParam[2] = __HammerTet_24[2][iG];
291 w = __HammerTet_24[3][iG];
292 return;
293 default:
294 w = 1e100;
295 }
296 }
297
298 // -----------------------------------------------------------------
299 // Quadrilateral Space - Product of 1D Gauss-Legendre
300 // -----------------------------------------------------------------
301 if (ps == QuadSpace)
302 {
303 const t_real *GLData = nullptr;
304 int GLSize = 0;
305 switch (scheme)
306 {
307 case INT_SCHEME_Quad_1:
308 GLData = &(__GaussLegendre_1[0][0]);
309 GLSize = 1;
310 break;
311 case INT_SCHEME_Quad_4:
312 GLData = &(__GaussLegendre_2[0][0]);
313 GLSize = 2;
314 break;
315 case INT_SCHEME_Quad_9:
316 GLData = &(__GaussLegendre_3[0][0]);
317 GLSize = 3;
318 break;
319 case INT_SCHEME_Quad_16:
320 GLData = &(__GaussLegendre_4[0][0]);
321 GLSize = 4;
322 break;
323 default:
324 w = 1e100;
325 return;
326 }
327 int iGi = iG % GLSize;
328 int iGj = iG / GLSize;
329 pParam[0] = GLData[iGi];
330 pParam[1] = GLData[iGj];
331 w = GLData[GLSize + iGi] * GLData[GLSize + iGj];
332 return;
333 }
334
335 // -----------------------------------------------------------------
336 // Hexahedron Space - Product of 1D Gauss-Legendre
337 // -----------------------------------------------------------------
338 if (ps == HexSpace)
339 {
340 const t_real *GLData = nullptr;
341 int GLSize = 0;
342 switch (scheme)
343 {
344 case INT_SCHEME_Hex_1:
345 GLData = &(__GaussLegendre_1[0][0]);
346 GLSize = 1;
347 break;
348 case INT_SCHEME_Hex_8:
349 GLData = &(__GaussLegendre_2[0][0]);
350 GLSize = 2;
351 break;
352 case INT_SCHEME_Hex_27:
353 GLData = &(__GaussLegendre_3[0][0]);
354 GLSize = 3;
355 break;
356 case INT_SCHEME_Hex_64:
357 GLData = &(__GaussLegendre_4[0][0]);
358 GLSize = 4;
359 break;
360 default:
361 w = 1e100;
362 return;
363 }
364 int iGi = iG % GLSize;
365 int iGj = (iG / GLSize) % GLSize;
366 int iGk = (iG / (GLSize * GLSize));
367
368 pParam[0] = GLData[iGi];
369 pParam[1] = GLData[iGj];
370 pParam[2] = GLData[iGk];
371 w = GLData[GLSize + iGi] * GLData[GLSize + iGj] * GLData[GLSize + iGk];
372 return;
373 }
374
375 // -----------------------------------------------------------------
376 // Pyramid Space - Product of Gauss-Legendre and Gauss-Jacobi
377 // -----------------------------------------------------------------
378 if (ps == PyramidSpace)
379 {
380 const t_real *GLData = nullptr;
381 const t_real *GJData = nullptr;
382 int GLSize = 0; // == GJSize
383 switch (scheme)
384 {
385 case INT_SCHEME_Pyramid_1:
386 GLData = &(__GaussLegendre_1[0][0]);
387 GJData = &(__GaussJacobi_01A2B0_1[0][0]);
388 GLSize = 1;
389 break;
390 case INT_SCHEME_Pyramid_8:
391 GLData = &(__GaussLegendre_2[0][0]);
392 GJData = &(__GaussJacobi_01A2B0_2[0][0]);
393 GLSize = 2;
394 break;
395 case INT_SCHEME_Pyramid_27:
396 GLData = &(__GaussLegendre_3[0][0]);
397 GJData = &(__GaussJacobi_01A2B0_3[0][0]);
398 GLSize = 3;
399 break;
400 case INT_SCHEME_Pyramid_64:
401 GLData = &(__GaussLegendre_4[0][0]);
402 GJData = &(__GaussJacobi_01A2B0_4[0][0]);
403 GLSize = 4;
404 break;
405 default:
406 w = 1e100;
407 return;
408 }
409 int iGi = iG % GLSize;
410 int iGj = (iG / GLSize) % GLSize;
411 int iGk = (iG / (GLSize * GLSize));
412
413 pParam[0] = GLData[iGi] * (1 - GJData[iGk]);
414 pParam[1] = GLData[iGj] * (1 - GJData[iGk]);
415 pParam[2] = GJData[iGk];
416 w = GLData[GLSize + iGi] * GLData[GLSize + iGj] * GJData[GLSize + iGk];
417 return;
418 }
419
420 // -----------------------------------------------------------------
421 // Prism Space - Product of Triangle and Line
422 // -----------------------------------------------------------------
423 if (ps == PrismSpace)
424 {
425 const t_real *GLData = nullptr;
426 const t_real *HammerData = nullptr;
427 int GLSize = 0;
428 int HammerSize = 0;
429 switch (scheme)
430 {
431 case INT_SCHEME_Prism_1:
432 GLData = &(__GaussLegendre_1[0][0]);
433 GLSize = 1;
434 HammerData = &(__HammerTri_1[0][0]);
435 HammerSize = 1;
436 break;
437 case INT_SCHEME_Prism_6:
438 GLData = &(__GaussLegendre_2[0][0]);
439 GLSize = 2;
440 HammerData = &(__HammerTri_3[0][0]);
441 HammerSize = 3;
442 break;
443 case INT_SCHEME_Prism_18:
444 GLData = &(__GaussLegendre_3[0][0]);
445 GLSize = 3;
446 HammerData = &(__HammerTri_6[0][0]);
447 HammerSize = 6;
448 break;
449 case INT_SCHEME_Prism_21:
450 GLData = &(__GaussLegendre_3[0][0]);
451 GLSize = 3;
452 HammerData = &(__HammerTri_7[0][0]);
453 HammerSize = 7;
454 break;
455 case INT_SCHEME_Prism_48:
456 GLData = &(__GaussLegendre_4[0][0]);
457 GLSize = 4;
458 HammerData = &(__HammerTri_12[0][0]);
459 HammerSize = 12;
460 break;
461 default:
462 w = 1e100;
463 return;
464 }
465 int iGi = iG % GLSize;
466 int iGj = iG / GLSize;
467
468 pParam[0] = HammerData[0 * HammerSize + iGj];
469 pParam[1] = HammerData[1 * HammerSize + iGj];
470 pParam[2] = GLData[iGi];
471
472 w = GLData[GLSize + iGi] * HammerData[2 * HammerSize + iGj];
473 return;
474 }
475
476 w = 1e100;
477 }
478
479} // namespace DNDS::Geom::Elem
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:108
DNDS::Geom::Elem::GetQuadratureScheme
constexpr t_index GetQuadratureScheme(ParamSpace ps, int int_order)
Definition QuadratureHub.hpp:37
DNDS::Geom::Elem::GetQuadraturePoint
void GetQuadraturePoint(ParamSpace ps, t_index scheme, int iG, TPoint &pParam, t_real &w)
Definition QuadratureHub.hpp:186
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::t_real
double t_real
Definition Geometric.hpp:8