Skip to content

Commit f82e2e0

Browse files
alxbilgerhugtalbot
alxbilger
and
hugtalbot
authored
[FEM.Elastic] Remove complex unproductive stiffness matrix data structure (#6133)
* [FEM.Elastic] Remove complex unproductive stiffness matrix data structure * rename ElementForce to ElementGradient * explicit namespace --------- Co-authored-by: Hugo <hugo.talbot@sofa-framework.org>
1 parent 65ccbee commit f82e2e0

11 files changed

Lines changed: 44 additions & 242 deletions

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f
5151

5252
private:
5353
using trait = sofa::component::solidmechanics::fem::elastic::trait<DataTypes, ElementType>;
54-
using ElementStiffness = typename trait::ElementStiffness;
54+
using ElementHessian = typename trait::ElementHessian;
5555
using StrainDisplacement = typename trait::StrainDisplacement;
5656
using Real = typename trait::Real;
5757

@@ -69,7 +69,7 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f
6969
/**
7070
* List of precomputed element stiffness matrices
7171
*/
72-
sofa::Data<sofa::type::vector<ElementStiffness> > d_elementStiffness;
72+
sofa::Data<sofa::type::vector<ElementHessian> > d_elementStiffness;
7373
};
7474

7575
#if !defined(ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP)

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323

2424
#include <sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h>
2525
#include <sofa/component/solidmechanics/fem/elastic/BaseLinearElasticityFEMForceField.inl>
26+
#include <sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h>
2627
#include <sofa/component/solidmechanics/fem/elastic/impl/LameParameters.h>
2728
#include <sofa/component/solidmechanics/fem/elastic/impl/VectorTools.h>
2829

@@ -98,7 +99,7 @@ void BaseElementLinearFEMForceField<DataTypes, ElementType>::precomputeElementSt
9899
const auto elasticityTensor = makeIsotropicElasticityTensor<DataTypes::spatial_dimensions, Real>(mu, lambda);
99100

100101
const std::array<sofa::Coord_t<DataTypes>, trait::NumberOfNodesInElement> nodesCoordinates = extractNodesVectorFromGlobalVector(element, restPositionAccessor.ref());
101-
elementStiffness[elementId] = integrate<DataTypes, ElementType, trait::matrixVectorProductType>(nodesCoordinates, elasticityTensor);
102+
elementStiffness[elementId] = integrate<DataTypes, ElementType>(nodesCoordinates, elasticityTensor);
102103
});
103104
}
104105

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -138,7 +138,7 @@ class ElementCorotationalFEMForceField :
138138

139139
private:
140140
using trait = sofa::component::solidmechanics::fem::elastic::trait<DataTypes, ElementType>;
141-
using ElementForce = typename trait::ElementForce;
141+
using ElementGradient = typename trait::ElementGradient;
142142
using RotationMatrix = sofa::type::Mat<trait::spatial_dimensions, trait::spatial_dimensions, sofa::Real_t<DataTypes>>;
143143

144144

@@ -158,19 +158,19 @@ class ElementCorotationalFEMForceField :
158158
protected:
159159

160160
void beforeElementForce(const sofa::core::MechanicalParams* mparams,
161-
sofa::type::vector<ElementForce>& f,
161+
sofa::type::vector<ElementGradient>& f,
162162
const sofa::VecCoord_t<DataTypes>& x) override;
163163

164164
void computeElementsForces(
165165
const sofa::simulation::Range<std::size_t>& range,
166166
const sofa::core::MechanicalParams* mparams,
167-
sofa::type::vector<ElementForce>& f,
167+
sofa::type::vector<ElementGradient>& f,
168168
const sofa::VecCoord_t<DataTypes>& x) override;
169169

170170
void computeElementsForcesDeriv(
171171
const sofa::simulation::Range<std::size_t>& range,
172172
const sofa::core::MechanicalParams* mparams,
173-
sofa::type::vector<ElementForce>& elementForcesDeriv,
173+
sofa::type::vector<ElementGradient>& elementForcesDeriv,
174174
const sofa::VecDeriv_t<DataTypes>& nodeDx) override;
175175

