[FEM] Compute and draw von Mises stress

Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt

@@ -13,6 +13,7 @@ set(HEADER_FILES
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.inl
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CauchyStressEvaluator.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.h
@@ -45,6 +46,8 @@ set(HEADER_FILES
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceFieldOptim.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceFieldOptim.inl
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VonMisesStress.h
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VonMisesStress.inl
 
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/ComputeStrategy.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/ElementStiffnessMatrix.h
@@ -84,6 +87,7 @@ set(SOURCE_FILES
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularAnisotropicFEMForceField.cpp
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceField.cpp
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceFieldOptim.cpp
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VonMisesStress.cpp
 )
 
 sofa_find_package(Sofa.Simulation.Core REQUIRED)

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

@@ -0,0 +1,45 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* This program is free software; you can redistribute it and/or modify it     *
+* under the terms of the GNU Lesser General Public License as published by    *
+* the Free Software Foundation; either version 2.1 of the License, or (at     *
+* your option) any later version.                                             *
+*                                                                             *
+* This program is distributed in the hope that it will be useful, but WITHOUT *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
+* for more details.                                                           *
+*                                                                             *
+* You should have received a copy of the GNU Lesser General Public License    *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+
+#include <sofa/component/solidmechanics/fem/elastic/config.h>
+#include <sofa/core/objectmodel/BaseComponent.h>
+
+#include <sofa/type/MatSym.h>
+
+namespace sofa::component::solidmechanics::fem::elastic
+{
+
+template<class DataTypes>
+class CauchyStressEvaluator : public virtual core::objectmodel::BaseComponent
+{
+public:
+    SOFA_ABSTRACT_CLASS(CauchyStressEvaluator, BaseComponent)
+
+    static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions;
+    using StressVoigtVector = sofa::type::Vec<sofa::type::NumberOfIndependentElements<spatial_dimensions>, sofa::Real_t<DataTypes>>;
+    using DeformationGradient = sofa::type::Mat<spatial_dimensions, spatial_dimensions, Real_t<DataTypes>>;
+
+    virtual StressVoigtVector computeStress(const DeformationGradient& F, sofa::Size elementId) = 0;
+};
+
+}

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

@@ -24,7 +24,7 @@
 #include <sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h>
 #include <sofa/component/solidmechanics/fem/elastic/config.h>
 #include <sofa/component/solidmechanics/fem/elastic/impl/trait.h>
-
+#include <sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h>
 #include <sofa/component/solidmechanics/fem/elastic/FEMForceField.h>
 
 #if !defined(ELASTICITY_COMPONENT_ELEMENT_LINEAR_SMALL_STRAIN_FEM_FORCE_FIELD_CPP)
@@ -37,20 +37,24 @@ namespace sofa::component::solidmechanics::fem::elastic
 template <class DataTypes, class ElementType>
 class LinearSmallStrainFEMForceField :
     public BaseElementLinearFEMForceField<DataTypes, ElementType>,
-    public FEMForceField<DataTypes, ElementType>
+    public FEMForceField<DataTypes, ElementType>,
+    public CauchyStressEvaluator<DataTypes>
 {
 public:
-    SOFA_CLASS2(
+    SOFA_CLASS3(
         SOFA_TEMPLATE2(LinearSmallStrainFEMForceField, DataTypes, ElementType),
             SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType),
-            SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType));
+            SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType),
+            CauchyStressEvaluator<DataTypes>);
 
 private:
     using trait = typename FEMForceField<DataTypes, ElementType>::trait;
     using ElementHessian = typename trait::ElementHessian;
     using ElementDisplacement = typename trait::ElementDisplacement;
     using StrainDisplacement = typename trait::StrainDisplacement;
     using ElementGradient = typename trait::ElementGradient;
+    using DeformationGradient = typename trait::DeformationGradient;
+    using StressVoigtVector = typename trait::StressVoigtVector;
 
 public:
     void init() override;
@@ -64,6 +68,8 @@ class LinearSmallStrainFEMForceField :
     // almost deprecated, but here for compatibility with unit tests
     void addKToMatrix(sofa::linearalgebra::BaseMatrix* matrix, SReal kFact, unsigned& offset) override;
 
+    StressVoigtVector computeStress(const DeformationGradient& F, sofa::Size elementId) override;
+
 protected:
 
     void computeElementsForces(

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

@@ -180,4 +180,32 @@ void LinearSmallStrainFEMForceField<DataTypes, ElementType>::addKToMatrix(
     }
 }
 
