From 85042a682ddf70e2feb5de25c56d508c00731a83 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:07:19 +0200 Subject: [PATCH 01/17] [FEM] Introduce higher order quadrature rule switch for Triangle --- .../src/sofa/fem/FiniteElement[Triangle].h | 40 +++++++++++++++++-- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h index 0d9659040de..7495e802b26 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include namespace sofa::fem { @@ -64,11 +66,41 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - return { - std::make_pair(sofa::type::Vec(1./3., 1./3.), 1./2.) - }; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule (default). + return std::array{ + std::make_pair(ReferenceCoord(1./3., 1./3.), Real(1./2.)) + }; + } + else if constexpr (Degree <= 2) + { + // Degree 2: 3-point interior rule. + return std::array{ + std::make_pair(ReferenceCoord(1./6., 1./6.), Real(1./6.)), + std::make_pair(ReferenceCoord(2./3., 1./6.), Real(1./6.)), + std::make_pair(ReferenceCoord(1./6., 2./3.), Real(1./6.)) + }; + } + else + { + static_assert(Degree <= 2, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; From 3968e90dfa48226cfc254151aa54ae12088b8642 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:17:19 +0200 Subject: [PATCH 02/17] Edge --- .../FEM/src/sofa/fem/FiniteElement[Edge].h | 43 ++++++++++++++++--- 1 file changed, 38 insertions(+), 5 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h index fbffb5dbfed..492caa7a380 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_EDGE_CPP) #include @@ -55,12 +57,43 @@ struct FiniteElement return {{-static_cast(0.5)}, {static_cast(0.5)}}; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr sofa::type::Vec q0(static_cast(0)); - return { - std::make_pair(q0, static_cast(2)) - }; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point midpoint rule (default). + return std::array{ + std::make_pair(ReferenceCoord(static_cast(0)), static_cast(2)) + }; + } + else if constexpr (Degree <= 3) + { + // Degrees 2-3: 2-point Gauss-Legendre rule. + constexpr Real sqrt3 = 1.73205080757; + constexpr Real g = static_cast(1) / sqrt3; + return std::array{ + std::make_pair(ReferenceCoord(-g), static_cast(1)), + std::make_pair(ReferenceCoord( g), static_cast(1)) + }; + } + else + { + static_assert(Degree <= 3, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: + case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; From 1151c81e555d93bdd6eb224106ffb6aec1ad1ed0 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:23:24 +0200 Subject: [PATCH 03/17] Tetrahedra --- .../src/sofa/fem/FiniteElement[Tetrahedron].h | 44 +++++++++++++++++-- 1 file changed, 40 insertions(+), 4 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h index ac14bb9c59b..4ab10ae3ed5 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_TETAHEDRON_CPP) #include @@ -68,11 +70,45 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr sofa::type::Vec q0(1./4., 1./4., 1./4.); - constexpr std::array q { std::make_pair(q0, 1./6.) }; - return q; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule (default). + return std::array{ + std::make_pair(ReferenceCoord(1./4., 1./4., 1./4.), Real(1./6.)) + }; + } + else if constexpr (Degree <= 2) + { + // Degree 2: 4-point rule. + constexpr Real sqrt5 = 2.2360679774997896; + constexpr Real a = (5. - sqrt5) / 20.; + constexpr Real b = (5. + 3. * sqrt5) / 20.; + return std::array{ + std::make_pair(ReferenceCoord(a, a, a), Real(1./24.)), + std::make_pair(ReferenceCoord(b, a, a), Real(1./24.)), + std::make_pair(ReferenceCoord(a, b, a), Real(1./24.)), + std::make_pair(ReferenceCoord(a, a, b), Real(1./24.)) + }; + } + else + { + static_assert(Degree <= 2, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; From a1ec92950380db9a352d6591e3ac9b87edcb5d81 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:41:34 +0200 Subject: [PATCH 04/17] Quad --- .../FEM/src/sofa/fem/FiniteElement[Quad].h | 67 +++++++++++++++---- 1 file changed, 55 insertions(+), 12 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h index 4e37585eace..1485a881001 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_QUAD_CPP) #include @@ -67,21 +69,62 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.) - constexpr Real sqrt6 = 2.44948974278; //sqrt(6.) - constexpr Real sqrt2 = 1.41421356237; //sqrt(2.) + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule. + return std::array{ + std::make_pair(ReferenceCoord(static_cast(0), static_cast(0)), static_cast(4)) + }; + } + else if constexpr (Degree <= 2) + { + // Degree 2: 3-point rule (default). + constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.) + constexpr Real sqrt6 = 2.44948974278; //sqrt(6.) + constexpr Real sqrt2 = 1.41421356237; //sqrt(2.) - constexpr sofa::type::Vec q0(sqrt2_3, 0.); - constexpr sofa::type::Vec q1(-1/sqrt6, -1./sqrt2); - constexpr sofa::type::Vec q2(-1/sqrt6, 1./sqrt2); + constexpr ReferenceCoord q0(sqrt2_3, 0.); + constexpr ReferenceCoord q1(-1/sqrt6, -1./sqrt2); + constexpr ReferenceCoord q2(-1/sqrt6, 1./sqrt2); - return { - std::make_pair(q0, 4./3.), - std::make_pair(q1, 4./3.), - std::make_pair(q2, 4./3.), - }; + return std::array{ + std::make_pair(q0, 4./3.), + std::make_pair(q1, 4./3.), + std::make_pair(q2, 4./3.) + }; + } + else if constexpr (Degree <= 3) + { + // Degree 3: 2x2 Gauss-Legendre rule. + constexpr Real sqrt3 = 1.73205080757; //sqrt(3.) + constexpr Real g = static_cast(1) / sqrt3; + return std::array{ + std::make_pair(ReferenceCoord(-g, -g), static_cast(1)), + std::make_pair(ReferenceCoord( g, -g), static_cast(1)), + std::make_pair(ReferenceCoord( g, g), static_cast(1)), + std::make_pair(ReferenceCoord(-g, g), static_cast(1)) + }; + } + else + { + static_assert(Degree <= 3, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; } + case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; From e42d7e8b2da9f0bb71368bc79f4ea0ad573b6532 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:41:41 +0200 Subject: [PATCH 05/17] Hexa --- .../src/sofa/fem/FiniteElement[Hexahedron].h | 74 +++++++++++++++---- 1 file changed, 60 insertions(+), 14 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h index 38339624e47..e3ab476e0f8 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_HEXAHEDRON_CPP) #include @@ -96,24 +98,68 @@ struct FiniteElement return gradient; } + template static constexpr auto quadraturePoints() { - constexpr Real sqrt3 = 1.73205080757; //sqrt(3.) - constexpr Real sqrt3_1 = static_cast(1) / sqrt3; - constexpr Real one = static_cast(1); + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule. + return std::array{ + std::make_pair(ReferenceCoord(static_cast(0), static_cast(0), static_cast(0)), static_cast(8)) + }; + } + else if constexpr (Degree <= 3) + { + // Degrees 2-3: 2x2x2 Gauss-Legendre rule (default). + constexpr Real sqrt3 = 1.73205080757; //sqrt(3.) + constexpr Real sqrt3_1 = static_cast(1) / sqrt3; + constexpr Real one = static_cast(1); - constexpr std::array q { - std::pair{referenceElementNodes[0] * sqrt3_1, one}, - std::pair{referenceElementNodes[1] * sqrt3_1, one}, - std::pair{referenceElementNodes[2] * sqrt3_1, one}, - std::pair{referenceElementNodes[3] * sqrt3_1, one}, - std::pair{referenceElementNodes[4] * sqrt3_1, one}, - std::pair{referenceElementNodes[5] * sqrt3_1, one}, - std::pair{referenceElementNodes[6] * sqrt3_1, one}, - std::pair{referenceElementNodes[7] * sqrt3_1, one}, - }; + return std::array { + std::pair{referenceElementNodes[0] * sqrt3_1, one}, + std::pair{referenceElementNodes[1] * sqrt3_1, one}, + std::pair{referenceElementNodes[2] * sqrt3_1, one}, + std::pair{referenceElementNodes[3] * sqrt3_1, one}, + std::pair{referenceElementNodes[4] * sqrt3_1, one}, + std::pair{referenceElementNodes[5] * sqrt3_1, one}, + std::pair{referenceElementNodes[6] * sqrt3_1, one}, + std::pair{referenceElementNodes[7] * sqrt3_1, one}, + }; + } + else if constexpr (Degree <= 5) + { + // Degrees 4-5: 3x3x3 Gauss-Legendre rule. + constexpr Real g = 0.77459666924; //sqrt(3./5.) + constexpr std::array node{ -g, static_cast(0), g }; + constexpr std::array weight{ static_cast(5./9.), static_cast(8./9.), static_cast(5./9.) }; + + std::array q{}; + sofa::Size k = 0; + for (sofa::Size i = 0; i < 3; ++i) + for (sofa::Size j = 0; j < 3; ++j) + for (sofa::Size l = 0; l < 3; ++l) + q[k++] = std::make_pair(ReferenceCoord(node[i], node[j], node[l]), weight[i] * weight[j] * weight[l]); + return q; + } + else + { + static_assert(Degree <= 5, "FiniteElement: no quadrature rule for the requested degree"); + } + } - return q; + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: + case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; } + case 4: + case 5: { static constexpr auto rule = quadraturePoints<5>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; From 9c42dbf2373723b6129a51d30518fe643c2cbf33 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 11 Aug 2026 15:28:18 +0200 Subject: [PATCH 06/17] Use uint for degree --- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h | 4 ++-- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h | 4 ++-- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h | 4 ++-- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h | 4 ++-- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h | 4 ++-- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h index 492caa7a380..53570c30d05 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h @@ -57,7 +57,7 @@ struct FiniteElement return {{-static_cast(0.5)}, {static_cast(0.5)}}; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -84,7 +84,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h index e3ab476e0f8..d855c9f8c72 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h @@ -98,7 +98,7 @@ struct FiniteElement return gradient; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -148,7 +148,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h index 1485a881001..261eeae7036 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h @@ -69,7 +69,7 @@ struct FiniteElement }; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -115,7 +115,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h index 4ab10ae3ed5..bca61671e0a 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h @@ -70,7 +70,7 @@ struct FiniteElement }; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -100,7 +100,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h index 7495e802b26..84c9ec19a74 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h @@ -66,7 +66,7 @@ struct FiniteElement }; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -92,7 +92,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { From 2780ef1257c41645062db1de98abaab6e10adaf4 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 21 Jul 2026 14:41:34 +0200 Subject: [PATCH 07/17] [SolidMechanics] Add FEMSourceTerm component --- .../SolidMechanics/FEM/Elastic/CMakeLists.txt | 3 + .../fem/elastic/FEMSourceTerm.cpp | 58 ++++++ .../fem/elastic/FEMSourceTerm.h | 158 ++++++++++++++++ .../fem/elastic/FEMSourceTerm.inl | 169 ++++++++++++++++++ .../solidmechanics/fem/elastic/init.cpp | 2 + 5 files changed, 390 insertions(+) create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.cpp create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index fddd6033a0d..829fcc8475a 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt @@ -15,6 +15,8 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTerm.h + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTerm.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.h @@ -69,6 +71,7 @@ set(SOURCE_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTerm.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMForceField.cpp diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.cpp new file mode 100644 index 00000000000..31935713335 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.cpp @@ -0,0 +1,58 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_CPP + +#include + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerFEMSourceTerm(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Consistent nodal load of a source term, integrated with the finite-element quadrature from a nodal source field") + .add< FEMSourceTerm >() + .add< FEMSourceTerm >() + .add< FEMSourceTerm >() + .add< FEMSourceTerm >() + .add< FEMSourceTerm >() + .add< FEMSourceTerm >() + .add< FEMSourceTerm >() + .add< FEMSourceTerm >() + .add< FEMSourceTerm >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h new file mode 100644 index 00000000000..b99a9c42173 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h @@ -0,0 +1,158 @@ +/****************************************************************************** +* 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 +#include +#include +#include + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class FEMSourceTerm + * @brief Computes nodal source terms by integrating a given source density field stored at the nodes. + * + * This class assembles and stores a geometric matrix M using a quadrature rule over the element domain. + * The matrix is then used to multiply the source density b to compute a nodal source term F that + * contributes to the RHS of the weak form through the addForce function. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + * @tparam TElementType The type of finite element (e.g., sofa::geometry::Tetrahedron). + */ +template +class FEMSourceTerm : + public sofa::core::behavior::ForceField, + public virtual sofa::core::behavior::TopologyAccessor +{ +public: + using DataTypes = TDataTypes; + using ElementType = TElementType; + SOFA_CLASS2(SOFA_TEMPLATE2(FEMSourceTerm, DataTypes, ElementType), + sofa::core::behavior::ForceField, + sofa::core::behavior::TopologyAccessor); + +protected: + using FiniteElement = sofa::fem::FiniteElement; + + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + + using GlobalMatrix = sofa::linearalgebra::CompressedRowSparseMatrixMechanical>; + +public: + + /** + * @brief Initializes the component. + * + * This method performs several initialization steps: + * 1. Initializes the base force field. + * 2. Initializes the topology accessor. + * 3. Resizes the nodal source density. + * 4. Assembles the global matrix M. + */ + void init() override; + + /** + * @brief Adds the force (f = M * b) to the RHS vector. + * + * This method computes the product of the geometric matrix and the source density vector, + * adding the result to the force vector `f`. + * + * @param mparams Mechanical parameters for the computation. + * @param f The force vector to which the source term will be added. + * @param x The current positions (unused in this implementation). + * @param v The current velocities (unused in this implementation). + */ + void addForce( + const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& f, + const sofa::DataVecCoord_t& x, + const sofa::DataVecDeriv_t& v) override; + + /** + * @brief A no-op as the load is prescribed on the rest configuration + */ + void addDForce(const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& df, + const sofa::DataVecDeriv_t& dx) override; + + /** + * @brief A no-op as the load is prescribed on the rest configuration + */ + void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; + + using sofa::core::behavior::ForceField::getPotentialEnergy; + /** + * @brief Not implemented, returns 0. + */ + SReal getPotentialEnergy(const sofa::core::MechanicalParams* mparams, + const sofa::DataVecCoord_t& x) const override; + + /** + * @brief Source term (per unit volume) sampled at each node. + */ + sofa::Data > d_nodalSourceDensity; + +protected: + + /** + * @brief Default constructor. + */ + FEMSourceTerm(); + + /** + * @brief Resizes the source density to the size of the mechanical state + */ + void resizeNodalSourceDensity(const std::size_t size); + + /** + * @brief Assembles and stores the geometry-only matrix \f$ M_{ij} = \int_{\Omega} N_i N_j \, d\Omega \f$ over each element on the rest configuration. + */ + void assembleGlobalMatrix(); + + /** + * @brief Geometry-only matrix \f$ M_{ij} = \int_{\Omega} N_i N_j \, d\Omega \f$ of the system. + * + * Stored in compressed sparse row format. Assembled once in init on the rest configuration. + */ + GlobalMatrix m_globalMatrix; +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl new file mode 100644 index 00000000000..d55cf8bdf7e --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl @@ -0,0 +1,169 @@ +/****************************************************************************** +* 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 +FEMSourceTerm::FEMSourceTerm() + : d_nodalSourceDensity(initData(&d_nodalSourceDensity, "nodalSourceDensity", + "Source term (per unit volume) sampled at each node. Interpolated inside the " + "element with the shape functions and integrated on the reference configuration.")) +{ +} + +template +void FEMSourceTerm::init() +{ + sofa::core::behavior::ForceField::init(); + + if (!this->isComponentStateInvalid()) + { + sofa::core::behavior::TopologyAccessor::init(); + } + + if (!this->isComponentStateInvalid() && this->mstate) + { + this->resizeNodalSourceDensity(this->mstate->getSize()); + } + + if (!this->isComponentStateInvalid() && this->l_topology && this->mstate) + { + this->assembleGlobalMatrix(); + } + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + } +} + +template +void FEMSourceTerm::resizeNodalSourceDensity(const std::size_t size) +{ + sofa::helper::WriteAccessor nodalSourceDensity = sofa::helper::getWriteAccessor(d_nodalSourceDensity); + + if (nodalSourceDensity.size() < size) + { + nodalSourceDensity.resize(size, sofa::Deriv_t{}); + } +} + +template +void FEMSourceTerm::assembleGlobalMatrix() +{ + // M_ij = integral of N_i N_j dV, evaluated on the rest configuration (geometry only). + const auto restPositionsAccessor = this->mstate->readRestPositions(); + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + + m_globalMatrix.clear(); + const auto size = this->mstate->getSize(); + m_globalMatrix.resize(size, size); + + for (const auto& element : elements) + { + const std::array, NumberOfNodesInElement> elementNodesRestCoordinates = + extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); + + sofa::type::Mat> elementMatrix; + + for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) + { + const auto N = FiniteElement::shapeFunctions(quadraturePoint); + const auto dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + + const auto jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + elementNodesRestCoordinates, dN_dq_ref); + const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); + + elementMatrix += (static_cast>(weight) * static_cast>(detJ)) * sofa::type::dyad(N, N); + } + + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + m_globalMatrix.add(element[i], element[j], elementMatrix(i, j)); + } + } + } + + m_globalMatrix.compress(); +} + +template +void FEMSourceTerm::addForce(const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& f, + const sofa::DataVecCoord_t& x, + const sofa::DataVecDeriv_t& v) +{ + SOFA_UNUSED(mparams); + SOFA_UNUSED(x); + SOFA_UNUSED(v); + + const sofa::helper::ReadAccessor nodalSourceDensity = sofa::helper::getReadAccessor(d_nodalSourceDensity); + auto forceAccessor = sofa::helper::getWriteAccessor(f); + + // f_i = sum_j M_ij b_j : apply the cached matrix to the nodal source density. + for (std::size_t xi = 0; xi < m_globalMatrix.rowIndex.size(); ++xi) + { + const auto rowId = m_globalMatrix.rowIndex[xi]; + typename GlobalMatrix::Range rowRange(m_globalMatrix.rowBegin[xi], m_globalMatrix.rowBegin[xi + 1]); + for (typename GlobalMatrix::Index xj = rowRange.begin(); xj < rowRange.end(); ++xj) + { + const auto columnId = m_globalMatrix.colsIndex[xj]; + const auto& value = m_globalMatrix.colsValue[xj]; + + forceAccessor[rowId] += nodalSourceDensity[columnId] * value; + } + } +} + +template +void FEMSourceTerm::addDForce(const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& df, + const sofa::DataVecDeriv_t& dx) +{ + SOFA_UNUSED(mparams); + SOFA_UNUSED(df); + SOFA_UNUSED(dx); +} + +template +void FEMSourceTerm::buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) +{ + SOFA_UNUSED(matrix); +} + +template +SReal FEMSourceTerm::getPotentialEnergy(const sofa::core::MechanicalParams* mparams, + const sofa::DataVecCoord_t& x) const +{ + SOFA_UNUSED(mparams); + SOFA_UNUSED(x); + return 0.0; +} + +} // 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..95e84fb5d1b 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 @@ -28,6 +28,7 @@ namespace sofa::component::solidmechanics::fem::elastic extern void registerBeamFEMForceField(sofa::core::ObjectFactory* factory); extern void registerCorotationalFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerFEMSourceTerm(sofa::core::ObjectFactory* factory); extern void registerFastTetrahedralCorotationalForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedralFEMForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedralFEMForceFieldAndMass(sofa::core::ObjectFactory* factory); @@ -68,6 +69,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) { registerBeamFEMForceField(factory); registerCorotationalFEMForceField(factory); + registerFEMSourceTerm(factory); registerFastTetrahedralCorotationalForceField(factory); registerHexahedralFEMForceField(factory); registerHexahedralFEMForceFieldAndMass(factory); From e8b080deee5c593f988b46a3960ac793365a870c Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 21 Jul 2026 18:19:01 +0200 Subject: [PATCH 08/17] Split assembly in two steps --- .../fem/elastic/FEMSourceTerm.h | 11 +++++ .../fem/elastic/FEMSourceTerm.inl | 48 +++++++++++++++---- 2 files changed, 49 insertions(+), 10 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h index b99a9c42173..724881a4bd2 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h @@ -62,6 +62,7 @@ class FEMSourceTerm : static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + using ElementMatrix = sofa::type::Mat>; using GlobalMatrix = sofa::linearalgebra::CompressedRowSparseMatrixMechanical>; public: @@ -135,6 +136,16 @@ class FEMSourceTerm : */ void assembleGlobalMatrix(); + /** + * @brief Computes the geometry-only matrix of each element. + */ + void calculateElementMatrix(const auto& elements, sofa::type::vector& elementMatrices); + + /** + * @brief Scatters the element matrices into the global matrix. + */ + void initializeGlobalMatrix(const auto& elements, const sofa::type::vector& elementMatrices); + /** * @brief Geometry-only matrix \f$ M_{ij} = \int_{\Omega} N_i N_j \, d\Omega \f$ of the system. * diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl index d55cf8bdf7e..528e0079750 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl @@ -74,21 +74,32 @@ void FEMSourceTerm::resizeNodalSourceDensity(const std:: template void FEMSourceTerm::assembleGlobalMatrix() { - // M_ij = integral of N_i N_j dV, evaluated on the rest configuration (geometry only). - const auto restPositionsAccessor = this->mstate->readRestPositions(); const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + sofa::type::vector elementMatrices; - m_globalMatrix.clear(); - const auto size = this->mstate->getSize(); - m_globalMatrix.resize(size, size); + // 1. compute the geometry-only matrix of each element + calculateElementMatrix(elements, elementMatrices); + + // 2. scatter the element matrices into the global matrix + initializeGlobalMatrix(elements, elementMatrices); +} - for (const auto& element : elements) +template +void FEMSourceTerm::calculateElementMatrix( + const auto& elements, sofa::type::vector& elementMatrices) +{ + const auto restPositionsAccessor = this->mstate->readRestPositions(); + elementMatrices.resize(elements.size()); + + for (std::size_t elementId = 0; elementId < elements.size(); ++elementId) { + const auto& element = elements[elementId]; + auto& elementMatrix = elementMatrices[elementId]; + const std::array, NumberOfNodesInElement> elementNodesRestCoordinates = extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); - sofa::type::Mat> elementMatrix; - + // M_ij = integral of N_i N_j dV, evaluated on the rest configuration (geometry only). for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) { const auto N = FiniteElement::shapeFunctions(quadraturePoint); @@ -98,8 +109,25 @@ void FEMSourceTerm::assembleGlobalMatrix() elementNodesRestCoordinates, dN_dq_ref); const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); - elementMatrix += (static_cast>(weight) * static_cast>(detJ)) * sofa::type::dyad(N, N); + const auto NT_N = sofa::type::dyad(N, N); + + elementMatrix += (weight * detJ) * NT_N; } + } +} + +template +void FEMSourceTerm::initializeGlobalMatrix( + const auto& elements, const sofa::type::vector& elementMatrices) +{ + m_globalMatrix.clear(); + const auto size = this->mstate->getSize(); + m_globalMatrix.resize(size, size); + + for (std::size_t elementId = 0; elementId < elements.size(); ++elementId) + { + const auto& element = elements[elementId]; + const auto& elementMatrix = elementMatrices[elementId]; for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) { @@ -126,7 +154,7 @@ void FEMSourceTerm::addForce(const sofa::core::Mechanica const sofa::helper::ReadAccessor nodalSourceDensity = sofa::helper::getReadAccessor(d_nodalSourceDensity); auto forceAccessor = sofa::helper::getWriteAccessor(f); - // f_i = sum_j M_ij b_j : apply the cached matrix to the nodal source density. + // f_i = sum_j M_ij b_j : apply the global matrix to the nodal source density. for (std::size_t xi = 0; xi < m_globalMatrix.rowIndex.size(); ++xi) { const auto rowId = m_globalMatrix.rowIndex[xi]; From c3bcc7ffdad3fd128b7db6ea16ae7d12d5360702 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 13 Aug 2026 11:17:32 +0200 Subject: [PATCH 09/17] [SolidMechanics] Enable choosing quadrature order in FEMSourceTerm --- .../solidmechanics/fem/elastic/FEMSourceTerm.h | 5 +++++ .../solidmechanics/fem/elastic/FEMSourceTerm.inl | 16 ++++++++++++++-- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h index 724881a4bd2..9b4e8f358d6 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h @@ -119,6 +119,11 @@ class FEMSourceTerm : */ sofa::Data > d_nodalSourceDensity; + /** + * @brief Degree of the quadrature rule integrating the element matrix M. + */ + sofa::Data d_quadratureDegree; + protected: /** diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl index 528e0079750..6e029f0bb98 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl @@ -28,10 +28,20 @@ namespace sofa::component::solidmechanics::fem::elastic template FEMSourceTerm::FEMSourceTerm() - : d_nodalSourceDensity(initData(&d_nodalSourceDensity, "nodalSourceDensity", + : d_nodalSourceDensity(initData(&d_nodalSourceDensity, "nodalSourceDensity", "Source term (per unit volume) sampled at each node. Interpolated inside the " "element with the shape functions and integrated on the reference configuration.")) + , d_quadratureDegree(initData(&d_quadratureDegree, static_cast(1), "quadratureDegree", + "Degree of the quadrature rule integrating the element matrix M.")) { + this->addUpdateCallback("reassembleSourceMatrix", {&d_quadratureDegree}, + [this](const sofa::core::DataTracker&) + { + if (!this->isComponentStateInvalid() && this->l_topology && this->mstate) + assembleGlobalMatrix(); + + return this->getComponentState(); + }, {}); } template @@ -91,6 +101,8 @@ void FEMSourceTerm::calculateElementMatrix( const auto restPositionsAccessor = this->mstate->readRestPositions(); elementMatrices.resize(elements.size()); + const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue()); + for (std::size_t elementId = 0; elementId < elements.size(); ++elementId) { const auto& element = elements[elementId]; @@ -100,7 +112,7 @@ void FEMSourceTerm::calculateElementMatrix( extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); // M_ij = integral of N_i N_j dV, evaluated on the rest configuration (geometry only). - for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) + for (const auto& [quadraturePoint, weight] : quadratureRule) { const auto N = FiniteElement::shapeFunctions(quadraturePoint); const auto dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); From 203b95fe6e38555c474e6d9146c078b91aa11027 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 1 Sep 2026 13:33:47 +0200 Subject: [PATCH 10/17] [FEM] Add constant source term integration --- .../SolidMechanics/FEM/Elastic/CMakeLists.txt | 8 +- .../fem/elastic/ConstantSourceTerm.cpp | 45 ++++++ .../fem/elastic/ConstantSourceTerm.h | 66 +++++++++ ...ceTerm.cpp => FEMSourceTermIntegrator.cpp} | 44 +++--- ...SourceTerm.h => FEMSourceTermIntegrator.h} | 102 +++++++++----- ...ceTerm.inl => FEMSourceTermIntegrator.inl} | 131 +++++++++++++----- .../solidmechanics/fem/elastic/init.cpp | 6 +- .../FEM/FEMSourceTermIntegrator.scn | 110 +++++++++++++++ 8 files changed, 413 insertions(+), 99 deletions(-) create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ConstantSourceTerm.cpp create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ConstantSourceTerm.h rename Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/{FEMSourceTerm.cpp => FEMSourceTermIntegrator.cpp} (61%) rename Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/{FEMSourceTerm.h => FEMSourceTermIntegrator.h} (63%) rename Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/{FEMSourceTerm.inl => FEMSourceTermIntegrator.inl} (65%) create mode 100644 examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index 829fcc8475a..4f1356b6330 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt @@ -13,10 +13,11 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ConstantSourceTerm.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.inl - ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTerm.h - ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTerm.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.h + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.h @@ -70,8 +71,9 @@ set(SOURCE_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseElementLinearFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ConstantSourceTerm.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.cpp - ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTerm.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMForceField.cpp diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ConstantSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ConstantSourceTerm.cpp new file mode 100644 index 00000000000..569b3d7e86f --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ConstantSourceTerm.cpp @@ -0,0 +1,45 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_CONSTANT_SOURCE_TERM_CPP + +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerConstantSourceTerm(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Constant source term (per unit volume) prescribed at the nodes") + .add< ConstantSourceTerm >() + .add< ConstantSourceTerm >() + .add< ConstantSourceTerm >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ConstantSourceTerm; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ConstantSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ConstantSourceTerm.h new file mode 100644 index 00000000000..e03bf47622a --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ConstantSourceTerm.h @@ -0,0 +1,66 @@ +/****************************************************************************** +* 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 +#include + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_CONSTANT_SOURCE_TERM_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class ConstantSourceTerm + * @brief A source density prescribed at the nodes, independent of the current displacement. + * + * The density is the inherited "property" Data (see BaseNodalProperty): a vector shorter than the + * mechanical state broadcasts its last value to the remaining nodes, so a uniform density is + * written with a single value. Link it to a FEMSourceTermIntegrator through l_constantSources. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + */ +template +class ConstantSourceTerm : public sofa::core::BaseNodalProperty> +{ +public: + using DataTypes = TDataTypes; + using Deriv = sofa::Deriv_t; + + SOFA_CLASS(SOFA_TEMPLATE(ConstantSourceTerm, DataTypes), + SOFA_TEMPLATE(sofa::core::BaseNodalProperty, Deriv)); + +protected: + + ConstantSourceTerm() : sofa::core::BaseNodalProperty(Deriv{}) {} +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_CONSTANT_SOURCE_TERM_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ConstantSourceTerm; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.cpp similarity index 61% rename from Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.cpp rename to Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.cpp index 31935713335..f7a908ed703 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.cpp @@ -19,9 +19,9 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#define SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_CPP +#define SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP -#include +#include #include #include @@ -30,29 +30,29 @@ namespace sofa::component::solidmechanics::fem::elastic { -void registerFEMSourceTerm(sofa::core::ObjectFactory* factory) +void registerFEMSourceTermIntegrator(sofa::core::ObjectFactory* factory) { - factory->registerObjects(sofa::core::ObjectRegistrationData("Consistent nodal load of a source term, integrated with the finite-element quadrature from a nodal source field") - .add< FEMSourceTerm >() - .add< FEMSourceTerm >() - .add< FEMSourceTerm >() - .add< FEMSourceTerm >() - .add< FEMSourceTerm >() - .add< FEMSourceTerm >() - .add< FEMSourceTerm >() - .add< FEMSourceTerm >() - .add< FEMSourceTerm >() + factory->registerObjects(sofa::core::ObjectRegistrationData("Consistent nodal load of a source term, integrated with the finite-element quadrature from a nodal source density field") + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() ); } -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h similarity index 63% rename from Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h rename to Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h index 9b4e8f358d6..7873aec5cdf 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h @@ -22,12 +22,14 @@ #pragma once #include +#include #include #include +#include #include #include -#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_CPP) +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP) #include #endif @@ -35,58 +37,69 @@ namespace sofa::component::solidmechanics::fem::elastic { /** - * @class FEMSourceTerm - * @brief Computes nodal source terms by integrating a given source density field stored at the nodes. + * @class FEMSourceTermIntegrator + * @brief Integrates source terms into consistent nodal loads. * - * This class assembles and stores a geometric matrix M using a quadrature rule over the element domain. - * The matrix is then used to multiply the source density b to compute a nodal source term F that - * contributes to the RHS of the weak form through the addForce function. + * A source term contributes \f$ \int_{\Omega} N_a \, r \, d\Omega \f$ to the right-hand side, where + * r is the per-node density carried by a linked ConstantSourceTerm (through l_constantSources) and + * does not depend on the displacement. Every term is thus summed and integrated once in init(); + * addForce merely accumulates the result. * * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). * @tparam TElementType The type of finite element (e.g., sofa::geometry::Tetrahedron). */ template -class FEMSourceTerm : +class FEMSourceTermIntegrator : public sofa::core::behavior::ForceField, public virtual sofa::core::behavior::TopologyAccessor { public: using DataTypes = TDataTypes; using ElementType = TElementType; - SOFA_CLASS2(SOFA_TEMPLATE2(FEMSourceTerm, DataTypes, ElementType), + SOFA_CLASS2(SOFA_TEMPLATE2(FEMSourceTermIntegrator, DataTypes, ElementType), sofa::core::behavior::ForceField, sofa::core::behavior::TopologyAccessor); protected: using FiniteElement = sofa::fem::FiniteElement; + using Real = sofa::Real_t; static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; using ElementMatrix = sofa::type::Mat>; using GlobalMatrix = sofa::linearalgebra::CompressedRowSparseMatrixMechanical>; public: + /** + * @brief Source terms integrated by this component. + * + * If left empty, the ConstantSourceTerm components found in the current context are used. + */ + sofa::MultiLink, ConstantSourceTerm, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_constantSources; + /** * @brief Initializes the component. * * This method performs several initialization steps: * 1. Initializes the base force field. * 2. Initializes the topology accessor. - * 3. Resizes the nodal source density. + * 3. Validates the linked source terms. * 4. Assembles the global matrix M. + * 5. Integrates the source terms into the constant nodal force. */ void init() override; /** - * @brief Adds the force (f = M * b) to the RHS vector. + * @brief Adds the nodal source term to the RHS vector. * - * This method computes the product of the geometric matrix and the source density vector, - * adding the result to the force vector `f`. + * The source terms were integrated once in init and are only accumulated here. * * @param mparams Mechanical parameters for the computation. * @param f The force vector to which the source term will be added. - * @param x The current positions (unused in this implementation). + * @param x The current positions (unused: the load is prescribed on the rest configuration). * @param v The current velocities (unused in this implementation). */ void addForce( @@ -96,29 +109,24 @@ class FEMSourceTerm : const sofa::DataVecDeriv_t& v) override; /** - * @brief A no-op as the load is prescribed on the rest configuration - */ + * @brief No-op: TODO for non-const source terms + */ void addDForce(const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& df, const sofa::DataVecDeriv_t& dx) override; /** - * @brief A no-op as the load is prescribed on the rest configuration - */ + * @brief No-op: TODO for non-const source terms + */ void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; using sofa::core::behavior::ForceField::getPotentialEnergy; /** - * @brief Not implemented, returns 0. + * @brief Potential energy of the constant nodal load, \f$ V = -\sum_a F_a \cdot (x_a - x_{0,a}) \f$. */ SReal getPotentialEnergy(const sofa::core::MechanicalParams* mparams, const sofa::DataVecCoord_t& x) const override; - /** - * @brief Source term (per unit volume) sampled at each node. - */ - sofa::Data > d_nodalSourceDensity; - /** * @brief Degree of the quadrature rule integrating the element matrix M. */ @@ -129,18 +137,32 @@ class FEMSourceTerm : /** * @brief Default constructor. */ - FEMSourceTerm(); + FEMSourceTermIntegrator(); /** - * @brief Resizes the source density to the size of the mechanical state + * @brief Ensures that valid source terms are linked, falling back to the current context. */ - void resizeNodalSourceDensity(const std::size_t size); + void validateSources(); /** * @brief Assembles and stores the geometry-only matrix \f$ M_{ij} = \int_{\Omega} N_i N_j \, d\Omega \f$ over each element on the rest configuration. */ void assembleGlobalMatrix(); + /** + * @brief Sums every displacement-independent source density and integrates it once into m_constantForce. + * + * Integration is linear, so the sum of the terms integrates to the sum of their contributions: + * a single matrix-vector product covers all of them. + */ + void assembleConstantForce(); + + /** + * @brief Applies the geometry-only matrix M to a nodal source term. + */ + void applyGlobalMatrix(const sofa::VecDeriv_t& nodalSourceTerm, + sofa::VecDeriv_t& result) const; + /** * @brief Computes the geometry-only matrix of each element. */ @@ -157,18 +179,26 @@ class FEMSourceTerm : * Stored in compressed sparse row format. Assembled once in init on the rest configuration. */ GlobalMatrix m_globalMatrix; + + /** + * @brief Nodal load of every term in l_constantSources, integrated once in init. + * + * @note Their contribution is integrated once, so editing the property of a linked term at run + * time has no effect until the scene is reinitialised. + */ + sofa::VecDeriv_t m_constantForce; }; -#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_CPP) -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTerm; +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; #endif } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl similarity index 65% rename from Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl rename to Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl index 6e029f0bb98..6e69510ed09 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTerm.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl @@ -20,32 +20,35 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once -#include +#include #include namespace sofa::component::solidmechanics::fem::elastic { template -FEMSourceTerm::FEMSourceTerm() - : d_nodalSourceDensity(initData(&d_nodalSourceDensity, "nodalSourceDensity", - "Source term (per unit volume) sampled at each node. Interpolated inside the " - "element with the shape functions and integrated on the reference configuration.")) +FEMSourceTermIntegrator::FEMSourceTermIntegrator() + : l_constantSources(initLink("constantSources", "Source terms of the weak form integrated by " + "this component. If empty, the ones found in the current context are used.")) , d_quadratureDegree(initData(&d_quadratureDegree, static_cast(1), "quadratureDegree", "Degree of the quadrature rule integrating the element matrix M.")) { + // Re-compute global matrix and constant forces in case of quadrature degree change this->addUpdateCallback("reassembleSourceMatrix", {&d_quadratureDegree}, [this](const sofa::core::DataTracker&) { if (!this->isComponentStateInvalid() && this->l_topology && this->mstate) + { assembleGlobalMatrix(); + assembleConstantForce(); + } return this->getComponentState(); }, {}); } template -void FEMSourceTerm::init() +void FEMSourceTermIntegrator::init() { sofa::core::behavior::ForceField::init(); @@ -54,14 +57,15 @@ void FEMSourceTerm::init() sofa::core::behavior::TopologyAccessor::init(); } - if (!this->isComponentStateInvalid() && this->mstate) + if (!this->isComponentStateInvalid()) { - this->resizeNodalSourceDensity(this->mstate->getSize()); + this->validateSources(); } if (!this->isComponentStateInvalid() && this->l_topology && this->mstate) { this->assembleGlobalMatrix(); + this->assembleConstantForce(); } if (!this->isComponentStateInvalid()) @@ -71,18 +75,28 @@ void FEMSourceTerm::init() } template -void FEMSourceTerm::resizeNodalSourceDensity(const std::size_t size) +void FEMSourceTermIntegrator::validateSources() { - sofa::helper::WriteAccessor nodalSourceDensity = sofa::helper::getWriteAccessor(d_nodalSourceDensity); - - if (nodalSourceDensity.size() < size) + // Gather all ConstantSourceTerm components in Context if empty + if (l_constantSources.empty()) { - nodalSourceDensity.resize(size, sofa::Deriv_t{}); + const auto sourcesInContext = this->getContext()->template getObjects >( + sofa::core::objectmodel::BaseContext::Local); + + for (const auto& source : sourcesInContext) + l_constantSources.add(source); + + msg_info_when(!sourcesInContext.empty(), this) << "No source term linked: the " + << sourcesInContext.size() << " one(s) found in the current context are used."; } + + msg_warning_when(l_constantSources.empty(), this) + << "No source term linked, and none found in the current context '" + << this->getContext()->getName() << "'. This component has zero force contribution."; } template -void FEMSourceTerm::assembleGlobalMatrix() +void FEMSourceTermIntegrator::assembleGlobalMatrix() { const auto& elements = FiniteElement::getElementSequence(*this->l_topology); sofa::type::vector elementMatrices; @@ -95,7 +109,7 @@ void FEMSourceTerm::assembleGlobalMatrix() } template -void FEMSourceTerm::calculateElementMatrix( +void FEMSourceTermIntegrator::calculateElementMatrix( const auto& elements, sofa::type::vector& elementMatrices) { const auto restPositionsAccessor = this->mstate->readRestPositions(); @@ -103,7 +117,7 @@ void FEMSourceTerm::calculateElementMatrix( const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue()); - for (std::size_t elementId = 0; elementId < elements.size(); ++elementId) + for (sofa::Index elementId = 0; elementId < elements.size(); ++elementId) { const auto& element = elements[elementId]; auto& elementMatrix = elementMatrices[elementId]; @@ -129,14 +143,14 @@ void FEMSourceTerm::calculateElementMatrix( } template -void FEMSourceTerm::initializeGlobalMatrix( +void FEMSourceTermIntegrator::initializeGlobalMatrix( const auto& elements, const sofa::type::vector& elementMatrices) { m_globalMatrix.clear(); const auto size = this->mstate->getSize(); m_globalMatrix.resize(size, size); - for (std::size_t elementId = 0; elementId < elements.size(); ++elementId) + for (sofa::Index elementId = 0; elementId < elements.size(); ++elementId) { const auto& element = elements[elementId]; const auto& elementMatrix = elementMatrices[elementId]; @@ -154,7 +168,44 @@ void FEMSourceTerm::initializeGlobalMatrix( } template -void FEMSourceTerm::addForce(const sofa::core::MechanicalParams* mparams, +void FEMSourceTermIntegrator::applyGlobalMatrix( + const sofa::VecDeriv_t& nodalSourceTerm, sofa::VecDeriv_t& result) const +{ + // f_i = sum_j M_ij b_j : apply the global matrix to the nodal source term. + for (sofa::Index xi = 0; xi < m_globalMatrix.rowIndex.size(); ++xi) + { + const auto rowId = m_globalMatrix.rowIndex[xi]; + typename GlobalMatrix::Range rowRange(m_globalMatrix.rowBegin[xi], m_globalMatrix.rowBegin[xi + 1]); + for (typename GlobalMatrix::Index xj = rowRange.begin(); xj < rowRange.end(); ++xj) + { + const auto columnId = m_globalMatrix.colsIndex[xj]; + const auto& value = m_globalMatrix.colsValue[xj]; + + result[rowId] += nodalSourceTerm[columnId] * value; + } + } +} + +template +void FEMSourceTermIntegrator::assembleConstantForce() +{ + const auto size = this->mstate->getSize(); + + // Aggregate all contributions to one vector before applying the global matrix + sofa::VecDeriv_t sourceTerms(size, sofa::Deriv_t{}); + + for (const auto& source : l_constantSources) + { + for (sofa::Index i = 0; i < size; ++i) + sourceTerms[i] += source->getNodeProperty(i); + } + + m_constantForce.assign(size, sofa::Deriv_t{}); + applyGlobalMatrix(sourceTerms, m_constantForce); +} + +template +void FEMSourceTermIntegrator::addForce(const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& f, const sofa::DataVecCoord_t& x, const sofa::DataVecDeriv_t& v) @@ -163,26 +214,21 @@ void FEMSourceTerm::addForce(const sofa::core::Mechanica SOFA_UNUSED(x); SOFA_UNUSED(v); - const sofa::helper::ReadAccessor nodalSourceDensity = sofa::helper::getReadAccessor(d_nodalSourceDensity); + if (this->isComponentStateInvalid()) + { + return; + } + auto forceAccessor = sofa::helper::getWriteAccessor(f); - // f_i = sum_j M_ij b_j : apply the global matrix to the nodal source density. - for (std::size_t xi = 0; xi < m_globalMatrix.rowIndex.size(); ++xi) + for (sofa::Index i = 0; i < m_constantForce.size(); ++i) { - const auto rowId = m_globalMatrix.rowIndex[xi]; - typename GlobalMatrix::Range rowRange(m_globalMatrix.rowBegin[xi], m_globalMatrix.rowBegin[xi + 1]); - for (typename GlobalMatrix::Index xj = rowRange.begin(); xj < rowRange.end(); ++xj) - { - const auto columnId = m_globalMatrix.colsIndex[xj]; - const auto& value = m_globalMatrix.colsValue[xj]; - - forceAccessor[rowId] += nodalSourceDensity[columnId] * value; - } + forceAccessor[i] += m_constantForce[i]; } } template -void FEMSourceTerm::addDForce(const sofa::core::MechanicalParams* mparams, +void FEMSourceTermIntegrator::addDForce(const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& df, const sofa::DataVecDeriv_t& dx) { @@ -192,18 +238,31 @@ void FEMSourceTerm::addDForce(const sofa::core::Mechanic } template -void FEMSourceTerm::buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) +void FEMSourceTermIntegrator::buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) { SOFA_UNUSED(matrix); } template -SReal FEMSourceTerm::getPotentialEnergy(const sofa::core::MechanicalParams* mparams, +SReal FEMSourceTermIntegrator::getPotentialEnergy(const sofa::core::MechanicalParams* mparams, const sofa::DataVecCoord_t& x) const { SOFA_UNUSED(mparams); - SOFA_UNUSED(x); - return 0.0; + + if (this->isComponentStateInvalid()) + { + return 0.0; + } + + const sofa::helper::ReadAccessor positionAccessor = sofa::helper::getReadAccessor(x); + const auto restPositionAccessor = this->mstate->readRestPositions(); + + SReal energy = 0.0; + for (sofa::Index i = 0; i < m_constantForce.size(); ++i) + { + energy -= dot(m_constantForce[i], positionAccessor[i] - restPositionAccessor.ref()[i]); + } + return energy; } } // 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 95e84fb5d1b..806da813432 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 @@ -27,8 +27,9 @@ namespace sofa::component::solidmechanics::fem::elastic { extern void registerBeamFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerConstantSourceTerm(sofa::core::ObjectFactory* factory); extern void registerCorotationalFEMForceField(sofa::core::ObjectFactory* factory); -extern void registerFEMSourceTerm(sofa::core::ObjectFactory* factory); +extern void registerFEMSourceTermIntegrator(sofa::core::ObjectFactory* factory); extern void registerFastTetrahedralCorotationalForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedralFEMForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedralFEMForceFieldAndMass(sofa::core::ObjectFactory* factory); @@ -68,8 +69,9 @@ const char* getModuleVersion() void registerObjects(sofa::core::ObjectFactory* factory) { registerBeamFEMForceField(factory); + registerConstantSourceTerm(factory); registerCorotationalFEMForceField(factory); - registerFEMSourceTerm(factory); + registerFEMSourceTermIntegrator(factory); registerFastTetrahedralCorotationalForceField(factory); registerHexahedralFEMForceField(factory); registerHexahedralFEMForceFieldAndMass(factory); diff --git a/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn b/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn new file mode 100644 index 00000000000..172a575b79c --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn @@ -0,0 +1,110 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From b750712907d2f4ca0d2598a45394d5b6fbbe39e3 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 1 Sep 2026 13:58:25 +0200 Subject: [PATCH 11/17] Demonstrate multiple sources --- .../FEM/FEMSourceTermIntegrator.scn | 229 ++++++++++++------ 1 file changed, 155 insertions(+), 74 deletions(-) diff --git a/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn b/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn index 172a575b79c..964c373b4b0 100644 --- a/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn +++ b/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn @@ -30,81 +30,162 @@ - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 233b23e5cdebdfb74a3da6dab843c738bb590ec5 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 1 Sep 2026 14:41:19 +0200 Subject: [PATCH 12/17] Add scene to show gravity replication --- .../SolidMechanics/FEM/ConstantSourceTerm.scn | 73 +++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100644 examples/Component/SolidMechanics/FEM/ConstantSourceTerm.scn diff --git a/examples/Component/SolidMechanics/FEM/ConstantSourceTerm.scn b/examples/Component/SolidMechanics/FEM/ConstantSourceTerm.scn new file mode 100644 index 00000000000..8fe1c5fa6eb --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/ConstantSourceTerm.scn @@ -0,0 +1,73 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 789e9c1ad246d83007736242bef8e451038c01da Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 1 Sep 2026 15:54:04 +0200 Subject: [PATCH 13/17] Transform example to unit test --- .../FEM/Elastic/tests/CMakeLists.txt | 1 + .../tests/FEMSourceTermIntegrator_test.cpp | 149 ++++++++++++++++++ .../FEM/FEMSourceTermIntegrator.scn | 79 ---------- 3 files changed, 150 insertions(+), 79 deletions(-) create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt index d5e6a1bcb25..b8a38e68eb4 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt @@ -6,6 +6,7 @@ set(SOURCE_FILES BaseTetrahedronFEMForceField_test.h BeamFEMForceField_test.cpp FastTetrahedralCorotationalForceField_test.cpp + FEMSourceTermIntegrator_test.cpp HexahedronFEMForceField_test.cpp StrainDisplacement_test.cpp TetrahedralCorotationalFEMForceField_test.cpp diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp new file mode 100644 index 00000000000..99d5598fdd3 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp @@ -0,0 +1,149 @@ +/****************************************************************************** +* 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 +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace sofa +{ +using sofa::simpleapi::createObject; +using sofa::simpleapi::createRootNode; +using DataTypes = defaulttype::Vec2Types; +using Integrator = component::solidmechanics::fem::elastic::FEMSourceTermIntegrator; +using DOF = component::statecontainer::MechanicalObject; +using VecCoord = DataTypes::VecCoord; +using VecDeriv = DataTypes::VecDeriv; + +class FEMSourceTermIntegrator_test : public testing::BaseTest +{ +protected: + simulation::Simulation* m_simulation = nullptr; + simulation::Node::SPtr m_root; + + void doSetUp() override + { + m_simulation = sofa::simulation::getSimulation(); + } + + void doTearDown() override + { + if (m_root != nullptr) + sofa::simulation::node::unload(m_root); + } + + simulation::Node::SPtr makeMesh() + { + this->loadPlugins({"Sofa.Component.StateContainer", + "Sofa.Component.Topology.Container.Constant", "Sofa.Component.SolidMechanics.FEM.Elastic"}); + + auto root = createRootNode(m_simulation, "root"); + createObject(root, "MechanicalObject", {{"template", "Vec2"}, {"position", "0 0 1 0 1 1 0 1"}}); + createObject(root, "MeshTopology", {{"name", "mesh"}, {"triangles", "0 1 2 0 2 3"}}); + return root; + } + + static VecDeriv addForce(Integrator* integrator) + { + core::MechanicalParams mparams; + Data f; + f.setValue(VecDeriv(4)); + Data x; + Data v; + integrator->addForce(&mparams, f, x, v); + return f.getValue(); + } +}; + +// Splitting one ConstantSourceTerm into several must not change the integrated force. +TEST_F(FEMSourceTermIntegrator_test, MultipleSourcesSumToOne) +{ + m_root = makeMesh(); + + createObject(m_root, "ConstantSourceTerm", {{"name", "full"}, {"template", "Vec2"}, {"property", "300 -600"}}); + auto* one = dynamic_cast(createObject(m_root, "FEMSourceTermIntegrator", + {{"name", "one"}, {"template", "Vec2,Triangle"}, {"topology", "@mesh"}, {"constantSources", "@full"}}).get()); + + createObject(m_root, "ConstantSourceTerm", {{"name", "a"}, {"template", "Vec2"}, {"property", "100 -200"}}); + createObject(m_root, "ConstantSourceTerm", {{"name", "b"}, {"template", "Vec2"}, {"property", "100 -200"}}); + createObject(m_root, "ConstantSourceTerm", {{"name", "c"}, {"template", "Vec2"}, {"property", "100 -200"}}); + auto* three = dynamic_cast(createObject(m_root, "FEMSourceTermIntegrator", + {{"name", "three"}, {"template", "Vec2,Triangle"}, {"topology", "@mesh"}, {"constantSources", "@a @b @c"}}).get()); + + simulation::node::initRoot(m_root.get()); + ASSERT_NE(one, nullptr); + ASSERT_NE(three, nullptr); + + const VecDeriv fOne = addForce(one); + const VecDeriv fThree = addForce(three); + ASSERT_EQ(fOne.size(), fThree.size()); + for (std::size_t i = 0; i < fOne.size(); ++i) + EXPECT_EQ(fOne[i], fThree[i]); +} + +// dE = -dX . F, the conservative-force identity getPotentialEnergy relies on. +TEST_F(FEMSourceTermIntegrator_test, PotentialEnergyMatchesWork) +{ + m_root = makeMesh(); + + createObject(m_root, "ConstantSourceTerm", {{"name", "bodyForce"}, {"template", "Vec2"}, {"property", "300 -600"}}); + auto* integrator = dynamic_cast(createObject(m_root, "FEMSourceTermIntegrator", + {{"name", "source"}, {"template", "Vec2,Triangle"}, {"topology", "@mesh"}, {"constantSources", "@bodyForce"}}).get()); + + simulation::node::initRoot(m_root.get()); + ASSERT_NE(integrator, nullptr); + + const VecDeriv f = addForce(integrator); + + VecCoord x0(4); + testing::copyFromData(x0, m_root->get()->readPositions()); + + VecCoord x1 = x0; + for (auto& xi : x1) + xi += DataTypes::Coord(0.01, -0.02); + + core::MechanicalParams mparams; + Data x0Data; + x0Data.setValue(x0); + Data x1Data; + x1Data.setValue(x1); + + const SReal e0 = integrator->getPotentialEnergy(&mparams, x0Data); + const SReal e1 = integrator->getPotentialEnergy(&mparams, x1Data); + + SReal work = 0; + for (std::size_t i = 0; i < f.size(); ++i) + work += dot(f[i], x1[i] - x0[i]); + + EXPECT_NEAR(e1 - e0, -work, 1e-9); +} + +} // namespace sofa diff --git a/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn b/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn index 964c373b4b0..f473fd8b62e 100644 --- a/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn +++ b/examples/Component/SolidMechanics/FEM/FEMSourceTermIntegrator.scn @@ -109,83 +109,4 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - From bfe4b7e8e1b143b3c366b6688a09abdfdf832af4 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 1 Sep 2026 16:29:48 +0200 Subject: [PATCH 14/17] [FEM] Add NonConstantSourceTerm and enable it in the FEMSourceTermIntegrator --- .../SolidMechanics/FEM/Elastic/CMakeLists.txt | 4 + .../fem/elastic/FEMSourceTermIntegrator.h | 77 +++++-- .../fem/elastic/FEMSourceTermIntegrator.inl | 203 +++++++++++++++--- .../fem/elastic/NonConstantSourceTerm.cpp | 63 ++++++ .../fem/elastic/NonConstantSourceTerm.h | 117 ++++++++++ .../fem/elastic/TractionSourceTerm.cpp | 45 ++++ .../fem/elastic/TractionSourceTerm.h | 114 ++++++++++ .../solidmechanics/fem/elastic/init.cpp | 4 + .../FEM/src/sofa/fem/FiniteElement.h | 2 +- .../Component/SolidMechanics/FEM/Traction.scn | 122 +++++++++++ 10 files changed, 706 insertions(+), 45 deletions(-) create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/NonConstantSourceTerm.cpp create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/NonConstantSourceTerm.h create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TractionSourceTerm.cpp create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TractionSourceTerm.h create mode 100644 examples/Component/SolidMechanics/FEM/Traction.scn diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index 4f1356b6330..eeda2d5edeb 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt @@ -14,6 +14,7 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ConstantSourceTerm.h + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/NonConstantSourceTerm.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.h @@ -38,6 +39,7 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TetrahedralCorotationalFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TetrahedronFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TetrahedronFEMForceField.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TractionSourceTerm.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangleFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangleFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangleFEMUtils.h @@ -72,6 +74,7 @@ set(SOURCE_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ConstantSourceTerm.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/NonConstantSourceTerm.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.cpp @@ -84,6 +87,7 @@ set(SOURCE_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/QuadBendingFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TetrahedralCorotationalFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TetrahedronFEMForceField.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TractionSourceTerm.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangleFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangleFEMUtils.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularAnisotropicFEMForceField.cpp diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h index 7873aec5cdf..aabdc15d7eb 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -40,10 +41,11 @@ namespace sofa::component::solidmechanics::fem::elastic * @class FEMSourceTermIntegrator * @brief Integrates source terms into consistent nodal loads. * - * A source term contributes \f$ \int_{\Omega} N_a \, r \, d\Omega \f$ to the right-hand side, where - * r is the per-node density carried by a linked ConstantSourceTerm (through l_constantSources) and - * does not depend on the displacement. Every term is thus summed and integrated once in init(); - * addForce merely accumulates the result. + * A source term contributes \f$ \int_{\Omega} N_a \, r \, d\Omega \f$ to the right-hand side. r is + * the per-node density carried by a linked ConstantSourceTerm (through l_constantSources) or + * NonConstantSourceTerm (through l_nonConstantSources). Constant terms do not depend on the + * displacement, so they are summed and integrated once in init(); non-constant ones are + * re-integrated at every call. * * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). * @tparam TElementType The type of finite element (e.g., sofa::geometry::Tetrahedron). @@ -63,12 +65,23 @@ class FEMSourceTermIntegrator : protected: using FiniteElement = sofa::fem::FiniteElement; using Real = sofa::Real_t; + using Coord = sofa::Coord_t; + using Deriv = sofa::Deriv_t; + using VecCoord = sofa::VecCoord_t; + using VecDeriv = sofa::VecDeriv_t; + using DataVecCoord = sofa::DataVecCoord_t; + using DataVecDeriv = sofa::DataVecDeriv_t; static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; using ElementMatrix = sofa::type::Mat>; using GlobalMatrix = sofa::linearalgebra::CompressedRowSparseMatrixMechanical>; + using Element = typename FiniteElement::TopologyElement; + using ShapeFunctions = sofa::type::Vec; + using Jacobian = sofa::type::Mat; + using SourceDerivative = sofa::type::Mat; public: @@ -80,6 +93,14 @@ class FEMSourceTermIntegrator : sofa::MultiLink, ConstantSourceTerm, sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_constantSources; + /** + * @brief Displacement-dependent source terms linked to this component. + * + * If left empty, the NonConstantSourceTerm components found in the current context are used. + */ + sofa::MultiLink, NonConstantSourceTerm, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_nonConstantSources; + /** * @brief Initializes the component. * @@ -95,28 +116,33 @@ class FEMSourceTermIntegrator : /** * @brief Adds the nodal source term to the RHS vector. * - * The source terms were integrated once in init and are only accumulated here. + * The constant terms were integrated once in init and are only accumulated here. The + * displacement-dependent ones are integrated at the current position. * * @param mparams Mechanical parameters for the computation. * @param f The force vector to which the source term will be added. - * @param x The current positions (unused: the load is prescribed on the rest configuration). + * @param x The current positions, used only by displacement-dependent terms. * @param v The current velocities (unused in this implementation). */ void addForce( const sofa::core::MechanicalParams* mparams, - sofa::DataVecDeriv_t& f, - const sofa::DataVecCoord_t& x, - const sofa::DataVecDeriv_t& v) override; + DataVecDeriv& f, + const DataVecCoord& x, + const DataVecDeriv& v) override; /** - * @brief No-op: TODO for non-const source terms + * @brief Applies the tangent of the displacement-dependent terms to dx. + * + * A no-op when l_nonConstantSources is empty. */ void addDForce(const sofa::core::MechanicalParams* mparams, - sofa::DataVecDeriv_t& df, - const sofa::DataVecDeriv_t& dx) override; + DataVecDeriv& df, + const DataVecDeriv& dx) override; /** - * @brief No-op: TODO for non-const source terms + * @brief Assembles the stiffness contribution of the displacement-dependent terms. + * + * A no-op when l_nonConstantSources is empty. */ void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; @@ -125,13 +151,18 @@ class FEMSourceTermIntegrator : * @brief Potential energy of the constant nodal load, \f$ V = -\sum_a F_a \cdot (x_a - x_{0,a}) \f$. */ SReal getPotentialEnergy(const sofa::core::MechanicalParams* mparams, - const sofa::DataVecCoord_t& x) const override; + const DataVecCoord& x) const override; /** * @brief Degree of the quadrature rule integrating the element matrix M. */ sofa::Data d_quadratureDegree; + /** + * @brief Whether to assemble/apply the stiffness of the displacement-dependent terms. + */ + sofa::Data d_useTangentStiffness; + protected: /** @@ -141,6 +172,8 @@ class FEMSourceTermIntegrator : /** * @brief Ensures that valid source terms are linked, falling back to the current context. + * + * Applies to both l_constantSources and l_nonConstantSources. */ void validateSources(); @@ -160,11 +193,12 @@ class FEMSourceTermIntegrator : /** * @brief Applies the geometry-only matrix M to a nodal source term. */ - void applyGlobalMatrix(const sofa::VecDeriv_t& nodalSourceTerm, - sofa::VecDeriv_t& result) const; + void applyGlobalMatrix(const VecDeriv& nodalSourceTerm, + VecDeriv& result) const; /** - * @brief Computes the geometry-only matrix of each element. + * @brief Computes the geometry-only matrix of each element, caching the Jacobian of the + * reference-to-physical mapping of every quadrature point along the way in m_referenceJacobian. */ void calculateElementMatrix(const auto& elements, sofa::type::vector& elementMatrices); @@ -186,7 +220,14 @@ class FEMSourceTermIntegrator : * @note Their contribution is integrated once, so editing the property of a linked term at run * time has no effect until the scene is reinitialised. */ - sofa::VecDeriv_t m_constantForce; + VecDeriv m_constantForce; + + /** + * @brief Jacobian of the reference-to-physical mapping, evaluated on the rest configuration. + * + * Assembled once in init on the rest configuration. + */ + sofa::type::vector m_referenceJacobian; }; #if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl index 6e69510ed09..900f85efcdb 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl @@ -22,6 +22,8 @@ #pragma once #include #include +#include +#include namespace sofa::component::solidmechanics::fem::elastic { @@ -30,8 +32,13 @@ template FEMSourceTermIntegrator::FEMSourceTermIntegrator() : l_constantSources(initLink("constantSources", "Source terms of the weak form integrated by " "this component. If empty, the ones found in the current context are used.")) + , l_nonConstantSources(initLink("nonConstantSources", "Displacement-dependent source terms " + "linked to this component. If empty, the ones found in the current context are " + "used.")) , d_quadratureDegree(initData(&d_quadratureDegree, static_cast(1), "quadratureDegree", "Degree of the quadrature rule integrating the element matrix M.")) + , d_useTangentStiffness(initData(&d_useTangentStiffness, true, "useTangentStiffness", + "Whether to assemble/apply the stiffness of the displacement-dependent terms.")) { // Re-compute global matrix and constant forces in case of quadrature degree change this->addUpdateCallback("reassembleSourceMatrix", {&d_quadratureDegree}, @@ -77,20 +84,28 @@ void FEMSourceTermIntegrator::init() template void FEMSourceTermIntegrator::validateSources() { - // Gather all ConstantSourceTerm components in Context if empty - if (l_constantSources.empty()) + // Gather all matching source components in Context if a link is left empty + auto fallbackToContext = [this](auto& link) { - const auto sourcesInContext = this->getContext()->template getObjects >( - sofa::core::objectmodel::BaseContext::Local); + using SourceType = typename std::remove_reference_t::DestType; - for (const auto& source : sourcesInContext) - l_constantSources.add(source); + if (link.empty()) + { + const auto sourcesInContext = this->getContext()->template getObjects( + sofa::core::objectmodel::BaseContext::Local); - msg_info_when(!sourcesInContext.empty(), this) << "No source term linked: the " - << sourcesInContext.size() << " one(s) found in the current context are used."; - } + for (const auto& source : sourcesInContext) + link.add(source); + + msg_info_when(!sourcesInContext.empty(), this) << "No source term linked: the " + << sourcesInContext.size() << " one(s) found in the current context are used."; + } + }; + + fallbackToContext(l_constantSources); + fallbackToContext(l_nonConstantSources); - msg_warning_when(l_constantSources.empty(), this) + msg_warning_when(l_constantSources.empty() && l_nonConstantSources.empty(), this) << "No source term linked, and none found in the current context '" << this->getContext()->getName() << "'. This component has zero force contribution."; } @@ -116,28 +131,35 @@ void FEMSourceTermIntegrator::calculateElementMatrix( elementMatrices.resize(elements.size()); const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue()); + const auto quadraturePointsPerElement = quadratureRule.size(); + + m_referenceJacobian.resize(elements.size() * quadraturePointsPerElement); for (sofa::Index elementId = 0; elementId < elements.size(); ++elementId) { const auto& element = elements[elementId]; auto& elementMatrix = elementMatrices[elementId]; - const std::array, NumberOfNodesInElement> elementNodesRestCoordinates = + const std::array elementNodesRestCoordinates = extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); // M_ij = integral of N_i N_j dV, evaluated on the rest configuration (geometry only). + sofa::Index quadraturePointIndex = 0; for (const auto& [quadraturePoint, weight] : quadratureRule) { const auto N = FiniteElement::shapeFunctions(quadraturePoint); const auto dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); - const auto jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + auto& jacobian = m_referenceJacobian[elementId * quadraturePointsPerElement + quadraturePointIndex]; + jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( elementNodesRestCoordinates, dN_dq_ref); const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); const auto NT_N = sofa::type::dyad(N, N); elementMatrix += (weight * detJ) * NT_N; + + ++quadraturePointIndex; } } } @@ -169,7 +191,7 @@ void FEMSourceTermIntegrator::initializeGlobalMatrix( template void FEMSourceTermIntegrator::applyGlobalMatrix( - const sofa::VecDeriv_t& nodalSourceTerm, sofa::VecDeriv_t& result) const + const VecDeriv& nodalSourceTerm, VecDeriv& result) const { // f_i = sum_j M_ij b_j : apply the global matrix to the nodal source term. for (sofa::Index xi = 0; xi < m_globalMatrix.rowIndex.size(); ++xi) @@ -192,7 +214,7 @@ void FEMSourceTermIntegrator::assembleConstantForce() const auto size = this->mstate->getSize(); // Aggregate all contributions to one vector before applying the global matrix - sofa::VecDeriv_t sourceTerms(size, sofa::Deriv_t{}); + VecDeriv sourceTerms(size, Deriv{}); for (const auto& source : l_constantSources) { @@ -200,18 +222,17 @@ void FEMSourceTermIntegrator::assembleConstantForce() sourceTerms[i] += source->getNodeProperty(i); } - m_constantForce.assign(size, sofa::Deriv_t{}); + m_constantForce.assign(size, Deriv{}); applyGlobalMatrix(sourceTerms, m_constantForce); } template void FEMSourceTermIntegrator::addForce(const sofa::core::MechanicalParams* mparams, - sofa::DataVecDeriv_t& f, - const sofa::DataVecCoord_t& x, - const sofa::DataVecDeriv_t& v) + DataVecDeriv& f, + const DataVecCoord& x, + const DataVecDeriv& v) { SOFA_UNUSED(mparams); - SOFA_UNUSED(x); SOFA_UNUSED(v); if (this->isComponentStateInvalid()) @@ -225,27 +246,157 @@ void FEMSourceTermIntegrator::addForce(const sofa::core: { forceAccessor[i] += m_constantForce[i]; } + + if (!l_nonConstantSources.empty()) + { + const sofa::helper::ReadAccessor positionAccessor = sofa::helper::getReadAccessor(x); + VecDeriv& nonConstantForce = forceAccessor.wref(); + + const auto restPositionsAccessor = this->mstate->readRestPositions(); + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue()); + + for (const auto& element : elements) + { + const std::array elementNodesRestCoordinates = + extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); + const std::array elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, positionAccessor.ref()); + + std::array elementNodesDisplacement; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + elementNodesDisplacement[i] = elementNodesCoordinates[i] - elementNodesRestCoordinates[i]; + } + + for (const auto& [quadraturePoint, weight] : quadratureRule) + { + const auto N = FiniteElement::shapeFunctions(quadraturePoint); + const auto dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + + const auto jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + elementNodesCoordinates, dN_dq_ref); + const auto weightTimesDetJ = static_cast(weight * sofa::type::absGeneralizedDeterminant(jacobian)); + + const auto restPosition = FiniteElement::Helper::evaluateValueInElement(elementNodesRestCoordinates, N); + const auto displacement = FiniteElement::Helper::evaluateValueInElement(elementNodesDisplacement, N); + + for (const auto& source : l_nonConstantSources) + { + const auto sourceDensity = source->evaluate(restPosition, displacement, jacobian); + + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + nonConstantForce[element[i]] += sourceDensity * (weightTimesDetJ * N[i]); + } + } + } + } + } } template void FEMSourceTermIntegrator::addDForce(const sofa::core::MechanicalParams* mparams, - sofa::DataVecDeriv_t& df, - const sofa::DataVecDeriv_t& dx) + DataVecDeriv& df, + const DataVecDeriv& dx) { - SOFA_UNUSED(mparams); - SOFA_UNUSED(df); - SOFA_UNUSED(dx); + if (this->isComponentStateInvalid() || l_nonConstantSources.empty() || !d_useTangentStiffness.getValue()) + { + return; + } + + // never mparams->kFactor() directly, so that Rayleigh stiffness damping is folded in + const auto kFactor = static_cast(sofa::core::mechanicalparams::kFactorIncludingRayleighDamping( + mparams, this->rayleighStiffness.getValue())); + + auto forceDerivAccessor = sofa::helper::getWriteAccessor(df); + const sofa::helper::ReadAccessor positionDerivAccessor = sofa::helper::getReadAccessor(dx); + const auto positionsAccessor = this->mstate->readPositions(); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue()); + + for (const auto& element : elements) + { + const std::array elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, positionsAccessor.ref()); + + for (const auto& [quadraturePoint, weight] : quadratureRule) + { + const auto N = FiniteElement::shapeFunctions(quadraturePoint); + const auto dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + + const auto jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + elementNodesCoordinates, dN_dq_ref); + + Deriv contraction{}; + for (const auto& source : l_nonConstantSources) + { + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + contraction += source->evaluateStiffness(jacobian, dN_dq_ref[j]) * positionDerivAccessor[element[j]]; + } + } + + for (sofa::Size a = 0; a < NumberOfNodesInElement; ++a) + { + forceDerivAccessor[element[a]] += contraction * (kFactor * static_cast(weight) * N[a]); + } + } + } } template void FEMSourceTermIntegrator::buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) { - SOFA_UNUSED(matrix); + if (this->isComponentStateInvalid() || l_nonConstantSources.empty() || !d_useTangentStiffness.getValue()) + { + return; + } + + auto dfdx = matrix->getForceDerivativeIn(this->mstate).withRespectToPositionsIn(this->mstate); + + const auto positionsAccessor = this->mstate->readPositions(); + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue()); + + for (const auto& element : elements) + { + const std::array elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, positionsAccessor.ref()); + + for (const auto& [quadraturePoint, weight] : quadratureRule) + { + const auto N = FiniteElement::shapeFunctions(quadraturePoint); + const auto dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + + const auto jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + elementNodesCoordinates, dN_dq_ref); + + std::array D{}; + for (const auto& source : l_nonConstantSources) + { + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + D[j] += source->evaluateStiffness(jacobian, dN_dq_ref[j]); + } + } + + for (sofa::Size a = 0; a < NumberOfNodesInElement; ++a) + { + for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) + { + dfdx(element[a] * spatial_dimensions, element[j] * spatial_dimensions) + += (static_cast(weight) * N[a]) * D[j]; + } + } + } + } } template SReal FEMSourceTermIntegrator::getPotentialEnergy(const sofa::core::MechanicalParams* mparams, - const sofa::DataVecCoord_t& x) const + const DataVecCoord& x) const { SOFA_UNUSED(mparams); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/NonConstantSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/NonConstantSourceTerm.cpp new file mode 100644 index 00000000000..18c1f475685 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/NonConstantSourceTerm.cpp @@ -0,0 +1,63 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#define SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_NON_CONSTANT_SOURCE_TERM_CPP + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerNonConstantSourceTerm(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Source term (per unit volume) depending on the current displacement; zero unless subclassed") + .add< NonConstantSourceTerm >() + .add< NonConstantSourceTerm >() + .add< NonConstantSourceTerm >() + .add< NonConstantSourceTerm >() + .add< NonConstantSourceTerm >() + .add< NonConstantSourceTerm >() + .add< NonConstantSourceTerm >() + .add< NonConstantSourceTerm >() + .add< NonConstantSourceTerm >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/NonConstantSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/NonConstantSourceTerm.h new file mode 100644 index 00000000000..1b6317326be --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/NonConstantSourceTerm.h @@ -0,0 +1,117 @@ +/****************************************************************************** +* 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 +#include +#include +#include + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_NON_CONSTANT_SOURCE_TERM_CPP) +#include +#include +#include +#include +#include +#include +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class NonConstantSourceTerm + * @brief A source density prescribed at the nodes, depending on the current displacement. + * + * evaluate() defaults to zero. This class contributes nothing until subclassed. Link it (or a + * subclass) to a FEMSourceTermIntegrator through l_nonConstantSources. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + */ +template +class NonConstantSourceTerm : public sofa::core::BaseNodalProperty> +{ +public: + using DataTypes = TDataTypes; + using ElementType = TElementType; + + SOFA_CLASS(SOFA_TEMPLATE2(NonConstantSourceTerm, DataTypes, ElementType), + SOFA_TEMPLATE(sofa::core::BaseNodalProperty, sofa::Deriv_t)); + + using Real = sofa::Real_t; + using Coord = sofa::Coord_t; + using Deriv = sofa::Deriv_t; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + using FiniteElement = sofa::fem::FiniteElement; + static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; + + /// Jacobian of the reference-to-physical mapping, evaluated where evaluate() is called. + using Jacobian = sofa::type::Mat; + + /// d(nodal force)/d(node j position), for one test node, at one integration point. + using SourceDerivative = sofa::type::Mat; + + /** + * @brief Source density at one integration point. Defaults to zero. + */ + virtual Deriv evaluate(const Coord& restPosition, const Deriv& displacement, const Jacobian& jacobian) const + { + SOFA_UNUSED(restPosition); + SOFA_UNUSED(displacement); + SOFA_UNUSED(jacobian); + return Deriv{}; + } + + /** + * @brief Stiffness contribution of node j at one integration point + * Defaults to zero. + */ + virtual SourceDerivative evaluateStiffness(const Jacobian& jacobian, + const sofa::type::Vec& gradientOfShapeFunction) const + { + SOFA_UNUSED(jacobian); + SOFA_UNUSED(gradientOfShapeFunction); + return SourceDerivative{}; + } + +protected: + + NonConstantSourceTerm() : sofa::core::BaseNodalProperty(Deriv{}) {} +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_NON_CONSTANT_SOURCE_TERM_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NonConstantSourceTerm; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TractionSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TractionSourceTerm.cpp new file mode 100644 index 00000000000..6fd2595ebd0 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TractionSourceTerm.cpp @@ -0,0 +1,45 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_TRACTION_SOURCE_TERM_CPP + +#include + +#include +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerTractionSourceTerm(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Pressure load following the current-configuration normal direction, on a Triangle or Quad boundary") + .add< TractionSourceTerm >() + .add< TractionSourceTerm >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API TractionSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API TractionSourceTerm; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TractionSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TractionSourceTerm.h new file mode 100644 index 00000000000..1d649530845 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TractionSourceTerm.h @@ -0,0 +1,114 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_TRACTION_SOURCE_TERM_CPP) +#include +#include +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class TractionSourceTerm + * @brief A pressure load, following the current-configuration normal direction. + * + * Restricted to Triangle and Quad, the same element types SurfacePressureForceField supports. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + * @tparam TElementType The boundary element type (Triangle or Quad). + */ +template +class TractionSourceTerm : public NonConstantSourceTerm +{ +public: + using DataTypes = TDataTypes; + using ElementType = TElementType; + + SOFA_CLASS(SOFA_TEMPLATE2(TractionSourceTerm, DataTypes, ElementType), + SOFA_TEMPLATE2(NonConstantSourceTerm, DataTypes, ElementType)); + + using Real = sofa::Real_t; + using Deriv = sofa::Deriv_t; + using Coord = sofa::Coord_t; + using Jacobian = typename NonConstantSourceTerm::Jacobian; + using SourceDerivative = typename NonConstantSourceTerm::SourceDerivative; + static constexpr sofa::Size TopologicalDimension = NonConstantSourceTerm::TopologicalDimension; + + /** + * @brief Pressure per unit area, following the current-configuration normal direction. + */ + sofa::Data d_pressure; + + /** + * @brief pressure * unit normal, jacobian.col(0) x jacobian.col(1) normalized — same tangent + * convention as SurfacePressureForceField, generalized through the Jacobian so Triangle and + * Quad share one formula. + */ + Deriv evaluate(const Coord& restPosition, const Deriv& displacement, const Jacobian& jacobian) const override + { + SOFA_UNUSED(restPosition); + SOFA_UNUSED(displacement); + return jacobian.col(0).cross(jacobian.col(1)).normalized() * d_pressure.getValue(); + } + + /** + * @brief d(pressure * jacobian.col(0) x jacobian.col(1))/d(node j position). + */ + SourceDerivative evaluateStiffness(const Jacobian& jacobian, + const sofa::type::Vec& gradientOfShapeFunction) const override + { + const auto t0 = jacobian.col(0); + const auto t1 = jacobian.col(1); + + SourceDerivative skewT0{}; + skewT0(0,1) = -t0[2]; skewT0(0,2) = t0[1]; + skewT0(1,0) = t0[2]; skewT0(1,2) = -t0[0]; + skewT0(2,0) = -t0[1]; skewT0(2,1) = t0[0]; + + SourceDerivative skewT1{}; + skewT1(0,1) = -t1[2]; skewT1(0,2) = t1[1]; + skewT1(1,0) = t1[2]; skewT1(1,2) = -t1[0]; + skewT1(2,0) = -t1[1]; skewT1(2,1) = t1[0]; + + return (skewT0 * gradientOfShapeFunction[1] - skewT1 * gradientOfShapeFunction[0]) * d_pressure.getValue(); + } + +protected: + + TractionSourceTerm() + : d_pressure(initData(&d_pressure, Real{0}, "pressure", + "Pressure per unit area, following the current-configuration normal direction.")) + {} +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_TRACTION_SOURCE_TERM_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API TractionSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API TractionSourceTerm; +#endif + +} // 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 806da813432..27449e0bf2c 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 @@ -36,9 +36,11 @@ extern void registerHexahedralFEMForceFieldAndMass(sofa::core::ObjectFactory* fa extern void registerHexahedronFEMForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedronFEMForceFieldAndMass(sofa::core::ObjectFactory* factory); extern void registerLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerNonConstantSourceTerm(sofa::core::ObjectFactory* factory); extern void registerQuadBendingFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTetrahedralCorotationalFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTetrahedronFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerTractionSourceTerm(sofa::core::ObjectFactory* factory); extern void registerTriangleFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangularAnisotropicFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangularFEMForceField(sofa::core::ObjectFactory* factory); @@ -78,9 +80,11 @@ void registerObjects(sofa::core::ObjectFactory* factory) registerHexahedronFEMForceField(factory); registerHexahedronFEMForceFieldAndMass(factory); registerLinearSmallStrainFEMForceField(factory); + registerNonConstantSourceTerm(factory); registerQuadBendingFEMForceField(factory); registerTetrahedralCorotationalFEMForceField(factory); registerTetrahedronFEMForceField(factory); + registerTractionSourceTerm(factory); registerTriangleFEMForceField(factory); registerTriangularAnisotropicFEMForceField(factory); registerTriangularFEMForceField(factory); diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h index a490fddd6a3..9c389914f9b 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h @@ -86,7 +86,7 @@ struct FiniteElementHelper const std::array& valuesAtNodes, const sofa::type::Vec& shapeFunctions) { - return std::inner_product(valuesAtNodes.begin(), valuesAtNodes.end(), shapeFunctions.begin(), T(0)); + return std::inner_product(valuesAtNodes.begin(), valuesAtNodes.end(), shapeFunctions.begin(), T{}); } }; diff --git a/examples/Component/SolidMechanics/FEM/Traction.scn b/examples/Component/SolidMechanics/FEM/Traction.scn new file mode 100644 index 00000000000..f2fcdf05795 --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/Traction.scn @@ -0,0 +1,122 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 2b3d3cb27daae44497c40360bdc519d7b3e98433 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Wed, 2 Sep 2026 18:33:49 +0200 Subject: [PATCH 15/17] Bypass assembly when no constant sources are present --- .../fem/elastic/FEMSourceTermIntegrator.h | 10 +--------- .../fem/elastic/FEMSourceTermIntegrator.inl | 13 +++---------- 2 files changed, 4 insertions(+), 19 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h index aabdc15d7eb..2cf86b8e808 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h @@ -197,8 +197,7 @@ class FEMSourceTermIntegrator : VecDeriv& result) const; /** - * @brief Computes the geometry-only matrix of each element, caching the Jacobian of the - * reference-to-physical mapping of every quadrature point along the way in m_referenceJacobian. + * @brief Computes the geometry-only matrix of each element. */ void calculateElementMatrix(const auto& elements, sofa::type::vector& elementMatrices); @@ -221,13 +220,6 @@ class FEMSourceTermIntegrator : * time has no effect until the scene is reinitialised. */ VecDeriv m_constantForce; - - /** - * @brief Jacobian of the reference-to-physical mapping, evaluated on the rest configuration. - * - * Assembled once in init on the rest configuration. - */ - sofa::type::vector m_referenceJacobian; }; #if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl index 900f85efcdb..d15a730aa02 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl @@ -44,7 +44,7 @@ FEMSourceTermIntegrator::FEMSourceTermIntegrator() this->addUpdateCallback("reassembleSourceMatrix", {&d_quadratureDegree}, [this](const sofa::core::DataTracker&) { - if (!this->isComponentStateInvalid() && this->l_topology && this->mstate) + if (!this->isComponentStateInvalid() && this->l_topology && this->mstate && !l_constantSources.empty()) { assembleGlobalMatrix(); assembleConstantForce(); @@ -69,7 +69,7 @@ void FEMSourceTermIntegrator::init() this->validateSources(); } - if (!this->isComponentStateInvalid() && this->l_topology && this->mstate) + if (!this->isComponentStateInvalid() && this->l_topology && this->mstate && !l_constantSources.empty()) { this->assembleGlobalMatrix(); this->assembleConstantForce(); @@ -131,9 +131,6 @@ void FEMSourceTermIntegrator::calculateElementMatrix( elementMatrices.resize(elements.size()); const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue()); - const auto quadraturePointsPerElement = quadratureRule.size(); - - m_referenceJacobian.resize(elements.size() * quadraturePointsPerElement); for (sofa::Index elementId = 0; elementId < elements.size(); ++elementId) { @@ -144,22 +141,18 @@ void FEMSourceTermIntegrator::calculateElementMatrix( extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); // M_ij = integral of N_i N_j dV, evaluated on the rest configuration (geometry only). - sofa::Index quadraturePointIndex = 0; for (const auto& [quadraturePoint, weight] : quadratureRule) { const auto N = FiniteElement::shapeFunctions(quadraturePoint); const auto dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); - auto& jacobian = m_referenceJacobian[elementId * quadraturePointsPerElement + quadraturePointIndex]; - jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + const auto jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( elementNodesRestCoordinates, dN_dq_ref); const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); const auto NT_N = sofa::type::dyad(N, N); elementMatrix += (weight * detJ) * NT_N; - - ++quadraturePointIndex; } } } From 2708665118baae61af0404613802b3042cd137f9 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Wed, 2 Sep 2026 19:39:19 +0200 Subject: [PATCH 16/17] Add tests for addDForce and buildStiffnessMatrix --- .../tests/FEMSourceTermIntegrator_test.cpp | 163 ++++++++++++++++++ 1 file changed, 163 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp index 99d5598fdd3..b038c7bb24d 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp @@ -29,7 +29,9 @@ #include #include +#include #include +#include #include #include @@ -43,6 +45,17 @@ using DOF = component::statecontainer::MechanicalObject; using VecCoord = DataTypes::VecCoord; using VecDeriv = DataTypes::VecDeriv; +// The non-constant (traction) path needs Vec3: the pressure/traction load is expressed +// through the triangle's normal, which only exists in 3D. +using DataTypes3 = defaulttype::Vec3Types; +using Integrator3 = component::solidmechanics::fem::elastic::FEMSourceTermIntegrator; +using Traction = component::solidmechanics::fem::elastic::TractionSourceTerm; +using DOF3 = component::statecontainer::MechanicalObject; +using VecCoord3 = DataTypes3::VecCoord; +using VecDeriv3 = DataTypes3::VecDeriv; +using Coord3 = DataTypes3::Coord; +using Deriv3 = DataTypes3::Deriv; + class FEMSourceTermIntegrator_test : public testing::BaseTest { protected: @@ -81,6 +94,39 @@ class FEMSourceTermIntegrator_test : public testing::BaseTest integrator->addForce(&mparams, f, x, v); return f.getValue(); } + + simulation::Node::SPtr makeTractionMesh(Integrator3*& integrator, Traction*& load) + { + this->loadPlugins({"Sofa.Component.StateContainer", + "Sofa.Component.Topology.Container.Constant", "Sofa.Component.SolidMechanics.FEM.Elastic"}); + + auto root = createRootNode(m_simulation, "root"); + createObject(root, "MechanicalObject", {{"template", "Vec3"}, {"position", "0 0 0 1 0 0 1 1 0 0 1 0"}}); + createObject(root, "MeshTopology", {{"name", "mesh"}, {"triangles", "0 1 2 0 2 3"}}); + + load = dynamic_cast(createObject(root, "TractionSourceTerm", + {{"name", "load"}, {"template", "Vec3,Triangle"}, {"pressure", "1000"}}).get()); + integrator = dynamic_cast(createObject(root, "FEMSourceTermIntegrator", + {{"name", "traction"}, {"template", "Vec3,Triangle"}, {"topology", "@mesh"}, + {"quadratureDegree", "1"}, {"nonConstantSources", "@load"}}).get()); + + return root; + } + + // Minimal dense accumulator: enough for the tiny meshes these tests use. + struct DenseStiffnessAccumulator : public core::behavior::StiffnessMatrixAccumulator + { + explicit DenseStiffnessAccumulator(sofa::Size size) : K(size, sofa::type::vector(size, 0.0)) {} + + void add(sofa::SignedIndex row, sofa::SignedIndex col, const sofa::type::Mat<3, 3, double>& value) override + { + for (sofa::Size i = 0; i < 3; ++i) + for (sofa::Size j = 0; j < 3; ++j) + K[row + i][col + j] += value(i, j); + } + + sofa::type::vector> K; + }; }; // Splitting one ConstantSourceTerm into several must not change the integrated force. @@ -146,4 +192,121 @@ TEST_F(FEMSourceTermIntegrator_test, PotentialEnergyMatchesWork) EXPECT_NEAR(e1 - e0, -work, 1e-9); } +// addDForce is the analytic tangent of a non-constant source (TractionSourceTerm). Check it +// against a central finite difference of addForce itself, taken at the same configuration. +TEST_F(FEMSourceTermIntegrator_test, TractionAddDForceMatchesFiniteDifference) +{ + Integrator3* integrator = nullptr; + Traction* load = nullptr; + m_root = makeTractionMesh(integrator, load); + simulation::node::initRoot(m_root.get()); + ASSERT_NE(integrator, nullptr); + ASSERT_NE(load, nullptr); + + const std::size_t n = m_root->get()->getSize(); + VecCoord3 x0(n); + testing::copyFromData(x0, m_root->get()->readPositions()); + + VecDeriv3 dx(n); + dx[0] = Deriv3(0.3, -0.2, 0.1); + dx[1] = Deriv3(-0.1, 0.4, -0.3); + dx[2] = Deriv3(0.2, 0.2, 0.2); + dx[3] = Deriv3(-0.2, 0.1, 0.3); + + const SReal h = 1e-6; + core::MechanicalParams mparams; + mparams.setKFactor(1.0); + + const auto evaluateForceAt = [&](SReal sign) + { + VecCoord3 x(n); + for (std::size_t i = 0; i < x.size(); ++i) + x[i] = x0[i] + dx[i] * (sign * h); + + Data xData; + xData.setValue(x); + Data vData; + Data fData; + fData.setValue(VecDeriv3(n)); + integrator->addForce(&mparams, fData, xData, vData); + return fData.getValue(); + }; + + const VecDeriv3 fPlus = evaluateForceAt(1.0); + const VecDeriv3 fMinus = evaluateForceAt(-1.0); + + VecDeriv3 finiteDifference(n); + for (std::size_t i = 0; i < finiteDifference.size(); ++i) + finiteDifference[i] = (fPlus[i] - fMinus[i]) / (2 * h); + + Data dfData; + dfData.setValue(VecDeriv3(n)); + Data dxData; + dxData.setValue(dx); + integrator->addDForce(&mparams, dfData, dxData); + const VecDeriv3 df = dfData.getValue(); + + for (std::size_t i = 0; i < df.size(); ++i) + for (unsigned d = 0; d < 3; ++d) + EXPECT_NEAR(df[i][d], finiteDifference[i][d], 1e-7) + << "node " << i << ", component " << d; +} + +// buildStiffnessMatrix and addDForce are two independently-written paths over the same +// per-element tangent (NonConstantSourceTerm::evaluateStiffness); they must agree exactly. +TEST_F(FEMSourceTermIntegrator_test, TractionBuildStiffnessMatrixMatchesAddDForce) +{ + Integrator3* integrator = nullptr; + Traction* load = nullptr; + m_root = makeTractionMesh(integrator, load); + simulation::node::initRoot(m_root.get()); + ASSERT_NE(integrator, nullptr); + ASSERT_NE(load, nullptr); + + auto* dof = m_root->get(); + const std::size_t n = dof->getSize(); + + VecDeriv3 dx(n); + dx[0] = Deriv3(0.3, -0.2, 0.1); + dx[1] = Deriv3(-0.1, 0.4, -0.3); + dx[2] = Deriv3(0.2, 0.2, 0.2); + dx[3] = Deriv3(-0.2, 0.1, 0.3); + + core::MechanicalParams mparams; + mparams.setKFactor(1.0); + + Data dfData; + dfData.setValue(VecDeriv3(n)); + Data dxData; + dxData.setValue(dx); + integrator->addDForce(&mparams, dfData, dxData); + const VecDeriv3 dfFromAddDForce = dfData.getValue(); + + DenseStiffnessAccumulator accumulator(n * 3); + core::behavior::StiffnessMatrix stiffnessMatrix; + stiffnessMatrix.setMatrixAccumulator(&accumulator, dof); + stiffnessMatrix.setMechanicalParams(&mparams); + integrator->buildStiffnessMatrix(&stiffnessMatrix); + + VecDeriv3 dfFromK(n); + for (std::size_t a = 0; a < n; ++a) + { + for (std::size_t j = 0; j < n; ++j) + { + for (unsigned row = 0; row < 3; ++row) + { + SReal sum = 0; + for (unsigned col = 0; col < 3; ++col) + sum += accumulator.K[a * 3 + row][j * 3 + col] * dx[j][col]; + dfFromK[a][row] += sum; + } + } + } + + for (std::size_t i = 0; i < dfFromK.size(); ++i) + for (unsigned d = 0; d < 3; ++d) + EXPECT_NEAR(dfFromK[i][d], dfFromAddDForce[i][d], 1e-13) + << "node " << i << ", component " << d; +} + } // namespace sofa From 14c9d850e352cc18c7b704213e9fbd8ace956401 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Wed, 2 Sep 2026 19:39:49 +0200 Subject: [PATCH 17/17] Rename example scene --- .../FEM/{Traction.scn => NonConstantSourceTerm.scn} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/Component/SolidMechanics/FEM/{Traction.scn => NonConstantSourceTerm.scn} (100%) diff --git a/examples/Component/SolidMechanics/FEM/Traction.scn b/examples/Component/SolidMechanics/FEM/NonConstantSourceTerm.scn similarity index 100% rename from examples/Component/SolidMechanics/FEM/Traction.scn rename to examples/Component/SolidMechanics/FEM/NonConstantSourceTerm.scn