Caribou
BaseHexahedron.h
1 #pragma once
2 
3 #include <Caribou/config.h>
4 #include <Caribou/constants.h>
5 #include <Caribou/Geometry/Element.h>
6 #include <Eigen/Core>
7 
8 
9 namespace caribou::geometry {
10 
15 template<typename Derived>
16 struct BaseHexahedron : public Element<Derived> {
17  // Types
18  using Base = Element<Derived>;
19 
20  using LocalCoordinates = typename Base::LocalCoordinates;
21  using WorldCoordinates = typename Base::WorldCoordinates;
22 
23  using GaussNode = typename Base::GaussNode;
24 
25  template <UNSIGNED_INTEGER_TYPE Dim>
26  using Vector = typename Base::template Vector<Dim>;
27 
28  template <UNSIGNED_INTEGER_TYPE Rows, UNSIGNED_INTEGER_TYPE Cols>
29  using Matrix = typename Base::template Matrix<Rows, Cols>;
30 
31  // Constants
32  static constexpr auto CanonicalDimension = Base::CanonicalDimension;
33  static constexpr auto Dimension = Base::Dimension;
34  static constexpr auto NumberOfNodesAtCompileTime = Base::NumberOfNodesAtCompileTime;
35  static constexpr auto NumberOfGaussNodesAtCompileTime = Base::NumberOfGaussNodesAtCompileTime;
36 
37  static_assert(Dimension == 3, "Hexahedrons can only be of dimension 3.");
38 
40  BaseHexahedron() = default;
41 
43  template<typename EigenType, REQUIRES(EigenType::RowsAtCompileTime == NumberOfNodesAtCompileTime)>
44  explicit BaseHexahedron(Eigen::EigenBase<EigenType> & nodes) :p_nodes(nodes.derived().template cast<typename Base::Scalar>()) {}
45 
47  template<typename EigenType, REQUIRES(EigenType::RowsAtCompileTime == NumberOfNodesAtCompileTime)>
48  explicit BaseHexahedron(const Eigen::EigenBase<EigenType> & nodes) :p_nodes(nodes.derived().template cast<typename Base::Scalar>()) {}
49 
51  template <
52  typename ...Nodes,
53  REQUIRES(NumberOfNodesAtCompileTime == sizeof...(Nodes)+1)
54  >
55  explicit BaseHexahedron(const WorldCoordinates & first_node, Nodes&&...remaining_nodes)
56  {
57  construct_from_nodes<0>(first_node, std::forward<Nodes>(remaining_nodes)...);
58  }
59 
60  // Public methods common to all hexahedral types
61 
66  inline auto faces() const {
67  return self().get_boundary_elements_nodes();
68  }
69 
80  inline auto frame() const -> Matrix<3, 3>
81  {
82  return frame(LocalCoordinates::Zeros());
83  }
84 
103  inline auto frame(const LocalCoordinates & local_point) const -> Matrix<3, 3>
104  {
105  // Position of the point inside the hexahedron where the frame should be computed
106  const auto p = this->world_coordinates( local_point );
107 
108  // Project of the point on the quad facing the u axis
109  const auto projected_on_u = this->world_coordinates({1, local_point[1],local_point[2]});
110 
111  // Project of the point on the quad facing the v axis
112  const auto projected_on_v = this->world_coordinates({local_point[0], 1,local_point[2]});
113 
114  /* @todo(jnbrunet2000@gmail.com): select between the pairs of axis (center-to-u, center-to-v),
115  (center-to-u, center-to-w) and (center-to-v, center-to-w) to find the best match (closer to orthogonal) */
116 
117  // Vector from the point to its projection on the quad facing the u axis
118  const auto point_to_u = projected_on_u - p;
119 
120  // Vector from the point to its projection on the quad facing the v axis
121  const auto point_to_v = projected_on_v - p;
122 
123  // The u-axis of the computed frame
124  const auto u = point_to_u.normalized();
125 
126  // The v-axis of the computed frame
127  auto v = point_to_v.normalized();
128 
129  // The w-axis of the computed frame
130  const WorldCoordinates w = u.cross(v).normalized();
131 
132  // v-axis (recompute the v-axis in case u and v aren't orthogonal
133  v = w.cross(u).normalized();
134 
135  Matrix<3, 3> m;
136  m << u, v, w;
137 
138  return m;
139  }
140 
141 private:
142  // Implementations
143  friend struct Element<Derived>;
144  [[nodiscard]]
145  inline auto get_number_of_nodes() const {return NumberOfNodesAtCompileTime;}
146  inline auto get_number_of_gauss_nodes() const {return NumberOfGaussNodesAtCompileTime;}
147  inline auto get_node(const UNSIGNED_INTEGER_TYPE & index) const {return WorldCoordinates(p_nodes.row(index));};
148  inline auto get_nodes() const -> const auto & {return p_nodes;};
149  inline auto get_center() const {return Base::world_coordinates(LocalCoordinates(0, 0, 0));};
150  inline auto get_number_of_boundary_elements() const -> UNSIGNED_INTEGER_TYPE {return 6;};
151  inline auto get_contains_local(const LocalCoordinates & xi, const FLOATING_POINT_TYPE & eps) const -> bool {
152  const auto & u = xi[0];
153  const auto & v = xi[1];
154  const auto & w = xi[2];
155  return IN_CLOSED_INTERVAL(-1-eps, u, 1+eps) and
156  IN_CLOSED_INTERVAL(-1-eps, v, 1+eps) and
157  IN_CLOSED_INTERVAL(-1-eps, w, 1+eps);
158  }
159 
160  auto self() -> Derived& { return *static_cast<Derived*>(this); }
161  auto self() const -> const Derived& { return *static_cast<const Derived*>(this); }
162 
163  template <size_t index, typename ...Nodes, REQUIRES(sizeof...(Nodes) >= 1)>
164  inline
165  void construct_from_nodes(const WorldCoordinates & first_node, Nodes&&...remaining_nodes) {
166  p_nodes.row(index) = first_node;
167  construct_from_nodes<index+1>(std::forward<Nodes>(remaining_nodes)...);
168  }
169 
170  template <size_t index>
171  inline
172  void construct_from_nodes(const WorldCoordinates & last_node) {
173  p_nodes.row(index) = last_node;
174  }
175 protected:
176  Matrix<NumberOfNodesAtCompileTime, Dimension> p_nodes;
177 };
178 
179 }
caribou::geometry::BaseHexahedron::frame
auto frame() const -> Matrix< 3, 3 >
Extract the orthogonal frame of the element by computing the cross product of the unit vectors from t...
Definition: BaseHexahedron.h:80
caribou::geometry::Element< Derived >::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
Definition: Element.h:28
caribou::geometry::BaseHexahedron::BaseHexahedron
BaseHexahedron()=default
Default empty constructor.
caribou::geometry::BaseHexahedron::faces
auto faces() const
Get the list of node indices of the faces.
Definition: BaseHexahedron.h:66
caribou::geometry::BaseHexahedron::BaseHexahedron
BaseHexahedron(const WorldCoordinates &first_node, Nodes &&...remaining_nodes)
Constructor from a serie of nodes.
Definition: BaseHexahedron.h:55
caribou::geometry::BaseHexahedron::frame
auto frame(const LocalCoordinates &local_point) const -> Matrix< 3, 3 >
Extract the frame positioned at the given position (in local coordinates) on the hexa by computing th...
Definition: BaseHexahedron.h:103
caribou::geometry::Element< Derived >::nodes
auto nodes() const
Get the set of nodes.
Definition: Element.h:67
caribou::geometry::BaseHexahedron::BaseHexahedron
BaseHexahedron(const Eigen::EigenBase< EigenType > &nodes)
Constructor from an Eigen matrix containing the positions of the hexahedron's nodes.
Definition: BaseHexahedron.h:48
caribou::geometry::BaseHexahedron::BaseHexahedron
BaseHexahedron(Eigen::EigenBase< EigenType > &nodes)
Constructor from an Eigen matrix containing the positions of the hexahedron's nodes.
Definition: BaseHexahedron.h:44
caribou::geometry::BaseHexahedron
Base class for hexahedral elements.
Definition: BaseHexahedron.h:16