Caribou
HyperelasticForcefield.h
1 #pragma once
2 
3 #include <functional>
4 #include <array>
5 
6 #include <sofa/version.h>
7 #include <sofa/core/behavior/ForceField.h>
8 #include <sofa/core/topology/BaseMeshTopology.h>
9 
10 #include <Caribou/config.h>
11 #include <Caribou/constants.h>
12 #include <Caribou/Geometry/Element.h>
13 #include <Caribou/Geometry/Triangle.h>
14 #include <Caribou/Geometry/Quad.h>
15 #include <Caribou/Geometry/Tetrahedron.h>
16 #include <Caribou/Geometry/Hexahedron.h>
17 
18 #include <SofaCaribou/Material/HyperelasticMaterial.h>
19 
20 #include <Eigen/Sparse>
21 #include <Eigen/Dense>
22 
23 #if (defined(SOFA_VERSION) && SOFA_VERSION < 201299)
24 namespace sofa { using Index = unsigned int; }
25 #endif
26 
27 namespace SofaCaribou::forcefield {
28 
29 using namespace sofa::core;
30 using namespace sofa::core::objectmodel;
31 using namespace sofa::core::behavior;
32 
33 // Traits to get the Sofa vector type from the dimension
34 template <std::size_t Dim> struct SofaVecType {};
35 template <> struct SofaVecType<1> { using Type = sofa::defaulttype::Vec1Types; };
36 template <> struct SofaVecType<2> { using Type = sofa::defaulttype::Vec2Types; };
37 template <> struct SofaVecType<3> { using Type = sofa::defaulttype::Vec3Types; };
38 
39 template <typename GaussNode, INTEGER_TYPE NumberOfGaussNodesAtCompileTime>
41  using Type = std::array<GaussNode, static_cast<std::size_t>(NumberOfGaussNodesAtCompileTime)>;
42 };
43 
44 template <typename GaussNode>
45 struct GaussContainer<GaussNode, caribou::Dynamic> {
46  using Type = std::vector<GaussNode>;
47 };
48 
49 template <typename Element>
50 class HyperelasticForcefield : public ForceField<typename SofaVecType<caribou::geometry::traits<Element>::Dimension>::Type> {
51 public:
52  SOFA_CLASS(SOFA_TEMPLATE(HyperelasticForcefield, Element), SOFA_TEMPLATE(ForceField, typename SofaVecType<caribou::geometry::traits<Element>::Dimension>::Type));
53 
54  // Type definitions
55  using DataTypes = typename SofaVecType<caribou::geometry::traits<Element>::Dimension>::Type;
56  using Inherit = ForceField<DataTypes>;
57  using VecCoord = typename DataTypes::VecCoord;
58  using VecDeriv = typename DataTypes::VecDeriv;
59  using Coord = typename DataTypes::Coord;
60  using Deriv = typename DataTypes::Deriv;
61  using Real = typename Coord::value_type;
62 
63  using LocalCoordinates = typename caribou::geometry::Element<Element>::LocalCoordinates;
64  using WorldCoordinates = typename caribou::geometry::Element<Element>::WorldCoordinates;
65 
66  static constexpr INTEGER_TYPE Dimension = caribou::geometry::traits<Element>::Dimension;
67  static constexpr INTEGER_TYPE NumberOfNodes = caribou::geometry::traits<Element>::NumberOfNodesAtCompileTime;
68  static constexpr INTEGER_TYPE NumberOfGaussNodes = caribou::geometry::traits<Element>::NumberOfGaussNodesAtCompileTime;
69 
70  template<int nRows, int nColumns>
71  using Matrix = typename caribou::geometry::Element<Element>::template Matrix<nRows, nColumns>;
72 
73  template<int nRows, int nColumns, int Options=0>
74  using MatrixI = typename caribou::geometry::Element<Element>::template MatrixI<nRows, nColumns>;
75 
76  template<int nRows, int nColumns>
77  using Map = Eigen::Map<const Matrix<nRows, nColumns>>;
78 
79  template<int nRows>
80  using Vector = typename caribou::geometry::Element<Element>::template Vector<nRows>;
81 
82  template<int nRows>
83  using MapVector = Eigen::Map<const Vector<nRows>>;
84 
85  using Mat33 = Matrix<3, 3>;
86  using Vec3 = Vector<3>;
87 
88  template <typename ObjectType>
89  using Link = SingleLink<HyperelasticForcefield<Element>, ObjectType, BaseLink::FLAG_STRONGLINK>;
90 
91  // Data structures
92  struct GaussNode {
93  Real weight;
94  Real jacobian_determinant;
95  Matrix<NumberOfNodes, Dimension> dN_dx;
96  Mat33 F = Mat33::Identity(); // Deformation gradient
97  };
98 
99  // Container of Gauss integration nodes for one Element
100  using GaussContainer = typename GaussContainer<GaussNode, NumberOfGaussNodes>::Type;
101 
102  // Public methods
103 
105 
106  [[nodiscard]] auto
107  getTemplateName() const -> std::string override {
108  return templateName(this);
109  }
110 
111  static auto templateName(const HyperelasticForcefield<Element>* = nullptr) -> std::string {
112  return "Unknown";
113  }
114  static auto canCreate(HyperelasticForcefield<Element>* o, BaseContext* context, BaseObjectDescription* arg) -> bool;
115 
116  void init() override;
117 
118  void addForce(
119  const MechanicalParams* mparams,
120  Data<VecDeriv>& d_f,
121  const Data<VecCoord>& d_x,
122  const Data<VecDeriv>& d_v) override;
123 
124  void addDForce(
125  const MechanicalParams* /*mparams*/,
126  Data<VecDeriv>& /*d_df*/,
127  const Data<VecDeriv>& /*d_dx*/) override;
128 
129  SReal getPotentialEnergy(
130  const MechanicalParams* /* mparams */,
131  const Data<VecCoord>& /* d_x */) const override;
132 
133  void addKToMatrix(sofa::defaulttype::BaseMatrix * /*matrix*/, SReal /*kFact*/, unsigned int & /*offset*/) override;
134 
135  void computeBBox(const sofa::core::ExecParams* params, bool onlyVisible) override;
136 
137  void draw(const sofa::core::visual::VisualParams* vparams) override;
138 
140  [[nodiscard]] inline
141  virtual auto number_of_elements() const -> std::size_t {
142  return 0;
143  }
144 
146  auto gauss_nodes_of(std::size_t element_id) const -> const auto & {
147  return p_elements_quadrature_nodes[element_id];
148  }
149 
151  auto K() -> Eigen::SparseMatrix<Real> {
152 
153  // K is symmetric, so we only stored "one side" of the matrix.
154  // But to accelerate the computation, coefficients were not
155  // stored only in the upper or lower triangular part, but instead
156  // in whatever triangular part (upper or lower) the first node
157  // index of the element was. This means that a coefficient (i,j)
158  // might be in the lower triangular part, while (k,l) is in the
159  // upper triangular part. But no coefficient will be both in the
160  // lower AND the upper part.
161 
162  std::vector<Eigen::Triplet<Real>> triplets;
163  triplets.reserve(p_K.size()*2);
164  for (int k = 0; k < p_K.outerSize(); ++k) {
165  for (typename Eigen::SparseMatrix<Real>::InnerIterator it(p_K, k); it; ++it) {
166  const auto i = it.row();
167  const auto j = it.col();
168  const auto v = it.value();
169  if (i != j) {
170  triplets.emplace_back(i, j, v);
171  triplets.emplace_back(j, i, v);
172  } else {
173  triplets.emplace_back(i, i, v);
174  }
175  }
176  }
177 
178  Eigen::SparseMatrix<Real> K;
179  K.resize(p_K.rows(), p_K.cols());
180  K.setFromTriplets(triplets.begin(), triplets.end());
181 
182  return K;
183  }
184 
186  auto eigenvalues() -> const Vector<Eigen::Dynamic> &;
187 
189  auto cond() -> Real;
190 
191 protected:
192  // These protected methods are implemented but can be overridden
201  inline
202  static auto mesh_is_compatible(const sofa::core::topology::BaseMeshTopology *) -> bool {
203  return false;
204  }
205 
206 private:
207 
208  // These private methods are implemented but can be overridden
209 
211  virtual auto get_element_nodes_indices(const std::size_t &) const -> const sofa::Index * {
212  return nullptr;
213  }
214 
216  virtual void initialize_elements();
217 
219  virtual void update_stiffness();
220 
222  virtual auto get_gauss_nodes(const std::size_t & element_id, const Element & element) const -> GaussContainer;
223 
224  // Data members
225  Link<sofa::core::topology::BaseMeshTopology> d_topology_container;
226  Link<material::HyperelasticMaterial<DataTypes>> d_material;
227  Data<bool> d_enable_multithreading;
228  Data<double> d_drawScale;
229 
230  // Private variables
231  std::vector<GaussContainer> p_elements_quadrature_nodes;
232  Eigen::SparseMatrix<Real> p_K;
233  Eigen::Matrix<Real, Eigen::Dynamic, 1> p_eigenvalues;
234  bool K_is_up_to_date = false;
235  bool eigenvalues_are_up_to_date = false;
236 };
237 
238 // Tetrahedron specialization
239 template <> auto HyperelasticForcefield<caribou::geometry::Tetrahedron < caribou::Linear>>::number_of_elements() const -> std::size_t;
240 template <> auto HyperelasticForcefield<caribou::geometry::Tetrahedron < caribou::Linear>>::mesh_is_compatible(const sofa::core::topology::BaseMeshTopology * topology) -> bool;
241 template <> auto HyperelasticForcefield<caribou::geometry::Tetrahedron < caribou::Linear>>::get_element_nodes_indices(const std::size_t & element_id) const -> const sofa::Index *;
242 template <> auto HyperelasticForcefield<caribou::geometry::Tetrahedron < caribou::Linear>>::templateName(const HyperelasticForcefield<caribou::geometry::Tetrahedron < caribou::Linear>> *) -> std::string;
243 extern template class HyperelasticForcefield<caribou::geometry::Tetrahedron < caribou::Linear>>;
244 
245 // Hexahedron specialization
246 template <> auto HyperelasticForcefield<caribou::geometry::Hexahedron < caribou::Linear>>::number_of_elements() const -> std::size_t;
247 template <> auto HyperelasticForcefield<caribou::geometry::Hexahedron < caribou::Linear>>::mesh_is_compatible(const sofa::core::topology::BaseMeshTopology * topology) -> bool;
248 template <> auto HyperelasticForcefield<caribou::geometry::Hexahedron < caribou::Linear>>::get_element_nodes_indices(const std::size_t & element_id) const -> const sofa::Index *;
249 template <> auto HyperelasticForcefield<caribou::geometry::Hexahedron < caribou::Linear>>::templateName(const HyperelasticForcefield<caribou::geometry::Hexahedron < caribou::Linear>> *) -> std::string;
250 extern template class HyperelasticForcefield<caribou::geometry::Hexahedron < caribou::Linear>>;
251 
252 } // namespace SofaCaribou::forcefield
SofaCaribou::forcefield::HyperelasticForcefield::gauss_nodes_of
auto gauss_nodes_of(std::size_t element_id) const -> const auto &
Get the set of Gauss integration nodes of an element.
Definition: HyperelasticForcefield.h:146
SofaCaribou::forcefield::HyperelasticForcefield::mesh_is_compatible
static auto mesh_is_compatible(const sofa::core::topology::BaseMeshTopology *) -> bool
Return true if the mesh topology is compatible with the type Element.
Definition: HyperelasticForcefield.h:202
SofaCaribou::forcefield::HyperelasticForcefield
Definition: HyperelasticForcefield.h:50
caribou::geometry::Element
Definition: Element.h:28
SofaCaribou::forcefield::GaussContainer
Definition: HyperelasticForcefield.h:40
caribou::geometry::traits
Definition: Element.h:14
SofaCaribou::forcefield::HyperelasticForcefield::number_of_elements
virtual auto number_of_elements() const -> std::size_t
Get the number of elements contained in this field.
Definition: HyperelasticForcefield.h:141
SofaCaribou::forcefield::HyperelasticForcefield::GaussNode
Definition: HyperelasticForcefield.h:92
SofaCaribou::forcefield::HyperelasticForcefield::K
auto K() -> Eigen::SparseMatrix< Real >
Get the complete tangent stiffness matrix as a compressed sparse matrix.
Definition: HyperelasticForcefield.h:151