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