Caribou
Quad.h
1 #pragma once
2 
3 #include <Caribou/config.h>
4 #include <Caribou/Geometry/BaseQuad.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 Quad;
13 
14 template<UNSIGNED_INTEGER_TYPE _Dimension>
15 struct traits<Quad <_Dimension, Linear>> {
16  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 2;
17  static constexpr UNSIGNED_INTEGER_TYPE Dimension = _Dimension;
18  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 4;
19  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 4;
20 
21  using BoundaryElementType = Segment<_Dimension, Linear>;
22  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 4;
23 };
24 
43 template<UNSIGNED_INTEGER_TYPE _Dimension>
44 struct Quad <_Dimension, Linear> : public BaseQuad<Quad <_Dimension, Linear>> {
45  // Types
46  using Base = BaseQuad<Quad <_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  static constexpr FLOATING_POINT_TYPE canonical_nodes [NumberOfNodesAtCompileTime][CanonicalDimension] = {
65  // u, v
66  {-1, -1}, // Node 0
67  {+1, -1}, // Node 1
68  {+1, +1}, // Node 2
69  {-1, +1} // Node 3
70  };
71 
72  // Constructors
73  using Base::Base;
74  Quad() : Base() {
75  if constexpr (Dimension == 2) {
76  this->p_nodes.row(0) = WorldCoordinates(canonical_nodes[0][0], canonical_nodes[0][1]); // Node 0
77  this->p_nodes.row(1) = WorldCoordinates(canonical_nodes[1][0], canonical_nodes[1][1]); // Node 1
78  this->p_nodes.row(2) = WorldCoordinates(canonical_nodes[2][0], canonical_nodes[2][1]); // Node 2
79  this->p_nodes.row(3) = WorldCoordinates(canonical_nodes[3][0], canonical_nodes[3][1]); // Node 3
80  } else {
81  this->p_nodes.row(0) = WorldCoordinates(canonical_nodes[0][0], canonical_nodes[0][1], 0); // Node 0
82  this->p_nodes.row(1) = WorldCoordinates(canonical_nodes[1][0], canonical_nodes[1][1], 0); // Node 1
83  this->p_nodes.row(2) = WorldCoordinates(canonical_nodes[2][0], canonical_nodes[2][1], 0); // Node 2
84  this->p_nodes.row(3) = WorldCoordinates(canonical_nodes[3][0], canonical_nodes[3][1], 0); // Node 3
85  }
86  };
87 
88 
89 private:
90  // Implementations
91  friend struct Element<Quad <_Dimension, Linear>>;
92  friend struct BaseQuad<Quad <_Dimension, Linear>>;
93  inline auto get_L(const LocalCoordinates & xi) const -> Vector<NumberOfNodesAtCompileTime> {
94  const auto & u = xi[0];
95  const auto & v = xi[1];
96  return {
97  static_cast<FLOATING_POINT_TYPE>(1 / 4. * (1 - u) * (1 - v)),
98  static_cast<FLOATING_POINT_TYPE>(1 / 4. * (1 + u) * (1 - v)),
99  static_cast<FLOATING_POINT_TYPE>(1 / 4. * (1 + u) * (1 + v)),
100  static_cast<FLOATING_POINT_TYPE>(1 / 4. * (1 - u) * (1 + v))
101  };
102  };
103 
104  inline auto get_dL(const LocalCoordinates & xi) const {
105  const auto & u = xi[0];
106  const auto & v = xi[1];
107 
108  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
109  // dL/du dL/dv
110  m << -1 / 4. * (1 - v), -1 / 4. * (1 - u), // Node 0
111  +1 / 4. * (1 - v), -1 / 4. * (1 + u), // Node 1
112  +1 / 4. * (1 + v), +1 / 4. * (1 + u), // Node 2
113  -1 / 4. * (1 + v), +1 / 4. * (1 - u); // Node 3
114  return m;
115  };
116 
117  inline auto get_gauss_nodes() const -> const auto & {
118  using Weight = FLOATING_POINT_TYPE;
119  static const std::vector<GaussNode> gauss_nodes {
120  GaussNode {LocalCoordinates(-1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 0
121  GaussNode {LocalCoordinates(+1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 1
122  GaussNode {LocalCoordinates(+1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 2
123  GaussNode {LocalCoordinates(-1/sqrt(3), +1/sqrt(3)), Weight(1)} // Node 3
124  };
125  return gauss_nodes;
126  }
127 
128  inline auto get_boundary_elements_nodes() const -> const auto & {
129  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 2>, 4> indices = {{
130  {0, 1}, // Edge 0
131  {1, 2}, // Edge 1
132  {2, 3}, // Edge 2
133  {3, 0} // Edge 3
134  }};
135  return indices;
136  }
137 };
138 
139 
140 template<UNSIGNED_INTEGER_TYPE _Dimension>
141 struct traits<Quad <_Dimension, Quadratic>> {
142  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 2;
143  static constexpr UNSIGNED_INTEGER_TYPE Dimension = _Dimension;
144  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 8;
145  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 4;
146 
147  using BoundaryElementType = Segment<_Dimension, Quadratic>;
148  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 4;
149 };
150 
169 template<UNSIGNED_INTEGER_TYPE _Dimension>
170 struct Quad <_Dimension, Quadratic> : public BaseQuad<Quad <_Dimension, Quadratic>> {
171  // Types
172  using Base = BaseQuad<Quad <_Dimension, Quadratic>>;
173  using LocalCoordinates = typename Base::LocalCoordinates;
174  using WorldCoordinates = typename Base::WorldCoordinates;
175 
176  using GaussNode = typename Base::GaussNode;
177 
178  template <UNSIGNED_INTEGER_TYPE Dim>
179  using Vector = typename Base::template Vector<Dim>;
180 
181  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
182  using Matrix = typename Base::template Matrix<Rows, Cols>;
183 
184  // Constants
185  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
186  static constexpr auto Dimension = Base::Dimension;
187  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
188  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
189 
190  static constexpr FLOATING_POINT_TYPE canonical_nodes [NumberOfNodesAtCompileTime][CanonicalDimension] = {
191  // u, v
192  {-1, -1}, // Node 0
193  {+1, -1}, // Node 1
194  {+1, +1}, // Node 2
195  {-1, +1}, // Node 3
196  { 0, -1}, // Node 4
197  {+1, 0}, // Node 5
198  { 0, +1}, // Node 6
199  {-1, 0} // Node 7
200  };
201 
202  // Constructors
203  using Base::Base;
204  Quad() : Base() {
205  if constexpr (Dimension == 2) {
206  this->p_nodes.row(0) = WorldCoordinates(canonical_nodes[0][0], canonical_nodes[0][1]); // Node 0
207  this->p_nodes.row(1) = WorldCoordinates(canonical_nodes[1][0], canonical_nodes[1][1]); // Node 1
208  this->p_nodes.row(2) = WorldCoordinates(canonical_nodes[2][0], canonical_nodes[2][1]); // Node 2
209  this->p_nodes.row(3) = WorldCoordinates(canonical_nodes[3][0], canonical_nodes[3][1]); // Node 3
210  this->p_nodes.row(4) = WorldCoordinates(canonical_nodes[4][0], canonical_nodes[4][1]); // Node 4
211  this->p_nodes.row(5) = WorldCoordinates(canonical_nodes[5][0], canonical_nodes[5][1]); // Node 5
212  this->p_nodes.row(6) = WorldCoordinates(canonical_nodes[6][0], canonical_nodes[6][1]); // Node 6
213  this->p_nodes.row(7) = WorldCoordinates(canonical_nodes[7][0], canonical_nodes[7][1]); // Node 7
214  } else {
215  this->p_nodes.row(0) = WorldCoordinates(canonical_nodes[0][0], canonical_nodes[0][1], 0); // Node 0
216  this->p_nodes.row(1) = WorldCoordinates(canonical_nodes[1][0], canonical_nodes[1][1], 0); // Node 1
217  this->p_nodes.row(2) = WorldCoordinates(canonical_nodes[2][0], canonical_nodes[2][1], 0); // Node 2
218  this->p_nodes.row(3) = WorldCoordinates(canonical_nodes[3][0], canonical_nodes[3][1], 0); // Node 3
219  this->p_nodes.row(4) = WorldCoordinates(canonical_nodes[4][0], canonical_nodes[4][1], 0); // Node 4
220  this->p_nodes.row(5) = WorldCoordinates(canonical_nodes[5][0], canonical_nodes[5][1], 0); // Node 5
221  this->p_nodes.row(6) = WorldCoordinates(canonical_nodes[6][0], canonical_nodes[6][1], 0); // Node 6
222  this->p_nodes.row(7) = WorldCoordinates(canonical_nodes[7][0], canonical_nodes[7][1], 0); // Node 7
223  }
224  };
225 
227  explicit Quad(const Quad<Dimension, Linear> & linear_Quad) : Base() {
228  this->p_nodes.row(0) = linear_Quad.node(0); // Node 0
229  this->p_nodes.row(1) = linear_Quad.node(1); // Node 1
230  this->p_nodes.row(2) = linear_Quad.node(2); // Node 2
231  this->p_nodes.row(3) = linear_Quad.node(3); // Node 3
232  this->p_nodes.row(4) = linear_Quad.world_coordinates(LocalCoordinates(canonical_nodes[4][0], canonical_nodes[4][1])); // Node 4
233  this->p_nodes.row(5) = linear_Quad.world_coordinates(LocalCoordinates(canonical_nodes[5][0], canonical_nodes[5][1])); // Node 5
234  this->p_nodes.row(6) = linear_Quad.world_coordinates(LocalCoordinates(canonical_nodes[6][0], canonical_nodes[6][1])); // Node 6
235  this->p_nodes.row(7) = linear_Quad.world_coordinates(LocalCoordinates(canonical_nodes[7][0], canonical_nodes[7][1])); // Node 7
236  }
237 
239  Quad(WorldCoordinates & p0, WorldCoordinates & p1, WorldCoordinates & p2, WorldCoordinates & p3)
240  : Quad(Quad<Dimension, Linear>(p0, p1, p2, p3)) {}
241 
243  Quad(const WorldCoordinates & p0, const WorldCoordinates & p1, const WorldCoordinates & p2, const WorldCoordinates & p3)
244  : Quad(Quad<Dimension, Linear>(p0, p1, p2, p3)) {}
245 
246 
247 private:
248  // Implementations
249  friend struct Element<Quad <_Dimension, Quadratic>>;
250  friend struct BaseQuad<Quad <_Dimension, Quadratic>>;
251  inline auto get_L(const LocalCoordinates & xi) const {
252  const auto & u = xi[0];
253  const auto & v = xi[1];
254 
255  Vector<NumberOfNodesAtCompileTime> m;
256  m[4] = 1/2.*(1 - u*u)*(1 - v); // Node 4
257  m[5] = 1/2.*(1 - v*v)*(1 + u); // Node 5
258  m[6] = 1/2.*(1 - u*u)*(1 + v); // Node 6
259  m[7] = 1/2.*(1 - v*v)*(1 - u); // Node 7
260 
261  m[0] = 1/4.*(1 - u)*(1 - v) - 1/2.*(m[4] + m[7]); // Node 0
262  m[1] = 1/4.*(1 + u)*(1 - v) - 1/2.*(m[4] + m[5]); // Node 1
263  m[2] = 1/4.*(1 + u)*(1 + v) - 1/2.*(m[5] + m[6]); // Node 2
264  m[3] = 1/4.*(1 - u)*(1 + v) - 1/2.*(m[6] + m[7]); // Node 3
265 
266  return m;
267  };
268 
269  inline auto get_dL(const LocalCoordinates & xi) const {
270  const auto & u = xi[0];
271  const auto & v = xi[1];
272 
273  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
274  // dL/du dL/dv
275  m << -1/4.*(1 - v) - 1/2.*(-u*(1 - v) + -1/2.*(1 - v*v)) , -1/4.*(1 - u) - 1/2.*(-1/2.*(1 - u*u) + -v*(1 - u)), // Node 0
276  1/4.*(1 - v) - 1/2.*(-u*(1 - v) + 1/2.*(1 - v*v)) , -1/4.*(1 + u) - 1/2.*(-1/2.*(1 - u*u) + -v*(1 + u)), // Node 1
277  1/4.*(1 + v) - 1/2.*(-u*(1 + v) + 1/2.*(1 - v*v)) , 1/4.*(1 + u) - 1/2.*( 1/2.*(1 - u*u) + -v*(1 + u)), // Node 2
278  -1/4.*(1 + v) - 1/2.*(-u*(1 + v) + -1/2.*(1 - v*v)) , 1/4.*(1 - u) - 1/2.*( 1/2.*(1 - u*u) + -v*(1 - u)), // Node 3
279  -u*(1 - v) , -1/2.*(1 - u*u) , // Node 4
280  1/2.*(1 - v*v) , -v*(1 + u) , // Node 5
281  -u*(1 + v) , 1/2.*(1 - u*u) , // Node 6
282  -1/2.*(1 - v*v) , -v*(1 - u); // Node 7
283 
284  return m;
285  };
286 
287  inline auto get_gauss_nodes() const -> const auto & {
288  using Weight = FLOATING_POINT_TYPE;
289  static const std::vector<GaussNode> gauss_nodes {
290  GaussNode {LocalCoordinates(-1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 0
291  GaussNode {LocalCoordinates(+1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 1
292  GaussNode {LocalCoordinates(+1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 2
293  GaussNode {LocalCoordinates(-1/sqrt(3), +1/sqrt(3)), Weight(1)} // Node 3
294  };
295  return gauss_nodes;
296  }
297 
298  inline auto get_boundary_elements_nodes() const -> const auto & {
299  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 3>, 4> indices = {{
300  {0, 1, 4}, // Edge 0
301  {1, 2, 5}, // Edge 1
302  {2, 3, 6}, // Edge 2
303  {3, 0, 7} // Edge 3
304  }};
305  return indices;
306  }
307 };
308 
309 }
caribou::geometry::Quad< _Dimension, Quadratic >::Quad
Quad(const Quad< Dimension, Linear > &linear_Quad)
Construct a quadratic Quad from a linear one.
Definition: Quad.h:227
caribou::geometry::Element< Derived >
caribou::geometry::Quad< _Dimension, Quadratic >::Quad
Quad(WorldCoordinates &p0, WorldCoordinates &p1, WorldCoordinates &p2, WorldCoordinates &p3)
Construct a quadratic Quad from four nodes.
Definition: Quad.h:239
caribou::geometry::BaseQuad
Definition: BaseQuad.h:11
caribou::geometry::Quad< _Dimension, Quadratic >::Quad
Quad(const WorldCoordinates &p0, const WorldCoordinates &p1, const WorldCoordinates &p2, const WorldCoordinates &p3)
Construct a quadratic Quad from four nodes.
Definition: Quad.h:243
caribou::geometry::Quad
Definition: Quad.h:12
caribou::geometry::traits
Definition: Element.h:14
caribou::geometry::Segment< _Dimension, Linear >
Linear segment.
Definition: Segment.h:30
caribou::geometry::Segment< _Dimension, Quadratic >
Quadratic segment.
Definition: Segment.h:112
Generated by doxygen 1.8.20