GooseFEM 1.4.1.dev2+g78f16df
Loading...
Searching...
No Matches
ElementHex8.h
Go to the documentation of this file.
1
10#ifndef GOOSEFEM_ELEMENTHEX8_H
11#define GOOSEFEM_ELEMENTHEX8_H
12
13#include "config.h"
14
15namespace GooseFEM {
16namespace Element {
17
21namespace Hex8 {
22
26namespace Gauss {
27
35inline size_t nip()
36{
37 return 8;
38}
39
45inline array_type::tensor<double, 2> xi()
46{
47 size_t nip = 8;
48 size_t ndim = 3;
49
50 array_type::tensor<double, 2> xi = xt::empty<double>({nip, ndim});
51
52 xi(0, 0) = -1.0 / std::sqrt(3.0);
53 xi(0, 1) = -1.0 / std::sqrt(3.0);
54 xi(0, 2) = -1.0 / std::sqrt(3.0);
55
56 xi(1, 0) = +1.0 / std::sqrt(3.0);
57 xi(1, 1) = -1.0 / std::sqrt(3.0);
58 xi(1, 2) = -1.0 / std::sqrt(3.0);
59
60 xi(2, 0) = +1.0 / std::sqrt(3.0);
61 xi(2, 1) = +1.0 / std::sqrt(3.0);
62 xi(2, 2) = -1.0 / std::sqrt(3.0);
63
64 xi(3, 0) = -1.0 / std::sqrt(3.0);
65 xi(3, 1) = +1.0 / std::sqrt(3.0);
66 xi(3, 2) = -1.0 / std::sqrt(3.0);
67
68 xi(4, 0) = -1.0 / std::sqrt(3.0);
69 xi(4, 1) = -1.0 / std::sqrt(3.0);
70 xi(4, 2) = +1.0 / std::sqrt(3.0);
71
72 xi(5, 0) = +1.0 / std::sqrt(3.0);
73 xi(5, 1) = -1.0 / std::sqrt(3.0);
74 xi(5, 2) = +1.0 / std::sqrt(3.0);
75
76 xi(6, 0) = +1.0 / std::sqrt(3.0);
77 xi(6, 1) = +1.0 / std::sqrt(3.0);
78 xi(6, 2) = +1.0 / std::sqrt(3.0);
79
80 xi(7, 0) = -1.0 / std::sqrt(3.0);
81 xi(7, 1) = +1.0 / std::sqrt(3.0);
82 xi(7, 2) = +1.0 / std::sqrt(3.0);
83
84 return xi;
85}
86
92inline array_type::tensor<double, 1> w()
93{
94 size_t nip = 8;
95
96 array_type::tensor<double, 1> w = xt::empty<double>({nip});
97
98 w(0) = 1.0;
99 w(1) = 1.0;
100 w(2) = 1.0;
101 w(3) = 1.0;
102 w(4) = 1.0;
103 w(5) = 1.0;
104 w(6) = 1.0;
105 w(7) = 1.0;
106
107 return w;
108}
109
110} // namespace Gauss
111
116namespace Nodal {
117
125inline size_t nip()
126{
127 return 8;
128}
129
135inline array_type::tensor<double, 2> xi()
136{
137 size_t nip = 8;
138 size_t ndim = 3;
139
140 array_type::tensor<double, 2> xi = xt::empty<double>({nip, ndim});
141
142 xi(0, 0) = -1.0;
143 xi(0, 1) = -1.0;
144 xi(0, 2) = -1.0;
145
146 xi(1, 0) = +1.0;
147 xi(1, 1) = -1.0;
148 xi(1, 2) = -1.0;
149
150 xi(2, 0) = +1.0;
151 xi(2, 1) = +1.0;
152 xi(2, 2) = -1.0;
153
154 xi(3, 0) = -1.0;
155 xi(3, 1) = +1.0;
156 xi(3, 2) = -1.0;
157
158 xi(4, 0) = -1.0;
159 xi(4, 1) = -1.0;
160 xi(4, 2) = +1.0;
161
162 xi(5, 0) = +1.0;
163 xi(5, 1) = -1.0;
164 xi(5, 2) = +1.0;
165
166 xi(6, 0) = +1.0;
167 xi(6, 1) = +1.0;
168 xi(6, 2) = +1.0;
169
170 xi(7, 0) = -1.0;
171 xi(7, 1) = +1.0;
172 xi(7, 2) = +1.0;
173
174 return xi;
175}
176
182inline array_type::tensor<double, 1> w()
183{
184 size_t nip = 8;
185
186 array_type::tensor<double, 1> w = xt::empty<double>({nip});
187
188 w(0) = 1.0;
189 w(1) = 1.0;
190 w(2) = 1.0;
191 w(3) = 1.0;
192 w(4) = 1.0;
193 w(5) = 1.0;
194 w(6) = 1.0;
195 w(7) = 1.0;
196
197 return w;
198}
199
200} // namespace Nodal
201
215class Quadrature : public QuadratureBaseCartesian<Quadrature> {
216public:
217 Quadrature() = default;
218
233 template <class T>
234 Quadrature(const T& x) : Quadrature(x, Gauss::xi(), Gauss::w())
235 {
236 }
237
254 template <class T, class X, class W>
255 Quadrature(const T& x, const X& xi, const W& w)
256 {
257 m_x = x;
258 m_w = w;
259 m_xi = xi;
260 m_nip = w.size();
261 m_nelem = m_x.shape(0);
262 m_N = xt::empty<double>({m_nip, s_nne});
263 m_dNxi = xt::empty<double>({m_nip, s_nne, s_ndim});
264
265 for (size_t q = 0; q < m_nip; ++q) {
266 m_N(q, 0) = 0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1)) * (1.0 - xi(q, 2));
267 m_N(q, 1) = 0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1)) * (1.0 - xi(q, 2));
268 m_N(q, 2) = 0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1)) * (1.0 - xi(q, 2));
269 m_N(q, 3) = 0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1)) * (1.0 - xi(q, 2));
270 m_N(q, 4) = 0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1)) * (1.0 + xi(q, 2));
271 m_N(q, 5) = 0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1)) * (1.0 + xi(q, 2));
272 m_N(q, 6) = 0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1)) * (1.0 + xi(q, 2));
273 m_N(q, 7) = 0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1)) * (1.0 + xi(q, 2));
274 }
275
276 for (size_t q = 0; q < m_nip; ++q) {
277 // - dN / dxi_0
278 m_dNxi(q, 0, 0) = -0.125 * (1.0 - xi(q, 1)) * (1.0 - xi(q, 2));
279 m_dNxi(q, 1, 0) = +0.125 * (1.0 - xi(q, 1)) * (1.0 - xi(q, 2));
280 m_dNxi(q, 2, 0) = +0.125 * (1.0 + xi(q, 1)) * (1.0 - xi(q, 2));
281 m_dNxi(q, 3, 0) = -0.125 * (1.0 + xi(q, 1)) * (1.0 - xi(q, 2));
282 m_dNxi(q, 4, 0) = -0.125 * (1.0 - xi(q, 1)) * (1.0 + xi(q, 2));
283 m_dNxi(q, 5, 0) = +0.125 * (1.0 - xi(q, 1)) * (1.0 + xi(q, 2));
284 m_dNxi(q, 6, 0) = +0.125 * (1.0 + xi(q, 1)) * (1.0 + xi(q, 2));
285 m_dNxi(q, 7, 0) = -0.125 * (1.0 + xi(q, 1)) * (1.0 + xi(q, 2));
286 // - dN / dxi_1
287 m_dNxi(q, 0, 1) = -0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 2));
288 m_dNxi(q, 1, 1) = -0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 2));
289 m_dNxi(q, 2, 1) = +0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 2));
290 m_dNxi(q, 3, 1) = +0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 2));
291 m_dNxi(q, 4, 1) = -0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 2));
292 m_dNxi(q, 5, 1) = -0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 2));
293 m_dNxi(q, 6, 1) = +0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 2));
294 m_dNxi(q, 7, 1) = +0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 2));
295 // - dN / dxi_2
296 m_dNxi(q, 0, 2) = -0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1));
297 m_dNxi(q, 1, 2) = -0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1));
298 m_dNxi(q, 2, 2) = -0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1));
299 m_dNxi(q, 3, 2) = -0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1));
300 m_dNxi(q, 4, 2) = +0.125 * (1.0 - xi(q, 0)) * (1.0 - xi(q, 1));
301 m_dNxi(q, 5, 2) = +0.125 * (1.0 + xi(q, 0)) * (1.0 - xi(q, 1));
302 m_dNxi(q, 6, 2) = +0.125 * (1.0 + xi(q, 0)) * (1.0 + xi(q, 1));
303 m_dNxi(q, 7, 2) = +0.125 * (1.0 - xi(q, 0)) * (1.0 + xi(q, 1));
304 }
305
306 GOOSEFEM_ASSERT(m_x.shape(1) == s_nne);
307 GOOSEFEM_ASSERT(m_x.shape(2) == s_ndim);
308 GOOSEFEM_ASSERT(xt::has_shape(m_xi, {m_nip, s_ndim}));
309 GOOSEFEM_ASSERT(xt::has_shape(m_w, {m_nip}));
310 GOOSEFEM_ASSERT(xt::has_shape(m_N, {m_nip, s_nne}));
311 GOOSEFEM_ASSERT(xt::has_shape(m_dNxi, {m_nip, s_nne, s_ndim}));
312
313 m_dNx = xt::empty<double>({m_nelem, m_nip, s_nne, s_ndim});
314 m_vol = xt::empty<double>(this->shape_qscalar());
315
316 this->compute_dN();
317 }
318
319private:
320 friend QuadratureBase<Quadrature>;
321 friend QuadratureBaseCartesian<Quadrature>;
322
323 template <class T, class R>
324 void int_N_scalar_NT_dV_impl(const T& qscalar, R& elemmat) const
325 {
326 GOOSEFEM_ASSERT(xt::has_shape(qscalar, this->shape_qscalar()));
327 GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat()));
328
329 elemmat.fill(0.0);
330
331#pragma omp parallel for
332 for (size_t e = 0; e < m_nelem; ++e) {
333
334 auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape<s_nne * s_ndim, s_nne * s_ndim>());
335
336 for (size_t q = 0; q < m_nip; ++q) {
337
338 auto N = xt::adapt(&m_N(q, 0), xt::xshape<s_nne>());
339 auto& vol = m_vol(e, q);
340 auto& rho = qscalar(e, q);
341
342 // M(m * ndim + i, n * ndim + i) += N(m) * scalar * N(n) * dV
343 for (size_t m = 0; m < s_nne; ++m) {
344 for (size_t n = 0; n < s_nne; ++n) {
345 M(m * s_ndim + 0, n * s_ndim + 0) += N(m) * rho * N(n) * vol;
346 M(m * s_ndim + 1, n * s_ndim + 1) += N(m) * rho * N(n) * vol;
347 M(m * s_ndim + 2, n * s_ndim + 2) += N(m) * rho * N(n) * vol;
348 }
349 }
350 }
351 }
352 }
353
354 template <class T, class R>
355 void int_gradN_dot_tensor2_dV_impl(const T& qtensor, R& elemvec) const
356 {
357 GOOSEFEM_ASSERT(xt::has_shape(qtensor, this->shape_qtensor<2>()));
358 GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec()));
359
360 elemvec.fill(0.0);
361
362#pragma omp parallel for
363 for (size_t e = 0; e < m_nelem; ++e) {
364
365 auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape<s_nne, s_ndim>());
366
367 for (size_t q = 0; q < m_nip; ++q) {
368
369 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<s_nne, s_ndim>());
370 auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<s_ndim, s_ndim>());
371 auto& v = m_vol(e, q);
372
373 for (size_t m = 0; m < s_nne; ++m) {
374 f(m, 0) +=
375 (dNx(m, 0) * sig(0, 0) + dNx(m, 1) * sig(1, 0) + dNx(m, 2) * sig(2, 0)) * v;
376 f(m, 1) +=
377 (dNx(m, 0) * sig(0, 1) + dNx(m, 1) * sig(1, 1) + dNx(m, 2) * sig(2, 1)) * v;
378 f(m, 2) +=
379 (dNx(m, 0) * sig(0, 2) + dNx(m, 1) * sig(1, 2) + dNx(m, 2) * sig(2, 2)) * v;
380 }
381 }
382 }
383 }
384
385 constexpr static size_t s_nne = 8;
386 constexpr static size_t s_ndim = 3;
387 constexpr static size_t s_tdim = 3;
388 size_t m_tdim = 3;
389 size_t m_nelem;
390 size_t m_nip;
391 array_type::tensor<double, 3> m_x;
392 array_type::tensor<double, 1> m_w;
393 array_type::tensor<double, 2> m_xi;
394 array_type::tensor<double, 2> m_N;
395 array_type::tensor<double, 3> m_dNxi;
396 array_type::tensor<double, 4> m_dNx;
397 array_type::tensor<double, 2> m_vol;
398};
399
400} // namespace Hex8
401} // namespace Element
402} // namespace GooseFEM
403
404#endif
GooseFEM::Element::Hex8::Quadrature
Interpolation and quadrature.
Definition ElementHex8.h:215
GooseFEM::Element::Hex8::Quadrature::Quadrature
Quadrature(const T &x, const X &xi, const W &w)
Constructor with custom integration.
Definition ElementHex8.h:255
GooseFEM::Element::Hex8::Quadrature::Quadrature
Quadrature(const T &x)
Constructor: use default Gauss integration.
Definition ElementHex8.h:234
GooseFEM::Element::QuadratureBaseCartesian
CRTP base class for interpolation and quadrature for a generic element in Cartesian coordinates.
Definition Element.h:546
GooseFEM::Element::QuadratureBaseCartesian< Quadrature >::compute_dN
void compute_dN()
Update the shape function gradients (called when the nodal positions are updated).
Definition Element.h:874
GooseFEM::Element::QuadratureBase
CRTP base class for quadrature.
Definition Element.h:151
GooseFEM::Element::QuadratureBase::shape_qscalar
auto shape_qscalar() const -> std::array< size_t, 2 >
Get the shape of a "qscalar" (a "qtensor" of rank 0)
Definition Element.h:354
GooseFEM::Element::QuadratureBase::shape_elemmat
auto shape_elemmat() const -> std::array< size_t, 3 >
Get the shape of an "elemmat".
Definition Element.h:273
GooseFEM::Element::QuadratureBase::shape_elemvec
auto shape_elemvec() const -> std::array< size_t, 3 >
Get the shape of an "elemvec".
Definition Element.h:252
GooseFEM::Element::QuadratureBase::AsTensor
auto AsTensor(const T &arg) const
Convert "qscalar" to "qtensor" of certain rank.
Definition Element.h:229
config.h
Basic configuration:
GOOSEFEM_ASSERT
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
Definition config.h:97
GooseFEM::Element::Hex8::Gauss::nip
size_t nip()
Number of integration points:
Definition ElementHex8.h:35
GooseFEM::Element::Hex8::Gauss::xi
array_type::tensor< double, 2 > xi()
Integration point coordinates (local coordinates).
Definition ElementHex8.h:45
GooseFEM::Element::Hex8::Gauss::w
array_type::tensor< double, 1 > w()
Integration point weights.
Definition ElementHex8.h:92
GooseFEM::Element::Hex8::Nodal::xi
array_type::tensor< double, 2 > xi()
Integration point coordinates (local coordinates).
Definition ElementHex8.h:135
GooseFEM::Element::Hex8::Nodal::w
array_type::tensor< double, 1 > w()
Integration point weights.
Definition ElementHex8.h:182
GooseFEM::Element::Hex8::Nodal::nip
size_t nip()
Number of integration points:
Definition ElementHex8.h:125
GooseFEM::array_type::tensor
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Definition config.h:177
GooseFEM
Toolbox to perform finite element computations.
Definition Allocate.h:14