+template <class DataTypes, class ElementType>
+auto LinearSmallStrainFEMForceField<DataTypes, ElementType>::computeStress(
+    const DeformationGradient& F, sofa::Size elementId) -> StressVoigtVector
+{
+    const auto strainTensor = static_cast<Real_t<DataTypes>>(1)/2 * (F + F.transposed())  - DeformationGradient::Identity();
+
+    sofa::type::Vec<type::NumberOfIndependentElements<trait::spatial_dimensions>, Real_t<DataTypes> > strainVoigt;
+    for (sofa::Size i = 0; i < type::NumberOfIndependentElements<trait::spatial_dimensions>; ++i)
+    {
+        const auto [p, q] = type::toTensorIndices<trait::spatial_dimensions>(i);
+        strainVoigt[i] = strainTensor(p, q);
+    }
+
+    const auto youngModulus = this->getYoungModulusInElement(elementId);
+    const auto poissonRatio = this->getPoissonRatioInElement(elementId);
+
+    LameLambda<Real_t<DataTypes>> lambda { 0 };
+    LameMu<Real_t<DataTypes>> mu { 0 };
+
+    sofa::component::solidmechanics::fem::elastic::toLameParameters<DataTypes::spatial_dimensions, Real_t<DataTypes>>(
+        YoungModulus<Real_t<DataTypes>>(youngModulus), PoissonRatio<Real_t<DataTypes>>(poissonRatio),
+        lambda, mu);
+
+    const auto elasticityTensor = makeIsotropicElasticityTensor<DataTypes::spatial_dimensions, Real_t<DataTypes>>(mu, lambda);
+
+    return elasticityTensor.toVoigtMatSym().toMat() * strainVoigt;
+}
+
 }  // namespace sofa::component::solidmechanics::fem::elastic

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp

@@ -0,0 +1,58 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* This program is free software; you can redistribute it and/or modify it     *
+* under the terms of the GNU Lesser General Public License as published by    *
+* the Free Software Foundation; either version 2.1 of the License, or (at     *
+* your option) any later version.                                             *
+*                                                                             *
+* This program is distributed in the hope that it will be useful, but WITHOUT *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
+* for more details.                                                           *
+*                                                                             *
+* You should have received a copy of the GNU Lesser General Public License    *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#include <sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl>
+#include <sofa/core/ObjectFactory.h>
+#include <sofa/fem/FiniteElement[all].h>
+
+namespace sofa::component::solidmechanics::fem::elastic
+{
+
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec2Types, sofa::geometry::Edge>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Edge>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec2Types, sofa::geometry::Triangle>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Triangle>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec2Types, sofa::geometry::Quad>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Pyramid>;
+
+void registerVonMisesStress(sofa::core::ObjectFactory* factory)
+{
+    factory->registerObjects(sofa::core::ObjectRegistrationData("Compute and draw von Mises stress based on a local least-square projection in each element")
+        .add< VonMisesStress<sofa::defaulttype::Vec1Types, sofa::geometry::Edge> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec2Types, sofa::geometry::Edge> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Edge> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec2Types, sofa::geometry::Triangle> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Triangle> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec2Types, sofa::geometry::Quad> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Quad> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Prism> >()
+        .add< VonMisesStress<sofa::defaulttype::Vec3Types, sofa::geometry::Pyramid> >()
+    );
+}
+
+}

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

