Caribou
Tetrahedron.h
1 #pragma once
2 
3 #include <Caribou/config.h>
4 #include <Caribou/Geometry/BaseTetrahedron.h>
5 #include <Caribou/Geometry/Triangle.h>
6 #include <Eigen/Core>
7 #include <array>
8 
9 namespace caribou::geometry {
10 
11 template<UNSIGNED_INTEGER_TYPE _Order = Linear>
12 struct Tetrahedron;
13 
14 template<>
15 struct traits<Tetrahedron <Linear>> {
16  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 3;
17  static constexpr UNSIGNED_INTEGER_TYPE Dimension = 3;
18  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 4;
19  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 1;
20 
21  using BoundaryElementType = Triangle<3, Linear>;
22  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 4;
23 };
24 
50 template<>
51 struct Tetrahedron <Linear> : public BaseTetrahedron<Tetrahedron <Linear>> {
52  // Types
53  using Base = BaseTetrahedron<Tetrahedron <Linear>>;
54  using LocalCoordinates = typename Base::LocalCoordinates;
55  using WorldCoordinates = typename Base::WorldCoordinates;
56 
57  using GaussNode = typename Base::GaussNode;
58 
59  template <UNSIGNED_INTEGER_TYPE Dim>
60  using Vector = typename Base::template Vector<Dim>;
61 
62  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
63  using Matrix = typename Base::template Matrix<Rows, Cols>;
64 
65  // Constants
66  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
67  static constexpr auto Dimension = Base::Dimension;
68  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
69  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
70 
71  // Constructors
72  using Base::Base;
73  Tetrahedron() : Base() {
74  this->p_nodes.row(0) = WorldCoordinates(0, 0, 0);
75  this->p_nodes.row(1) = WorldCoordinates(1, 0, 0);
76  this->p_nodes.row(2) = WorldCoordinates(0, 1, 0);
77  this->p_nodes.row(3) = WorldCoordinates(0, 0, 1);
78  };
79 
80  // Public methods
85  inline auto edges() const {
86  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 2>, 6> indices = {{
87  {0, 1}, // Edge 0
88  {1, 2}, // Edge 1
89  {2, 0}, // Edge 2
90  {0, 3}, // Edge 3
91  {3, 1}, // Edge 4
92  {3, 2} // Edge 5
93  }};
94  return indices;
95  }
96 
97 
98 private:
99  // Implementations
100  friend struct Element<Tetrahedron <Linear>>;
101  friend struct BaseTetrahedron<Tetrahedron <Linear>>;
102  inline auto get_L(const LocalCoordinates & xi) const -> Vector<NumberOfNodesAtCompileTime> {
103  const auto & u = xi[0];
104  const auto & v = xi[1];
105  const auto & w = xi[2];
106  return {
107  static_cast<FLOATING_POINT_TYPE>(1 - u - v - w),
108  static_cast<FLOATING_POINT_TYPE>(u),
109  static_cast<FLOATING_POINT_TYPE>(v),
110  static_cast<FLOATING_POINT_TYPE>(w)
111  };
112  };
113 
114  inline auto get_dL(const LocalCoordinates & /*xi*/) const {
115  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
116  // dL/du dL/dv dL/dw
117  m << -1, -1, -1, // Node 0
118  1, 0, 0, // Node 1
119  0, 1, 0, // Node 1
120  0, 0, 1; // Node 3
121  return m;
122  };
123 
124  inline auto get_gauss_nodes() const -> const auto & {
125  using Weight = FLOATING_POINT_TYPE;
126  static const std::vector<GaussNode> gauss_nodes {
127  GaussNode {LocalCoordinates(1/4., 1/4., 1/4.), Weight(1/6.)} // Node 0
128  };
129  return gauss_nodes;
130  }
131 
132  inline auto get_boundary_elements_nodes() const -> const auto & {
133  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 3>, 4> indices = {{
134  {0, 2, 1}, // Face 0
135  {0, 1, 3}, // Face 1
136  {0, 3, 2}, // Face 2
137  {3, 1, 2} // Face 3
138  }};
139  return indices;
140  }
141 };
142 
143 
144 template<>
145 struct traits<Tetrahedron <Quadratic>> {
146  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 3;
147  static constexpr UNSIGNED_INTEGER_TYPE Dimension = 3;
148  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 10;
149  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 4;
150 
151  using BoundaryElementType = Triangle<3, Quadratic>;
152  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 4;
153 };
154 
180 template<>
181 struct Tetrahedron <Quadratic> : public BaseTetrahedron<Tetrahedron <Quadratic>> {
182  // Types
183  using Base = BaseTetrahedron<Tetrahedron <Quadratic>>;
184  using LocalCoordinates = typename Base::LocalCoordinates;
185  using WorldCoordinates = typename Base::WorldCoordinates;
186 
187  using GaussNode = typename Base::GaussNode;
188 
189  template <UNSIGNED_INTEGER_TYPE Dim>
190  using Vector = typename Base::template Vector<Dim>;
191 
192  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
193  using Matrix = typename Base::template Matrix<Rows, Cols>;
194 
195  // Constants
196  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
197  static constexpr auto Dimension = Base::Dimension;
198  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
199  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
200 
201  // Constructors
202  using Base::Base;
203  Tetrahedron() : Base() {
204  this->p_nodes.row(0) = WorldCoordinates(0.0, 0.0, 0.0); // Node 0
205  this->p_nodes.row(1) = WorldCoordinates(1.0, 0.0, 0.0); // Node 1
206  this->p_nodes.row(2) = WorldCoordinates(0.0, 1.0, 0.0); // Node 2
207  this->p_nodes.row(3) = WorldCoordinates(0.0, 0.0, 1.0); // Node 3
208  this->p_nodes.row(4) = WorldCoordinates(0.5, 0.0, 0.0); // Node 4
209  this->p_nodes.row(5) = WorldCoordinates(0.5, 0.5, 0.0); // Node 5
210  this->p_nodes.row(6) = WorldCoordinates(0.0, 0.5, 0.0); // Node 6
211  this->p_nodes.row(7) = WorldCoordinates(0.0, 0.0, 0.5); // Node 7
212  this->p_nodes.row(8) = WorldCoordinates(0.5, 0.0, 0.5); // Node 8
213  this->p_nodes.row(9) = WorldCoordinates(0.0, 0.5, 0.5); // Node 9
214  };
215 
216  // Construct a quadratic Tetrahedron from a linear one
217  explicit Tetrahedron(const Tetrahedron<Linear> & linear_Tetrahedron) : Base() {
218  this->p_nodes.row(0) = linear_Tetrahedron.node(0); // Node 0
219  this->p_nodes.row(1) = linear_Tetrahedron.node(1); // Node 1
220  this->p_nodes.row(2) = linear_Tetrahedron.node(2); // Node 2
221  this->p_nodes.row(3) = linear_Tetrahedron.node(3); // Node 3
222  this->p_nodes.row(4) = linear_Tetrahedron.world_coordinates(LocalCoordinates(0.5, 0.0, 0.0)); // Node 4
223  this->p_nodes.row(5) = linear_Tetrahedron.world_coordinates(LocalCoordinates(0.5, 0.5, 0.0)); // Node 5
224  this->p_nodes.row(6) = linear_Tetrahedron.world_coordinates(LocalCoordinates(0.0, 0.5, 0.0)); // Node 6
225  this->p_nodes.row(7) = linear_Tetrahedron.world_coordinates(LocalCoordinates(0.0, 0.0, 0.5)); // Node 7
226  this->p_nodes.row(8) = linear_Tetrahedron.world_coordinates(LocalCoordinates(0.5, 0.0, 0.5)); // Node 8
227  this->p_nodes.row(9) = linear_Tetrahedron.world_coordinates(LocalCoordinates(0.0, 0.5, 0.5)); // Node 9
228  }
229 
231  Tetrahedron(WorldCoordinates & p0, WorldCoordinates & p1, WorldCoordinates & p2, WorldCoordinates & p3)
232  : Tetrahedron(Tetrahedron<Linear>(p0, p1, p2, p3)) {}
233 
235  Tetrahedron(const WorldCoordinates & p0, const WorldCoordinates & p1, const WorldCoordinates & p2, const WorldCoordinates & p3)
236  : Tetrahedron(Tetrahedron<Linear>(p0, p1, p2, p3)) {}
237 
238  // Public methods
243  inline auto edges() const {
244  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 3>, 6> indices = {{
245  {0, 1, 4}, // Edge 0
246  {1, 2, 5}, // Edge 1
247  {2, 0, 6}, // Edge 2
248  {0, 3, 7}, // Edge 3
249  {3, 1, 8}, // Edge 4
250  {3, 2, 9} // Edge 5
251  }};
252  return indices;
253  }
254 
255 private:
256  // Implementations
257  friend struct Element<Tetrahedron <Quadratic>>;
258  friend struct BaseTetrahedron<Tetrahedron <Quadratic>>;
259  inline auto get_L(const LocalCoordinates & xi) const {
260  const auto & u = xi[0];
261  const auto & v = xi[1];
262  const auto & w = xi[2];
263  const auto l = 1 - u - v - w;
264 
265  Vector<NumberOfNodesAtCompileTime> m;
266  m << l * (2*l - 1), // Node 0
267  u * (2*u - 1), // Node 1
268  v * (2*v - 1), // Node 2
269  w * (2*w - 1), // Node 3
270  4 * l * u, // Node 4
271  4 * u * v, // Node 5
272  4 * v * l, // Node 6
273  4 * l * w, // Node 7
274  4 * u * w, // Node 8
275  4 * v * w; // Node 9
276 
277  return m;
278  };
279 
280  inline auto get_dL(const LocalCoordinates & xi) const {
281  const auto & u = xi[0];
282  const auto & v = xi[1];
283  const auto & w = xi[2];
284  const auto l = 1 - u - v - w;
285 
286  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
287  // dL/du dL/dv dL/dw
288  m << 1 - (4 * l), 1 - (4 * l), 1 - (4 * l), // Node 0
289  (4 * u) - 1 , 0 , 0 , // Node 1
290  0 , (4 * v) - 1 , 0 , // Node 2
291  0 , 0 , (4 * w) - 1 , // Node 3
292  4 * (l - u), -4 * u , -4 * u , // Node 4
293  4 * v , 4 * u , 0 , // Node 5
294  -4 * v , 4 * (l - v), -4 * v , // Node 6
295  -4 * w , -4 * w , 4 * (l - w), // Node 7
296  4 * w , 0 , 4 * u , // Node 8
297  0 , 4 * w , 4 * v ; // Node 9
298  return m;
299  };
300 
301  inline auto get_gauss_nodes() const -> const auto & {
302  using Weight = FLOATING_POINT_TYPE;
303  static const std::vector<GaussNode> gauss_nodes {
304  GaussNode {LocalCoordinates(0.1381966011250105, 0.1381966011250105, 0.1381966011250105), Weight(1/24.)}, // Node 0
305  GaussNode {LocalCoordinates(0.1381966011250105, 0.1381966011250105, 0.5854101966249685), Weight(1/24.)}, // Node 1
306  GaussNode {LocalCoordinates(0.1381966011250105, 0.5854101966249685, 0.1381966011250105), Weight(1/24.)}, // Node 2
307  GaussNode {LocalCoordinates(0.5854101966249685, 0.1381966011250105, 0.1381966011250105), Weight(1/24.)} // Node 3
308  };
309  return gauss_nodes;
310  }
311 
312  inline auto get_boundary_elements_nodes() const -> const auto & {
313  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 6>, 4> indices = {{
314  {0, 2, 1, 6, 5, 4}, // Face 0
315  {0, 1, 3, 4, 8, 7}, // Face 1
316  {0, 3, 2, 7, 9, 6}, // Face 2
317  {3, 1, 2, 8, 5, 9} // Face 3
318  }};
319  return indices;
320  }
321 };
322 
323 }
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::BaseTetrahedron
Definition: BaseTetrahedron.h:12
caribou::geometry::Tetrahedron
Definition: Tetrahedron.h:12
caribou::geometry::Tetrahedron< Linear >
Linear Tetrahedron.
Definition: Tetrahedron.h:51
caribou::geometry::Tetrahedron< Quadratic >::Tetrahedron
Tetrahedron(const WorldCoordinates &p0, const WorldCoordinates &p1, const WorldCoordinates &p2, const WorldCoordinates &p3)
Construct a quadratic Tetrahedron from four nodes.
Definition: Tetrahedron.h:235
caribou::geometry::traits
Definition: Element.h:14
caribou::geometry::Tetrahedron< Quadratic >::edges
auto edges() const
Get the list of node indices of the edges.
Definition: Tetrahedron.h:243
caribou::geometry::Tetrahedron< Quadratic >::Tetrahedron
Tetrahedron(WorldCoordinates &p0, WorldCoordinates &p1, WorldCoordinates &p2, WorldCoordinates &p3)
Construct a quadratic Tetrahedron from four nodes.
Definition: Tetrahedron.h:231
caribou::geometry::Element::node
auto node(const UNSIGNED_INTEGER_TYPE &index) const
Get the Node at given index.
Definition: Element.h:64
caribou::geometry::Triangle
Definition: Triangle.h:12
caribou::geometry::Tetrahedron< Linear >::edges
auto edges() const
Get the list of node indices of the edges.
Definition: Tetrahedron.h:85
Generated by doxygen 1.8.20