Caribou
Hexahedron.h
1 #pragma once
2 
3 #include <Caribou/config.h>
4 #include <Caribou/Geometry/BaseHexahedron.h>
5 #include <Caribou/Geometry/Quad.h>
6 #include <Eigen/Core>
7 
8 namespace caribou::geometry {
9 
10 template<UNSIGNED_INTEGER_TYPE Order = Linear>
11 struct Hexahedron;
12 
13 template<>
14 struct traits<Hexahedron <Linear>> {
15  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 3;
16  static constexpr UNSIGNED_INTEGER_TYPE Dimension = 3;
17  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 8;
18  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 8;
19 
20  using BoundaryElementType = Quad<3, Linear>;
21  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 6;
22 };
23 
43 template<>
44 struct Hexahedron <Linear> : public BaseHexahedron<Hexahedron <Linear>> {
45  // Types
46  using Base = BaseHexahedron<Hexahedron <Linear>>;
47  using LocalCoordinates = typename Base::LocalCoordinates;
48  using WorldCoordinates = typename Base::WorldCoordinates;
49 
50  using GaussNode = typename Base::GaussNode;
51 
52  template <UNSIGNED_INTEGER_TYPE Dim>
53  using Vector = typename Base::template Vector<Dim>;
54 
55  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
56  using Matrix = typename Base::template Matrix<Rows, Cols>;
57 
58  // Constants
59  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
60  static constexpr auto Dimension = Base::Dimension;
61  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
62  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
63 
64  static constexpr FLOATING_POINT_TYPE canonical_nodes [NumberOfNodesAtCompileTime][CanonicalDimension] = {
65  // u, v, w
66  {-1, -1, -1}, // Node 0
67  {+1, -1, -1}, // Node 1
68  {+1, +1, -1}, // Node 2
69  {-1, +1, -1}, // Node 3
70  {-1, -1, +1}, // Node 4
71  {+1, -1, +1}, // Node 5
72  {+1, +1, +1}, // Node 6
73  {-1, +1, +1}, // Node 7
74  };
75 
76  // Constructors
77  using Base::Base;
78  Hexahedron() : Base() {
79  this->p_nodes.row(0) = WorldCoordinates(canonical_nodes[0][0], canonical_nodes[0][1], canonical_nodes[0][2]);
80  this->p_nodes.row(1) = WorldCoordinates(canonical_nodes[1][0], canonical_nodes[1][1], canonical_nodes[1][2]);
81  this->p_nodes.row(2) = WorldCoordinates(canonical_nodes[2][0], canonical_nodes[2][1], canonical_nodes[2][2]);
82  this->p_nodes.row(3) = WorldCoordinates(canonical_nodes[3][0], canonical_nodes[3][1], canonical_nodes[3][2]);
83  this->p_nodes.row(4) = WorldCoordinates(canonical_nodes[4][0], canonical_nodes[4][1], canonical_nodes[4][2]);
84  this->p_nodes.row(5) = WorldCoordinates(canonical_nodes[5][0], canonical_nodes[5][1], canonical_nodes[5][2]);
85  this->p_nodes.row(6) = WorldCoordinates(canonical_nodes[6][0], canonical_nodes[6][1], canonical_nodes[6][2]);
86  this->p_nodes.row(7) = WorldCoordinates(canonical_nodes[7][0], canonical_nodes[7][1], canonical_nodes[7][2]);
87  };
88 
89  // Public methods
94  inline auto edges() const {
95  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 2>, 12> indices = {{
96  {0, 1}, // Edge 0
97  {1, 2}, // Edge 1
98  {2, 3}, // Edge 2
99  {3, 0}, // Edge 3
100  {4, 5}, // Edge 4
101  {5, 6}, // Edge 5
102  {6, 7}, // Edge 6
103  {7, 4}, // Edge 7
104  {0, 4}, // Edge 8
105  {1, 5}, // Edge 9
106  {2, 6}, // Edge 10
107  {3, 7} // Edge 11
108  }};
109  return indices;
110  }
111 
112 
113 private:
114  // Implementations
115  friend struct Element<Hexahedron <Linear>>;
116  friend struct BaseHexahedron<Hexahedron <Linear>>;
117  inline auto get_L(const LocalCoordinates & xi) const {
118  const auto & u = xi[0];
119  const auto & v = xi[1];
120  const auto & w = xi[2];
121 
122  Vector<NumberOfNodesAtCompileTime> m;
123  m << (1/8.) * (1 - u) * (1 - v) * (1 - w),
124  (1/8.) * (1 + u) * (1 - v) * (1 - w),
125  (1/8.) * (1 + u) * (1 + v) * (1 - w),
126  (1/8.) * (1 - u) * (1 + v) * (1 - w),
127  (1/8.) * (1 - u) * (1 - v) * (1 + w),
128  (1/8.) * (1 + u) * (1 - v) * (1 + w),
129  (1/8.) * (1 + u) * (1 + v) * (1 + w),
130  (1/8.) * (1 - u) * (1 + v) * (1 + w);
131  return m;
132  };
133 
134  inline auto get_dL(const LocalCoordinates & xi) const {
135  const auto & u = xi[0];
136  const auto & v = xi[1];
137  const auto & w = xi[2];
138 
139  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
140  // dL/du dL/dv dL/dw
141  m << -1/8. * (1 - v) * (1 - w), -1/8. * (1 - u) * (1 - w), -1/8. * (1 - u) * (1 - v), // Node 0
142  +1/8. * (1 - v) * (1 - w), -1/8. * (1 + u) * (1 - w), -1/8. * (1 + u) * (1 - v), // Node 1
143  +1/8. * (1 + v) * (1 - w), +1/8. * (1 + u) * (1 - w), -1/8. * (1 + u) * (1 + v), // Node 2
144  -1/8. * (1 + v) * (1 - w), +1/8. * (1 - u) * (1 - w), -1/8. * (1 - u) * (1 + v), // Node 3
145  -1/8. * (1 - v) * (1 + w), -1/8. * (1 - u) * (1 + w), +1/8. * (1 - u) * (1 - v), // Node 4
146  +1/8. * (1 - v) * (1 + w), -1/8. * (1 + u) * (1 + w), +1/8. * (1 + u) * (1 - v), // Node 5
147  +1/8. * (1 + v) * (1 + w), +1/8. * (1 + u) * (1 + w), +1/8. * (1 + u) * (1 + v), // Node 6
148  -1/8. * (1 + v) * (1 + w), +1/8. * (1 - u) * (1 + w), +1/8. * (1 - u) * (1 + v); // Node 7
149  return m;
150  };
151 
152  inline auto get_gauss_nodes() const -> const auto & {
153  using Weight = FLOATING_POINT_TYPE;
154  static const std::vector<GaussNode> gauss_nodes {
155  GaussNode {LocalCoordinates(-1/sqrt(3), -1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 0
156  GaussNode {LocalCoordinates(+1/sqrt(3), -1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 1
157  GaussNode {LocalCoordinates(-1/sqrt(3), +1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 2
158  GaussNode {LocalCoordinates(+1/sqrt(3), +1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 3
159  GaussNode {LocalCoordinates(-1/sqrt(3), -1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 4
160  GaussNode {LocalCoordinates(+1/sqrt(3), -1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 5
161  GaussNode {LocalCoordinates(-1/sqrt(3), +1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 6
162  GaussNode {LocalCoordinates(+1/sqrt(3), +1/sqrt(3), +1/sqrt(3)), Weight(1)} // Node 7
163  };
164  return gauss_nodes;
165  }
166 
167  inline auto get_boundary_elements_nodes() const -> const auto & {
168  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 4>, 6> indices = {{
169  {0, 3, 2, 1}, // Face 0
170  {0, 4, 7, 3}, // Face 1
171  {1, 2, 6, 5}, // Face 2
172  {0, 1, 5, 4}, // Face 3
173  {2, 3, 7, 6}, // Face 4
174  {4, 5, 6, 7} // Face 5
175  }};
176  return indices;
177  }
178 };
179 
180 
181 template<>
182 struct traits<Hexahedron <Quadratic>> {
183  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 3;
184  static constexpr UNSIGNED_INTEGER_TYPE Dimension = 3;
185  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 20;
186  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 4;
187 
188  using BoundaryElementType = Quad<3, Quadratic>;
189  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 6;
190 };
191 
213 template<>
214 struct Hexahedron <Quadratic> : public BaseHexahedron<Hexahedron <Quadratic>> {
215  // Types
216  using Base = BaseHexahedron<Hexahedron <Quadratic>>;
217  using LocalCoordinates = typename Base::LocalCoordinates;
218  using WorldCoordinates = typename Base::WorldCoordinates;
219 
220  using GaussNode = typename Base::GaussNode;
221 
222  template <UNSIGNED_INTEGER_TYPE Dim>
223  using Vector = typename Base::template Vector<Dim>;
224 
225  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
226  using Matrix = typename Base::template Matrix<Rows, Cols>;
227 
228  // Constants
229  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
230  static constexpr auto Dimension = Base::Dimension;
231  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
232  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
233 
234  static constexpr FLOATING_POINT_TYPE canonical_nodes [NumberOfNodesAtCompileTime][CanonicalDimension] = {
235  // u, v, w
236  {-1, -1, -1}, // Node 0
237  {+1, -1, -1}, // Node 1
238  {+1, +1, -1}, // Node 2
239  {-1, +1, -1}, // Node 3
240  {-1, -1, +1}, // Node 4
241  {+1, -1, +1}, // Node 5
242  {+1, +1, +1}, // Node 6
243  {-1, +1, +1}, // Node 7
244  { 0, -1, -1}, // Node 8
245  {+1, 0, -1}, // Node 9
246  { 0, +1, -1}, // Node 10
247  {-1, 0, -1}, // Node 11
248  { 0, -1, +1}, // Node 12
249  {+1, 0, +1}, // Node 13
250  { 0, +1, +1}, // Node 14
251  {-1, 0, +1}, // Node 15
252  {-1, -1, 0}, // Node 16
253  {+1, -1, 0}, // Node 17
254  {+1, +1, 0}, // Node 18
255  {-1, +1, 0}, // Node 19
256  };
257 
258  // Constructors
259  using Base::Base;
260  Hexahedron() : Base() {
261  this->p_nodes.row(0) = WorldCoordinates(canonical_nodes[0][0], canonical_nodes[0][1], canonical_nodes[0][2]); // Node 0
262  this->p_nodes.row(1) = WorldCoordinates(canonical_nodes[1][0], canonical_nodes[1][1], canonical_nodes[1][2]); // Node 1
263  this->p_nodes.row(2) = WorldCoordinates(canonical_nodes[2][0], canonical_nodes[2][1], canonical_nodes[2][2]); // Node 2
264  this->p_nodes.row(3) = WorldCoordinates(canonical_nodes[3][0], canonical_nodes[3][1], canonical_nodes[3][2]); // Node 3
265  this->p_nodes.row(4) = WorldCoordinates(canonical_nodes[4][0], canonical_nodes[4][1], canonical_nodes[4][2]); // Node 4
266  this->p_nodes.row(5) = WorldCoordinates(canonical_nodes[5][0], canonical_nodes[5][1], canonical_nodes[5][2]); // Node 5
267  this->p_nodes.row(6) = WorldCoordinates(canonical_nodes[6][0], canonical_nodes[6][1], canonical_nodes[6][2]); // Node 6
268  this->p_nodes.row(7) = WorldCoordinates(canonical_nodes[7][0], canonical_nodes[7][1], canonical_nodes[7][2]); // Node 7
269  this->p_nodes.row(8) = WorldCoordinates(canonical_nodes[8][0], canonical_nodes[8][1], canonical_nodes[8][2]); // Node 8
270  this->p_nodes.row(9) = WorldCoordinates(canonical_nodes[9][0], canonical_nodes[9][1], canonical_nodes[9][2]); // Node 9
271  this->p_nodes.row(10) = WorldCoordinates(canonical_nodes[10][0], canonical_nodes[10][1], canonical_nodes[10][2]); // Node 10
272  this->p_nodes.row(11) = WorldCoordinates(canonical_nodes[11][0], canonical_nodes[11][1], canonical_nodes[11][2]); // Node 11
273  this->p_nodes.row(12) = WorldCoordinates(canonical_nodes[12][0], canonical_nodes[12][1], canonical_nodes[12][2]); // Node 12
274  this->p_nodes.row(13) = WorldCoordinates(canonical_nodes[13][0], canonical_nodes[13][1], canonical_nodes[13][2]); // Node 13
275  this->p_nodes.row(14) = WorldCoordinates(canonical_nodes[14][0], canonical_nodes[14][1], canonical_nodes[14][2]); // Node 14
276  this->p_nodes.row(15) = WorldCoordinates(canonical_nodes[15][0], canonical_nodes[15][1], canonical_nodes[15][2]); // Node 15
277  this->p_nodes.row(16) = WorldCoordinates(canonical_nodes[16][0], canonical_nodes[16][1], canonical_nodes[16][2]); // Node 16
278  this->p_nodes.row(17) = WorldCoordinates(canonical_nodes[17][0], canonical_nodes[17][1], canonical_nodes[17][2]); // Node 17
279  this->p_nodes.row(18) = WorldCoordinates(canonical_nodes[18][0], canonical_nodes[18][1], canonical_nodes[18][2]); // Node 18
280  this->p_nodes.row(19) = WorldCoordinates(canonical_nodes[19][0], canonical_nodes[19][1], canonical_nodes[19][2]); // Node 19
281 
282  };
283 
284  // Construct a quadratic Hexahedron from a linear one
285  explicit Hexahedron(const Hexahedron<Linear> & linear_Hexahedron) : Base() {
286  this->p_nodes.row(0) = linear_Hexahedron.node(0); // Node 0
287  this->p_nodes.row(1) = linear_Hexahedron.node(1); // Node 1
288  this->p_nodes.row(2) = linear_Hexahedron.node(2); // Node 2
289  this->p_nodes.row(3) = linear_Hexahedron.node(3); // Node 3
290  this->p_nodes.row(4) = linear_Hexahedron.node(4); // Node 4
291  this->p_nodes.row(5) = linear_Hexahedron.node(5); // Node 5
292  this->p_nodes.row(6) = linear_Hexahedron.node(6); // Node 6
293  this->p_nodes.row(7) = linear_Hexahedron.node(7); // Node 7
294  this->p_nodes.row(8) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[8][0], canonical_nodes[8][1], canonical_nodes[8][2])); // Node 8
295  this->p_nodes.row(9) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[9][0], canonical_nodes[9][1], canonical_nodes[9][2])); // Node 9
296  this->p_nodes.row(10) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[10][0], canonical_nodes[10][1], canonical_nodes[10][2])); // Node 10
297  this->p_nodes.row(11) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[11][0], canonical_nodes[11][1], canonical_nodes[11][2])); // Node 11
298  this->p_nodes.row(12) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[12][0], canonical_nodes[12][1], canonical_nodes[12][2])); // Node 12
299  this->p_nodes.row(13) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[13][0], canonical_nodes[13][1], canonical_nodes[13][2])); // Node 13
300  this->p_nodes.row(14) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[14][0], canonical_nodes[14][1], canonical_nodes[14][2])); // Node 14
301  this->p_nodes.row(15) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[15][0], canonical_nodes[15][1], canonical_nodes[15][2])); // Node 15
302  this->p_nodes.row(16) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[16][0], canonical_nodes[16][1], canonical_nodes[16][2])); // Node 16
303  this->p_nodes.row(17) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[17][0], canonical_nodes[17][1], canonical_nodes[17][2])); // Node 17
304  this->p_nodes.row(18) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[18][0], canonical_nodes[18][1], canonical_nodes[18][2])); // Node 18
305  this->p_nodes.row(19) = linear_Hexahedron.world_coordinates(LocalCoordinates(canonical_nodes[19][0], canonical_nodes[19][1], canonical_nodes[19][2])); // Node 19
306  }
307 
309  Hexahedron(WorldCoordinates & p0, WorldCoordinates & p1, WorldCoordinates & p2, WorldCoordinates & p3, WorldCoordinates & p4, WorldCoordinates & p5, WorldCoordinates & p6, WorldCoordinates & p7)
310  : Hexahedron(Hexahedron<Linear>(p0, p1, p2, p3, p4, p5, p6, p7)) {}
311 
313  Hexahedron(const WorldCoordinates & p0, const WorldCoordinates & p1, const WorldCoordinates & p2, const WorldCoordinates & p3, const WorldCoordinates & p4, const WorldCoordinates & p5, const WorldCoordinates & p6, const WorldCoordinates & p7)
314  : Hexahedron(Hexahedron<Linear>(p0, p1, p2, p3, p4, p5, p6, p7)) {}
315 
316  // Public methods
321  inline auto edges() const {
322  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 3>, 12> indices = {{
323  {0, 1, 8}, // Edge 0
324  {1, 2, 9}, // Edge 1
325  {2, 3, 10}, // Edge 2
326  {3, 0, 11}, // Edge 3
327  {4, 5, 12}, // Edge 4
328  {5, 6, 13}, // Edge 5
329  {6, 7, 14}, // Edge 6
330  {7, 4, 15}, // Edge 7
331  {0, 4, 16}, // Edge 8
332  {1, 5, 17}, // Edge 9
333  {2, 6, 18}, // Edge 10
334  {3, 7, 19} // Edge 11
335  }};
336  return indices;
337  }
338 
339 private:
340  // Implementations
341  friend struct Element<Hexahedron <Quadratic>>;
342  friend struct BaseHexahedron<Hexahedron <Quadratic>>;
343  inline auto get_L(const LocalCoordinates & xi) const {
344  const auto & u = xi[0];
345  const auto & v = xi[1];
346  const auto & w = xi[2];
347 
348  Vector<NumberOfNodesAtCompileTime> m;
349  m[ 8] = 1/4.*(1 - u*u)*(1 - v)*(1 - w);
350  m[ 9] = 1/4.*(1 - v*v)*(1 + u)*(1 - w);
351  m[10] = 1/4.*(1 - u*u)*(1 + v)*(1 - w);
352  m[11] = 1/4.*(1 - v*v)*(1 - u)*(1 - w);
353  m[12] = 1/4.*(1 - u*u)*(1 - v)*(1 + w);
354  m[13] = 1/4.*(1 - v*v)*(1 + u)*(1 + w);
355  m[14] = 1/4.*(1 - u*u)*(1 + v)*(1 + w);
356  m[15] = 1/4.*(1 - v*v)*(1 - u)*(1 + w);
357  m[16] = 1/4.*(1 - w*w)*(1 - u)*(1 - v);
358  m[17] = 1/4.*(1 - w*w)*(1 + u)*(1 - v);
359  m[18] = 1/4.*(1 - w*w)*(1 + u)*(1 + v);
360  m[19] = 1/4.*(1 - w*w)*(1 - u)*(1 + v);
361 
362  m[0] = 1/8.*(1 - u)*(1 - v)*(1 - w) - 1/2.*(m[ 8] + m[11] + m[16]);
363  m[1] = 1/8.*(1 + u)*(1 - v)*(1 - w) - 1/2.*(m[ 8] + m[ 9] + m[17]);
364  m[2] = 1/8.*(1 + u)*(1 + v)*(1 - w) - 1/2.*(m[ 9] + m[10] + m[18]);
365  m[3] = 1/8.*(1 - u)*(1 + v)*(1 - w) - 1/2.*(m[10] + m[11] + m[19]);
366  m[4] = 1/8.*(1 - u)*(1 - v)*(1 + w) - 1/2.*(m[12] + m[15] + m[16]);
367  m[5] = 1/8.*(1 + u)*(1 - v)*(1 + w) - 1/2.*(m[12] + m[13] + m[17]);
368  m[6] = 1/8.*(1 + u)*(1 + v)*(1 + w) - 1/2.*(m[13] + m[14] + m[18]);
369  m[7] = 1/8.*(1 - u)*(1 + v)*(1 + w) - 1/2.*(m[14] + m[15] + m[19]);
370  return m;
371  };
372 
373  inline auto get_dL(const LocalCoordinates & xi) const {
374  using ShapeDerivative = Vector<3>;
375  const auto & u = xi[0];
376  const auto & v = xi[1];
377  const auto & w = xi[2];
378 
379  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
380 
381  // dL/du dL/dv dL/dw
382  m.row( 8) = ShapeDerivative({-1/2.*u*(1 - v)*(1 - w), -1/4.*(1 - u*u)*(1 - w), -1/4.*(1 - u*u)*(1 - v)});
383  m.row( 9) = ShapeDerivative({ 1/4.*(1 - v*v)*(1 - w), -1/2.*v*(1 + u)*(1 - w), -1/4.*(1 - v*v)*(1 + u)});
384  m.row(10) = ShapeDerivative({-1/2.*u*(1 + v)*(1 - w), 1/4.*(1 - u*u)*(1 - w), -1/4.*(1 - u*u)*(1 + v)});
385  m.row(11) = ShapeDerivative({-1/4.*(1 - v*v)*(1 - w), -1/2.*v*(1 - u)*(1 - w), -1/4.*(1 - v*v)*(1 - u)});
386  m.row(12) = ShapeDerivative({-1/2.*u*(1 - v)*(1 + w), -1/4.*(1 - u*u)*(1 + w), 1/4.*(1 - u*u)*(1 - v)});
387  m.row(13) = ShapeDerivative({ 1/4.*(1 - v*v)*(1 + w), -1/2.*v*(1 + u)*(1 + w), 1/4.*(1 - v*v)*(1 + u)});
388  m.row(14) = ShapeDerivative({-1/2.*u*(1 + v)*(1 + w), 1/4.*(1 - u*u)*(1 + w), 1/4.*(1 - u*u)*(1 + v)});
389  m.row(15) = ShapeDerivative({-1/4.*(1 - v*v)*(1 + w), -1/2.*v*(1 - u)*(1 + w), 1/4.*(1 - v*v)*(1 - u)});
390  m.row(16) = ShapeDerivative({-1/4.*(1 - w*w)*(1 - v), -1/4.*(1 - w*w)*(1 - u), -1/2.*w*(1 - u)*(1 - v)});
391  m.row(17) = ShapeDerivative({ 1/4.*(1 - w*w)*(1 - v), -1/4.*(1 - w*w)*(1 + u), -1/2.*w*(1 + u)*(1 - v)});
392  m.row(18) = ShapeDerivative({ 1/4.*(1 - w*w)*(1 + v), 1/4.*(1 - w*w)*(1 + u), -1/2.*w*(1 + u)*(1 + v)});
393  m.row(19) = ShapeDerivative({-1/4.*(1 - w*w)*(1 + v), 1/4.*(1 - w*w)*(1 - u), -1/2.*w*(1 - u)*(1 + v)});
394 
395  const auto du = m.col(0); // dL/du
396  const auto dv = m.col(1); // dL/dv
397  const auto dw = m.col(2); // dL/dw
398  // dL/du dL/dv dL/dw
399  m.row(0) = ShapeDerivative({-1/8.*(1 - v)*(1 - w) - 1/2.*(du[8 ] + du[11] + du[16]), -1/8.*(1 - u)*(1 - w) - 1/2.*(dv[ 8] + dv[11] + dv[16]), -1/8.*(1 - u)*(1 - v) - 1/2.*(dw[ 8] + dw[11] + dw[16])});
400  m.row(1) = ShapeDerivative({ 1/8.*(1 - v)*(1 - w) - 1/2.*(du[8 ] + du[ 9] + du[17]), -1/8.*(1 + u)*(1 - w) - 1/2.*(dv[ 8] + dv[ 9] + dv[17]), -1/8.*(1 + u)*(1 - v) - 1/2.*(dw[ 8] + dw[ 9] + dw[17])});
401  m.row(2) = ShapeDerivative({ 1/8.*(1 + v)*(1 - w) - 1/2.*(du[9 ] + du[10] + du[18]), 1/8.*(1 + u)*(1 - w) - 1/2.*(dv[ 9] + dv[10] + dv[18]), -1/8.*(1 + u)*(1 + v) - 1/2.*(dw[ 9] + dw[10] + dw[18])});
402  m.row(3) = ShapeDerivative({-1/8.*(1 + v)*(1 - w) - 1/2.*(du[10] + du[11] + du[19]), 1/8.*(1 - u)*(1 - w) - 1/2.*(dv[10] + dv[11] + dv[19]), -1/8.*(1 - u)*(1 + v) - 1/2.*(dw[10] + dw[11] + dw[19])});
403  m.row(4) = ShapeDerivative({-1/8.*(1 - v)*(1 + w) - 1/2.*(du[12] + du[15] + du[16]), -1/8.*(1 - u)*(1 + w) - 1/2.*(dv[12] + dv[15] + dv[16]), 1/8.*(1 - u)*(1 - v) - 1/2.*(dw[12] + dw[15] + dw[16])});
404  m.row(5) = ShapeDerivative({ 1/8.*(1 - v)*(1 + w) - 1/2.*(du[12] + du[13] + du[17]), -1/8.*(1 + u)*(1 + w) - 1/2.*(dv[12] + dv[13] + dv[17]), 1/8.*(1 + u)*(1 - v) - 1/2.*(dw[12] + dw[13] + dw[17])});
405  m.row(6) = ShapeDerivative({ 1/8.*(1 + v)*(1 + w) - 1/2.*(du[13] + du[14] + du[18]), 1/8.*(1 + u)*(1 + w) - 1/2.*(dv[13] + dv[14] + dv[18]), 1/8.*(1 + u)*(1 + v) - 1/2.*(dw[13] + dw[14] + dw[18])});
406  m.row(7) = ShapeDerivative({-1/8.*(1 + v)*(1 + w) - 1/2.*(du[14] + du[15] + du[19]), 1/8.*(1 - u)*(1 + w) - 1/2.*(dv[14] + dv[15] + dv[19]), 1/8.*(1 - u)*(1 + v) - 1/2.*(dw[14] + dw[15] + dw[19])});
407 
408 
409  return m;
410  };
411 
412  inline auto get_gauss_nodes() const -> const auto & {
413  using Weight = FLOATING_POINT_TYPE;
414  static const std::vector<GaussNode> gauss_nodes {
415  GaussNode {LocalCoordinates(-1/sqrt(3), -1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 0
416  GaussNode {LocalCoordinates(+1/sqrt(3), -1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 1
417  GaussNode {LocalCoordinates(-1/sqrt(3), +1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 2
418  GaussNode {LocalCoordinates(+1/sqrt(3), +1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 3
419  GaussNode {LocalCoordinates(-1/sqrt(3), -1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 4
420  GaussNode {LocalCoordinates(+1/sqrt(3), -1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 5
421  GaussNode {LocalCoordinates(-1/sqrt(3), +1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 6
422  GaussNode {LocalCoordinates(+1/sqrt(3), +1/sqrt(3), +1/sqrt(3)), Weight(1)} // Node 7
423  };
424  return gauss_nodes;
425  }
426 
427  inline auto get_boundary_elements_nodes() const -> const auto & {
428  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 8>, 6> indices = {{
429  {0, 3, 2, 1, 11, 10, 9, 8}, // Face 0
430  {0, 4, 7, 3, 16, 15, 19, 11}, // Face 1
431  {1, 2, 6, 5, 9, 18, 13, 17}, // Face 2
432  {0, 1, 5, 4, 8, 17, 12, 16}, // Face 3
433  {2, 3, 7, 6, 10, 19, 14, 18}, // Face 4
434  {4, 5, 6, 7, 12, 13, 14, 15} // Face 5
435  }};
436  return indices;
437  }
438 };
439 
440 }
caribou::geometry::Hexahedron
Definition: Hexahedron.h:11
caribou::geometry::Element::world_coordinates
auto world_coordinates(const LocalCoordinates &coordinates) const -> WorldCoordinates
Get the world coordinates of a point from its local coordinates.
Definition: Element.h:157
caribou::geometry::Element< Derived >
caribou::geometry::Hexahedron< Linear >::edges
auto edges() const
Get the list of node indices of the edges.
Definition: Hexahedron.h:94
caribou::geometry::Quad
Definition: Quad.h:12
caribou::geometry::Hexahedron< Quadratic >::Hexahedron
Hexahedron(const WorldCoordinates &p0, const WorldCoordinates &p1, const WorldCoordinates &p2, const WorldCoordinates &p3, const WorldCoordinates &p4, const WorldCoordinates &p5, const WorldCoordinates &p6, const WorldCoordinates &p7)
Construct a quadratic Hexahedron from eight nodes.
Definition: Hexahedron.h:313
caribou::geometry::traits
Definition: Element.h:14
caribou::geometry::Element::node
auto node(const UNSIGNED_INTEGER_TYPE &index) const
Get the Node at given index.
Definition: Element.h:64
caribou::geometry::Hexahedron< Quadratic >::Hexahedron
Hexahedron(WorldCoordinates &p0, WorldCoordinates &p1, WorldCoordinates &p2, WorldCoordinates &p3, WorldCoordinates &p4, WorldCoordinates &p5, WorldCoordinates &p6, WorldCoordinates &p7)
Construct a quadratic Hexahedron from eight nodes.
Definition: Hexahedron.h:309
caribou::geometry::Hexahedron< Quadratic >::edges
auto edges() const
Get the list of node indices of the edges.
Definition: Hexahedron.h:321
caribou::geometry::Hexahedron< Linear >
Linear Hexahedron.
Definition: Hexahedron.h:44
caribou::geometry::BaseHexahedron
Base class for hexahedral elements.
Definition: BaseHexahedron.h:16
Generated by doxygen 1.8.20