176176
sofa::type::vector<RotationMatrix> m_rotations;

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -72,7 +72,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::init()
7272

7373
template <class DataTypes, class ElementType>
7474
void ElementCorotationalFEMForceField<DataTypes, ElementType>::beforeElementForce(
75-
const sofa::core::MechanicalParams* mparams, sofa::type::vector<ElementForce>& f,
75+
const sofa::core::MechanicalParams* mparams, sofa::type::vector<ElementGradient>& f,
7676
const sofa::VecCoord_t<DataTypes>& x)
7777
{
7878
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
@@ -82,7 +82,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::beforeElementForc
8282
template <class DataTypes, class ElementType>
8383
void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsForces(
8484
const sofa::simulation::Range<std::size_t>& range, const sofa::core::MechanicalParams* mparams,
85-
sofa::type::vector<ElementForce>& elementForces, const sofa::VecCoord_t<DataTypes>& nodePositions)
85+
sofa::type::vector<ElementGradient>& elementForces, const sofa::VecCoord_t<DataTypes>& nodePositions)
8686
{
8787
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
8888
auto restPositionAccessor = this->mstate->readRestPositions();
@@ -130,7 +130,7 @@ template <class DataTypes, class ElementType>
130130
void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsForcesDeriv(
131131
const sofa::simulation::Range<std::size_t>& range,
132132
const sofa::core::MechanicalParams* mparams,
133-
sofa::type::vector<ElementForce>& elementForcesDeriv,
133+
sofa::type::vector<ElementGradient>& elementForcesDeriv,
134134
const sofa::VecDeriv_t<DataTypes>& nodeDx)
135135
{
136136
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
@@ -193,7 +193,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::buildStiffnessMat
193193
{
194194
for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2)
195195
{
196-
stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); // extract the submatrix corresponding to the coupling of nodes n1 and n2
196+
stiffnessMatrix.getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); // extract the submatrix corresponding to the coupling of nodes n1 and n2
197197
dfdx(element[n1] * trait::spatial_dimensions, element[n2] * trait::spatial_dimensions) += -elementRotation * localMatrix * elementRotation_T;
198198
}
199199
}

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -62,10 +62,10 @@ class ElementLinearSmallStrainFEMForceField :
6262

6363
private:
6464
using trait = sofa::component::solidmechanics::fem::elastic::trait<DataTypes, ElementType>;
65-
using ElementStiffness = typename trait::ElementStiffness;
65+
using ElementHessian = typename trait::ElementHessian;
6666
using ElementDisplacement = typename trait::ElementDisplacement;
6767
using StrainDisplacement = typename trait::StrainDisplacement;
68-
using ElementForce = typename trait::ElementForce;
68+
using ElementGradient = typename trait::ElementGradient;
6969

7070
public:
7171
void init() override;
@@ -84,13 +84,13 @@ class ElementLinearSmallStrainFEMForceField :
8484
void computeElementsForces(
8585
const sofa::simulation::Range<std::size_t>& range,
8686
const sofa::core::MechanicalParams* mparams,
87-
sofa::type::vector<ElementForce>& elementForces,
87+
sofa::type::vector<ElementGradient>& elementForces,
8888
const sofa::VecCoord_t<DataTypes>& nodePositions) override;
8989

9090
void computeElementsForcesDeriv(
9191
const sofa::simulation::Range<std::size_t>& range,
9292
const sofa::core::MechanicalParams* mparams,
93-
sofa::type::vector<ElementForce>& elementForcesDeriv,
93+
sofa::type::vector<ElementGradient>& elementForcesDeriv,
9494
const sofa::VecDeriv_t<DataTypes>& nodeDx) override;
9595

