From a08f8e939d8f24f0e16d79168ba6fe44c711bc94 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 27 Jul 2026 15:30:19 +0200 Subject: [PATCH 01/28] initialize a skeleton for a component computing von Mises stress --- .../SolidMechanics/FEM/Elastic/CMakeLists.txt | 3 + .../fem/elastic/VonMisesStress.cpp | 57 +++++++++++++++++++ .../fem/elastic/VonMisesStress.h | 44 ++++++++++++++ .../fem/elastic/VonMisesStress.inl | 44 ++++++++++++++ .../solidmechanics/fem/elastic/init.cpp | 2 + 5 files changed, 150 insertions(+) create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index fddd6033a0d..77329f6a861 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt @@ -45,6 +45,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 +86,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) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp new file mode 100644 index 00000000000..6ce240cabf1 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp @@ -0,0 +1,57 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; + +void registerVonMisesStress(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Compute von Mises stress") + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + ); +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h new file mode 100644 index 00000000000..32dc0bd07fd --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -0,0 +1,44 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +class VonMisesStress : public core::objectmodel::BaseComponent +{ +public: + SOFA_CLASS( + SOFA_TEMPLATE2(VonMisesStress, DataTypes, ElementType), + core::objectmodel::BaseComponent); + +protected: + VonMisesStress(); + + void handleEvent(core::objectmodel::Event*) override; +}; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl new file mode 100644 index 00000000000..145d007fb6b --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -0,0 +1,44 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +VonMisesStress::VonMisesStress() +{ + // This component must receive events + f_listening.setValue(true); +} +template +void VonMisesStress::handleEvent(core::objectmodel::Event* event) +{ + if (simulation::AnimateEndEvent::checkEventType(event)) + { + msg_info() << "Computing von Mises stress"; + } +} + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp index 2c37e45f155..87d14c740ff 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp @@ -41,6 +41,7 @@ extern void registerTriangleFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangularAnisotropicFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangularFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangularFEMForceFieldOptim(sofa::core::ObjectFactory* factory); +extern void registerVonMisesStress(sofa::core::ObjectFactory* factory); extern "C" { SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule(); @@ -81,6 +82,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) registerTriangularAnisotropicFEMForceField(factory); registerTriangularFEMForceField(factory); registerTriangularFEMForceFieldOptim(factory); + registerVonMisesStress(factory); } void init() From 1852a6329e4b7f217e5995ffd162e43dd4c595ad Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 27 Jul 2026 15:30:40 +0200 Subject: [PATCH 02/28] add an example scene that will show von Mises stress --- examples/Demos/elongation.scn | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 examples/Demos/elongation.scn diff --git a/examples/Demos/elongation.scn b/examples/Demos/elongation.scn new file mode 100644 index 00000000000..d9dc53ebc8d --- /dev/null +++ b/examples/Demos/elongation.scn @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + From a9ca8894533b35ec3edfb71ec89b3788620f5eac Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 27 Jul 2026 16:09:15 +0200 Subject: [PATCH 03/28] draw hexahedra in DrawColoredMesh --- .../src/sofa/core/visual/DrawColoredMesh.h | 57 ++++++++++++++++++- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/Sofa/framework/Core/src/sofa/core/visual/DrawColoredMesh.h b/Sofa/framework/Core/src/sofa/core/visual/DrawColoredMesh.h index 692a64eb2d3..f0ac63f6a89 100644 --- a/Sofa/framework/Core/src/sofa/core/visual/DrawColoredMesh.h +++ b/Sofa/framework/Core/src/sofa/core/visual/DrawColoredMesh.h @@ -149,6 +149,58 @@ struct SOFA_CORE_API DrawElementColoredMesh }; +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Hexahedron; + friend BaseDrawColoredMesh; + static constexpr std::size_t NumberQuadsInHexahedron = 6; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + const auto& elements = topology->getHexahedra(); + const auto& facets = topology->getQuads(); + + if(facets.empty()) + { + msg_error_once("DrawElementColoredMesh") << "Drawing hexahedra needs the associated quads in the topology."; + return; + } + + renderedPoints.clear(); + colors.clear(); + + for (auto i : elementIndices) + { + const auto& element = elements[i]; + const auto& facetsInElement = topology->getQuadsInHexahedron(i); + const auto center = this->elementCenter(position, element); + + for (std::size_t j = 0; j < NumberQuadsInHexahedron; ++j) + { + const auto faceId = facetsInElement[j]; + for (const auto vertexId : facets[faceId]) + { + const auto p = this->applyElementSpace(position[vertexId], center); + + renderedPoints.push_back(p); + colors.push_back(nodesColors[vertexId]); + } + } + } + + drawTool->drawQuads(renderedPoints, colors); + } + +}; + class SOFA_CORE_API DrawColoredMesh { public: @@ -159,6 +211,7 @@ class SOFA_CORE_API DrawColoredMesh void draw(sofa::helper::visual::DrawTool* drawTool, const PositionContainer& position, const ColorContainer& nodesColors, sofa::core::topology::BaseMeshTopology* topology) { std::get>(m_meshes).drawAllElements(drawTool, position, topology, nodesColors); + std::get>(m_meshes).drawAllElements(drawTool, position, topology, nodesColors); } private: @@ -166,8 +219,8 @@ class SOFA_CORE_API DrawColoredMesh // DrawElementColoredMesh, // DrawElementColoredMesh, // DrawElementColoredMesh, - DrawElementColoredMesh - // DrawElementColoredMesh, + DrawElementColoredMesh, + DrawElementColoredMesh // DrawElementColoredMesh, // DrawElementColoredMesh > m_meshes; From a86f607cb452f538a2af2037e60cc7982488b70d Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 27 Jul 2026 16:09:28 +0200 Subject: [PATCH 04/28] don't average in DrawToolGL --- Sofa/GL/src/sofa/gl/DrawToolGL.cpp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/Sofa/GL/src/sofa/gl/DrawToolGL.cpp b/Sofa/GL/src/sofa/gl/DrawToolGL.cpp index fbbe4087f13..b339d3b1b95 100644 --- a/Sofa/GL/src/sofa/gl/DrawToolGL.cpp +++ b/Sofa/GL/src/sofa/gl/DrawToolGL.cpp @@ -993,15 +993,9 @@ void DrawToolGL::drawQuads(const std::vector &points, const std::vector Date: Mon, 27 Jul 2026 16:09:51 +0200 Subject: [PATCH 05/28] create data and accessors in VonMisesStress --- .../fem/elastic/VonMisesStress.cpp | 1 + .../fem/elastic/VonMisesStress.h | 19 +++++++++--- .../fem/elastic/VonMisesStress.inl | 31 ++++++++++++++++++- examples/Demos/elongation.scn | 7 +++-- 4 files changed, 50 insertions(+), 8 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp index 6ce240cabf1..b252f6c4e0c 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp @@ -21,6 +21,7 @@ ******************************************************************************/ #include #include +#include namespace sofa::component::solidmechanics::fem::elastic { diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index 32dc0bd07fd..17403b1b32a 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -22,23 +22,34 @@ #pragma once #include -#include +#include +#include +#include namespace sofa::component::solidmechanics::fem::elastic { template -class VonMisesStress : public core::objectmodel::BaseComponent +class VonMisesStress : public core::behavior::TopologyAccessor, public core::behavior::SingleStateAccessor { public: - SOFA_CLASS( + SOFA_CLASS2( SOFA_TEMPLATE2(VonMisesStress, DataTypes, ElementType), - core::objectmodel::BaseComponent); + core::behavior::TopologyAccessor, + core::behavior::SingleStateAccessor); + + void init() override; + + Data > > d_nodalStress; protected: + + using trait = sofa::component::solidmechanics::fem::elastic::trait; + VonMisesStress(); void handleEvent(core::objectmodel::Event*) override; + }; } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 145d007fb6b..e986279792d 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -28,16 +28,45 @@ namespace sofa::component::solidmechanics::fem::elastic template VonMisesStress::VonMisesStress() + : d_nodalStress(initData(&d_nodalStress, sofa::type::vector >{}, "nodalStress", "Nodal von Mises stress values")) { // This component must receive events f_listening.setValue(true); } + +template +void VonMisesStress::init() +{ + core::behavior::SingleStateAccessor::init(); + + if (!this->isComponentStateInvalid()) + { + this->validateTopology(); + } + + if (!this->isComponentStateInvalid()) + { + auto nodalStress = sofa::helper::getWriteOnlyAccessor(d_nodalStress); + nodalStress->resize(this->mstate->getSize()); + + // std::iota(nodalStress->begin(), nodalStress->end(), 0); + } +} + template void VonMisesStress::handleEvent(core::objectmodel::Event* event) { if (simulation::AnimateEndEvent::checkEventType(event)) { - msg_info() << "Computing von Mises stress"; + auto nodalStress = sofa::helper::getWriteOnlyAccessor(d_nodalStress); + nodalStress->resize(this->mstate->getSize()); + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + + for (const auto& element : elements) + { + + } } } diff --git a/examples/Demos/elongation.scn b/examples/Demos/elongation.scn index d9dc53ebc8d..911aad743ea 100644 --- a/examples/Demos/elongation.scn +++ b/examples/Demos/elongation.scn @@ -3,7 +3,7 @@ - + @@ -24,8 +24,9 @@ - + - + From ee9dbca828f557f9ae2bf21901d9299180a48489 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Tue, 28 Jul 2026 10:08:01 +0200 Subject: [PATCH 06/28] start loop on quadrature points --- .../component/solidmechanics/fem/elastic/VonMisesStress.inl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index e986279792d..89b87a11248 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -65,7 +65,9 @@ void VonMisesStress::handleEvent(core::objectmodel::Even for (const auto& element : elements) { - + for (const auto& [quadraturePoint, weight] : trait::FiniteElement::quadraturePoints()) + { + } } } } From 1d0476aa1ddb67910b54d28f4531fc80340c6383 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 30 Jul 2026 16:04:42 +0200 Subject: [PATCH 07/28] global architecture --- .../fem/elastic/VonMisesStress.h | 29 +++- .../fem/elastic/VonMisesStress.inl | 146 +++++++++++++++++- examples/Demos/elongation.scn | 15 +- 3 files changed, 172 insertions(+), 18 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index 17403b1b32a..d6d6bc416ef 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -38,18 +38,39 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh core::behavior::TopologyAccessor, core::behavior::SingleStateAccessor); - void init() override; +protected: + using FiniteElement = sofa::fem::FiniteElement; - Data > > d_nodalStress; -protected: + 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; + + // a stress tensor represented as a vector using the Voigt mapping + using StressVoigtVector = std::array, sofa::type::NumberOfIndependentElements>; + +public: - using trait = sofa::component::solidmechanics::fem::elastic::trait; + void init() override; + + // A stress value for each node in an element + using LocalStressValues = std::array, NumberOfNodesInElement>; + Data> d_nodalStress; + + Data d_continuousField; VonMisesStress(); void handleEvent(core::objectmodel::Event*) override; +protected: + using ElementMassMatrix = sofa::type::Mat>; + sofa::type::vector m_elementMassMatrices; + + void calculateElementMassMatrix(const auto& elements, sofa::type::vector &elementMassMatrices); + + static StressVoigtVector deviatoricStress(const StressVoigtVector& sigma); }; } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 89b87a11248..9b681fe4982 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -21,6 +21,7 @@ ******************************************************************************/ #pragma once #include +#include #include namespace sofa::component::solidmechanics::fem::elastic @@ -28,7 +29,12 @@ namespace sofa::component::solidmechanics::fem::elastic template VonMisesStress::VonMisesStress() - : d_nodalStress(initData(&d_nodalStress, sofa::type::vector >{}, "nodalStress", "Nodal von Mises stress values")) + : d_nodalStress(initData(&d_nodalStress, sofa::type::vector{}, "nodalStress", + "Local nodal von Mises stress values")) + , d_continuousField(initData(&d_continuousField, false, "continuousField", + "Compute von Mises stress as a continuous field across the elements. Necessitate the solve " + "of a sparse linear system. Otherwise, the von Mises stress is computed locally. A local " + "stress value may indicate discretization errors if the field does not appear continuous.")) { // This component must receive events f_listening.setValue(true); @@ -47,9 +53,10 @@ void VonMisesStress::init() if (!this->isComponentStateInvalid()) { auto nodalStress = sofa::helper::getWriteOnlyAccessor(d_nodalStress); - nodalStress->resize(this->mstate->getSize()); + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + nodalStress->resize(elements.size()); - // std::iota(nodalStress->begin(), nodalStress->end(), 0); + this->calculateElementMassMatrix(elements, m_elementMassMatrices); } } @@ -59,17 +66,142 @@ void VonMisesStress::handleEvent(core::objectmodel::Even if (simulation::AnimateEndEvent::checkEventType(event)) { auto nodalStress = sofa::helper::getWriteOnlyAccessor(d_nodalStress); + nodalStress->clear(); nodalStress->resize(this->mstate->getSize()); - const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + auto restPositionAccessor = this->mstate->readRestPositions(); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + const auto nbElements = elements.size(); + helper::IotaView indices{static_cast(0ul), nbElements}; + + std::for_each( + indices.begin(), indices.end(), + [&](const auto elementId) + { + const auto& element = elements[elementId]; + + std::array, NumberOfNodesInElement> nodeCoordinatesInElement; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + nodeCoordinatesInElement[i] = restPositionAccessor[element[i]]; + } + + std::array nodalStressInElement; + + for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) + { + sofa::type::Vec> b; + + for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) + { + // gradient of shape functions in the reference element evaluated at the quadrature + // point + const sofa::type::Mat> + dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + + // jacobian of the mapping from the reference space to the physical space, evaluated + // at the quadrature point + sofa::type::Mat> + jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + nodeCoordinatesInElement, dN_dq_ref); + + const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); + + // shape functions in the reference element evaluated at the quadrature point + const auto N = FiniteElement::shapeFunctions(quadraturePoint); + + StressVoigtVector stress; + + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + b[j] += N[j] * weight * detJ * stress[i]; + } + } + + const auto stressCoordinate = m_elementMassMatrices[elementId] * b; + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + nodalStressInElement[j][i] = stressCoordinate[j]; + // nodalStress[elementId][j] = stressCoordinate[j]; + } + } + + // nodalStressInElement + }); + } +} + +template +void VonMisesStress::calculateElementMassMatrix( + const auto& elements, sofa::type::vector& elementMassMatrices) +{ + const auto nbElements = elements.size(); + elementMassMatrices.resize(nbElements); + + auto restPositionAccessor = this->mstate->readRestPositions(); - for (const auto& element : elements) + SCOPED_TIMER("elementMassMatrix"); + helper::IotaView indices{static_cast(0ul), nbElements}; + std::for_each( + indices.begin(), indices.end(), + [&](const auto elementId) { - for (const auto& [quadraturePoint, weight] : trait::FiniteElement::quadraturePoints()) + const auto& element = elements[elementId]; + auto& elementMassMatrix = elementMassMatrices[elementId]; + + std::array, NumberOfNodesInElement> nodeCoordinatesInElement; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) { + nodeCoordinatesInElement[i] = restPositionAccessor[element[i]]; } - } + + for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) + { + // gradient of shape functions in the reference element evaluated at the quadrature + // point + const sofa::type::Mat> + dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + + // jacobian of the mapping from the reference space to the physical space, evaluated + // at the quadrature point + sofa::type::Mat> + jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + nodeCoordinatesInElement, dN_dq_ref); + + const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); + + // shape functions in the reference element evaluated at the quadrature point + const auto N = FiniteElement::shapeFunctions(quadraturePoint); + + const auto NT_N = sofa::type::dyad(N, N); + + elementMassMatrix += (weight * detJ) * NT_N; + } + + sofa::type::invertMatrix(elementMassMatrix, elementMassMatrix); + }); +} + +template +auto VonMisesStress::deviatoricStress(const StressVoigtVector& sigma) -> StressVoigtVector +{ + StressVoigtVector s = sigma; + + Real_t trace {}; + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + trace += sigma[type::tensorToVoigtIndex(i, i)]; } + + static constexpr Real_t dim_inv = static_cast>(1.0 / spatial_dimensions); + + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + s[type::tensorToVoigtIndex(i, i)] -= dim_inv * trace; + } + + return s; } } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/examples/Demos/elongation.scn b/examples/Demos/elongation.scn index 911aad743ea..c0d5deb631b 100644 --- a/examples/Demos/elongation.scn +++ b/examples/Demos/elongation.scn @@ -1,5 +1,5 @@ - + @@ -13,20 +13,21 @@ + keyTimes="0 2 4 100000" movements="0 0 0 -0.1 0 0 0.1 0 0 0.1 0 0"/> + + keyTimes="0 2 4 100000" movements="0 0 0 0.1 0 0 -0.1 0 0 -0.1 0 0"/> - + - - + + + From 21e6de1889c6e240c1b907aac2986d72c3825826 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 30 Jul 2026 16:10:31 +0200 Subject: [PATCH 08/28] add function to compute von Mises stress --- .../solidmechanics/fem/elastic/VonMisesStress.h | 1 + .../fem/elastic/VonMisesStress.inl | 16 ++++++++++++++++ 2 files changed, 17 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index d6d6bc416ef..6a851cfdd5e 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -71,6 +71,7 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh void calculateElementMassMatrix(const auto& elements, sofa::type::vector &elementMassMatrices); static StressVoigtVector deviatoricStress(const StressVoigtVector& sigma); + static Real_t vonMisesStress(const StressVoigtVector& deviatoricStress); }; } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 9b681fe4982..aec78201c59 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -204,4 +204,20 @@ auto VonMisesStress::deviatoricStress(const StressVoigtV return s; } +template +Real_t +VonMisesStress::vonMisesStress(const StressVoigtVector& deviatoricStress) +{ + sofa::type::Mat> s; + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + for (sofa::Size j = 0; j < spatial_dimensions; ++j) + { + s[i][j] = deviatoricStress[type::tensorToVoigtIndex(i, j)]; + } + } + + return std::sqrt(static_cast>(3) / 2 * sofa::type::trace(s.multTranspose(s))); +} + } // namespace sofa::component::solidmechanics::fem::elastic From 00c7ac61fe2639bee9d256b4f2615506907a7b0e Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 30 Jul 2026 16:14:50 +0200 Subject: [PATCH 09/28] compute von Mises on nodes --- .../solidmechanics/fem/elastic/VonMisesStress.inl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index aec78201c59..7ebdbac03db 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -123,11 +123,14 @@ void VonMisesStress::handleEvent(core::objectmodel::Even for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) { nodalStressInElement[j][i] = stressCoordinate[j]; - // nodalStress[elementId][j] = stressCoordinate[j]; } } - // nodalStressInElement + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + nodalStress[elementId][i] = vonMisesStress(deviatoricStress(nodalStressInElement[i])); + } + }); } } From 576cc1f90b6f8db31f3d55e4f47cbb4249db3b59 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 30 Jul 2026 16:21:43 +0200 Subject: [PATCH 10/28] rename --- .../solidmechanics/fem/elastic/VonMisesStress.inl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 7ebdbac03db..1f5321d9ddc 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -102,11 +102,10 @@ void VonMisesStress::handleEvent(core::objectmodel::Even // jacobian of the mapping from the reference space to the physical space, evaluated // at the quadrature point - sofa::type::Mat> - jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( - nodeCoordinatesInElement, dN_dq_ref); + const sofa::type::Mat> J_q = + FiniteElement::Helper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); - const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); + const auto detJ = sofa::type::absGeneralizedDeterminant(J_q); // shape functions in the reference element evaluated at the quadrature point const auto N = FiniteElement::shapeFunctions(quadraturePoint); From a1c539ff555de60ab4841aacf0b1bd60cfb8b1dd Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 30 Jul 2026 16:27:43 +0200 Subject: [PATCH 11/28] compute deformation gradient --- .../fem/elastic/VonMisesStress.h | 15 +++++ .../fem/elastic/VonMisesStress.inl | 62 ++++++++++++++++--- 2 files changed, 67 insertions(+), 10 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index 6a851cfdd5e..6d4497d4b59 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -46,10 +46,13 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh 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 = std::array, sofa::type::NumberOfIndependentElements>; + using DeformationGradient = sofa::type::Mat>; + public: void init() override; @@ -68,6 +71,18 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh using ElementMassMatrix = sofa::type::Mat>; sofa::type::vector m_elementMassMatrices; + struct PrecomputedData + { + sofa::type::Mat> jacobian { sofa::type::NOINIT }; + sofa::type::Mat> jacobianInv { sofa::type::NOINIT }; + Real_t detJacobian {}; + sofa::type::Mat> dN_dQ { sofa::type::NOINIT }; + }; + + sofa::type::vector> m_precomputedData; + + void precomputeData(); + void calculateElementMassMatrix(const auto& elements, sofa::type::vector &elementMassMatrices); static StressVoigtVector deviatoricStress(const StressVoigtVector& sigma); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 1f5321d9ddc..2ef09b3b31d 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -56,6 +56,7 @@ void VonMisesStress::init() const auto& elements = FiniteElement::getElementSequence(*this->l_topology); nodalStress->resize(elements.size()); + this->precomputeData(); this->calculateElementMassMatrix(elements, m_elementMassMatrices); } } @@ -88,27 +89,32 @@ void VonMisesStress::handleEvent(core::objectmodel::Even } std::array nodalStressInElement; + static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); + static constexpr auto quadraturePoints = FiniteElement::quadraturePoints(); + for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) { sofa::type::Vec> b; - for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) + for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) { - // gradient of shape functions in the reference element evaluated at the quadrature - // point - const sofa::type::Mat> - dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + const auto& weight = quadraturePoints[q].second; + const auto& precomputed = m_precomputedData[elementId][q]; + + // gradient of shape functions in the reference element evaluated at the quadrature point + const auto& dN_dq_ref = gradients[q]; + + // jacobian of the mapping from the reference space to the CURRENT physical space + const auto J_q = FiniteElement::Helper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); - // jacobian of the mapping from the reference space to the physical space, evaluated - // at the quadrature point - const sofa::type::Mat> J_q = - FiniteElement::Helper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); + // Deformation Gradient F = J_curr * J_rest_inv + const DeformationGradient F = J_q * precomputed.jacobianInv; const auto detJ = sofa::type::absGeneralizedDeterminant(J_q); // shape functions in the reference element evaluated at the quadrature point - const auto N = FiniteElement::shapeFunctions(quadraturePoint); + const auto N = FiniteElement::shapeFunctions(quadraturePoints[q].first); StressVoigtVector stress; @@ -134,6 +140,42 @@ void VonMisesStress::handleEvent(core::objectmodel::Even } } +template +void VonMisesStress::precomputeData() +{ + if (this->l_topology == nullptr) return; + + auto restPositionAccessor = this->mstate->readRestPositions(); + const auto& restPosition = restPositionAccessor.ref(); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + m_precomputedData.resize(elements.size()); + + static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); + + for (std::size_t i = 0; i < elements.size(); ++i) + { + const auto& element = elements[i]; + std::array, NumberOfNodesInElement> nodeCoordinatesInElement; + for (sofa::Size n = 0; n < NumberOfNodesInElement; ++n) + nodeCoordinatesInElement[n] = restPosition[element[n]]; + + for (std::size_t j = 0; j < NumberOfQuadraturePoints; ++j) + { + const auto& dN_dq_ref = gradients[j]; + PrecomputedData& data = m_precomputedData[i][j]; + data.jacobian = sofa::fem::FiniteElementHelper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); + data.jacobianInv = sofa::type::inverse(data.jacobian); + data.detJacobian = sofa::type::absGeneralizedDeterminant(data.jacobian); + + for (sofa::Size n = 0; n < NumberOfNodesInElement; ++n) + { + data.dN_dQ[n] = data.jacobianInv.multTranspose(dN_dq_ref[n]); + } + } + } +} + template void VonMisesStress::calculateElementMassMatrix( const auto& elements, sofa::type::vector& elementMassMatrices) From 7e43f1fb62a5b0789f0cff2225aceb4dd43b1df2 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 30 Jul 2026 16:30:34 +0200 Subject: [PATCH 12/28] reorganize to optimize --- .../fem/elastic/VonMisesStress.inl | 28 +++++++++++-------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 2ef09b3b31d..f7f1f8dc761 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -92,12 +92,10 @@ void VonMisesStress::handleEvent(core::objectmodel::Even static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); static constexpr auto quadraturePoints = FiniteElement::quadraturePoints(); + std::array>, sofa::type::NumberOfIndependentElements> b; + for (auto& vec : b) vec.clear(); - for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) - { - sofa::type::Vec> b; - - for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) + for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) { const auto& weight = quadraturePoints[q].second; const auto& precomputed = m_precomputedData[elementId][q]; @@ -116,20 +114,28 @@ void VonMisesStress::handleEvent(core::objectmodel::Even // shape functions in the reference element evaluated at the quadrature point const auto N = FiniteElement::shapeFunctions(quadraturePoints[q].first); + // TODO: Compute the actual stress at the quadrature point based on F or other data StressVoigtVector stress; - for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + const auto commonFactor = weight * detJ; + for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) { - b[j] += N[j] * weight * detJ * stress[i]; + const auto stressFactor = commonFactor * stress[i]; + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + b[i][j] += N[j] * stressFactor; + } } } - const auto stressCoordinate = m_elementMassMatrices[elementId] * b; - for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) { - nodalStressInElement[j][i] = stressCoordinate[j]; + const auto stressCoordinate = m_elementMassMatrices[elementId] * b[i]; + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + nodalStressInElement[j][i] = stressCoordinate[j]; + } } - } for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) { From 46477658e88d6f92d8bb686c991fef5fa1b99379 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 30 Jul 2026 16:53:10 +0200 Subject: [PATCH 13/28] a data for a color map --- .../component/solidmechanics/fem/elastic/VonMisesStress.h | 6 +++++- .../solidmechanics/fem/elastic/VonMisesStress.inl | 7 +++++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index 6d4497d4b59..ff9aa2a4de7 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -23,8 +23,9 @@ #include #include -#include #include +#include +#include namespace sofa::component::solidmechanics::fem::elastic { @@ -56,6 +57,7 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh public: void init() override; + void draw(const core::visual::VisualParams* vparams) override; // A stress value for each node in an element using LocalStressValues = std::array, NumberOfNodesInElement>; @@ -63,6 +65,8 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh Data d_continuousField; + Data d_colorMap; + VonMisesStress(); void handleEvent(core::objectmodel::Event*) override; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index f7f1f8dc761..821215010bb 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -35,6 +35,7 @@ VonMisesStress::VonMisesStress() "Compute von Mises stress as a continuous field across the elements. Necessitate the solve " "of a sparse linear system. Otherwise, the von Mises stress is computed locally. A local " "stress value may indicate discretization errors if the field does not appear continuous.")) + , d_colorMap(initData(&d_colorMap, sofa::helper::ColorMap(), "colorMap", "Color map")) { // This component must receive events f_listening.setValue(true); @@ -270,4 +271,10 @@ VonMisesStress::vonMisesStress(const StressVoigtVector& return std::sqrt(static_cast>(3) / 2 * sofa::type::trace(s.multTranspose(s))); } +template +void VonMisesStress::draw(const core::visual::VisualParams* vparams) +{ + +} + } // namespace sofa::component::solidmechanics::fem::elastic From 806f9e6b8a9f191b1e6d54c7580e4effb43b9e91 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 09:52:53 +0200 Subject: [PATCH 14/28] support of more elements (empty for now) in DrawColoredMesh --- .../src/sofa/core/visual/DrawColoredMesh.h | 119 +++++++++++++++++- 1 file changed, 117 insertions(+), 2 deletions(-) diff --git a/Sofa/framework/Core/src/sofa/core/visual/DrawColoredMesh.h b/Sofa/framework/Core/src/sofa/core/visual/DrawColoredMesh.h index f0ac63f6a89..3bb49114094 100644 --- a/Sofa/framework/Core/src/sofa/core/visual/DrawColoredMesh.h +++ b/Sofa/framework/Core/src/sofa/core/visual/DrawColoredMesh.h @@ -67,6 +67,21 @@ struct BaseDrawColoredMesh protected: + template + static sofa::type::RGBAColor getColor(const ColorContainer& container, std::size_t elementId, std::size_t localId, std::size_t globalId) + { + if constexpr (std::is_same_v) + { + // Global approach + return container[globalId]; + } + else + { + // Local approach + return container[elementId][localId]; + } + } + template PositionType applyElementSpace(const PositionType& position, const PositionType& elementCenter) const { @@ -92,6 +107,60 @@ struct BaseDrawColoredMesh template struct DrawElementColoredMesh{}; +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Edge; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Triangle; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Quad; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + template<> struct SOFA_CORE_API DrawElementColoredMesh : public BaseDrawColoredMesh> @@ -139,7 +208,12 @@ struct SOFA_CORE_API DrawElementColoredMesh const auto p = this->applyElementSpace(position[vertexId], center); renderedPoints.push_back(p); - colors.push_back(nodesColors[vertexId]); + + // Find local ID of the vertex in the Tetrahedron to support local color mapping + std::size_t localId = 0; + for (; localId < 4; ++localId) if (element[localId] == vertexId) break; + + colors.push_back(this->getColor(nodesColors, i, localId, vertexId)); } } } @@ -191,7 +265,12 @@ struct SOFA_CORE_API DrawElementColoredMesh const auto p = this->applyElementSpace(position[vertexId], center); renderedPoints.push_back(p); - colors.push_back(nodesColors[vertexId]); + + // Find local ID of the vertex in the Hexahedron to support local color mapping + std::size_t localId = 0; + for (; localId < 8; ++localId) if (element[localId] == vertexId) break; + + colors.push_back(this->getColor(nodesColors, i, localId, vertexId)); } } } @@ -201,6 +280,42 @@ struct SOFA_CORE_API DrawElementColoredMesh }; +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Prism; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Pyramid; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + class SOFA_CORE_API DrawColoredMesh { public: From b5407768f397598c5a663faceb9511f4a93c267e Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 09:53:07 +0200 Subject: [PATCH 15/28] draw stresses in the component --- .../fem/elastic/VonMisesStress.inl | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 821215010bb..bffb364429b 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #include #include @@ -274,7 +276,42 @@ VonMisesStress::vonMisesStress(const StressVoigtVector& template void VonMisesStress::draw(const core::visual::VisualParams* vparams) { + const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); + const auto nodalStress = d_nodalStress.getValue(); + if (nodalStress.empty()) + return; + + const auto positions = this->mstate->readPositions(); + const auto& colorMap = d_colorMap.getValue(); + + Real_t minStress = std::numeric_limits>::max(); + Real_t maxStress = std::numeric_limits>::lowest(); + + for (const auto& elementStress : nodalStress) + { + for (const auto& stress : elementStress) + { + minStress = std::min(minStress, stress); + maxStress = std::max(maxStress, stress); + } + } + + const auto evaluator = colorMap.getEvaluator(minStress, maxStress); + + sofa::type::vector> nodesColors; + for (const auto& elementStress : nodalStress) + { + std::array nodesColorsInElement; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + nodesColorsInElement[i] = (evaluator(elementStress[i])); + } + nodesColors.push_back(nodesColorsInElement); + } + + core::visual::DrawElementColoredMesh renderer; + renderer.drawAllElements(vparams->drawTool(), positions.ref(), this->l_topology.get(), nodesColors); } } // namespace sofa::component::solidmechanics::fem::elastic From 901700eced29d1125651e871b91283d4ac95e51d Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 09:55:51 +0200 Subject: [PATCH 16/28] fix indentation --- .../fem/elastic/VonMisesStress.inl | 54 +++++++++---------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index bffb364429b..8f9bd701bbc 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -99,46 +99,46 @@ void VonMisesStress::handleEvent(core::objectmodel::Even for (auto& vec : b) vec.clear(); for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) - { - const auto& weight = quadraturePoints[q].second; - const auto& precomputed = m_precomputedData[elementId][q]; - - // gradient of shape functions in the reference element evaluated at the quadrature point - const auto& dN_dq_ref = gradients[q]; + { + const auto& weight = quadraturePoints[q].second; + const auto& precomputed = m_precomputedData[elementId][q]; - // jacobian of the mapping from the reference space to the CURRENT physical space - const auto J_q = FiniteElement::Helper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); + // gradient of shape functions in the reference element evaluated at the quadrature point + const auto& dN_dq_ref = gradients[q]; - // Deformation Gradient F = J_curr * J_rest_inv - const DeformationGradient F = J_q * precomputed.jacobianInv; + // jacobian of the mapping from the reference space to the CURRENT physical space + const auto J_q = FiniteElement::Helper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); - const auto detJ = sofa::type::absGeneralizedDeterminant(J_q); + // Deformation Gradient F = J_curr * J_rest_inv + const DeformationGradient F = J_q * precomputed.jacobianInv; - // shape functions in the reference element evaluated at the quadrature point - const auto N = FiniteElement::shapeFunctions(quadraturePoints[q].first); + const auto detJ = sofa::type::absGeneralizedDeterminant(J_q); - // TODO: Compute the actual stress at the quadrature point based on F or other data - StressVoigtVector stress; + // shape functions in the reference element evaluated at the quadrature point + const auto N = FiniteElement::shapeFunctions(quadraturePoints[q].first); - const auto commonFactor = weight * detJ; - for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) - { - const auto stressFactor = commonFactor * stress[i]; - for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) - { - b[i][j] += N[j] * stressFactor; - } - } - } + // TODO: Compute the actual stress at the quadrature point based on F or other data + StressVoigtVector stress; + const auto commonFactor = weight * detJ; for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) { - const auto stressCoordinate = m_elementMassMatrices[elementId] * b[i]; + const auto stressFactor = commonFactor * stress[i]; for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) { - nodalStressInElement[j][i] = stressCoordinate[j]; + b[i][j] += N[j] * stressFactor; } } + } + + for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) + { + const auto stressCoordinate = m_elementMassMatrices[elementId] * b[i]; + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + nodalStressInElement[j][i] = stressCoordinate[j]; + } + } for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) { From 047a2f46d8a79b99a959409d5ce11310dcb47dbf Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 10:19:44 +0200 Subject: [PATCH 17/28] introduce base class to compute stress --- .../SolidMechanics/FEM/Elastic/CMakeLists.txt | 1 + .../fem/elastic/CauchyStressEvaluator.h | 43 +++++++++++++++++++ .../fem/elastic/VonMisesStress.h | 7 +++ .../fem/elastic/VonMisesStress.inl | 29 +++++++++++++ 4 files changed, 80 insertions(+) create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index 77329f6a861..85447e27cdc 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/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 diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h new file mode 100644 index 00000000000..2d3645376ba --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h @@ -0,0 +1,43 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +class CauchyStressEvaluator : public core::objectmodel::BaseComponent +{ +public: + SOFA_ABSTRACT_CLASS(CauchyStressEvaluator, BaseComponent) + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + using StressVoigtVector = std::array, sofa::type::NumberOfIndependentElements>; + using DeformationGradient = sofa::type::Mat>; + + virtual StressVoigtVector computeStress(const DeformationGradient& F, sofa::Size elementId) = 0; +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index ff9aa2a4de7..39f40f50bbe 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -67,11 +68,17 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh Data d_colorMap; + sofa::SingleLink, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_stressEvaluator; + VonMisesStress(); void handleEvent(core::objectmodel::Event*) override; protected: + + void validateStressEvaluatorLink(); + using ElementMassMatrix = sofa::type::Mat>; sofa::type::vector m_elementMassMatrices; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 8f9bd701bbc..c7e8ededa04 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -21,6 +21,7 @@ ******************************************************************************/ #pragma once #include +#include #include #include #include @@ -38,6 +39,7 @@ VonMisesStress::VonMisesStress() "of a sparse linear system. Otherwise, the von Mises stress is computed locally. A local " "stress value may indicate discretization errors if the field does not appear continuous.")) , d_colorMap(initData(&d_colorMap, sofa::helper::ColorMap(), "colorMap", "Color map")) + , l_stressEvaluator(initLink("stressEvaluator", "The component in charge of evaluating the Cauchy stress.")) { // This component must receive events f_listening.setValue(true); @@ -53,6 +55,11 @@ void VonMisesStress::init() this->validateTopology(); } + if (!this->isComponentStateInvalid()) + { + validateStressEvaluatorLink(); + } + if (!this->isComponentStateInvalid()) { auto nodalStress = sofa::helper::getWriteOnlyAccessor(d_nodalStress); @@ -149,6 +156,28 @@ void VonMisesStress::handleEvent(core::objectmodel::Even } } +template +void VonMisesStress::validateStressEvaluatorLink() +{ + if (l_stressEvaluator.empty()) + { + msg_info() << "Link to a valid stress evaluator should be set to ensure right behavior. The first " + "stress evaluator found in current context will be used."; + l_stressEvaluator.set(this->getContext()->template get>()); + } + + if (l_stressEvaluator == nullptr) + { + msg_error() << "No stress evaluator component found at path: '" << this->l_stressEvaluator.getLinkedPath() + << "', nor in current context: " << this->getContext()->name + << ". Object must have a stress evaluator. " + << "The list of available stress evaluator components is: " + << sofa::core::ObjectFactory::getInstance() + ->listClassesDerivedFrom>(); + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } +} + template void VonMisesStress::precomputeData() { From 52c076f88101eef2857a5cad06510a20603bf1e9 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 10:30:33 +0200 Subject: [PATCH 18/28] LinearSmallStrainFEMForceField now derives from CauchyStressEvaluator and can be linked to VonMisesStress --- .../fem/elastic/CauchyStressEvaluator.h | 4 ++-- .../fem/elastic/LinearSmallStrainFEMForceField.h | 12 +++++++++--- .../fem/elastic/LinearSmallStrainFEMForceField.inl | 9 +++++++++ .../solidmechanics/fem/elastic/impl/trait.h | 4 ++++ examples/Demos/elongation.scn | 2 +- 5 files changed, 25 insertions(+), 6 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h index 2d3645376ba..c68d3c333c2 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h @@ -28,13 +28,13 @@ namespace sofa::component::solidmechanics::fem::elastic { template -class CauchyStressEvaluator : public core::objectmodel::BaseComponent +class CauchyStressEvaluator : public virtual core::objectmodel::BaseComponent { public: SOFA_ABSTRACT_CLASS(CauchyStressEvaluator, BaseComponent) static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; - using StressVoigtVector = std::array, sofa::type::NumberOfIndependentElements>; + using StressVoigtVector = sofa::type::Vec, sofa::Real_t>; using DeformationGradient = sofa::type::Mat>; virtual StressVoigtVector computeStress(const DeformationGradient& F, sofa::Size elementId) = 0; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h index 899edf4c084..50a833497ed 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h @@ -24,7 +24,7 @@ #include #include #include - +#include #include #if !defined(ELASTICITY_COMPONENT_ELEMENT_LINEAR_SMALL_STRAIN_FEM_FORCE_FIELD_CPP) @@ -37,13 +37,15 @@ namespace sofa::component::solidmechanics::fem::elastic template class LinearSmallStrainFEMForceField : public BaseElementLinearFEMForceField, - public FEMForceField + public FEMForceField, + public CauchyStressEvaluator { public: SOFA_CLASS2( SOFA_TEMPLATE2(LinearSmallStrainFEMForceField, DataTypes, ElementType), SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType), - SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType)); + SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType), + CauchyStressEvaluator); private: using trait = typename FEMForceField::trait; @@ -51,6 +53,8 @@ class LinearSmallStrainFEMForceField : 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( diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl index 6fea5667a99..64924cc48e7 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl @@ -180,4 +180,13 @@ void LinearSmallStrainFEMForceField::addKToMatrix( } } +template +auto LinearSmallStrainFEMForceField::computeStress( + const DeformationGradient& F, sofa::Size elementId) -> StressVoigtVector +{ + StressVoigtVector s; + + return s; +} + } // namespace sofa::component::solidmechanics::fem::elastic 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 ccc0b9f98a6..91b94764ba5 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 @@ -67,6 +67,10 @@ struct trait sofa::Real_t>; using ElementGradient = sofa::type::Vec>; + + using DeformationGradient = sofa::type::Mat>; + + using StressVoigtVector = sofa::type::Vec, sofa::Real_t>; }; } diff --git a/examples/Demos/elongation.scn b/examples/Demos/elongation.scn index c0d5deb631b..6ab03a433d8 100644 --- a/examples/Demos/elongation.scn +++ b/examples/Demos/elongation.scn @@ -25,7 +25,7 @@ - + From c918348ff3bf07c8f30b567bfb9732897397bb07 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 10:59:58 +0200 Subject: [PATCH 19/28] compute stress in linear elasticity --- .../LinearSmallStrainFEMForceField.inl | 23 +++++++++++++++++-- .../fem/elastic/VonMisesStress.h | 3 ++- .../fem/elastic/VonMisesStress.inl | 14 +++++++++-- examples/Demos/elongation.scn | 4 ++-- 4 files changed, 37 insertions(+), 7 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl index 64924cc48e7..f11f640cf75 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl @@ -184,9 +184,28 @@ template auto LinearSmallStrainFEMForceField::computeStress( const DeformationGradient& F, sofa::Size elementId) -> StressVoigtVector { - StressVoigtVector s; + const auto strainTensor = static_cast>(1)/2 * (F + F.transposed() - DeformationGradient::Identity()); - return s; + sofa::type::Vec, Real_t > strainVoigt; + for (sofa::Size i = 0; i < type::NumberOfIndependentElements; ++i) + { + const auto [p, q] = type::toTensorIndices(i); + strainVoigt[i] = strainTensor(p, q); + } + + const auto youngModulus = this->getYoungModulusInElement(elementId); + const auto poissonRatio = this->getPoissonRatioInElement(elementId); + + LameLambda> lambda { 0 }; + LameMu> mu { 0 }; + + sofa::component::solidmechanics::fem::elastic::toLameParameters>( + YoungModulus>(youngModulus), PoissonRatio>(poissonRatio), + lambda, mu); + + const auto elasticityTensor = makeIsotropicElasticityTensor>(mu, lambda); + + return elasticityTensor.toVoigtMatSym().toMat() * strainVoigt; } } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index 39f40f50bbe..c5a52a6060c 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -51,7 +51,7 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh static constexpr sofa::Size NumberOfQuadraturePoints = FiniteElement::quadraturePoints().size(); // a stress tensor represented as a vector using the Voigt mapping - using StressVoigtVector = std::array, sofa::type::NumberOfIndependentElements>; + using StressVoigtVector = sofa::type::Vec, sofa::Real_t>; using DeformationGradient = sofa::type::Mat>; @@ -59,6 +59,7 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh void init() override; void draw(const core::visual::VisualParams* vparams) override; + void computeBBox(const core::ExecParams* params, bool) override; // A stress value for each node in an element using LocalStressValues = std::array, NumberOfNodesInElement>; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index c7e8ededa04..946c71e2539 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -124,8 +124,7 @@ void VonMisesStress::handleEvent(core::objectmodel::Even // shape functions in the reference element evaluated at the quadrature point const auto N = FiniteElement::shapeFunctions(quadraturePoints[q].first); - // TODO: Compute the actual stress at the quadrature point based on F or other data - StressVoigtVector stress; + const StressVoigtVector stress = l_stressEvaluator->computeStress(F, elementId); const auto commonFactor = weight * detJ; for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) @@ -307,6 +306,8 @@ void VonMisesStress::draw(const core::visual::VisualPara { const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); + vparams->drawTool()->disableLighting(); + const auto nodalStress = d_nodalStress.getValue(); if (nodalStress.empty()) return; @@ -343,4 +344,13 @@ void VonMisesStress::draw(const core::visual::VisualPara renderer.drawAllElements(vparams->drawTool(), positions.ref(), this->l_topology.get(), nodesColors); } +template +void VonMisesStress::computeBBox(const core::ExecParams* params, bool) +{ + if (!this->mstate) return; + + const auto bbox = this->mstate->computeBBox(); //this may compute twice the mstate bbox, but there is no way to determine if the bbox has already been computed + this->f_bbox.setValue(std::move(bbox)); +} + } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/examples/Demos/elongation.scn b/examples/Demos/elongation.scn index 6ab03a433d8..a2ee63e61b0 100644 --- a/examples/Demos/elongation.scn +++ b/examples/Demos/elongation.scn @@ -1,5 +1,5 @@ - + @@ -29,5 +29,5 @@ - + From c5c02cea5372413cf51bfc7472f2420b179c425f Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 11:13:33 +0200 Subject: [PATCH 20/28] multiple fixes --- .../elastic/LinearSmallStrainFEMForceField.inl | 2 +- .../solidmechanics/fem/elastic/VonMisesStress.inl | 15 ++++++++------- examples/Demos/elongation.scn | 14 +++++++------- 3 files changed, 16 insertions(+), 15 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl index f11f640cf75..caeaaef83f8 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl @@ -184,7 +184,7 @@ template auto LinearSmallStrainFEMForceField::computeStress( const DeformationGradient& F, sofa::Size elementId) -> StressVoigtVector { - const auto strainTensor = static_cast>(1)/2 * (F + F.transposed() - DeformationGradient::Identity()); + const auto strainTensor = static_cast>(1)/2 * (F + F.transposed()) - DeformationGradient::Identity(); sofa::type::Vec, Real_t > strainVoigt; for (sofa::Size i = 0; i < type::NumberOfIndependentElements; ++i) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 946c71e2539..32c1c01317a 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -76,14 +76,15 @@ void VonMisesStress::handleEvent(core::objectmodel::Even { if (simulation::AnimateEndEvent::checkEventType(event)) { + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + const auto nbElements = elements.size(); + auto nodalStress = sofa::helper::getWriteOnlyAccessor(d_nodalStress); nodalStress->clear(); - nodalStress->resize(this->mstate->getSize()); + nodalStress->resize(nbElements); - auto restPositionAccessor = this->mstate->readRestPositions(); + auto positionAccessor = this->mstate->readPositions(); - const auto& elements = FiniteElement::getElementSequence(*this->l_topology); - const auto nbElements = elements.size(); helper::IotaView indices{static_cast(0ul), nbElements}; std::for_each( @@ -95,7 +96,7 @@ void VonMisesStress::handleEvent(core::objectmodel::Even std::array, NumberOfNodesInElement> nodeCoordinatesInElement; for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) { - nodeCoordinatesInElement[i] = restPositionAccessor[element[i]]; + nodeCoordinatesInElement[i] = positionAccessor[element[i]]; } std::array nodalStressInElement; @@ -220,7 +221,7 @@ void VonMisesStress::calculateElementMassMatrix( const auto nbElements = elements.size(); elementMassMatrices.resize(nbElements); - auto restPositionAccessor = this->mstate->readRestPositions(); + auto positionAccessor = this->mstate->readPositions(); SCOPED_TIMER("elementMassMatrix"); helper::IotaView indices{static_cast(0ul), nbElements}; @@ -234,7 +235,7 @@ void VonMisesStress::calculateElementMassMatrix( std::array, NumberOfNodesInElement> nodeCoordinatesInElement; for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) { - nodeCoordinatesInElement[i] = restPositionAccessor[element[i]]; + nodeCoordinatesInElement[i] = positionAccessor[element[i]]; } for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) diff --git a/examples/Demos/elongation.scn b/examples/Demos/elongation.scn index a2ee63e61b0..d466824ff85 100644 --- a/examples/Demos/elongation.scn +++ b/examples/Demos/elongation.scn @@ -1,5 +1,5 @@ - + @@ -12,20 +12,20 @@ - - + + + - + + - + From 9498cb56d69925348629db885abf4bdf07cd32f4 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 11:28:11 +0200 Subject: [PATCH 21/28] avoid reallocation --- .../solidmechanics/fem/elastic/CauchyStressEvaluator.h | 2 ++ .../sofa/component/solidmechanics/fem/elastic/VonMisesStress.h | 3 +++ .../component/solidmechanics/fem/elastic/VonMisesStress.inl | 3 +-- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h index c68d3c333c2..b890687e46b 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h @@ -24,6 +24,8 @@ #include #include +#include + namespace sofa::component::solidmechanics::fem::elastic { diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index c5a52a6060c..2ed857d6f73 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -26,6 +26,7 @@ #include #include #include +#include #include namespace sofa::component::solidmechanics::fem::elastic @@ -99,6 +100,8 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh static StressVoigtVector deviatoricStress(const StressVoigtVector& sigma); static Real_t vonMisesStress(const StressVoigtVector& deviatoricStress); + + core::visual::DrawElementColoredMesh m_renderer; }; } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 32c1c01317a..60286ec5e58 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -341,8 +341,7 @@ void VonMisesStress::draw(const core::visual::VisualPara nodesColors.push_back(nodesColorsInElement); } - core::visual::DrawElementColoredMesh renderer; - renderer.drawAllElements(vparams->drawTool(), positions.ref(), this->l_topology.get(), nodesColors); + m_renderer.drawAllElements(vparams->drawTool(), positions.ref(), this->l_topology.get(), nodesColors); } template From cd34f87a0eb074661103ddcb1879aa9cb2e718a1 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 11:53:33 +0200 Subject: [PATCH 22/28] finalize scene and a Data for lighting --- .../fem/elastic/VonMisesStress.h | 1 + .../fem/elastic/VonMisesStress.inl | 3 +- Sofa/GL/src/sofa/gl/DrawToolGL.cpp | 7 ++++ .../SolidMechanics/FEM/VonMisesStress.scn | 37 +++++++++++++++++++ examples/Demos/elongation.scn | 33 ----------------- 5 files changed, 47 insertions(+), 34 deletions(-) create mode 100644 examples/Component/SolidMechanics/FEM/VonMisesStress.scn delete mode 100644 examples/Demos/elongation.scn diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index 2ed857d6f73..d96847c53e0 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -69,6 +69,7 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh Data d_continuousField; Data d_colorMap; + Data d_lighting; sofa::SingleLink, sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_stressEvaluator; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 60286ec5e58..f5b310519ae 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -39,6 +39,7 @@ VonMisesStress::VonMisesStress() "of a sparse linear system. Otherwise, the von Mises stress is computed locally. A local " "stress value may indicate discretization errors if the field does not appear continuous.")) , d_colorMap(initData(&d_colorMap, sofa::helper::ColorMap(), "colorMap", "Color map")) + , d_lighting(initData(&d_lighting, true, "lighting", "If true, light is simulated on the mesh. Otherwise, no lighting effect.")) , l_stressEvaluator(initLink("stressEvaluator", "The component in charge of evaluating the Cauchy stress.")) { // This component must receive events @@ -307,7 +308,7 @@ void VonMisesStress::draw(const core::visual::VisualPara { const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); - vparams->drawTool()->disableLighting(); + vparams->drawTool()->setLightingEnabled(d_lighting.getValue()); const auto nodalStress = d_nodalStress.getValue(); if (nodalStress.empty()) diff --git a/Sofa/GL/src/sofa/gl/DrawToolGL.cpp b/Sofa/GL/src/sofa/gl/DrawToolGL.cpp index b339d3b1b95..7f2fbd705d4 100644 --- a/Sofa/GL/src/sofa/gl/DrawToolGL.cpp +++ b/Sofa/GL/src/sofa/gl/DrawToolGL.cpp @@ -979,6 +979,10 @@ void DrawToolGL::drawQuads(const std::vector &points, const type::RGBAColo void DrawToolGL::drawQuads(const std::vector &points, const std::vector& colors) { + glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE); + glEnable(GL_COLOR_MATERIAL); + setMaterial(colors[0]); + glBegin(GL_QUADS); { for (std::size_t i=0; i &points, const std::vector + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Demos/elongation.scn b/examples/Demos/elongation.scn deleted file mode 100644 index d466824ff85..00000000000 --- a/examples/Demos/elongation.scn +++ /dev/null @@ -1,33 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - From e9a55ddecf389efe3f4f13d7fa48691d1a7d5a84 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 11:55:27 +0200 Subject: [PATCH 23/28] remove unused Data --- .../component/solidmechanics/fem/elastic/VonMisesStress.h | 3 --- .../component/solidmechanics/fem/elastic/VonMisesStress.inl | 4 ---- 2 files changed, 7 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index d96847c53e0..ca44c646559 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -44,7 +44,6 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh protected: using FiniteElement = sofa::fem::FiniteElement; - static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; static constexpr sofa::Size NumberOfDofsInElement = NumberOfNodesInElement * spatial_dimensions; @@ -66,8 +65,6 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh using LocalStressValues = std::array, NumberOfNodesInElement>; Data> d_nodalStress; - Data d_continuousField; - Data d_colorMap; Data d_lighting; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index f5b310519ae..838112a74a8 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -34,10 +34,6 @@ template VonMisesStress::VonMisesStress() : d_nodalStress(initData(&d_nodalStress, sofa::type::vector{}, "nodalStress", "Local nodal von Mises stress values")) - , d_continuousField(initData(&d_continuousField, false, "continuousField", - "Compute von Mises stress as a continuous field across the elements. Necessitate the solve " - "of a sparse linear system. Otherwise, the von Mises stress is computed locally. A local " - "stress value may indicate discretization errors if the field does not appear continuous.")) , d_colorMap(initData(&d_colorMap, sofa::helper::ColorMap(), "colorMap", "Color map")) , d_lighting(initData(&d_lighting, true, "lighting", "If true, light is simulated on the mesh. Otherwise, no lighting effect.")) , l_stressEvaluator(initLink("stressEvaluator", "The component in charge of evaluating the Cauchy stress.")) From ea48d8941af82d2ba9a61301b96dc1b38ce389e2 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 12:13:08 +0200 Subject: [PATCH 24/28] refactoring to make it clearer --- .../fem/elastic/VonMisesStress.h | 20 ++++- .../fem/elastic/VonMisesStress.inl | 90 ++++++++++--------- 2 files changed, 67 insertions(+), 43 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index ca44c646559..ec9ca0f50d6 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -79,8 +79,8 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh void validateStressEvaluatorLink(); - using ElementMassMatrix = sofa::type::Mat>; - sofa::type::vector m_elementMassMatrices; + using ElementGramMatrix = sofa::type::Mat>; + sofa::type::vector m_elementInverseGramMatrices; struct PrecomputedData { @@ -94,7 +94,21 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh void precomputeData(); - void calculateElementMassMatrix(const auto& elements, sofa::type::vector &elementMassMatrices); + /** + * @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& 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 projectQuadraturePointValuesToNodes( + sofa::Size elementId, + const std::array& valuesAtQuadraturePoints) const; static StressVoigtVector deviatoricStress(const StressVoigtVector& sigma); static Real_t vonMisesStress(const StressVoigtVector& deviatoricStress); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 838112a74a8..fc516c77674 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -64,7 +64,7 @@ void VonMisesStress::init() nodalStress->resize(elements.size()); this->precomputeData(); - this->calculateElementMassMatrix(elements, m_elementMassMatrices); + this->calculateElementInverseGramMatrices(elements, m_elementInverseGramMatrices); } } @@ -96,53 +96,21 @@ void VonMisesStress::handleEvent(core::objectmodel::Even nodeCoordinatesInElement[i] = positionAccessor[element[i]]; } - std::array nodalStressInElement; static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); - static constexpr auto quadraturePoints = FiniteElement::quadraturePoints(); - - std::array>, sofa::type::NumberOfIndependentElements> b; - for (auto& vec : b) vec.clear(); + std::array stressAtQuadraturePoints; for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) { - const auto& weight = quadraturePoints[q].second; - const auto& precomputed = m_precomputedData[elementId][q]; - - // gradient of shape functions in the reference element evaluated at the quadrature point const auto& dN_dq_ref = gradients[q]; - - // jacobian of the mapping from the reference space to the CURRENT physical space const auto J_q = FiniteElement::Helper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); // Deformation Gradient F = J_curr * J_rest_inv - const DeformationGradient F = J_q * precomputed.jacobianInv; - - const auto detJ = sofa::type::absGeneralizedDeterminant(J_q); - - // shape functions in the reference element evaluated at the quadrature point - const auto N = FiniteElement::shapeFunctions(quadraturePoints[q].first); - - const StressVoigtVector stress = l_stressEvaluator->computeStress(F, elementId); - - const auto commonFactor = weight * detJ; - for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) - { - const auto stressFactor = commonFactor * stress[i]; - for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) - { - b[i][j] += N[j] * stressFactor; - } - } + const DeformationGradient F = J_q * m_precomputedData[elementId][q].jacobianInv; + stressAtQuadraturePoints[q] = l_stressEvaluator->computeStress(F, elementId); } - for (sofa::Size i = 0; i < sofa::type::NumberOfIndependentElements; ++i) - { - const auto stressCoordinate = m_elementMassMatrices[elementId] * b[i]; - for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) - { - nodalStressInElement[j][i] = stressCoordinate[j]; - } - } + // Least-square projection from quadrature points to nodes + const auto nodalStressInElement = projectQuadraturePointValuesToNodes(elementId, stressAtQuadraturePoints); for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) { @@ -212,8 +180,8 @@ void VonMisesStress::precomputeData() } template -void VonMisesStress::calculateElementMassMatrix( - const auto& elements, sofa::type::vector& elementMassMatrices) +void VonMisesStress::calculateElementInverseGramMatrices( + const auto& elements, sofa::type::vector& elementMassMatrices) { const auto nbElements = elements.size(); elementMassMatrices.resize(nbElements); @@ -262,6 +230,48 @@ void VonMisesStress::calculateElementMassMatrix( }); } +template +auto VonMisesStress::projectQuadraturePointValuesToNodes( + sofa::Size elementId, + const std::array& valuesAtQuadraturePoints) const +-> std::array +{ + using Real = sofa::Real_t; + static constexpr auto quadraturePoints = FiniteElement::quadraturePoints(); + static constexpr auto voigtSize = sofa::type::NumberOfIndependentElements; + + std::array projectedNodalValues; + for (auto& val : projectedNodalValues) val.clear(); + + // Project each Voigt component independently: M * x_i = b_i + for (sofa::Size i = 0; i < voigtSize; ++i) + { + sofa::type::Vec b; + b.clear(); + + for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) + { + const auto& weight = quadraturePoints[q].second; + const auto detJ = m_precomputedData[elementId][q].detJacobian; + const auto N = FiniteElement::shapeFunctions(quadraturePoints[q].first); + + const Real val = valuesAtQuadraturePoints[q][i] * weight * detJ; + for (sofa::Size node = 0; node < NumberOfNodesInElement; ++node) + { + b[node] += N[node] * val; + } + } + + const auto x = m_elementInverseGramMatrices[elementId] * b; + for (sofa::Size node = 0; node < NumberOfNodesInElement; ++node) + { + projectedNodalValues[node][i] = x[node]; + } + } + + return projectedNodalValues; +} + template auto VonMisesStress::deviatoricStress(const StressVoigtVector& sigma) -> StressVoigtVector { From 1a31cb6f24659db112e7714ff974f6d82dca7883 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 12:17:58 +0200 Subject: [PATCH 25/28] cleaning --- .../component/solidmechanics/fem/elastic/VonMisesStress.cpp | 2 +- .../component/solidmechanics/fem/elastic/VonMisesStress.h | 2 +- .../component/solidmechanics/fem/elastic/VonMisesStress.inl | 5 ++++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp index b252f6c4e0c..6eb5fe5af53 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp @@ -40,7 +40,7 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStressregisterObjects(sofa::core::ObjectRegistrationData("Compute von Mises stress") + factory->registerObjects(sofa::core::ObjectRegistrationData("Compute and draw von Mises stress based on a local least-square projection in each element") .add< VonMisesStress >() .add< VonMisesStress >() .add< VonMisesStress >() diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index ec9ca0f50d6..ae58da53a69 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -59,7 +59,7 @@ class VonMisesStress : public core::behavior::TopologyAccessor, public core::beh void init() override; void draw(const core::visual::VisualParams* vparams) override; - void computeBBox(const core::ExecParams* params, bool) override; + void computeBBox(const core::ExecParams* params, bool onlyVisible) override; // A stress value for each node in an element using LocalStressValues = std::array, NumberOfNodesInElement>; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index fc516c77674..5c8515b588c 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -352,8 +352,11 @@ void VonMisesStress::draw(const core::visual::VisualPara } template -void VonMisesStress::computeBBox(const core::ExecParams* params, bool) +void VonMisesStress::computeBBox(const core::ExecParams* params, bool onlyVisible) { + SOFA_UNUSED(params); + SOFA_UNUSED(onlyVisible); + if (!this->mstate) return; const auto bbox = this->mstate->computeBBox(); //this may compute twice the mstate bbox, but there is no way to determine if the bbox has already been computed From 86e0fbff9b23d95df727eb79e60bf0b1e71e751b Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 13:11:43 +0200 Subject: [PATCH 26/28] fix --- .../solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h index 50a833497ed..45bc796d644 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h @@ -41,7 +41,7 @@ class LinearSmallStrainFEMForceField : public CauchyStressEvaluator { public: - SOFA_CLASS2( + SOFA_CLASS3( SOFA_TEMPLATE2(LinearSmallStrainFEMForceField, DataTypes, ElementType), SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType), SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType), From c7b8aae680ed84edc7e4d5c53911e720e361cdc9 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 13:38:56 +0200 Subject: [PATCH 27/28] add documentation --- .../fem/elastic/VonMisesStress.h | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h index ae58da53a69..6ddcb89e844 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -32,6 +32,24 @@ 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 VonMisesStress : public core::behavior::TopologyAccessor, public core::behavior::SingleStateAccessor { From f204250bacde5dcdee08374b9221d7e96eb0a1ca Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 31 Jul 2026 14:17:24 +0200 Subject: [PATCH 28/28] add timer --- .../component/solidmechanics/fem/elastic/VonMisesStress.inl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl index 5c8515b588c..1070d6758d4 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -73,6 +73,8 @@ void VonMisesStress::handleEvent(core::objectmodel::Even { if (simulation::AnimateEndEvent::checkEventType(event)) { + SCOPED_TIMER("vonMisesStress"); + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); const auto nbElements = elements.size(); @@ -147,6 +149,7 @@ template void VonMisesStress::precomputeData() { if (this->l_topology == nullptr) return; + if (this->mstate == nullptr) return; auto restPositionAccessor = this->mstate->readRestPositions(); const auto& restPosition = restPositionAccessor.ref();