@@ -0,0 +1,137 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* This program is free software; you can redistribute it and/or modify it     *
+* under the terms of the GNU Lesser General Public License as published by    *
+* the Free Software Foundation; either version 2.1 of the License, or (at     *
+* your option) any later version.                                             *
+*                                                                             *
+* This program is distributed in the hope that it will be useful, but WITHOUT *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
+* for more details.                                                           *
+*                                                                             *
+* You should have received a copy of the GNU Lesser General Public License    *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+
+#include <sofa/component/solidmechanics/fem/elastic/config.h>
+#include <sofa/component/solidmechanics/fem/elastic/impl/trait.h>
+#include <sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h>
+#include <sofa/core/behavior/SingleStateAccessor.h>
+#include <sofa/core/behavior/TopologyAccessor.h>
+#include <sofa/core/visual/DrawColoredMesh.h>
+#include <sofa/helper/ColorMap.h>
+
+namespace sofa::component::solidmechanics::fem::elastic
+{
+
+/**
+ * @brief Component that computes and visualizes the von Mises stress field on a FEM mesh.
+ *
+ * This component evaluates the Cauchy stress at each quadrature point of the elements
+ * using a linked @see CauchyStressEvaluator. Since stress values are naturally calculated
+ * at quadrature points (e.g., Gauss points), but visualization typically requires nodal
+ * values, this component performs a least-squares projection.
+ *
+ * The projection maps the stress tensor components from the quadrature points to the
+ * element nodes by solving a local system involving the Gram matrix (the mass matrix
+ * of the shape functions).
+ *
+ * Once the stress tensor is projected to the nodes, the deviatoric part and the
+ * resulting von Mises scalar are computed for visualization via a color map.
+ *
+ * @tparam DataTypes The SOFA DataTypes (e.g. Vec3d)
+ * @tparam ElementType The geometric element type (e.g. Tetrahedron, Hexahedron)
+ */
+template <class DataTypes, class ElementType>
+class VonMisesStress : public core::behavior::TopologyAccessor, public core::behavior::SingleStateAccessor<DataTypes>
+{
+public:
+    SOFA_CLASS2(
+        SOFA_TEMPLATE2(VonMisesStress, DataTypes, ElementType),
+        core::behavior::TopologyAccessor,
+        core::behavior::SingleStateAccessor<DataTypes>);
+
+protected:
+    using FiniteElement = sofa::fem::FiniteElement<ElementType, DataTypes>;
+
+    static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions;
+    static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes;
+    static constexpr sofa::Size NumberOfDofsInElement = NumberOfNodesInElement * spatial_dimensions;
+    static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension;
+    static constexpr sofa::Size NumberOfQuadraturePoints = FiniteElement::quadraturePoints().size();
+
+    // a stress tensor represented as a vector using the Voigt mapping
+    using StressVoigtVector = sofa::type::Vec<sofa::type::NumberOfIndependentElements<spatial_dimensions>, sofa::Real_t<DataTypes>>;
+
+    using DeformationGradient = sofa::type::Mat<spatial_dimensions, spatial_dimensions, Real_t<DataTypes>>;
+
+public:
+
+    void init() override;
+    void draw(const core::visual::VisualParams* vparams) override;
+    void computeBBox(const core::ExecParams* params, bool onlyVisible) override;
+
+    // A stress value for each node in an element
+    using LocalStressValues = std::array<sofa::Real_t<DataTypes>, NumberOfNodesInElement>;
+    Data<sofa::type::vector<LocalStressValues>> d_nodalStress;
+
+    Data<helper::ColorMap> d_colorMap;
+    Data<bool> d_lighting;
+
+    sofa::SingleLink<MyType, CauchyStressEvaluator<DataTypes>,
+        sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_stressEvaluator;
+
+    VonMisesStress();
+
+    void handleEvent(core::objectmodel::Event*) override;
+
+protected:
+
+    void validateStressEvaluatorLink();
+
+    using ElementGramMatrix = sofa::type::Mat<NumberOfNodesInElement, NumberOfNodesInElement, sofa::Real_t<DataTypes>>;
+    sofa::type::vector<ElementGramMatrix> m_elementInverseGramMatrices;
+
+    struct PrecomputedData
+    {
+        sofa::type::Mat<spatial_dimensions, TopologicalDimension, Real_t<DataTypes>> jacobian { sofa::type::NOINIT };
+        sofa::type::Mat<TopologicalDimension, spatial_dimensions, Real_t<DataTypes>> jacobianInv { sofa::type::NOINIT };
+        Real_t<DataTypes> detJacobian {};
+        sofa::type::Mat<NumberOfNodesInElement, spatial_dimensions, Real_t<DataTypes>> dN_dQ { sofa::type::NOINIT };
+    };
+
+    sofa::type::vector<std::array<PrecomputedData, NumberOfQuadraturePoints>> m_precomputedData;
+
+    void precomputeData();
+
+    /**
+     * @brief Computes the inverse of the Gram matrix (integrated N^T * N) for each element.
+     * This matrix is used for the least-square projection of values from quadrature points to nodes.
+     */
+    void calculateElementInverseGramMatrices(const auto& elements, sofa::type::vector<ElementGramMatrix>& inverseGramMatrices);
+
+    /**
+     * @brief Projects values evaluated at quadrature points to nodal values using least-squares.
+     * @param elementId Index of the element
+     * @param valuesAtQuadraturePoints Array of values evaluated at each quadrature point of the element
+     * @return Array of projected values at each node of the element
+     */
+    std::array<StressVoigtVector, NumberOfNodesInElement> projectQuadraturePointValuesToNodes(
+        sofa::Size elementId,
+        const std::array<StressVoigtVector, NumberOfQuadraturePoints>& valuesAtQuadraturePoints) const;
+
+    static StressVoigtVector deviatoricStress(const StressVoigtVector& sigma);
+    static Real_t<DataTypes> vonMisesStress(const StressVoigtVector& deviatoricStress);
+
+    core::visual::DrawElementColoredMesh<ElementType> m_renderer;
+};
+
+}  // namespace sofa::component::solidmechanics::fem::elastic