PhysicsBasedAnimationToolkit 0.0.10
Cross-platform C++20 library of algorithms and data structures commonly used in computer graphics research on physically-based simulation.
Loading...
Searching...
No Matches
Tetrahedron.h
Go to the documentation of this file.
1
11#ifndef PBAT_FEM_TETRAHEDRON_H
12#define PBAT_FEM_TETRAHEDRON_H
13
14#include "pbat/Aliases.h"
15#include "pbat/common/Concepts.h"
16#include "QuadratureRules.h"
17
18#include <array>
19
20namespace pbat {
21namespace fem {
22
23namespace detail {
24
25template <int Order>
26struct Tetrahedron;
27
28} // namespace detail
29
37template <int Order>
38using Tetrahedron = typename detail::Tetrahedron<Order>;
39
40namespace detail {
41
42template <>
43struct Tetrahedron<1>
44{
45 using AffineBaseType = Tetrahedron<1>;
46
47 static int constexpr kOrder = 1;
48 static int constexpr kDims = 3;
49 static int constexpr kNodes = 4;
50 static std::array<int, kNodes * kDims> constexpr Coordinates =
51 {0,0,0,1,0,0,0,1,0,0,0,1};
52 static std::array<int, AffineBaseType::kNodes> constexpr Vertices = {0,1,2,3};
53 static bool constexpr bHasConstantJacobian = true;
54
55 template <int PolynomialOrder, common::CFloatingPoint TScalar = Scalar>
56 using QuadratureType = math::SymmetricSimplexPolynomialQuadratureRule<kDims, PolynomialOrder, TScalar>;
57
58 template <class TDerived, class TScalar = typename TDerived::Scalar>
59 [[maybe_unused]] static Eigen::Vector<TScalar, kNodes> N([[maybe_unused]] Eigen::DenseBase<TDerived> const& X_)
60 {
61#include "pbat/warning/Push.h"
62#include "pbat/warning/SignConversion.h"
63 using namespace pbat::math;
64 Eigen::Vector<TScalar, kNodes> Nm;
65 auto const X = X_.reshaped();
66 Nm[0] = -X[0] - X[1] - X[2] + 1;
67 Nm[1] = X[0];
68 Nm[2] = X[1];
69 Nm[3] = X[2];
70 return Nm;
71#include "pbat/warning/Pop.h"
72 }
73
74 template <class TDerived, class TScalar = typename TDerived::Scalar>
75 [[maybe_unused]] static Eigen::Matrix<TScalar, kNodes, kDims> GradN([[maybe_unused]] Eigen::DenseBase<TDerived> const& X_)
76 {
77#include "pbat/warning/Push.h"
78#include "pbat/warning/SignConversion.h"
79 Eigen::Matrix<TScalar, kNodes, kDims> GNm;
80 TScalar* GNp = GNm.data();
81 [[maybe_unused]] auto const X = X_.reshaped();
82 GNp[0] = -1;
83 GNp[1] = 1;
84 GNp[2] = 0;
85 GNp[3] = 0;
86 GNp[4] = -1;
87 GNp[5] = 0;
88 GNp[6] = 1;
89 GNp[7] = 0;
90 GNp[8] = -1;
91 GNp[9] = 0;
92 GNp[10] = 0;
93 GNp[11] = 1;
94 return GNm;
95#include "pbat/warning/Pop.h"
96 }
97};
98
99template <>
100struct Tetrahedron<2>
101{
102 using AffineBaseType = Tetrahedron<1>;
103
104 static int constexpr kOrder = 2;
105 static int constexpr kDims = 3;
106 static int constexpr kNodes = 10;
107 static std::array<int, kNodes * kDims> constexpr Coordinates =
108 {0,0,0,1,0,0,2,0,0,0,1,0,1,1,0,0,2,0,0,0,1,1,0,1,0,1,1,0,0,2};
109 static std::array<int, AffineBaseType::kNodes> constexpr Vertices = {0,2,5,9};
110 static bool constexpr bHasConstantJacobian = false;
111
112 template <int PolynomialOrder, common::CFloatingPoint TScalar = Scalar>
113 using QuadratureType = math::SymmetricSimplexPolynomialQuadratureRule<kDims, PolynomialOrder, TScalar>;
114
115 template <class TDerived, class TScalar = typename TDerived::Scalar>
116 [[maybe_unused]] static Eigen::Vector<TScalar, kNodes> N([[maybe_unused]] Eigen::DenseBase<TDerived> const& X_)
117 {
118#include "pbat/warning/Push.h"
119#include "pbat/warning/SignConversion.h"
120 using namespace pbat::math;
121 Eigen::Vector<TScalar, kNodes> Nm;
122 auto const X = X_.reshaped();
123 auto const a0 = X[0] + X[1] + X[2] - 1;
124 auto const a1 = 2*X[1];
125 auto const a2 = 2*X[2];
126 auto const a3 = 2*X[0] - 1;
127 auto const a4 = 4*a0;
128 auto const a5 = 4*X[0];
129 Nm[0] = a0*(a1 + a2 + a3);
130 Nm[1] = -a4*X[0];
131 Nm[2] = a3*X[0];
132 Nm[3] = -a4*X[1];
133 Nm[4] = a5*X[1];
134 Nm[5] = (a1 - 1)*X[1];
135 Nm[6] = -a4*X[2];
136 Nm[7] = a5*X[2];
137 Nm[8] = 4*X[1]*X[2];
138 Nm[9] = (a2 - 1)*X[2];
139 return Nm;
140#include "pbat/warning/Pop.h"
141 }
142
143 template <class TDerived, class TScalar = typename TDerived::Scalar>
144 [[maybe_unused]] static Eigen::Matrix<TScalar, kNodes, kDims> GradN([[maybe_unused]] Eigen::DenseBase<TDerived> const& X_)
145 {
146#include "pbat/warning/Push.h"
147#include "pbat/warning/SignConversion.h"
148 Eigen::Matrix<TScalar, kNodes, kDims> GNm;
149 TScalar* GNp = GNm.data();
150 [[maybe_unused]] auto const X = X_.reshaped();
151 auto const a0 = 4*X[0];
152 auto const a1 = 4*X[1];
153 auto const a2 = 4*X[2];
154 auto const a3 = a1 + a2;
155 auto const a4 = a0 + a3 - 3;
156 auto const a5 = -a1;
157 auto const a6 = -a2;
158 auto const a7 = -a0;
159 auto const a8 = a0 - 4;
160 GNp[0] = a4;
161 GNp[1] = -a3 - 8*X[0] + 4;
162 GNp[2] = a0 - 1;
163 GNp[3] = a5;
164 GNp[4] = a1;
165 GNp[5] = 0;
166 GNp[6] = a6;
167 GNp[7] = a2;
168 GNp[8] = 0;
169 GNp[9] = 0;
170 GNp[10] = a4;
171 GNp[11] = a7;
172 GNp[12] = 0;
173 GNp[13] = -a2 - a8 - 8*X[1];
174 GNp[14] = a0;
175 GNp[15] = a1 - 1;
176 GNp[16] = a6;
177 GNp[17] = 0;
178 GNp[18] = a2;
179 GNp[19] = 0;
180 GNp[20] = a4;
181 GNp[21] = a7;
182 GNp[22] = 0;
183 GNp[23] = a5;
184 GNp[24] = 0;
185 GNp[25] = 0;
186 GNp[26] = -a1 - a8 - 8*X[2];
187 GNp[27] = a0;
188 GNp[28] = a1;
189 GNp[29] = a2 - 1;
190 return GNm;
191#include "pbat/warning/Pop.h"
192 }
193};
194
195template <>
196struct Tetrahedron<3>
197{
198 using AffineBaseType = Tetrahedron<1>;
199
200 static int constexpr kOrder = 3;
201 static int constexpr kDims = 3;
202 static int constexpr kNodes = 20;
203 static std::array<int, kNodes * kDims> constexpr Coordinates =
204 {0,0,0,1,0,0,2,0,0,3,0,0,0,1,0,1,1,0,2,1,0,0,2,0,1,2,0,0,3,0,0,0,1,1,0,1,2,0,1,0,1,1,1,1,1,0,2,1,0,0,2,1,0,2,0,1,2,0,0,3};
205 static std::array<int, AffineBaseType::kNodes> constexpr Vertices = {0,3,9,19};
206 static bool constexpr bHasConstantJacobian = false;
207
208 template <int PolynomialOrder, common::CFloatingPoint TScalar = Scalar>
209 using QuadratureType = math::SymmetricSimplexPolynomialQuadratureRule<kDims, PolynomialOrder, TScalar>;
210
211 template <class TDerived, class TScalar = typename TDerived::Scalar>
212 [[maybe_unused]] static Eigen::Vector<TScalar, kNodes> N([[maybe_unused]] Eigen::DenseBase<TDerived> const& X_)
213 {
214#include "pbat/warning/Push.h"
215#include "pbat/warning/SignConversion.h"
216 using namespace pbat::math;
217 Eigen::Vector<TScalar, kNodes> Nm;
218 auto const X = X_.reshaped();
219 auto const a0 = 3*X[0];
220 auto const a1 = a0 - 1;
221 auto const a2 = 3*X[1];
222 auto const a3 = 3*X[2];
223 auto const a4 = a2 + a3;
224 auto const a5 = X[0] + X[1] + X[2] - 1;
225 auto const a6 = a0 - 2;
226 auto const a7 = a5*(a4 + a6);
227 auto const a8 = (9.0/2.0)*X[0];
228 auto const a9 = a1*a8;
229 auto const a10 = (9.0/2.0)*X[1];
230 auto const a11 = 27*a5*X[0];
231 auto const a12 = a2 - 1;
232 auto const a13 = a10*a12;
233 auto const a14 = a12*X[1];
234 auto const a15 = (9.0/2.0)*X[2];
235 auto const a16 = 27*X[1]*X[2];
236 auto const a17 = a3 - 1;
237 auto const a18 = a17*X[2];
238 Nm[0] = -1.0/2.0*a7*(a1 + a4);
239 Nm[1] = a7*a8;
240 Nm[2] = -a5*a9;
241 Nm[3] = (1.0/2.0)*a1*a6*X[0];
242 Nm[4] = a10*a7;
243 Nm[5] = -a11*X[1];
244 Nm[6] = a9*X[1];
245 Nm[7] = -a13*a5;
246 Nm[8] = a14*a8;
247 Nm[9] = (1.0/2.0)*a14*(a2 - 2);
248 Nm[10] = a15*a7;
249 Nm[11] = -a11*X[2];
250 Nm[12] = a9*X[2];
251 Nm[13] = -a16*a5;
252 Nm[14] = a16*X[0];
253 Nm[15] = a13*X[2];
254 Nm[16] = -a15*a17*a5;
255 Nm[17] = a18*a8;
256 Nm[18] = a10*a18;
257 Nm[19] = (1.0/2.0)*a18*(a3 - 2);
258 return Nm;
259#include "pbat/warning/Pop.h"
260 }
261
262 template <class TDerived, class TScalar = typename TDerived::Scalar>
263 [[maybe_unused]] static Eigen::Matrix<TScalar, kNodes, kDims> GradN([[maybe_unused]] Eigen::DenseBase<TDerived> const& X_)
264 {
265#include "pbat/warning/Push.h"
266#include "pbat/warning/SignConversion.h"
267 Eigen::Matrix<TScalar, kNodes, kDims> GNm;
268 TScalar* GNp = GNm.data();
269 [[maybe_unused]] auto const X = X_.reshaped();
270 auto const a0 = X[0] + X[1] + X[2] - 1;
271 auto const a1 = (3.0/2.0)*X[0];
272 auto const a2 = a1 - 1.0/2.0;
273 auto const a3 = (3.0/2.0)*X[1];
274 auto const a4 = (3.0/2.0)*X[2];
275 auto const a5 = a3 + a4;
276 auto const a6 = -a2 - a5;
277 auto const a7 = 3*X[1];
278 auto const a8 = 3*X[2];
279 auto const a9 = 3*X[0] - 2;
280 auto const a10 = a7 + a8 + a9;
281 auto const a11 = 3*a0*a6 + a10*a6 + a10*(-a1 - a5 + 3.0/2.0);
282 auto const a12 = (9.0/2.0)*X[0];
283 auto const a13 = (9.0/2.0)*X[1];
284 auto const a14 = (9.0/2.0)*X[2];
285 auto const a15 = a10*(a12 + a13 + a14 - 9.0/2.0);
286 auto const a16 = (27.0/2.0)*X[0];
287 auto const a17 = (27.0/2.0)*X[1];
288 auto const a18 = (27.0/2.0)*X[2];
289 auto const a19 = a16 + a17 + a18;
290 auto const a20 = a19 - 27.0/2.0;
291 auto const a21 = a19 - 9;
292 auto const a22 = a20*X[0] + a21*X[0];
293 auto const a23 = a16 - 9.0/2.0;
294 auto const a24 = -a23;
295 auto const a25 = a24*X[0];
296 auto const a26 = -a20;
297 auto const a27 = a20*X[1] + a21*X[1];
298 auto const a28 = 27*X[0];
299 auto const a29 = a28*X[1];
300 auto const a30 = -a29;
301 auto const a31 = 27*X[1];
302 auto const a32 = -a28 - a31 - 27*X[2] + 27;
303 auto const a33 = a32*X[1];
304 auto const a34 = a16*X[1];
305 auto const a35 = a17 - 9.0/2.0;
306 auto const a36 = -a35;
307 auto const a37 = a36*X[1];
308 auto const a38 = a35*X[1];
309 auto const a39 = a20*X[2] + a21*X[2];
310 auto const a40 = a28*X[2];
311 auto const a41 = -a40;
312 auto const a42 = a32*X[2];
313 auto const a43 = a16*X[2];
314 auto const a44 = a31*X[2];
315 auto const a45 = -a44;
316 auto const a46 = a18 - 9.0/2.0;
317 auto const a47 = -a46;
318 auto const a48 = a47*X[2];
319 auto const a49 = a46*X[2];
320 auto const a50 = a32*X[0];
321 auto const a51 = a23*X[0];
322 auto const a52 = a17*X[2];
323 GNp[0] = a11;
324 GNp[1] = a15 + a22;
325 GNp[2] = a0*a24 + a25 + a26*X[0];
326 GNp[3] = a2*a9 + (a12 - 3)*X[0] + (a12 - 3.0/2.0)*X[0];
327 GNp[4] = a27;
328 GNp[5] = a30 + a33;
329 GNp[6] = a23*X[1] + a34;
330 GNp[7] = a37;
331 GNp[8] = a38;
332 GNp[9] = 0;
333 GNp[10] = a39;
334 GNp[11] = a41 + a42;
335 GNp[12] = a23*X[2] + a43;
336 GNp[13] = a45;
337 GNp[14] = a44;
338 GNp[15] = 0;
339 GNp[16] = a48;
340 GNp[17] = a49;
341 GNp[18] = 0;
342 GNp[19] = 0;
343 GNp[20] = a11;
344 GNp[21] = a22;
345 GNp[22] = a25;
346 GNp[23] = 0;
347 GNp[24] = a15 + a27;
348 GNp[25] = a30 + a50;
349 GNp[26] = a51;
350 GNp[27] = a0*a36 + a26*X[1] + a37;
351 GNp[28] = a34 + a35*X[0];
352 GNp[29] = (a13 - 3)*X[1] + (a13 - 3.0/2.0)*X[1] + (a3 - 1.0/2.0)*(a7 - 2);
353 GNp[30] = a39;
354 GNp[31] = a41;
355 GNp[32] = 0;
356 GNp[33] = a42 + a45;
357 GNp[34] = a40;
358 GNp[35] = a35*X[2] + a52;
359 GNp[36] = a48;
360 GNp[37] = 0;
361 GNp[38] = a49;
362 GNp[39] = 0;
363 GNp[40] = a11;
364 GNp[41] = a22;
365 GNp[42] = a25;
366 GNp[43] = 0;
367 GNp[44] = a27;
368 GNp[45] = a30;
369 GNp[46] = 0;
370 GNp[47] = a37;
371 GNp[48] = 0;
372 GNp[49] = 0;
373 GNp[50] = a15 + a39;
374 GNp[51] = a41 + a50;
375 GNp[52] = a51;
376 GNp[53] = a33 + a45;
377 GNp[54] = a29;
378 GNp[55] = a38;
379 GNp[56] = a0*a47 + a26*X[2] + a48;
380 GNp[57] = a43 + a46*X[0];
381 GNp[58] = a46*X[1] + a52;
382 GNp[59] = (a14 - 3)*X[2] + (a14 - 3.0/2.0)*X[2] + (a4 - 1.0/2.0)*(a8 - 2);
383 return GNm;
384#include "pbat/warning/Pop.h"
385 }
386};
387
388} // namespace detail
389} // namespace fem
390} // namespace pbat
391
392#endif // PBAT_FEM_TETRAHEDRON_H
Concepts.h
Concepts for common types.
pbat::fem
Finite Element Method (FEM)
Definition Concepts.h:19
pbat::fem::Tetrahedron
typename detail::Tetrahedron< Order > Tetrahedron
Tetrahedron finite element.
Definition Tetrahedron.h:38
pbat::math
Math related functionality.
Definition Concepts.h:19
pbat::math::SymmetricSimplexPolynomialQuadratureRule
typename detail::SymmetricSimplexPolynomialQuadratureRule< Dims, Order, TScalar > SymmetricSimplexPolynomialQuadratureRule
Symmetric quadrature rule for reference simplices in dimensions .
Definition SymmetricQuadratureRules.h:96
pbat
The main namespace of the library.
Definition Aliases.h:15