9696
};

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ template <class DataTypes, class ElementType>
4444
void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::computeElementsForces(
4545
const sofa::simulation::Range<std::size_t>& range,
4646
const sofa::core::MechanicalParams* mparams,
47-
sofa::type::vector<ElementForce>& elementForces,
47+
sofa::type::vector<ElementGradient>& elementForces,
4848
const sofa::VecCoord_t<DataTypes>& nodePositions)
4949
{
5050
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
@@ -75,7 +75,7 @@ template <class DataTypes, class ElementType>
7575
void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::computeElementsForcesDeriv(
7676
const sofa::simulation::Range<std::size_t>& range,
7777
const sofa::core::MechanicalParams* mparams,
78-
sofa::type::vector<ElementForce>& elementForcesDeriv,
78+
sofa::type::vector<ElementGradient>& elementForcesDeriv,
7979
const sofa::VecDeriv_t<DataTypes>& nodeDx)
8080
{
8181
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
@@ -133,7 +133,7 @@ void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::buildStiffne
133133
{
134134
for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2)
135135
{
136-
stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2
136+
stiffnessMatrix.getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2
137137
dfdx(element[n1] * trait::spatial_dimensions, element[n2] * trait::spatial_dimensions) += -localMatrix;
138138
}
139139
}
@@ -169,7 +169,7 @@ void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::addKToMatrix
169169
{
170170
for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2)
171171
{
172-
stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2
172+
stiffnessMatrix.getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2
173173

174174
const auto value = (-static_cast<sofa::Real_t<DataTypes>>(kFact)) * static_cast<sofa::type::ScalarOrMatrix<LocalMatType>>(localMatrix);
175175
matrix->add(

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ class FEMForceField :
5050

5151
private:
5252
using trait = sofa::component::solidmechanics::fem::elastic::trait<DataTypes, ElementType>;
53-
using ElementForce = typename trait::ElementForce;
53+
using ElementGradient = typename trait::ElementGradient;
5454

5555
public:
5656
void init() override;
@@ -81,17 +81,17 @@ class FEMForceField :
8181
/// Methods related to addForce
8282
/// @{
8383
void computeElementsForces(const sofa::core::MechanicalParams* mparams,
84-
sofa::type::vector<ElementForce>& f,
84+
sofa::type::vector<ElementGradient>& f,
8585
const sofa::VecCoord_t<DataTypes>& x);
8686

8787
virtual void beforeElementForce(const sofa::core::MechanicalParams* mparams,
88-
sofa::type::vector<ElementForce>& f,
88+
sofa::type::vector<ElementGradient>& f,
8989
const sofa::VecCoord_t<DataTypes>& x) {}
9090

9191
virtual void computeElementsForces(
9292
const sofa::simulation::Range<std::size_t>& range,
9393
const sofa::core::MechanicalParams* mparams,
94-
sofa::type::vector<ElementForce>& f,
94+
sofa::type::vector<ElementGradient>& f,
9595
const sofa::VecCoord_t<DataTypes>& x) = 0;
9696

9797
void dispatchElementForcesToNodes(
@@ -103,15 +103,15 @@ class FEMForceField :
103103
/// Methods related to addDForce
104104
/// @{
105105
void computeElementsForcesDeriv(const sofa::core::MechanicalParams* mparams,
106-
sofa::type::vector<ElementForce>& df,
106+
sofa::type::vector<ElementGradient>& df,
107107
const sofa::VecDeriv_t<DataTypes>& dx);
108108

109109
virtual void beforeElementForceDeriv(const sofa::core::MechanicalParams* mparams) {}
110110

111111
virtual void computeElementsForcesDeriv(
112112
const sofa::simulation::Range<std::size_t>& range,
113113
const sofa::core::MechanicalParams* mparams,
114-
sofa::type::vector<ElementForce>& df,
114+
sofa::type::vector<ElementGradient>& df,
115115
const sofa::VecDeriv_t<DataTypes>& dx) = 0;
116116

117117
/**

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ void FEMForceField<DataTypes, ElementType>::addForce(
8888
template <class DataTypes, class ElementType>
8989
void FEMForceField<DataTypes, ElementType>::computeElementsForces(
9090
const sofa::core::MechanicalParams* mparams,
91-
sofa::type::vector<ElementForce>& f,
91+
sofa::type::vector<ElementGradient>& f,
9292
const sofa::VecCoord_t<DataTypes>& x)
9393
{
9494
SCOPED_TIMER("ElementForces");
@@ -159,7 +159,7 @@ void FEMForceField<DataTypes, ElementType>::addDForce(
159159

160160
template <class DataTypes, class ElementType>
161161
void FEMForceField<DataTypes, ElementType>::computeElementsForcesDeriv(
162-
const sofa::core::MechanicalParams* mparams, sofa::type::vector<ElementForce>& df,
162+
const sofa::core::MechanicalParams* mparams, sofa::type::vector<ElementGradient>& df,
163163
const sofa::VecDeriv_t<DataTypes>& dx)
164164
{
165165
SCOPED_TIMER("ElementForcesDeriv");

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h‎

Lines changed: 7 additions & 115 deletions
Original file line numberDiff line numberDiff line change
@@ -25,125 +25,16 @@
2525
#include <sofa/type/FullySymmetric4Tensor.h>
2626
#include <sofa/component/solidmechanics/fem/elastic/impl/OrthotropicElasticityTensor.h>
2727
#include <sofa/component/solidmechanics/fem/elastic/impl/StrainDisplacement.h>
28+
#include <sofa/component/solidmechanics/fem/elastic/impl/trait.h>
2829
#include <sofa/type/Mat.h>
2930

3031
#include <ostream>
3132

3233
namespace sofa::component::solidmechanics::fem::elastic
3334
{
3435

35-
/**
36-
* Specifies the type of matrix-vector product to be used with the stiffness matrix.
37-
*/
38-
enum class MatrixVectorProductType
39-
{
40-
/**
41-
* The matrix-vector product is computed using the factorization of the matrix
42-
*/
43-
Factorization,
44-
45-
/**
46-
* The matrix-vector product is computed using the dense matrix representation
47-
*/
48-
Dense
49-
};
50-
51-
/**
52-
* Represents an element stiffness matrix. It contains the dense matrix representation, but also
53-
* its factorization. Both the factorization or the dense matrix can be used for the product of the
54-
* stiffness matrix with a vector. Although it gives the same result, it has an impact on the
55-
* performances.
56-
*/
57-
template <class DataTypes, class ElementType, MatrixVectorProductType matrixVectorProductType>
58-
struct FactorizedElementStiffness
59-
{
60-
private:
61-
using Real = sofa::Real_t<DataTypes>;
62-
using FiniteElement = sofa::fem::FiniteElement<ElementType, DataTypes>;
63-
static constexpr auto NbQuadraturePoints = FiniteElement::quadraturePoints().size();
64-
static constexpr sofa::Size NumberOfIndependentElements = sofa::type::NumberOfIndependentElements<DataTypes::spatial_dimensions>;
65-
66-
/// The factors of the stiffness matrix
67-
/// @{
68-
std::array<StrainDisplacement<DataTypes, ElementType>, NbQuadraturePoints> B;
69-
OrthotropicElasticityTensor<DataTypes::spatial_dimensions, sofa::Real_t<DataTypes>> elasticityTensor;
70-
std::array<Real, NbQuadraturePoints> factors;
71-
/// @}
72-
73-
/**
74-
* The assembled stiffness matrix
75-
*/
76-
sofa::type::Mat<
77-
ElementType::NumberOfNodes * DataTypes::spatial_dimensions,
78-
ElementType::NumberOfNodes * DataTypes::spatial_dimensions,
79-
sofa::Real_t<DataTypes>> stiffnessMatrix;
80-
81-
public:
82-
const StrainDisplacement<DataTypes, ElementType>& getB(std::size_t i) const { return B[i]; }
83-
Real getFactor(std::size_t i) const { return factors[i]; }
84-
const OrthotropicElasticityTensor<DataTypes::spatial_dimensions, sofa::Real_t<DataTypes>>& getElasticityTensor() const { return elasticityTensor; }
85-
86-
void setElasticityTensor(
87-
const sofa::type::FullySymmetric4Tensor<DataTypes::spatial_dimensions, sofa::Real_t<DataTypes>>& elasticityTensor_)
88-
{
89-
this->elasticityTensor.set(elasticityTensor_);
90-
}
91-
92-
void addFactor(std::size_t quadraturePointIndex,
93-
const Real factor,
94-
const StrainDisplacement<DataTypes, ElementType>& B_)
95-
{
96-
this->factors[quadraturePointIndex] = factor;
97-
this->B[quadraturePointIndex] = B_;
98-
99-
const auto K_i = factor * B_.multTranspose(elasticityTensor.toMat() * B_.B);
100-
stiffnessMatrix += K_i;
101-
}
102-
103-
const auto& getAssembledMatrix() const { return stiffnessMatrix; }
104-
105-
using Vec = sofa::type::Vec<
106-
ElementType::NumberOfNodes * DataTypes::spatial_dimensions,
107-
sofa::Real_t<DataTypes>>;
108-
109-
inline Vec operator*(const Vec& v) const
110-
{
111-
if constexpr (matrixVectorProductType == MatrixVectorProductType::Factorization)
112-
{
113-
if constexpr (NbQuadraturePoints > 1)
114-
{
115-
Vec result;
116-
for (std::size_t i = 0; i < NbQuadraturePoints; ++i)
117-
{
118-
const auto& B = this->B[i];
119-
result += B.multTranspose(this->factors[i] * (this->elasticityTensor * (B * v)));
120-
}
121-
return result;
122-
}
123-
else
124-
{
125-
return this->B[0].multTranspose(this->factors[0] * (this->elasticityTensor * (this->B[0] * v)));
126-
}
127-
}
128-
else
129-
{
130-
return this->stiffnessMatrix * v;
131-
}
132-
}
133-
134-
friend std::ostream& operator<<(std::ostream& os, const FactorizedElementStiffness& obj)
135-
{
136-
return os << obj.getAssembledMatrix();
137-
}
138-
139-
friend std::istream& operator>>(std::istream& in, FactorizedElementStiffness& obj)
140-
{
141-
return in;
142-
}
143-
};
144-
145-
template <class DataTypes, class ElementType, MatrixVectorProductType matrixVectorProductType = MatrixVectorProductType::Dense>
146-
FactorizedElementStiffness<DataTypes, ElementType, matrixVectorProductType> integrate(
36+
template <class DataTypes, class ElementType>
37+
auto integrate(
14738
const std::array<sofa::Coord_t<DataTypes>, ElementType::NumberOfNodes>& nodesCoordinates,
14839
const sofa::type::FullySymmetric4Tensor<DataTypes::spatial_dimensions, sofa::Real_t<DataTypes>>& elasticityTensor)
14940
{
@@ -154,8 +45,9 @@ FactorizedElementStiffness<DataTypes, ElementType, matrixVectorProductType> inte
15445
static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes;
15546
static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension;
15647

157-
FactorizedElementStiffness<DataTypes, ElementType, matrixVectorProductType> K;
158-
K.setElasticityTensor(elasticityTensor);
48+
typename trait<DataTypes, ElementType>::ElementHessian K;
49+
50+
auto C = elasticityTensor.toVoigtMatSym().toMat();
15951

16052
std::size_t quadraturePointIndex = 0;
16153
for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints())
@@ -180,7 +72,7 @@ FactorizedElementStiffness<DataTypes, ElementType, matrixVectorProductType> inte
18072
dN_dq[i] = J_inv.transposed() * dN_dq_ref[i];
18173

18274
const auto B = makeStrainDisplacement<DataTypes, ElementType>(dN_dq);
183-
K.addFactor(quadraturePointIndex, weight * detJ, B);
75+
K += (weight * detJ) * B.multTranspose(C * B);
18476

18577
++quadraturePointIndex;
18678
}

0 commit comments

Comments
 (0)