DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Quad9.hpp
Go to the documentation of this file.
1#pragma once
2// Auto-generated by tools/gen_shape_functions -- DO NOT EDIT
3// Element: Quad9
4// Regenerate: /usr/bin/python3 -m tools.gen_shape_functions.generate
5
6#include "DNDS/Defines.hpp"
7#include "Geom/Geometric.hpp"
8#include "Geom/ElemEnum.hpp"
9#include "Geom/ElementTraitsBase.hpp"
10
11namespace DNDS::Geom::Elem
12{
13
14 // Forward declaration (primary template is in ElementTraitsBase.hpp)
15 template <ElemType> struct ShapeFuncImpl;
16
17 // <GEN_SHAPE_FUNCS_BEGIN>
18 template <>
19 struct ShapeFuncImpl<Quad9>
20 {
21 template <class TPoint, class TArray>
22 DNDS_DEVICE_CALLABLE static inline void Diff0(const TPoint &p, TArray &&v)
23 {
24 t_real xi = p[0];
25 t_real et = p[1];
26 const t_real _t0 = xi*(xi - 1);
27 const t_real _t1 = et*(et - 1);
28 const t_real _t2 = ((0.25))*_t1;
29 const t_real _t3 = xi*(xi + 1);
30 const t_real _t4 = et*(et + 1);
31 const t_real _t5 = ((0.25))*_t4;
32 const t_real _t6 = ((xi) * (xi)) - 1;
33 const t_real _t7 = ((0.5))*_t6;
34 const t_real _t8 = ((et) * (et)) - 1;
35 const t_real _t9 = ((0.5))*_t8;
36 v(0, 0) = _t0*_t2;
37 v(0, 1) = _t2*_t3;
38 v(0, 2) = _t3*_t5;
39 v(0, 3) = _t0*_t5;
40 v(0, 4) = -_t1*_t7;
41 v(0, 5) = -_t3*_t9;
42 v(0, 6) = -_t4*_t7;
43 v(0, 7) = -_t0*_t9;
44 v(0, 8) = _t6*_t8;
45 }
46
47 template <class TPoint, class TArray>
48 DNDS_DEVICE_CALLABLE static inline void Diff1(const TPoint &p, TArray &&v)
49 {
50 t_real xi = p[0];
51 t_real et = p[1];
52 const t_real _t0 = 2*xi;
53 const t_real _t1 = _t0 - 1;
54 const t_real _t2 = et - 1;
55 const t_real _t3 = ((0.25))*et;
56 const t_real _t4 = _t2*_t3;
57 const t_real _t5 = _t0 + 1;
58 const t_real _t6 = et + 1;
59 const t_real _t7 = _t3*_t6;
60 const t_real _t8 = et*xi;
61 const t_real _t9 = ((et) * (et)) - 1;
62 const t_real _t10 = ((0.5))*_t9;
63 const t_real _t11 = xi - 1;
64 const t_real _t12 = 2*et;
65 const t_real _t13 = _t12 - 1;
66 const t_real _t14 = ((0.25))*xi;
67 const t_real _t15 = _t13*_t14;
68 const t_real _t16 = xi + 1;
69 const t_real _t17 = _t12 + 1;
70 const t_real _t18 = _t14*_t17;
71 const t_real _t19 = ((xi) * (xi)) - 1;
72 const t_real _t20 = ((0.5))*_t19;
73 v(0, 0) = _t1*_t4;
74 v(0, 1) = _t4*_t5;
75 v(0, 2) = _t5*_t7;
76 v(0, 3) = _t1*_t7;
77 v(0, 4) = -_t2*_t8;
78 v(0, 5) = -_t10*_t5;
79 v(0, 6) = -_t6*_t8;
80 v(0, 7) = -_t1*_t10;
81 v(0, 8) = _t0*_t9;
82 v(1, 0) = _t11*_t15;
83 v(1, 1) = _t15*_t16;
84 v(1, 2) = _t16*_t18;
85 v(1, 3) = _t11*_t18;
86 v(1, 4) = -_t13*_t20;
87 v(1, 5) = -_t16*_t8;
88 v(1, 6) = -_t17*_t20;
89 v(1, 7) = -_t11*_t8;
90 v(1, 8) = _t12*_t19;
91 }
92
93 template <class TPoint, class TArray>
94 DNDS_DEVICE_CALLABLE static inline void Diff2(const TPoint &p, TArray &&v)
95 {
96 t_real xi = p[0];
97 t_real et = p[1];
98 const t_real _t0 = et - 1;
99 const t_real _t1 = ((0.5))*et;
100 const t_real _t2 = _t0*_t1;
101 const t_real _t3 = et*(et + 1);
102 const t_real _t4 = ((0.5))*_t3;
103 const t_real _t5 = ((et) * (et));
104 const t_real _t6 = 1 - _t5;
105 const t_real _t7 = -_t1;
106 const t_real _t8 = ((0.5))*xi;
107 const t_real _t9 = et*xi;
108 const t_real _t10 = -_t8 + _t9;
109 const t_real _t11 = _t8 + _t9;
110 const t_real _t12 = 2*et;
111 const t_real _t13 = 2*xi;
112 const t_real _t14 = xi - 1;
113 const t_real _t15 = _t14*_t8;
114 const t_real _t16 = xi + 1;
115 const t_real _t17 = _t16*_t8;
116 const t_real _t18 = ((xi) * (xi));
117 const t_real _t19 = 1 - _t18;
118 v(0, 0) = _t2;
119 v(0, 1) = _t2;
120 v(0, 2) = _t4;
121 v(0, 3) = _t4;
122 v(0, 4) = -_t0*et;
123 v(0, 5) = _t6;
124 v(0, 6) = -_t3;
125 v(0, 7) = _t6;
126 v(0, 8) = 2*_t5 - 2;
127 v(1, 0) = _t10 + _t7 + (0.25);
128 v(1, 1) = _t1 + _t10 + (-0.25);
129 v(1, 2) = _t1 + _t11 + (0.25);
130 v(1, 3) = _t11 + _t7 + (-0.25);
131 v(1, 4) = xi*(1 - _t12);
132 v(1, 5) = et*(-_t13 - 1);
133 v(1, 6) = xi*(-_t12 - 1);
134 v(1, 7) = et*(1 - _t13);
135 v(1, 8) = 4*_t9;
136 v(2, 0) = _t15;
137 v(2, 1) = _t17;
138 v(2, 2) = _t17;
139 v(2, 3) = _t15;
140 v(2, 4) = _t19;
141 v(2, 5) = -_t16*xi;
142 v(2, 6) = _t19;
143 v(2, 7) = -_t14*xi;
144 v(2, 8) = 2*_t18 - 2;
145 }
146
147 template <class TPoint, class TArray>
148 DNDS_DEVICE_CALLABLE static inline void Diff3(const TPoint &p, TArray &&v)
149 {
150 t_real xi = p[0];
151 t_real et = p[1];
152 const t_real _t0 = et + (-0.5);
153 const t_real _t1 = et + (0.5);
154 const t_real _t2 = 2*et;
155 const t_real _t3 = -_t2;
156 const t_real _t4 = xi + (-0.5);
157 const t_real _t5 = xi + (0.5);
158 const t_real _t6 = 2*xi;
159 const t_real _t7 = -_t6;
160 v(1, 0) = _t0;
161 v(1, 1) = _t0;
162 v(1, 2) = _t1;
163 v(1, 3) = _t1;
164 v(1, 4) = 1 - _t2;
165 v(1, 5) = _t3;
166 v(1, 6) = -_t2 - 1;
167 v(1, 7) = _t3;
168 v(1, 8) = 4*et;
169 v(2, 0) = _t4;
170 v(2, 1) = _t5;
171 v(2, 2) = _t5;
172 v(2, 3) = _t4;
173 v(2, 4) = _t7;
174 v(2, 5) = -_t6 - 1;
175 v(2, 6) = _t7;
176 v(2, 7) = 1 - _t6;
177 v(2, 8) = 4*xi;
178 }
179 };
180 // <GEN_SHAPE_FUNCS_END>
181
182
183 /**
184 * @brief Element traits for 9-node biquadratic quadrilateral (Quad9)
185 *
186 * Quad9 is a high-order 2D quadrilateral element with:
187 * - 4 corner nodes (vertices)
188 * - 4 edge mid-nodes
189 * - 1 interior node at the center
190 *
191 * Used for high-order finite element and spectral element methods.
192 */
193 template <>
194 struct ElementTraits<Quad9>
195 {
196 // ============================================================
197 // Core Element Identification
198 // ============================================================
199
200 static constexpr ElemType elemType = Quad9;
201 static constexpr int dim = 2;
202 static constexpr int order = 2;
203 static constexpr int numVertices = 4;
204 static constexpr int numNodes = 9;
205 static constexpr int numFaces = 4;
206 static constexpr ParamSpace paramSpace = QuadSpace;
207 static constexpr t_real paramSpaceVol = 4.0;
208
209 // ============================================================
210 // Geometry Definition
211 // ============================================================
212
213 /**
214 * @brief Standard coordinates of nodes in parametric space
215 *
216 * Reference square [-1,1] x [-1,1] with nodes:
217 * Nodes 0-3: corners (same as Quad4)
218 * Nodes 4-7: edge midpoints
219 * Node 8: center (0, 0)
220 */
221 static constexpr std::array<t_real, 3 * 9> standardCoords = {
222 -1, -1, 0, // Node 0: bottom-left corner
223 1, -1, 0, // Node 1: bottom-right corner
224 1, 1, 0, // Node 2: top-right corner
225 -1, 1, 0, // Node 3: top-left corner
226 0, -1, 0, // Node 4: bottom edge midpoint
227 1, 0, 0, // Node 5: right edge midpoint
228 0, 1, 0, // Node 6: top edge midpoint
229 -1, 0, 0, // Node 7: left edge midpoint
230 0, 0, 0}; // Node 8: center
231
232 // ============================================================
233 // Face/Edge Definitions
234 // ============================================================
235
236 /**
237 * @brief Get the element type of a face (edge)
238 * @return Line3 (all edges are quadratic lines with 3 nodes)
239 */
240 static constexpr ElemType GetFaceType(t_index /*iFace*/) { return Line3; }
241
242 /**
243 * @brief Node indices for each face (edge)
244 *
245 * Each edge has 3 nodes: 2 vertices + 1 midpoint
246 * Edge 0: nodes 0-1-4 (bottom)
247 * Edge 1: nodes 1-2-5 (right)
248 * Edge 2: nodes 2-3-6 (top)
249 * Edge 3: nodes 3-0-7 (left)
250 */
251 static constexpr std::array<std::array<t_index, 10>, 4> faceNodes = {{
252 {0, 1, 4}, // Edge 0: bottom
253 {1, 2, 5}, // Edge 1: right
254 {2, 3, 6}, // Edge 2: top
255 {3, 0, 7}}}; // Edge 3: left
256
257 // ============================================================
258 // Order Elevation (P-Refinement)
259 // ============================================================
260
261 /// @brief Element type after order elevation (O2 has no higher elevation defined)
262 static constexpr ElemType elevatedType = UnknownElem;
263
264 /// @brief Number of additional nodes created during elevation (none for O2)
265 static constexpr int numElevNodes = 0;
266
267 // ============================================================
268 // Bisection (Adaptive Refinement)
269 // ============================================================
270
271 /// @brief Number of sub-elements created when bisecting (4 Quad4 elements)
272 static constexpr int numBisect = 4;
273
274 /// @brief Number of bisection variants (only 1 way to uniformly bisect)
275 static constexpr int numBisectVariants = 1;
276
277 /**
278 * @brief Get the element type of a sub-element after bisection
279 * @return Quad4 (all sub-elements are bilinear quads)
280 */
281 static constexpr ElemType GetBisectElemType(t_index /*i*/) { return Quad4; }
282
283 /**
284 * @brief Node indices for each sub-element created by bisection
285 *
286 * Bisecting creates 4 sub-quads meeting at the center node 8:
287 * Sub-element 0: bottom-left quad
288 * Sub-element 1: bottom-right quad
289 * Sub-element 2: top-left quad
290 * Sub-element 3: top-right quad
291 */
292 static constexpr std::array<tBisectSub, 4> bisectElements = {{
293 {0, 4, 8, 7}, // Bottom-left sub-quad
294 {4, 1, 5, 8}, // Bottom-right sub-quad
295 {7, 8, 6, 3}, // Top-left sub-quad
296 {8, 5, 2, 6}}}; // Top-right sub-quad
297
298 // ============================================================
299 // VTK/Visualization Support
300 // ============================================================
301
302 /// @brief VTK cell type identifier (23 = VTK_QUADRATIC_QUAD)
303 static constexpr int vtkCellType = 23;
304
305 /**
306 * @brief VTK node ordering map
307 *
308 * VTK uses 8 nodes for quadratic quad (excludes interior node 8):
309 * VTK nodes 0-3 = corner nodes 0-3
310 * VTK nodes 4-7 = edge mid-nodes 4-7
311 */
312 static constexpr std::array<int, 8> vtkNodeOrder = {0, 1, 2, 3, 4, 5, 6, 7};
313 };
314
315} // namespace DNDS::Geom::Elem
Defines.hpp
Core type aliases, constants, and metaprogramming utilities for the DNDS framework.
DNDS_DEVICE_CALLABLE
#define DNDS_DEVICE_CALLABLE
Definition Defines.hpp:76
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::t_real
double t_real
Definition Geometric.hpp:8
DNDS::Geom::Elem::ElementTraits< Quad9 >::GetBisectElemType
static constexpr ElemType GetBisectElemType(t_index)
Get the element type of a sub-element after bisection.
Definition Quad9.hpp:281
DNDS::Geom::Elem::ElementTraits< Quad9 >::GetFaceType
static constexpr ElemType GetFaceType(t_index)
Get the element type of a face (edge)
Definition Quad9.hpp:240
DNDS::Geom::Elem::ShapeFuncImpl< Quad9 >::Diff3
static DNDS_DEVICE_CALLABLE void Diff3(const TPoint &p, TArray &&v)
Definition Quad9.hpp:148
DNDS::Geom::Elem::ShapeFuncImpl< Quad9 >::Diff1
static DNDS_DEVICE_CALLABLE void Diff1(const TPoint &p, TArray &&v)
Definition Quad9.hpp:48
DNDS::Geom::Elem::ShapeFuncImpl< Quad9 >::Diff0
static DNDS_DEVICE_CALLABLE void Diff0(const TPoint &p, TArray &&v)
Definition Quad9.hpp:22
DNDS::Geom::Elem::ShapeFuncImpl< Quad9 >::Diff2
static DNDS_DEVICE_CALLABLE void Diff2(const TPoint &p, TArray &&v)
Definition Quad9.hpp:94
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
order
double order
Definition test_ODE.cpp:257