Caribou
RectangularQuad.h
1 #pragma once
2 
3 #include <Caribou/config.h>
4 #include <Caribou/Geometry/BaseRectangularQuad.h>
5 #include <Caribou/Geometry/Quad.h>
6 #include <Caribou/Geometry/Segment.h>
7 #include <Eigen/Core>
8 #include <array>
9 
10 namespace caribou::geometry {
11 
12 template<UNSIGNED_INTEGER_TYPE Dimension, UNSIGNED_INTEGER_TYPE Order = Linear>
13 struct RectangularQuad;
14 
15 
16 template<UNSIGNED_INTEGER_TYPE _Dimension>
17 struct traits<RectangularQuad <_Dimension, Linear>> {
18  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 2;
19  static constexpr UNSIGNED_INTEGER_TYPE Dimension = _Dimension;
20  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 4;
21  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 4;
22 
23  using BoundaryElementType = Segment<_Dimension, Linear>;
24  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 4;
25 };
26 
45 template<UNSIGNED_INTEGER_TYPE _Dimension>
46 struct RectangularQuad <_Dimension, Linear> : public BaseRectangularQuad<RectangularQuad <_Dimension, Linear>> {
47  // Types
48  using Base = BaseRectangularQuad<RectangularQuad <_Dimension, Linear>>;
49  using LocalCoordinates = typename Base::LocalCoordinates;
50  using WorldCoordinates = typename Base::WorldCoordinates;
51 
52  using GaussNode = typename Base::GaussNode;
53 
54  template <UNSIGNED_INTEGER_TYPE Dim>
55  using Vector = typename Base::template Vector<Dim>;
56 
57  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
58  using Matrix = typename Base::template Matrix<Rows, Cols>;
59 
60  // Constants
61  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
62  static constexpr auto Dimension = Base::Dimension;
63  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
64  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
65 
66  static constexpr auto canonical_nodes = Quad<Dimension, Linear>::canonical_nodes;
67 
68  using Size = typename Base::Size;
69  using Rotation = typename Base::Rotation;
70 
71  // Constructors
72  using Base::Base;
73 
74  // Constructor from a regular quad
75  explicit RectangularQuad(const Quad<Dimension, Linear> & quad) {
76  this->p_center = quad.center();
77  this->p_H = Size((quad.node(1)-quad.node(0)).norm(), (quad.node(3)-quad.node(0)).norm());
78  this->p_R = quad.frame({0, 0});
79  };
80 
82  template<typename EigenType, REQUIRES(EigenType::RowsAtCompileTime == NumberOfNodesAtCompileTime)>
83  explicit RectangularQuad(Eigen::EigenBase<EigenType> & nodes)
84  : RectangularQuad(Quad<Dimension, Linear>(nodes)) {}
85 
87  template<typename EigenType, REQUIRES(EigenType::RowsAtCompileTime == NumberOfNodesAtCompileTime)>
88  explicit RectangularQuad(const Eigen::EigenBase<EigenType> & nodes)
89  : RectangularQuad(Quad<Dimension, Linear>(nodes)) {}
90 
91 private:
92  // Implementations
93  friend struct Element<RectangularQuad <_Dimension, Linear>>;
94  friend struct BaseQuad<RectangularQuad <_Dimension, Linear>>;
95  inline auto get_L(const LocalCoordinates & xi) const -> Vector<NumberOfNodesAtCompileTime> {
96  const auto & u = xi[0];
97  const auto & v = xi[1];
98  return {
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  static_cast<FLOATING_POINT_TYPE>(1 / 4. * (1 + u) * (1 + v)),
102  static_cast<FLOATING_POINT_TYPE>(1 / 4. * (1 - u) * (1 + v))
103  };
104  };
105 
106  inline auto get_dL(const LocalCoordinates & xi) const {
107  const auto & u = xi[0];
108  const auto & v = xi[1];
109 
110  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
111  // dL/du dL/dv
112  m << -1 / 4. * (1 - v), -1 / 4. * (1 - u), // Node 0
113  +1 / 4. * (1 - v), -1 / 4. * (1 + u), // Node 1
114  +1 / 4. * (1 + v), +1 / 4. * (1 + u), // Node 2
115  -1 / 4. * (1 + v), +1 / 4. * (1 - u); // Node 3
116  return m;
117  };
118 
119  inline auto get_node(const UNSIGNED_INTEGER_TYPE & index) const -> WorldCoordinates {
120  caribou_assert(index < NumberOfNodesAtCompileTime and "Trying to get a node from an invalid node index.");
121  return this->T(LocalCoordinates(canonical_nodes[index][0], canonical_nodes[index][1]));
122  }
123 
124  inline auto get_nodes() const {
125  Matrix<NumberOfNodesAtCompileTime, Dimension> nodes;
126  for (Eigen::Index node_id = 0; node_id < NumberOfNodesAtCompileTime; ++node_id) {
127  nodes.row(node_id) = get_node(node_id);
128  }
129  return nodes;
130  }
131 
132  inline auto get_gauss_nodes() const -> const auto & {
133  using Weight = FLOATING_POINT_TYPE;
134  static const std::vector<GaussNode> gauss_nodes {
135  GaussNode {LocalCoordinates(-1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 0
136  GaussNode {LocalCoordinates(+1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 1
137  GaussNode {LocalCoordinates(+1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 2
138  GaussNode {LocalCoordinates(-1/sqrt(3), +1/sqrt(3)), Weight(1)} // Node 3
139  };
140  return gauss_nodes;
141  }
142 
143  inline auto get_boundary_elements_nodes() const -> const auto & {
144  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 2>, 4> indices = {{
145  {0, 1}, // Edge 0
146  {1, 2}, // Edge 1
147  {2, 3}, // Edge 2
148  {3, 0} // Edge 3
149  }};
150  return indices;
151  }
152 };
153 
154 
155 template<UNSIGNED_INTEGER_TYPE _Dimension>
156 struct traits<RectangularQuad <_Dimension, Quadratic>> {
157  static constexpr UNSIGNED_INTEGER_TYPE CanonicalDimension = 2;
158  static constexpr UNSIGNED_INTEGER_TYPE Dimension = _Dimension;
159  static constexpr INTEGER_TYPE NumberOfNodesAtCompileTime = 8;
160  static constexpr INTEGER_TYPE NumberOfGaussNodesAtCompileTime = 8;
161 
162  using BoundaryElementType = Segment<_Dimension, Quadratic>;
163  static constexpr INTEGER_TYPE NumberOfBoundaryElementsAtCompileTime = 4;
164 };
165 
184 template<UNSIGNED_INTEGER_TYPE _Dimension>
185 struct RectangularQuad <_Dimension, Quadratic> : public BaseRectangularQuad<RectangularQuad <_Dimension, Quadratic>> {
186  // Types
187  using Base = BaseRectangularQuad<RectangularQuad <_Dimension, Quadratic>>;
188  using LocalCoordinates = typename Base::LocalCoordinates;
189  using WorldCoordinates = typename Base::WorldCoordinates;
190 
191  using GaussNode = typename Base::GaussNode;
192 
193  template <UNSIGNED_INTEGER_TYPE Dim>
194  using Vector = typename Base::template Vector<Dim>;
195 
196  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
197  using Matrix = typename Base::template Matrix<Rows, Cols>;
198 
199  // Constants
200  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
201  static constexpr auto Dimension = Base::Dimension;
202  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
203  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
204 
205  static constexpr auto canonical_nodes = Quad<Dimension, Quadratic>::canonical_nodes;
206 
207  using Size = typename Base::Size;
208  using Rotation = typename Base::Rotation;
209 
210  // Constructors
211  using Base::Base;
212 
213  // Constructor from a regular quad
214  template<UNSIGNED_INTEGER_TYPE Order>
215  explicit RectangularQuad(const Quad<Dimension, Order> & quad) {
216  this->p_center = quad.center();
217  this->p_H = Size((quad.node(1)-quad.node(0)).norm(), (quad.node(3)-quad.node(0)).norm());
218  this->p_R = quad.frame({0, 0});
219  }
220 
222  template<typename EigenType, REQUIRES(EigenType::RowsAtCompileTime == NumberOfNodesAtCompileTime)>
223  explicit RectangularQuad(Eigen::EigenBase<EigenType> & nodes)
224  : RectangularQuad(Quad<Dimension, Quadratic>(nodes)) {}
225 
227  template<typename EigenType, REQUIRES(EigenType::RowsAtCompileTime == NumberOfNodesAtCompileTime)>
228  explicit RectangularQuad(const Eigen::EigenBase<EigenType> & nodes)
229  : RectangularQuad(Quad<Dimension, Quadratic>(nodes)) {}
230 
231 private:
232  // Implementations
233  friend struct Element<RectangularQuad <_Dimension, Quadratic>>;
234  friend struct BaseQuad<RectangularQuad <_Dimension, Quadratic>>;
235  inline auto get_L(const LocalCoordinates & xi) const {
236  const auto & u = xi[0];
237  const auto & v = xi[1];
238 
239  Vector<NumberOfNodesAtCompileTime> m;
240  m[4] = 1/2.*(1 - u*u)*(1 - v); // Node 4
241  m[5] = 1/2.*(1 - v*v)*(1 + u); // Node 5
242  m[6] = 1/2.*(1 - u*u)*(1 + v); // Node 6
243  m[7] = 1/2.*(1 - v*v)*(1 - u); // Node 7
244 
245  m[0] = 1/4.*(1 - u)*(1 - v) - 1/2.*(m[4] + m[7]); // Node 0
246  m[1] = 1/4.*(1 + u)*(1 - v) - 1/2.*(m[4] + m[5]); // Node 1
247  m[2] = 1/4.*(1 + u)*(1 + v) - 1/2.*(m[5] + m[6]); // Node 2
248  m[3] = 1/4.*(1 - u)*(1 + v) - 1/2.*(m[6] + m[7]); // Node 3
249 
250  return m;
251  };
252 
253  inline auto get_dL(const LocalCoordinates & xi) const {
254  const auto & u = xi[0];
255  const auto & v = xi[1];
256 
257  Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> m;
258  // dL/du dL/dv
259  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
260  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
261  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
262  -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
263  -u*(1 - v) , -1/2.*(1 - u*u) , // Node 4
264  1/2.*(1 - v*v) , -v*(1 + u) , // Node 5
265  -u*(1 + v) , 1/2.*(1 - u*u) , // Node 6
266  -1/2.*(1 - v*v) , -v*(1 - u); // Node 7
267 
268  return m;
269  };
270 
271  inline auto get_node(const UNSIGNED_INTEGER_TYPE & index) const -> WorldCoordinates {
272  caribou_assert(index < NumberOfNodesAtCompileTime and "Trying to get a node from an invalid node index.");
273  return this->T(LocalCoordinates(canonical_nodes[index][0], canonical_nodes[index][1]));
274  }
275 
276  inline auto get_nodes() const {
277  Matrix<NumberOfNodesAtCompileTime, Dimension> nodes;
278  for (UNSIGNED_INTEGER_TYPE node_id = 0; node_id < NumberOfNodesAtCompileTime; ++node_id) {
279  nodes.row(node_id) = get_node(node_id);
280  }
281  return nodes;
282  }
283 
284  inline auto get_gauss_nodes() const -> const auto & {
285  using Weight = FLOATING_POINT_TYPE;
286  static const std::vector<GaussNode> gauss_nodes {
287  GaussNode {LocalCoordinates(-1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 0
288  GaussNode {LocalCoordinates(+1/sqrt(3), -1/sqrt(3)), Weight(1)}, // Node 1
289  GaussNode {LocalCoordinates(+1/sqrt(3), +1/sqrt(3)), Weight(1)}, // Node 2
290  GaussNode {LocalCoordinates(-1/sqrt(3), +1/sqrt(3)), Weight(1)} // Node 3
291  };
292  return gauss_nodes;
293  }
294 
295  inline auto get_boundary_elements_nodes() const -> const auto & {
296  static const std::array<std::array<UNSIGNED_INTEGER_TYPE, 3>, 4> indices = {{
297  {0, 1, 4}, // Edge 0
298  {1, 2, 5}, // Edge 1
299  {2, 3, 6}, // Edge 2
300  {3, 0, 7} // Edge 3
301  }};
302  return indices;
303  }
304 };
305 
306 }
caribou::geometry::Element< RectangularQuad< _Dimension, Linear > >
caribou::geometry::BaseQuad
Definition: BaseQuad.h:11
caribou::geometry::BaseRectangularQuad
Definition: BaseRectangularQuad.h:12
caribou::geometry::RectangularQuad< _Dimension, Quadratic >::RectangularQuad
RectangularQuad(Eigen::EigenBase< EigenType > &nodes)
Constructor from an Eigen matrix containing the positions of the quad's nodes.
Definition: RectangularQuad.h:223
caribou::geometry::RectangularQuad
Definition: RectangularQuad.h:13
caribou::geometry::Quad
Definition: Quad.h:12
caribou::geometry::traits
Definition: Element.h:14
caribou::geometry::RectangularQuad< _Dimension, Quadratic >::RectangularQuad
RectangularQuad(const Eigen::EigenBase< EigenType > &nodes)
Constructor from an Eigen matrix containing the positions of the quad's nodes.
Definition: RectangularQuad.h:228
caribou::geometry::RectangularQuad< _Dimension, Linear >::RectangularQuad
RectangularQuad(const Eigen::EigenBase< EigenType > &nodes)
Constructor from an Eigen matrix containing the positions of the quad's nodes.
Definition: RectangularQuad.h:88
caribou::geometry::Segment< _Dimension, Linear >
Linear segment.
Definition: Segment.h:30
caribou::geometry::Segment< _Dimension, Quadratic >
Quadratic segment.
Definition: Segment.h:112
caribou::geometry::RectangularQuad< _Dimension, Linear >::RectangularQuad
RectangularQuad(Eigen::EigenBase< EigenType > &nodes)
Constructor from an Eigen matrix containing the positions of the quad's nodes.
Definition: RectangularQuad.h:83
Generated by doxygen 1.8.20