diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h index 1a980373a5e..efbc9dda6db 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h @@ -51,7 +51,7 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f private: using trait = sofa::component::solidmechanics::fem::elastic::trait; - using ElementStiffness = typename trait::ElementStiffness; + using ElementHessian = typename trait::ElementHessian; using StrainDisplacement = typename trait::StrainDisplacement; using Real = typename trait::Real; @@ -69,7 +69,7 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f /** * List of precomputed element stiffness matrices */ - sofa::Data > d_elementStiffness; + sofa::Data > d_elementStiffness; }; #if !defined(ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl index de380c5e4cc..e2191f8485d 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -98,7 +99,7 @@ void BaseElementLinearFEMForceField::precomputeElementSt const auto elasticityTensor = makeIsotropicElasticityTensor(mu, lambda); const std::array, trait::NumberOfNodesInElement> nodesCoordinates = extractNodesVectorFromGlobalVector(element, restPositionAccessor.ref()); - elementStiffness[elementId] = integrate(nodesCoordinates, elasticityTensor); + elementStiffness[elementId] = integrate(nodesCoordinates, elasticityTensor); }); } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h index 484a5c89697..1425fdfa925 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h @@ -138,7 +138,7 @@ class ElementCorotationalFEMForceField : private: using trait = sofa::component::solidmechanics::fem::elastic::trait; - using ElementForce = typename trait::ElementForce; + using ElementGradient = typename trait::ElementGradient; using RotationMatrix = sofa::type::Mat>; @@ -158,19 +158,19 @@ class ElementCorotationalFEMForceField : protected: void beforeElementForce(const sofa::core::MechanicalParams* mparams, - sofa::type::vector& f, + sofa::type::vector& f, const sofa::VecCoord_t& x) override; void computeElementsForces( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& f, + sofa::type::vector& f, const sofa::VecCoord_t& x) override; void computeElementsForcesDeriv( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& elementForcesDeriv, + sofa::type::vector& elementForcesDeriv, const sofa::VecDeriv_t& nodeDx) override; sofa::type::vector m_rotations; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl index a3a74a36a81..91a3d89bec0 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl @@ -72,7 +72,7 @@ void ElementCorotationalFEMForceField::init() template void ElementCorotationalFEMForceField::beforeElementForce( - const sofa::core::MechanicalParams* mparams, sofa::type::vector& f, + const sofa::core::MechanicalParams* mparams, sofa::type::vector& f, const sofa::VecCoord_t& x) { const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); @@ -82,7 +82,7 @@ void ElementCorotationalFEMForceField::beforeElementForc template void ElementCorotationalFEMForceField::computeElementsForces( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& elementForces, const sofa::VecCoord_t& nodePositions) + sofa::type::vector& elementForces, const sofa::VecCoord_t& nodePositions) { const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); auto restPositionAccessor = this->mstate->readRestPositions(); @@ -130,7 +130,7 @@ template void ElementCorotationalFEMForceField::computeElementsForcesDeriv( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& elementForcesDeriv, + sofa::type::vector& elementForcesDeriv, const sofa::VecDeriv_t& nodeDx) { const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); @@ -193,7 +193,7 @@ void ElementCorotationalFEMForceField::buildStiffnessMat { for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2) { - stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); // extract the submatrix corresponding to the coupling of nodes n1 and n2 + stiffnessMatrix.getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); // extract the submatrix corresponding to the coupling of nodes n1 and n2 dfdx(element[n1] * trait::spatial_dimensions, element[n2] * trait::spatial_dimensions) += -elementRotation * localMatrix * elementRotation_T; } } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h index db8436d83fc..853867ff38e 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h @@ -62,10 +62,10 @@ class ElementLinearSmallStrainFEMForceField : private: using trait = sofa::component::solidmechanics::fem::elastic::trait; - using ElementStiffness = typename trait::ElementStiffness; + using ElementHessian = typename trait::ElementHessian; using ElementDisplacement = typename trait::ElementDisplacement; using StrainDisplacement = typename trait::StrainDisplacement; - using ElementForce = typename trait::ElementForce; + using ElementGradient = typename trait::ElementGradient; public: void init() override; @@ -84,13 +84,13 @@ class ElementLinearSmallStrainFEMForceField : void computeElementsForces( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& elementForces, + sofa::type::vector& elementForces, const sofa::VecCoord_t& nodePositions) override; void computeElementsForcesDeriv( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& elementForcesDeriv, + sofa::type::vector& elementForcesDeriv, const sofa::VecDeriv_t& nodeDx) override; }; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl index fc58db7abe4..8967dc3d4a2 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl @@ -44,7 +44,7 @@ template void ElementLinearSmallStrainFEMForceField::computeElementsForces( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& elementForces, + sofa::type::vector& elementForces, const sofa::VecCoord_t& nodePositions) { const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); @@ -75,7 +75,7 @@ template void ElementLinearSmallStrainFEMForceField::computeElementsForcesDeriv( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& elementForcesDeriv, + sofa::type::vector& elementForcesDeriv, const sofa::VecDeriv_t& nodeDx) { const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); @@ -133,7 +133,7 @@ void ElementLinearSmallStrainFEMForceField::buildStiffne { for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2) { - stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2 + stiffnessMatrix.getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2 dfdx(element[n1] * trait::spatial_dimensions, element[n2] * trait::spatial_dimensions) += -localMatrix; } } @@ -169,7 +169,7 @@ void ElementLinearSmallStrainFEMForceField::addKToMatrix { for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2) { - stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2 + stiffnessMatrix.getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2 const auto value = (-static_cast>(kFact)) * static_cast>(localMatrix); matrix->add( diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h index c5eb9f4a3ee..85043d0aa8e 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h @@ -50,7 +50,7 @@ class FEMForceField : private: using trait = sofa::component::solidmechanics::fem::elastic::trait; - using ElementForce = typename trait::ElementForce; + using ElementGradient = typename trait::ElementGradient; public: void init() override; @@ -81,17 +81,17 @@ class FEMForceField : /// Methods related to addForce /// @{ void computeElementsForces(const sofa::core::MechanicalParams* mparams, - sofa::type::vector& f, + sofa::type::vector& f, const sofa::VecCoord_t& x); virtual void beforeElementForce(const sofa::core::MechanicalParams* mparams, - sofa::type::vector& f, + sofa::type::vector& f, const sofa::VecCoord_t& x) {} virtual void computeElementsForces( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& f, + sofa::type::vector& f, const sofa::VecCoord_t& x) = 0; void dispatchElementForcesToNodes( @@ -103,13 +103,13 @@ class FEMForceField : /// Methods related to addDForce /// @{ void computeElementsForcesDeriv(const sofa::core::MechanicalParams* mparams, - sofa::type::vector& df, + sofa::type::vector& df, const sofa::VecDeriv_t& dx); virtual void computeElementsForcesDeriv( const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, - sofa::type::vector& df, + sofa::type::vector& df, const sofa::VecDeriv_t& dx) = 0; /** diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl index bdae03ead2b..7ec3cc4f106 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl @@ -88,7 +88,7 @@ void FEMForceField::addForce( template void FEMForceField::computeElementsForces( const sofa::core::MechanicalParams* mparams, - sofa::type::vector& f, + sofa::type::vector& f, const sofa::VecCoord_t& x) { SCOPED_TIMER("ElementForces"); @@ -159,7 +159,7 @@ void FEMForceField::addDForce( template void FEMForceField::computeElementsForcesDeriv( - const sofa::core::MechanicalParams* mparams, sofa::type::vector& df, + const sofa::core::MechanicalParams* mparams, sofa::type::vector& df, const sofa::VecDeriv_t& dx) { SCOPED_TIMER("ElementForcesDeriv"); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h index 53d4b6e3fb5..97fa1d979e4 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -32,118 +33,8 @@ namespace sofa::component::solidmechanics::fem::elastic { -/** - * Specifies the type of matrix-vector product to be used with the stiffness matrix. - */ -enum class MatrixVectorProductType -{ - /** - * The matrix-vector product is computed using the factorization of the matrix - */ - Factorization, - - /** - * The matrix-vector product is computed using the dense matrix representation - */ - Dense -}; - -/** - * Represents an element stiffness matrix. It contains the dense matrix representation, but also - * its factorization. Both the factorization or the dense matrix can be used for the product of the - * stiffness matrix with a vector. Although it gives the same result, it has an impact on the - * performances. - */ -template -struct FactorizedElementStiffness -{ -private: - using Real = sofa::Real_t; - using FiniteElement = sofa::fem::FiniteElement; - static constexpr auto NbQuadraturePoints = FiniteElement::quadraturePoints().size(); - static constexpr sofa::Size NumberOfIndependentElements = sofa::type::NumberOfIndependentElements; - - /// The factors of the stiffness matrix - /// @{ - std::array, NbQuadraturePoints> B; - OrthotropicElasticityTensor> elasticityTensor; - std::array factors; - /// @} - - /** - * The assembled stiffness matrix - */ - sofa::type::Mat< - ElementType::NumberOfNodes * DataTypes::spatial_dimensions, - ElementType::NumberOfNodes * DataTypes::spatial_dimensions, - sofa::Real_t> stiffnessMatrix; - -public: - const StrainDisplacement& getB(std::size_t i) const { return B[i]; } - Real getFactor(std::size_t i) const { return factors[i]; } - const OrthotropicElasticityTensor>& getElasticityTensor() const { return elasticityTensor; } - - void setElasticityTensor( - const sofa::type::FullySymmetric4Tensor>& elasticityTensor_) - { - this->elasticityTensor.set(elasticityTensor_); - } - - void addFactor(std::size_t quadraturePointIndex, - const Real factor, - const StrainDisplacement& B_) - { - this->factors[quadraturePointIndex] = factor; - this->B[quadraturePointIndex] = B_; - - const auto K_i = factor * B_.multTranspose(elasticityTensor.toMat() * B_.B); - stiffnessMatrix += K_i; - } - - const auto& getAssembledMatrix() const { return stiffnessMatrix; } - - using Vec = sofa::type::Vec< - ElementType::NumberOfNodes * DataTypes::spatial_dimensions, - sofa::Real_t>; - - inline Vec operator*(const Vec& v) const - { - if constexpr (matrixVectorProductType == MatrixVectorProductType::Factorization) - { - if constexpr (NbQuadraturePoints > 1) - { - Vec result; - for (std::size_t i = 0; i < NbQuadraturePoints; ++i) - { - const auto& B = this->B[i]; - result += B.multTranspose(this->factors[i] * (this->elasticityTensor * (B * v))); - } - return result; - } - else - { - return this->B[0].multTranspose(this->factors[0] * (this->elasticityTensor * (this->B[0] * v))); - } - } - else - { - return this->stiffnessMatrix * v; - } - } - - friend std::ostream& operator<<(std::ostream& os, const FactorizedElementStiffness& obj) - { - return os << obj.getAssembledMatrix(); - } - - friend std::istream& operator>>(std::istream& in, FactorizedElementStiffness& obj) - { - return in; - } -}; - -template -FactorizedElementStiffness integrate( +template +auto integrate( const std::array, ElementType::NumberOfNodes>& nodesCoordinates, const sofa::type::FullySymmetric4Tensor>& elasticityTensor) { @@ -154,8 +45,9 @@ FactorizedElementStiffness inte static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; - FactorizedElementStiffness K; - K.setElasticityTensor(elasticityTensor); + typename trait::ElementHessian K; + + auto C = elasticityTensor.toVoigtMatSym().toMat(); std::size_t quadraturePointIndex = 0; for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) @@ -180,7 +72,7 @@ FactorizedElementStiffness inte dN_dq[i] = J_inv.transposed() * dN_dq_ref[i]; const auto B = makeStrainDisplacement(dN_dq); - K.addFactor(quadraturePointIndex, weight * detJ, B); + K += (weight * detJ) * B.multTranspose(C * B); ++quadraturePointIndex; } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/OrthotropicElasticityTensor.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/OrthotropicElasticityTensor.h index 771af1650ba..792305d8042 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/OrthotropicElasticityTensor.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/OrthotropicElasticityTensor.h @@ -50,92 +50,4 @@ auto makeIsotropicElasticityTensor(LameMu mu, LameLambda lambda) }}; } - -template -constexpr sofa::type::Vec orthotropicElasticityTensorProduct(const sofa::type::Mat& tensor, const sofa::type::Vec& v) -{ - return tensor * v; -} - -/** - * Specialization for 3D where the operations containing known zeros in the elasticity tensor are - * omitted. - * - * WARNING: this specialization is only valid for a give Voigt index mapping. If this mapping happens - * to change, this specialization will need to be updated accordingly. - */ -template -constexpr sofa::type::Vec<6, real> orthotropicElasticityTensorProduct(const sofa::type::Mat<6, 6, real>& tensor, const sofa::type::Vec<6, real>& v) -{ - sofa::type::Vec<6, real> result { sofa::type::NOINIT }; - - result[0] = tensor(0, 0) * v[0] + tensor(0, 1) * v[1] + tensor(0, 2) * v[2]; - result[1] = tensor(1, 0) * v[0] + tensor(1, 1) * v[1] + tensor(1, 2) * v[2]; - result[2] = tensor(2, 0) * v[0] + tensor(2, 1) * v[1] + tensor(2, 2) * v[2]; - result[3] = tensor(3, 3) * v[3]; - result[4] = tensor(4, 4) * v[4]; - result[5] = tensor(5, 5) * v[5]; - - return result; -} - -/** - * @class OrthotropicElasticityTensor - * @brief Represents an elasticity tensor for orthotropic materials. - * - * This class is a wrapper on a matrix (6x6 in 3D) with a specialization for the matrix-vector - * product, speeding up computations compared to a regular product. It uses the known shape of an - * orthotropic material to optimize the matrix-vector product omitting the operations containing - * zero coefficients. - */ -template -struct OrthotropicElasticityTensor -{ - using Real = real; - static constexpr auto NbIndependentElements = sofa::type::NumberOfIndependentElements; - - explicit OrthotropicElasticityTensor(const sofa::type::Mat& mat) : C(mat) {} - OrthotropicElasticityTensor() = default; - - void set(const sofa::type::Mat& mat) - { - C = mat; - } - - void set(const sofa::type::FullySymmetric4Tensor& tensor) - { - C = tensor.toVoigtMatSym().toMat(); - } - - sofa::type::Vec operator*(const sofa::type::Vec& v) const - { - return orthotropicElasticityTensorProduct(C, v); - } - - template - sofa::type::Mat operator*(const sofa::type::Mat& v) const noexcept - { - return C * v; - } - - Real operator()(sofa::Size i, sofa::Size j) const - { - return C(i, j); - } - - Real& operator()(sofa::Size i, sofa::Size j) - { - return C(i, j); - } - - const sofa::type::Mat& toMat() const - { - return C; - } - -private: - sofa::type::Mat C; - -}; - } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h index 78c28124837..ccc0b9f98a6 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h @@ -21,10 +21,11 @@ ******************************************************************************/ #pragma once +#include #include -#include #include -#include + +#include namespace sofa::component::solidmechanics::fem::elastic { @@ -59,17 +60,13 @@ struct trait /// the concatenation of the displacement of the element nodes in a single vector using ElementDisplacement = sofa::type::Vec; - /// tells how to compute the matrix-vector product of the stiffness matrix with a displacement - /// vector. It does not change the result, but it can have an impact on performances. - static constexpr MatrixVectorProductType matrixVectorProductType = - NbQuadraturePoints > 1 - ? MatrixVectorProductType::Dense - : MatrixVectorProductType::Factorization; - - /// the type of the element stiffness matrix - using ElementStiffness = sofa::component::solidmechanics::fem::elastic::FactorizedElementStiffness; + /// the type of the element hessian matrix + using ElementHessian = sofa::type::Mat< + NumberOfDofsInElement, + NumberOfDofsInElement, + sofa::Real_t>; - using ElementForce = sofa::type::Vec>; + using ElementGradient = sofa::type::Vec>; }; }