Caribou
Triangle.h
1 #pragma once
2 
3 #include <Caribou/config.h>
4 #include <Caribou/Geometry/BaseTriangle.h>
5 #include <Caribou/Geometry/Segment.h>
6 #include <Eigen/Core>
7 #include <array>
8 
9 namespace caribou::geometry {
10 
11 template<UNSIGNED_INTEGER_TYPE _Dimension, UNSIGNED_INTEGER_TYPE _Order = Linear>
12 struct Triangle;
13 
14 template<UNSIGNED_INTEGER_TYPE _Dimension>
15 struct traits<Triangle <_Dimension, Linear>> {
16  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 2;
17  static constexpr UNSIGNED_INTEGER_TYPE Dimension = _Dimension;
18  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 3;
19  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 1;
20 
21  using BoundaryElementType = Segment<_Dimension, Linear>;
22  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 3;
23 };
24 
43 template<UNSIGNED_INTEGER_TYPE _Dimension>
44 struct Triangle <_Dimension, Linear> : public BaseTriangle<Triangle <_Dimension, Linear>> {
45  // Types
46  using Base = BaseTriangle<Triangle <_Dimension, 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  // Constructors
65  using Base::Base;
66  Triangle() : Base() {
67  if constexpr (Dimension == 2) {
68  this->p_nodes.row(0) = WorldCoordinates(0, 0);
69  this->p_nodes.row(1) = WorldCoordinates(1, 0);
70  this->p_nodes.row(2) = WorldCoordinates(0, 1);
71  } else {
72  this->p_nodes.row(0) = WorldCoordinates(0, 0, 0);
73  this->p_nodes.row(1) = WorldCoordinates(1, 0, 0);
74  this->p_nodes.row(2) = WorldCoordinates(0, 1, 0);
75  }
76  };
77 
78 
79 private:
80  // Implementations
81  friend struct Element<Triangle <_Dimension, Linear>>;
82  friend struct BaseTriangle<Triangle <_Dimension, Linear>>;
83  inline auto get_L(const LocalCoordinates & xi) const -> Vector<NumberOfNodesAtCompileTime> {
84  const auto & u = xi[0];
85  const auto & v = xi[1];
86  return {
87  static_cast<FLOATING_POINT_TYPE>(1 - u - v),
88  static_cast<FLOATING_POINT_TYPE>(u),
89  static_cast<FLOATING_POINT_TYPE>(v)
90  };
91  };
92 
93  inline auto get_dL(const LocalCoordinates & /*xi*/) const {
94  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
95  // dL/du dL/dv
96  m << -1, -1, // Node 0
97  1, 0, // Node 1
98  0, 1; // Node 2
99  return m;
100  };
101 
102  inline auto get_gauss_nodes() const -> const auto & {
103  using Weight = FLOATING_POINT_TYPE;
104  static const std::vector<GaussNode> gauss_nodes {
105  GaussNode {LocalCoordinates(1/3., 1/3.), Weight(1/2.)} // Node 0
106  };
107  return gauss_nodes;
108  }
109 
110  inline auto get_boundary_elements_nodes() const -> const auto & {
111  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 2>, 3> indices = {{
112  {0, 1}, // Edge 0
113  {1, 2}, // Edge 1
114  {2, 0} // Edge 2
115  }};
116  return indices;
117  }
118 };
119 
120 
121 template<UNSIGNED_INTEGER_TYPE _Dimension>
122 struct traits<Triangle <_Dimension, Quadratic>> {
123  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 2;
124  static constexpr UNSIGNED_INTEGER_TYPE Dimension = _Dimension;
125  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 6;
126  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 3;
127 
128  using BoundaryElementType = Segment<_Dimension, Quadratic>;
129  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 3;
130 };
131 
150 template<UNSIGNED_INTEGER_TYPE _Dimension>
151 struct Triangle <_Dimension, Quadratic> : public BaseTriangle<Triangle <_Dimension, Quadratic>> {
152  // Types
153  using Base = BaseTriangle<Triangle <_Dimension, Quadratic>>;
154  using LocalCoordinates = typename Base::LocalCoordinates;
155  using WorldCoordinates = typename Base::WorldCoordinates;
156 
157  using GaussNode = typename Base::GaussNode;
158 
159  template <UNSIGNED_INTEGER_TYPE Dim>
160  using Vector = typename Base::template Vector<Dim>;
161 
162  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
163  using Matrix = typename Base::template Matrix<Rows, Cols>;
164 
165  // Constants
166  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
167  static constexpr auto Dimension = Base::Dimension;
168  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
169  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
170 
171  // Constructors
172  using Base::Base;
173  Triangle() : Base() {
174  if constexpr (Dimension == 2) {
175  this->p_nodes.row(0) = WorldCoordinates(0.0, 0.0); // Node 0
176  this->p_nodes.row(1) = WorldCoordinates(1.0, 0.0); // Node 1
177  this->p_nodes.row(2) = WorldCoordinates(0.0, 1.0); // Node 2
178  this->p_nodes.row(3) = WorldCoordinates(0.5, 0.0); // Node 3
179  this->p_nodes.row(4) = WorldCoordinates(0.5, 0.5); // Node 4
180  this->p_nodes.row(5) = WorldCoordinates(0.0, 0.5); // Node 5
181  } else {
182  this->p_nodes.row(0) = WorldCoordinates(0.0, 0.0, 0.0); // Node 0
183  this->p_nodes.row(1) = WorldCoordinates(1.0, 0.0, 0.0); // Node 1
184  this->p_nodes.row(2) = WorldCoordinates(0.0, 1.0, 0.0); // Node 2
185  this->p_nodes.row(3) = WorldCoordinates(0.5, 0.0, 0.0); // Node 3
186  this->p_nodes.row(4) = WorldCoordinates(0.5, 0.5, 0.0); // Node 4
187  this->p_nodes.row(5) = WorldCoordinates(0.0, 0.5, 0.0); // Node 5
188  }
189  };
190 
191  // Construct a quadratic triangle from a linear one
192  explicit Triangle(const Triangle<Dimension, Linear> & linear_triangle) : Base() {
193  this->p_nodes.row(0) = linear_triangle.node(0); // Node 0
194  this->p_nodes.row(1) = linear_triangle.node(1); // Node 1
195  this->p_nodes.row(2) = linear_triangle.node(2); // Node 2
196  this->p_nodes.row(3) = linear_triangle.world_coordinates(LocalCoordinates(0.5, 0.0)); // Node 3
197  this->p_nodes.row(4) = linear_triangle.world_coordinates(LocalCoordinates(0.5, 0.5)); // Node 4
198  this->p_nodes.row(5) = linear_triangle.world_coordinates(LocalCoordinates(0.0, 0.5)); // Node 5
199  }
200 
202  Triangle(WorldCoordinates & p0, WorldCoordinates & p1, WorldCoordinates & p2)
203  : Triangle(Triangle<Dimension, Linear>(p0, p1, p2)) {}
204 
206  Triangle(const WorldCoordinates & p0, const WorldCoordinates & p1, const WorldCoordinates & p2)
207  : Triangle(Triangle<Dimension, Linear>(p0, p1, p2)) {}
208 
209 private:
210  // Implementations
211  friend struct Element<Triangle <_Dimension, Quadratic>>;
212  friend struct BaseTriangle<Triangle <_Dimension, Quadratic>>;
213  inline auto get_L(const LocalCoordinates & xi) const {
214  const auto & u = xi[0];
215  const auto & v = xi[1];
216  const auto l = 1 - u - v;
217 
218  Vector<NumberOfNodesAtCompileTime> m;
219  m << l * (2*l - 1), // Node 0
220  u * (2*u - 1), // Node 1
221  v * (2*v - 1), // Node 2
222  4 * u * l, // Node 3
223  4 * u * v, // Node 4
224  4 * v * l; // Node 5
225 
226  return m;
227  };
228 
229  inline auto get_dL(const LocalCoordinates & xi) const {
230  const auto & u = xi[0];
231  const auto & v = xi[1];
232  const auto l = 1 - u - v;
233 
234  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
235  // dL/du dL/dv
236  m << 1 - 4 * l , 1 - 4 * l , // Node 0
237  4 * u - 1 , 0 , // Node 1
238  0 , 4 * v - 1 , // Node 2
239  4 * (l - u), -4 * u , // Node 3
240  4 * v , 4 * u , // Node 4
241  - 4 * v , 4 * (l - v); // Node 5
242  return m;
243  };
244 
245  inline auto get_gauss_nodes() const -> const auto & {
246  using Weight = FLOATING_POINT_TYPE;
247  static const std::vector<GaussNode> gauss_nodes {
248  GaussNode {LocalCoordinates(2/3., 1/6.), Weight(1/6.)}, // Node 0
249  GaussNode {LocalCoordinates(1/6., 2/3.), Weight(1/6.)}, // Node 1
250  GaussNode {LocalCoordinates(1/6., 1/6.), Weight(1/6.)} // Node 2
251  };
252  return gauss_nodes;
253  }
254 
255  inline auto get_boundary_elements_nodes() const -> const auto & {
256  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 3>, 3> indices = {{
257  {0, 1, 3}, // Edge 0
258  {1, 2, 4}, // Edge 1
259  {2, 0, 5} // Edge 2
260  }};
261  return indices;
262  }
263 };
264 
265 }
caribou::geometry::Triangle< _Dimension, Quadratic >::Triangle
Triangle(WorldCoordinates &p0, WorldCoordinates &p1, WorldCoordinates &p2)
Construct a quadratic triangle from three nodes.
Definition: Triangle.h:202
caribou::geometry::Element< Derived >
caribou::geometry::BaseTriangle
Definition: BaseTriangle.h:11
caribou::geometry::traits
Definition: Element.h:14
caribou::geometry::Segment< _Dimension, Linear >
Linear segment.
Definition: Segment.h:30
caribou::geometry::Triangle
Definition: Triangle.h:12
caribou::geometry::Triangle< _Dimension, Quadratic >::Triangle
Triangle(const WorldCoordinates &p0, const WorldCoordinates &p1, const WorldCoordinates &p2)
Construct a quadratic triangle from three nodes.
Definition: Triangle.h:206
caribou::geometry::Segment< _Dimension, Quadratic >
Quadratic segment.
Definition: Segment.h:112
Generated by doxygen 1.8.20