From 10fddf98aa97183618b73aaad2b55c9b2df2dc0e Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 17:34:02 +0200 Subject: [PATCH 01/17] Add link for non-constant source terms --- .../fem/elastic/FEMSourceTermIntegrator.h | 40 ++++++++++++++----- 1 file changed, 30 insertions(+), 10 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 6f230e46aa5..2e06c493576 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 @@ -40,8 +40,12 @@ namespace sofa::component::solidmechanics::fem::elastic * @brief Integrates a source density 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 density evaluated by a linked BaseSourceTerm (through l_constantSources) at - * each quadrature point. + * r is the density evaluated by a linked BaseSourceTerm at each quadrature point. + * + * The link a term is attached to decides when it is integrated: + * l_constantSources are integrated once on the rest configuration and cached, while + * l_nonConstantSources are integrated again at every call to addForce, on the configuration it is + * given. The latter are follower loads. * * @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). @@ -68,13 +72,17 @@ class FEMSourceTermIntegrator : public: /** - * @brief Source terms integrated by this component. - * - * If left empty, the BaseSourceTerm components found in the current context are used. + * @brief Source terms integrated once on the rest configuration. */ sofa::MultiLink, BaseSourceTerm, sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_constantSources; + /** + * @brief Source terms integrated again at every step, on the current configuration. + */ + sofa::MultiLink, BaseSourceTerm, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_nonConstantSources; + /** * @brief Initializes the component. * @@ -90,6 +98,9 @@ class FEMSourceTermIntegrator : /** * @brief Adds the nodal source term to the RHS vector. * + * The cached load of the constant sources is added as is, and the non-constant ones are + * integrated on the current configuration x that this method is given. + * * @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. @@ -134,6 +145,9 @@ class FEMSourceTermIntegrator : /** * @brief Ensures that valid source terms are linked, falling back to the current context. + * + * If both constant and non-constant source links are left empty, the BaseSourceTerm components + * found in the current context are used as non-constant sources. */ void validateSources(); @@ -143,14 +157,20 @@ class FEMSourceTermIntegrator : void doUpdateInternal() override; /** - * @brief Runs the quadrature and accumulates every linked source term into m_constantForce. - * - * For each element and each quadrature point, a QuadratureContext is built and handed to every - * source term; the density it returns is weighted by \f$ w \, |\det J| \, N_a \f$ and scattered - * onto the element nodes. + * @brief Integrates the terms of l_constantSources on the rest configuration into + * m_constantForce, as consistent nodal forces. */ void assembleConstantForce(); + /** + * @brief Integrates the terms of l_nonConstantSources as consistent nodal forces. + * + * @param x The configuration the terms are integrated on. + * @param f The force vector the nodal forces are accumulated into. + */ + void assembleNonConstantForce(const sofa::VecCoord_t& x, + sofa::VecDeriv_t& f) const; + /** * @brief Nodal load of every term in l_constantSources. */ From 9970b051a25c8838aa59f3799f6070845e54a54a Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 17:35:35 +0200 Subject: [PATCH 02/17] initLink the non-constant source term --- .../solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) 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 264bea419e6..d283d9ace5f 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 @@ -28,8 +28,11 @@ namespace sofa::component::solidmechanics::fem::elastic 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_constantSources(initLink("constantSources", "Source terms of the weak form integrated once " + "on the rest configuration.")) + , l_nonConstantSources(initLink("nonConstantSources", "Source terms of the weak form " + "integrated at every step on the current configuration. If both source links are " + "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 source terms.")) { From 43bf5ad8d2fc193928140a5b5bf0144dc154160d Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 17:39:49 +0200 Subject: [PATCH 03/17] Add an assembler for non-constant forces --- .../fem/elastic/FEMSourceTermIntegrator.inl | 63 ++++++++++++++++++- 1 file changed, 62 insertions(+), 1 deletion(-) 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 d283d9ace5f..5f08e800495 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 @@ -171,6 +171,65 @@ void FEMSourceTermIntegrator::assembleConstantForce() } } +template +void FEMSourceTermIntegrator::assembleNonConstantForce( + const sofa::VecCoord_t& x, sofa::VecDeriv_t& f) const +{ + if (l_nonConstantSources.empty()) + { + return; + } + + 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, NumberOfNodesInElement> elementNodesRestCoordinates = + extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); + const std::array, NumberOfNodesInElement> elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, x); + + std::array, NumberOfNodesInElement> 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 measure = static_cast(sofa::type::absGeneralizedDeterminant(jacobian)); + + const auto restPosition = + FiniteElement::Helper::evaluateValueInElement(elementNodesRestCoordinates, N); + const auto displacement = + FiniteElement::Helper::evaluateValueInElement(elementNodesDisplacement, N); + + const QuadratureContext context{ + element, N, dN_dq_ref, jacobian, measure, restPosition, displacement}; + + const auto weightTimesMeasure = static_cast(weight) * measure; + + for (const auto& source : l_nonConstantSources) + { + const auto density = source->evaluate(context); + + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + f[element[i]] += density * (weightTimesMeasure * N[i]); + } + } + } + } +} + template void FEMSourceTermIntegrator::addForce(const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& f, @@ -178,7 +237,6 @@ void FEMSourceTermIntegrator::addForce(const sofa::core: const sofa::DataVecDeriv_t& v) { SOFA_UNUSED(mparams); - SOFA_UNUSED(x); SOFA_UNUSED(v); if (this->isComponentStateInvalid()) @@ -192,6 +250,9 @@ void FEMSourceTermIntegrator::addForce(const sofa::core: { forceAccessor[i] += m_constantForce[i]; } + + const sofa::helper::ReadAccessor positionAccessor = sofa::helper::getReadAccessor(x); + assembleNonConstantForce(positionAccessor.ref(), forceAccessor.wref()); } template From 00b7f37ed27f6ddb9d765335ee795f1dc912ebff Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 17:44:41 +0200 Subject: [PATCH 04/17] Gather sources from context to non-constant link if none found provided --- .../fem/elastic/FEMSourceTermIntegrator.inl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) 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 5f08e800495..f4debab6820 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 @@ -97,20 +97,21 @@ void FEMSourceTermIntegrator::doUpdateInternal() template void FEMSourceTermIntegrator::validateSources() { - // Gather all BaseSourceTerm components in Context if empty - if (l_constantSources.empty()) + // Gather all BaseSourceTerm components in Context if both links are empty + if (l_constantSources.empty() && l_nonConstantSources.empty()) { const auto sourcesInContext = this->getContext()->template getObjects >( sofa::core::objectmodel::BaseContext::Local); for (const auto& source : sourcesInContext) - l_constantSources.add(source); + l_nonConstantSources.add(source); msg_info_when(!sourcesInContext.empty(), this) << "No source term linked: the " - << sourcesInContext.size() << " one(s) found in the current context are used."; + << sourcesInContext.size() << " one(s) found in the current context are integrated on " + "the current configuration."; } - 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."; } From 12737136849e9a2b77fb68c000cf742f7b55a0d9 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 17:56:03 +0200 Subject: [PATCH 05/17] Demonstrate PressureSourceTerm --- .../SolidMechanics/FEM/PressureSourceTerm.scn | 184 ++++++++++++++++++ 1 file changed, 184 insertions(+) create mode 100644 examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn diff --git a/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn new file mode 100644 index 00000000000..d12d4e8608f --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn @@ -0,0 +1,184 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From e7bc9f40c8b67214fd005879fe2b9b983806f82e Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:28:20 +0200 Subject: [PATCH 06/17] Add pure virtual evaluateStiffness on the BaseSourceTerm --- .../fem/elastic/BaseSourceTerm.h | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h index 457b3cec534..110db1575cd 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h @@ -144,12 +144,20 @@ class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent SOFA_CLASS(SOFA_TEMPLATE2(BaseSourceTerm, DataTypes, ElementType), sofa::core::objectmodel::BaseComponent); + using Real = sofa::Real_t; using Deriv = sofa::Deriv_t; using QuadratureContext_t = QuadratureContext; + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + /// Derivative of a nodal force with respect to the position of one node of its element. + using SourceDerivative = sofa::type::Mat; + /** * @brief Source density at one quadrature point, per unit physical measure. * + * The integrator weights it by \f$ w \, |\det J| \, N_a \f$. + * * @param context Geometry of the quadrature point. */ virtual Deriv evaluate(const QuadratureContext_t& context) const = 0; @@ -165,6 +173,20 @@ class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent return {}; } + /** + * @brief Derivative of the measure-weighted density with respect to a node position. + * + * \f$ \partial (|\det J| \, r) / \partial x_b \f$, the measure included. The integrator weights + * this one by \f$ w \, N_a \f$ alone. + * + * A term that does not follow the configuration returns zero. + * + * @param context Geometry of the quadrature point. + * @param node Index in the element of the node the derivative is taken with respect to. + */ + virtual SourceDerivative evaluateStiffness(const QuadratureContext& context, + sofa::Size node) const = 0; + protected: BaseSourceTerm() = default; From 1e0cb8ed7fd03c45b3c48a2d5fb1b2d630e0da98 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:30:59 +0200 Subject: [PATCH 07/17] Add stiffness to PressureSourceTerm --- .../fem/elastic/BaseSourceTerm.h | 2 +- .../fem/elastic/PressureSourceTerm.h | 8 ++++ .../fem/elastic/PressureSourceTerm.inl | 37 +++++++++++++++++++ 3 files changed, 46 insertions(+), 1 deletion(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h index 110db1575cd..396adc0df01 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h @@ -184,7 +184,7 @@ class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent * @param context Geometry of the quadrature point. * @param node Index in the element of the node the derivative is taken with respect to. */ - virtual SourceDerivative evaluateStiffness(const QuadratureContext& context, + virtual SourceDerivative evaluateStiffness(const QuadratureContext_t& context, sofa::Size node) const = 0; protected: diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h index 6ba2126bff6..8ebdc3de796 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h @@ -78,8 +78,10 @@ class PressureSourceTerm : public BaseSourceTerm SOFA_CLASS(SOFA_TEMPLATE2(PressureSourceTerm, DataTypes, ElementType), SOFA_TEMPLATE2(BaseSourceTerm, DataTypes, ElementType)); + using Real = sofa::Real_t; using Deriv = sofa::Deriv_t; using QuadratureContext_t = QuadratureContext; + using SourceDerivative = typename BaseSourceTerm::SourceDerivative; using NodalPressure = ::sofa::component::solidmechanics::fem::elastic::NodalPressure; /** @@ -103,6 +105,12 @@ class PressureSourceTerm : public BaseSourceTerm */ sofa::type::vector integrandInputs() const override; + /** + * @brief The pressure times the derivative of the area normal of the element. + */ + SourceDerivative evaluateStiffness(const QuadratureContext_t& context, + sofa::Size node) const override; + protected: PressureSourceTerm(); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl index 4947aa792a8..94aa102dbde 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl @@ -76,4 +76,41 @@ PressureSourceTerm::integrandInputs() const return {&l_pressure->d_property}; } +template +auto PressureSourceTerm::evaluateStiffness( + const QuadratureContext_t& context, sofa::Size node) const -> SourceDerivative +{ + if (!l_pressure) + { + return SourceDerivative{}; + } + + const auto pressure = this->interpolateProperty(*l_pressure, context); + const auto& gradient = context.gradientShapeFunctions[node]; + + if constexpr (QuadratureContext_t::spatial_dimensions == 3) + { + const auto tangent0 = context.jacobian.col(0); + const auto tangent1 = context.jacobian.col(1); + + const SourceDerivative skewTangent0 { + {Real{0}, -tangent0[2], tangent0[1]}, + {tangent0[2], Real{0}, -tangent0[0]}, + {-tangent0[1], tangent0[0], Real{0}}}; + + const SourceDerivative skewTangent1 { + {Real{0}, -tangent1[2], tangent1[1]}, + {tangent1[2], Real{0}, -tangent1[0]}, + {-tangent1[1], tangent1[0], Real{0}}}; + + return (skewTangent0 * gradient[1] - skewTangent1 * gradient[0]) * pressure; + } + else + { + const SourceDerivative rotation {{Real{0}, Real{1}}, {Real{-1}, Real{0}}}; + + return rotation * (gradient[0] * pressure); + } +} + } // namespace sofa::component::solidmechanics::fem::elastic From f46eded32d752f7f6a6e992e9800e10c57ff555b Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:32:04 +0200 Subject: [PATCH 08/17] Move skew helper to Base --- .../fem/elastic/BaseSourceTerm.h | 55 +++++++++++++++++++ .../fem/elastic/PressureSourceTerm.inl | 28 +--------- 2 files changed, 57 insertions(+), 26 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h index 396adc0df01..84608175f8c 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h @@ -125,6 +125,61 @@ sofa::type::Vec elementNormal( return normal; } +/** + * @brief Derivative of the area normal of an element with respect to one of its node positions. + * + * The area normal is \f$ |\det J| \, n \f$: in 3D the cross product of the two columns of the + * jacobian, in 2D the tangent rotated by a quarter turn. Differentiating it with respect to the + * position of the node whose shape function is N gives + * + * \f[ [t_0]_\times \frac{\partial N}{\partial q_1} - [t_1]_\times \frac{\partial N}{\partial q_0} + * \quad \text{in 3D,} \qquad R \, \frac{\partial N}{\partial q_0} \quad \text{in 2D,} \f] + * + * where \f$ t_i \f$ is the i-th column of the jacobian, \f$ [\,]_\times \f$ the cross-product + * matrix and R the quarter turn. + * + * Defined only where the element spans one dimension less than the space it lives in, as + * elementNormal is. + * + * @param jacobian dx/dq of the reference-to-physical mapping at the point of interest. + * @param gradientOfShapeFunction dN/dq of the node the derivative is taken with respect to. + */ +template +sofa::type::Mat elementAreaNormalDerivative( + const sofa::type::Mat& jacobian, + const sofa::type::Vec& gradientOfShapeFunction) +{ + static_assert(TopologicalDimension + 1 == spatial_dimensions, + "An area normal is only defined for an element of codimension 1."); + + using Derivative = sofa::type::Mat; + + if constexpr (spatial_dimensions == 3) + { + const auto tangent0 = jacobian.col(0); + const auto tangent1 = jacobian.col(1); + + const Derivative skewTangent0 { + {Real{0}, -tangent0[2], tangent0[1]}, + {tangent0[2], Real{0}, -tangent0[0]}, + {-tangent0[1], tangent0[0], Real{0}}}; + + const Derivative skewTangent1 { + {Real{0}, -tangent1[2], tangent1[1]}, + {tangent1[2], Real{0}, -tangent1[0]}, + {-tangent1[1], tangent1[0], Real{0}}}; + + return skewTangent0 * gradientOfShapeFunction[1] + - skewTangent1 * gradientOfShapeFunction[0]; + } + else + { + const Derivative rotation {{Real{0}, Real{1}}, {Real{-1}, Real{0}}}; + + return rotation * gradientOfShapeFunction[0]; + } +} + /** * @class BaseSourceTerm * @brief A source density whose value is determined by the geometry, not by the solution. diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl index 94aa102dbde..9caffcbe7de 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl @@ -85,32 +85,8 @@ auto PressureSourceTerm::evaluateStiffness( return SourceDerivative{}; } - const auto pressure = this->interpolateProperty(*l_pressure, context); - const auto& gradient = context.gradientShapeFunctions[node]; - - if constexpr (QuadratureContext_t::spatial_dimensions == 3) - { - const auto tangent0 = context.jacobian.col(0); - const auto tangent1 = context.jacobian.col(1); - - const SourceDerivative skewTangent0 { - {Real{0}, -tangent0[2], tangent0[1]}, - {tangent0[2], Real{0}, -tangent0[0]}, - {-tangent0[1], tangent0[0], Real{0}}}; - - const SourceDerivative skewTangent1 { - {Real{0}, -tangent1[2], tangent1[1]}, - {tangent1[2], Real{0}, -tangent1[0]}, - {-tangent1[1], tangent1[0], Real{0}}}; - - return (skewTangent0 * gradient[1] - skewTangent1 * gradient[0]) * pressure; - } - else - { - const SourceDerivative rotation {{Real{0}, Real{1}}, {Real{-1}, Real{0}}}; - - return rotation * (gradient[0] * pressure); - } + return elementAreaNormalDerivative(context.jacobian, context.gradientShapeFunctions[node]) + * this->interpolateProperty(*l_pressure, context); } } // namespace sofa::component::solidmechanics::fem::elastic From 9ed123abb79b5b7515c3dcfc06c54e235eae1d43 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:36:56 +0200 Subject: [PATCH 09/17] Implement stiffness derivative for stress source term --- .../solidmechanics/fem/elastic/StressSourceTerm.h | 7 +++++++ .../solidmechanics/fem/elastic/StressSourceTerm.inl | 13 +++++++++++++ 2 files changed, 20 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.h index 5ea8c00fb91..28c363f6f41 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.h @@ -120,6 +120,7 @@ class BaseStressSourceTerm : public BaseSourceTerm using Deriv = sofa::Deriv_t; using QuadratureContext_t = QuadratureContext; + using SourceDerivative = typename BaseSourceTerm::SourceDerivative; /** * @brief Nodal values of the stress tensor this term integrates. @@ -142,6 +143,12 @@ class BaseStressSourceTerm : public BaseSourceTerm */ sofa::type::vector integrandInputs() const override; + /** + * @brief The stress tensor applied to the derivative of the area normal of the element. + */ + SourceDerivative evaluateStiffness(const QuadratureContext_t& context, + sofa::Size node) const override; + protected: BaseStressSourceTerm(); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.inl index d6a322d2c9c..aedb23fb3dc 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.inl @@ -77,4 +77,17 @@ BaseStressSourceTerm::integrandInputs() co return {&l_stress->d_property}; } +template +auto BaseStressSourceTerm::evaluateStiffness( + const QuadratureContext_t& context, sofa::Size node) const -> SourceDerivative +{ + if (!l_stress) + { + return SourceDerivative{}; + } + + return this->interpolateProperty(*l_stress, context) + * elementAreaNormalDerivative(context.jacobian, context.gradientShapeFunctions[node]); +} + } // namespace sofa::component::solidmechanics::fem::elastic From 8e763d448ee0d4c78b3e5c8368dc257ae9b851b6 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:41:29 +0200 Subject: [PATCH 10/17] Implement it in VectorSourceTerm --- .../solidmechanics/fem/elastic/VectorSourceTerm.h | 11 +++++++++++ .../fem/elastic/VectorSourceTerm.inl | 15 +++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h index 424b9f7e527..5df045e0d2c 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h @@ -78,6 +78,7 @@ class VectorSourceTerm : public BaseSourceTerm using Deriv = sofa::Deriv_t; using QuadratureContext_t = QuadratureContext; + using SourceDerivative = typename BaseSourceTerm::SourceDerivative; using NodalSourceDensity = ::sofa::component::solidmechanics::fem::elastic::NodalSourceDensity; @@ -102,6 +103,16 @@ class VectorSourceTerm : public BaseSourceTerm */ sofa::type::vector integrandInputs() const override; + /** + * @brief The source density times the derivative of the measure of the element. + * + * \f$ \partial (|\det J| \, r) / \partial x_b + * = |\det J| \; r \otimes \left( (J^{-1})^T \frac{\partial N_b}{\partial q} \right) \f$, + * the inverse being the left pseudo-inverse when J is not square. + */ + SourceDerivative evaluateStiffness(const QuadratureContext_t& context, + sofa::Size node) const override; + protected: VectorSourceTerm(); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.inl index 1cead57afb4..a8b3c2b3593 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.inl @@ -76,4 +76,19 @@ VectorSourceTerm::integrandInputs() const return {&l_sourceDensity->d_property}; } +template +auto VectorSourceTerm::evaluateStiffness( + const QuadratureContext_t& context, sofa::Size node) const -> SourceDerivative +{ + if (!l_sourceDensity) + { + return SourceDerivative{}; + } + + const auto measureDerivative = sofa::type::inverse(context.jacobian).transposed() + * context.gradientShapeFunctions[node] * context.measure; + + return sofa::type::dyad(this->interpolateProperty(*l_sourceDensity, context), measureDerivative); +} + } // namespace sofa::component::solidmechanics::fem::elastic From 746f9a24d7078b2b107436985fe2435a110bfc09 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:43:10 +0200 Subject: [PATCH 11/17] Add buildStiffness matrix to the integrator --- .../fem/elastic/FEMSourceTermIntegrator.h | 9 ++- .../fem/elastic/FEMSourceTermIntegrator.inl | 70 ++++++++++++++++++- 2 files changed, 77 insertions(+), 2 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 2e06c493576..7f58113f315 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 @@ -120,7 +120,9 @@ class FEMSourceTermIntegrator : const sofa::DataVecDeriv_t& dx) override; /** - * @brief No-op. + * @brief Assembles the tangent of the terms of l_nonConstantSources. + * + * The tangent of a term that follows the configuration is not symmetric. */ void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; @@ -136,6 +138,11 @@ class FEMSourceTermIntegrator : */ sofa::Data d_quadratureDegree; + /** + * @brief Whether the non-symmetric tangent of the non-constant source terms is assembled. + */ + sofa::Data d_useTangentStiffness; + protected: /** 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 f4debab6820..d2cb70710c1 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 @@ -35,6 +35,8 @@ FEMSourceTermIntegrator::FEMSourceTermIntegrator() "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 source terms.")) + , d_useTangentStiffness(initData(&d_useTangentStiffness, true, "useTangentStiffness", + "Whether the (non-symmetric) tangent of the non-constant source terms is assembled.")) { this->addUpdateCallback("reassembleConstantForce", {&d_quadratureDegree}, [this](const sofa::core::DataTracker&) @@ -269,7 +271,73 @@ void FEMSourceTermIntegrator::addDForce(const sofa::core 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 restPositionsAccessor = this->mstate->readRestPositions(); + 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, NumberOfNodesInElement> elementNodesRestCoordinates = + extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); + const std::array, NumberOfNodesInElement> elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, positionsAccessor.ref()); + + std::array, NumberOfNodesInElement> 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 measure = static_cast(sofa::type::absGeneralizedDeterminant(jacobian)); + + const auto restPosition = + FiniteElement::Helper::evaluateValueInElement(elementNodesRestCoordinates, N); + const auto displacement = + FiniteElement::Helper::evaluateValueInElement(elementNodesDisplacement, N); + + const QuadratureContext context{ + element, N, dN_dq_ref, jacobian, measure, restPosition, displacement}; + + std::array::SourceDerivative, + NumberOfNodesInElement> sourceDerivative{}; + + for (const auto& source : l_nonConstantSources) + { + for (sofa::Size b = 0; b < NumberOfNodesInElement; ++b) + { + sourceDerivative[b] += source->evaluateStiffness(context, b); + } + } + + for (sofa::Size a = 0; a < NumberOfNodesInElement; ++a) + { + const auto weightTimesShapeFunction = static_cast(weight) * N[a]; + + for (sofa::Size b = 0; b < NumberOfNodesInElement; ++b) + { + dfdx(element[a] * spatial_dimensions, element[b] * spatial_dimensions) + += sourceDerivative[b] * weightTimesShapeFunction; + } + } + } + } } template From ae43f8c194e232048a520f8e5f257aaec585624e Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:44:51 +0200 Subject: [PATCH 12/17] Add addDForce to the integrator --- .../fem/elastic/FEMSourceTermIntegrator.h | 2 +- .../fem/elastic/FEMSourceTermIntegrator.inl | 68 ++++++++++++++++++- 2 files changed, 66 insertions(+), 4 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 7f58113f315..b07abafab4d 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 @@ -113,7 +113,7 @@ class FEMSourceTermIntegrator : const sofa::DataVecDeriv_t& v) override; /** - * @brief No-op. + * @brief Applies the tangent of the terms of l_nonConstantSources to a displacement increment. */ void addDForce(const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& df, 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 d2cb70710c1..39f822cfb11 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 @@ -263,9 +263,71 @@ void FEMSourceTermIntegrator::addDForce(const sofa::core sofa::DataVecDeriv_t& df, const sofa::DataVecDeriv_t& dx) { - SOFA_UNUSED(mparams); - SOFA_UNUSED(df); - SOFA_UNUSED(dx); + if (this->isComponentStateInvalid() || l_nonConstantSources.empty() + || !d_useTangentStiffness.getValue()) + { + return; + } + + const auto kFactor = static_cast(sofa::core::mechanicalparams::kFactorIncludingRayleighDamping( + mparams, this->rayleighStiffness.getValue())); + + auto forceDerivativeAccessor = sofa::helper::getWriteAccessor(df); + const sofa::helper::ReadAccessor displacementDerivativeAccessor = sofa::helper::getReadAccessor(dx); + + const auto restPositionsAccessor = this->mstate->readRestPositions(); + 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, NumberOfNodesInElement> elementNodesRestCoordinates = + extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); + const std::array, NumberOfNodesInElement> elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, positionsAccessor.ref()); + + std::array, NumberOfNodesInElement> 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 measure = static_cast(sofa::type::absGeneralizedDeterminant(jacobian)); + + const auto restPosition = + FiniteElement::Helper::evaluateValueInElement(elementNodesRestCoordinates, N); + const auto displacement = + FiniteElement::Helper::evaluateValueInElement(elementNodesDisplacement, N); + + const QuadratureContext context{ + element, N, dN_dq_ref, jacobian, measure, restPosition, displacement}; + + sofa::Deriv_t contraction{}; + for (const auto& source : l_nonConstantSources) + { + for (sofa::Size b = 0; b < NumberOfNodesInElement; ++b) + { + contraction += source->evaluateStiffness(context, b) + * displacementDerivativeAccessor[element[b]]; + } + } + + for (sofa::Size a = 0; a < NumberOfNodesInElement; ++a) + { + forceDerivativeAccessor[element[a]] += + contraction * (kFactor * static_cast(weight) * N[a]); + } + } + } } template From 28a7f71f8a95b0fd455a7a7cbbfc211a51b9e8a6 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:46:20 +0200 Subject: [PATCH 13/17] Fix comments and reduce Newton iters --- .../SolidMechanics/FEM/PressureSourceTerm.scn | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn index d12d4e8608f..4c50ab40e7b 100644 --- a/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn +++ b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn @@ -8,15 +8,6 @@ Cyan: SurfacePressureForceField on a triangulated surface Orange: SurfacePressureForceField on the quads of that surface - This is a follower load: the traction is p n on the deformed surface, so both the normal and - the area are those of the current configuration. The source term is therefore attached to the - nonConstantSources link of the integrator, which integrates it again at every call to addForce. - Nothing in PressureSourceTerm itself says whether it is a dead or a follower load: the link it - hangs from decides. - - Neither component contributes a tangent stiffness here, so the Newton iteration is a fixed - point on the load and needs several iterations to converge. - The pressure is applied on the same surface twice, once through its triangles and once through its quads, which SurfacePressureForceField and FEMSourceTermIntegrator handle by separate paths. @@ -53,7 +44,7 @@ - + @@ -87,7 +78,7 @@ - + @@ -131,7 +122,7 @@ - + @@ -154,7 +145,7 @@ - + From 12206eb4eceb8a47db51b82ac5e6503b191149dc Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 18:48:39 +0200 Subject: [PATCH 14/17] Debug build --- .../solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl | 1 + 1 file changed, 1 insertion(+) 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 39f822cfb11..a25cd6876b5 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,7 @@ #pragma once #include #include +#include namespace sofa::component::solidmechanics::fem::elastic { From 7b3ab06a39d092d7eb9794d1f6d92fc0e4e40418 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 10 Sep 2026 19:05:58 +0200 Subject: [PATCH 15/17] Add example for StressSourceTerm --- .../SolidMechanics/FEM/StressSourceTerm.scn | 107 ++++++++++++++++-- 1 file changed, 95 insertions(+), 12 deletions(-) diff --git a/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn b/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn index fb3723ee104..f381dc7cacc 100644 --- a/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn +++ b/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn @@ -1,13 +1,19 @@ - @@ -30,20 +36,97 @@ - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - + @@ -62,7 +145,7 @@ - + - + @@ -95,7 +178,7 @@ - + From 611a182f624669fe724b4c53b3f544af35e25327 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 6 Nov 2026 13:30:55 +0100 Subject: [PATCH 16/17] Touch up comments --- .../fem/elastic/BaseSourceTerm.h | 77 ++++++------------- .../fem/elastic/FEMSourceTermIntegrator.h | 47 +++-------- .../fem/elastic/PressureSourceTerm.h | 12 +-- .../fem/elastic/VectorSourceTerm.h | 17 ++-- .../SolidMechanics/FEM/PressureSourceTerm.scn | 20 ++--- .../SolidMechanics/FEM/StressSourceTerm.scn | 23 ++---- 6 files changed, 59 insertions(+), 137 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h index 84608175f8c..63a516b3b31 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h @@ -44,10 +44,7 @@ namespace sofa::component::solidmechanics::fem::elastic /** * @struct QuadratureContext - * @brief Everything the integrator knows at one quadrature point. - * - * Built once per quadrature point and handed to every integrated term. - * A source term reads from it and returns an integrand. + * @brief Geometry of one quadrature point. * * @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). @@ -72,36 +69,32 @@ struct QuadratureContext using GradientShapeFunctions = sofa::type::Mat; using Jacobian = sofa::type::Mat; - /// Node indices of the element being integrated, with which a term gathers its own nodal - /// degrees of freedom. + /// Node indices of the element. const Element& element; - /// Shape function value at this quadrature point. + /// Shape function values. ShapeFunctions N; - /// Reference-space gradients of the shape functions at this quadrature point. + /// Reference-space gradients of the shape functions. GradientShapeFunctions gradientShapeFunctions; - /// dx/dq of the reference-to-physical mapping, on the configuration the integrator chose. + /// dx/dq of the reference-to-physical mapping. Jacobian jacobian; - /// \f$ |\det J| \f$, for information only: the integrator applies it, a term must not. + /// \f$ |\det J| \f$ Real measure; - /// Interpolated rest position at this quadrature point. + /// Interpolated rest position. Coord restPosition; - /// Interpolated displacement at this quadrature point. + /// Interpolated displacement. Deriv displacement; }; /** - * @brief Unit normal of an element, from the jacobian of its mapping. - * - * Defined only where the element spans one dimension less than the space it lives in: a surface - * element in 3D, an edge in 2D. Its orientation follows the node ordering of the element. + * @brief Unit normal of an element of codimension 1, oriented by its node ordering. * - * @param jacobian dx/dq of the reference-to-physical mapping at the point of interest. + * @param jacobian dx/dq of the reference-to-physical mapping. */ template sofa::type::Vec elementNormal( @@ -126,23 +119,15 @@ sofa::type::Vec elementNormal( } /** - * @brief Derivative of the area normal of an element with respect to one of its node positions. - * - * The area normal is \f$ |\det J| \, n \f$: in 3D the cross product of the two columns of the - * jacobian, in 2D the tangent rotated by a quarter turn. Differentiating it with respect to the - * position of the node whose shape function is N gives + * @brief Derivative of the area normal of an element with respect to a node position: * - * \f[ [t_0]_\times \frac{\partial N}{\partial q_1} - [t_1]_\times \frac{\partial N}{\partial q_0} - * \quad \text{in 3D,} \qquad R \, \frac{\partial N}{\partial q_0} \quad \text{in 2D,} \f] + * \f[ \frac{\partial (|\det J| \, n)}{\partial x} + * = \sum_i \frac{\partial (|\det J| \, n)}{\partial t_i} \frac{\partial N}{\partial q_i}, \f] * - * where \f$ t_i \f$ is the i-th column of the jacobian, \f$ [\,]_\times \f$ the cross-product - * matrix and R the quarter turn. + * with \f$ t_i \f$ the columns of J and N the shape function of the node. * - * Defined only where the element spans one dimension less than the space it lives in, as - * elementNormal is. - * - * @param jacobian dx/dq of the reference-to-physical mapping at the point of interest. - * @param gradientOfShapeFunction dN/dq of the node the derivative is taken with respect to. + * @param jacobian dx/dq of the reference-to-physical mapping. + * @param gradientOfShapeFunction dN/dq of the node. */ template sofa::type::Mat elementAreaNormalDerivative( @@ -182,9 +167,7 @@ sofa::type::Mat elementAreaNormalD /** * @class BaseSourceTerm - * @brief A source density whose value is determined by the geometry, not by the solution. - * - * The component calculates the integrand in evaluate() given a QuadratureContext. + * @brief A source term of the weak form. * * @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). @@ -205,23 +188,18 @@ class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; - /// Derivative of a nodal force with respect to the position of one node of its element. + /// Derivative of a nodal force with respect to a node position. using SourceDerivative = sofa::type::Mat; /** * @brief Source density at one quadrature point, per unit physical measure. * - * The integrator weights it by \f$ w \, |\det J| \, N_a \f$. - * * @param context Geometry of the quadrature point. */ virtual Deriv evaluate(const QuadratureContext_t& context) const = 0; /** - * @brief The Data this term's integrand reads, beyond the geometry it is handed. - * - * Access to this Data is enabled so that a callback on FEMSourceTermIntegrator can be activated. - * They are tracked and the cached integral is rebuilt whenever one of them is written. + * @brief Data the integrand depends on; used for callbacks. */ virtual sofa::type::vector integrandInputs() const { @@ -229,15 +207,11 @@ class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent } /** - * @brief Derivative of the measure-weighted density with respect to a node position. - * - * \f$ \partial (|\det J| \, r) / \partial x_b \f$, the measure included. The integrator weights - * this one by \f$ w \, N_a \f$ alone. - * - * A term that does not follow the configuration returns zero. + * @brief Derivative of the measure-weighted density with respect to a node position, + * \f$ \partial (|\det J| \, r) / \partial x_b \f$, with r the source density. * * @param context Geometry of the quadrature point. - * @param node Index in the element of the node the derivative is taken with respect to. + * @param node Index of the node in the element. */ virtual SourceDerivative evaluateStiffness(const QuadratureContext_t& context, sofa::Size node) const = 0; @@ -247,12 +221,9 @@ class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent BaseSourceTerm() = default; /** - * @brief Value of a nodal property interpolated at the quadrature point. - * - * The property is gathered at the nodes of the element being integrated and combined with the - * shape functions evaluated at that point. + * @brief Nodal property interpolated at the quadrature point. * - * @param property The component holding the nodal values. + * @param property The nodal values. * @param context Geometry of the quadrature point. */ template 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 b07abafab4d..dcee7fc3096 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 @@ -37,15 +37,9 @@ namespace sofa::component::solidmechanics::fem::elastic /** * @class FEMSourceTermIntegrator - * @brief Integrates a source density 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 density evaluated by a linked BaseSourceTerm at each quadrature point. - * - * The link a term is attached to decides when it is integrated: - * l_constantSources are integrated once on the rest configuration and cached, while - * l_nonConstantSources are integrated again at every call to addForce, on the configuration it is - * given. The latter are follower loads. + * @brief Integrates source terms into consistent nodal loads, + * \f$ \int_{\Omega} N_a \, r \, d\Omega \f$, with \f$ N_a \f$ the shape function of node a and r + * the source density. * * @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). @@ -78,31 +72,21 @@ class FEMSourceTermIntegrator : sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_constantSources; /** - * @brief Source terms integrated again at every step, on the current configuration. + * @brief Source terms integrated on the current configuration. */ sofa::MultiLink, BaseSourceTerm, sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_nonConstantSources; /** * @brief Initializes the component. - * - * This method performs several initialization steps: - * 1. Initializes the base force field. - * 2. Initializes the topology accessor. - * 3. Validates the linked source terms. - * 4. Tracks the Data the constant source terms read. - * 5. Integrates the source terms into the nodal force. */ void init() override; /** - * @brief Adds the nodal source term to the RHS vector. - * - * The cached load of the constant sources is added as is, and the non-constant ones are - * integrated on the current configuration x that this method is given. + * @brief Adds the nodal loads of the source terms. * * @param mparams Mechanical parameters for the computation. - * @param f The force vector to which the source term will be added. + * @param f The force vector, to which the source term will be added. * @param x The current positions. * @param v The current velocities. */ @@ -113,16 +97,14 @@ class FEMSourceTermIntegrator : const sofa::DataVecDeriv_t& v) override; /** - * @brief Applies the tangent of the terms of l_nonConstantSources to a displacement increment. + * @brief Applies the tangent of the non-constant source terms to a displacement increment. */ void addDForce(const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& df, const sofa::DataVecDeriv_t& dx) override; /** - * @brief Assembles the tangent of the terms of l_nonConstantSources. - * - * The tangent of a term that follows the configuration is not symmetric. + * @brief Assembles the tangent of the non-constant source terms. */ void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; @@ -145,16 +127,10 @@ class FEMSourceTermIntegrator : protected: - /** - * @brief Default constructor. - */ FEMSourceTermIntegrator(); /** * @brief Ensures that valid source terms are linked, falling back to the current context. - * - * If both constant and non-constant source links are left empty, the BaseSourceTerm components - * found in the current context are used as non-constant sources. */ void validateSources(); @@ -164,13 +140,12 @@ class FEMSourceTermIntegrator : void doUpdateInternal() override; /** - * @brief Integrates the terms of l_constantSources on the rest configuration into - * m_constantForce, as consistent nodal forces. + * @brief Integrates the constant source terms on the rest configuration. */ void assembleConstantForce(); /** - * @brief Integrates the terms of l_nonConstantSources as consistent nodal forces. + * @brief Integrates the non-constant source terms. * * @param x The configuration the terms are integrated on. * @param f The force vector the nodal forces are accumulated into. @@ -179,7 +154,7 @@ class FEMSourceTermIntegrator : sofa::VecDeriv_t& f) const; /** - * @brief Nodal load of every term in l_constantSources. + * @brief Nodal load of the constant source terms. */ sofa::VecDeriv_t m_constantForce; }; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h index 8ebdc3de796..f92a0808543 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h @@ -59,11 +59,7 @@ class NodalPressure : public sofa::core::BaseNodalProperty sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_pressure; /** - * @brief Initializes the component and checks that a pressure is linked. + * @brief Initializes the component. */ void init() override; /** - * @brief The linked pressure interpolated at the quadrature point, times the unit normal. + * @brief The pressure interpolated at the quadrature point, times the unit normal. */ Deriv evaluate(const QuadratureContext_t& context) const override; /** - * @brief The nodal values this term interpolates, for FEMSourceTermIntegrator to track. + * @brief The nodal values of the pressure. */ sofa::type::vector integrandInputs() const override; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h index 5df045e0d2c..c0eb99e2675 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h @@ -59,10 +59,6 @@ class NodalSourceDensity : public sofa::core::BaseNodalProperty sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_sourceDensity; /** - * @brief Initializes the component and checks that a source density is linked. + * @brief Initializes the component. */ void init() override; /** - * @brief The linked source density interpolated at the quadrature point. + * @brief The source density interpolated at the quadrature point. */ Deriv evaluate(const QuadratureContext_t& context) const override; /** - * @brief The nodal values this term interpolates, for FEMSourceTermIntegrator to track. + * @brief The nodal values of the source density. */ sofa::type::vector integrandInputs() const override; /** - * @brief The source density times the derivative of the measure of the element. + * @brief The source density times the derivative of the measure: * - * \f$ \partial (|\det J| \, r) / \partial x_b - * = |\det J| \; r \otimes \left( (J^{-1})^T \frac{\partial N_b}{\partial q} \right) \f$, - * the inverse being the left pseudo-inverse when J is not square. + * \f[ \frac{\partial (|\det J| \, r)}{\partial x_b} + * = |\det J| \; r \otimes \left( J^{-T} \frac{\partial N_b}{\partial q} \right) \f] */ SourceDerivative evaluateStiffness(const QuadratureContext_t& context, sofa::Size node) const override; diff --git a/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn index 4c50ab40e7b..3ba89a90079 100644 --- a/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn +++ b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn @@ -1,16 +1,10 @@ - @@ -39,7 +33,7 @@ - + @@ -117,7 +111,7 @@ - + diff --git a/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn b/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn index f381dc7cacc..61fa3584ca3 100644 --- a/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn +++ b/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn @@ -1,19 +1,10 @@ - @@ -41,7 +32,7 @@ - + @@ -117,7 +108,7 @@ - + From 5a37d53e6187e1daebe6b291289ce1ff040ecab2 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 6 Nov 2026 14:40:19 +0100 Subject: [PATCH 17/17] update scene for integration schemes --- .../SolidMechanics/FEM/PressureSourceTerm.scn | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn index 3ba89a90079..81d0718da1b 100644 --- a/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn +++ b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn @@ -11,9 +11,9 @@ + - @@ -38,8 +38,7 @@ - - + @@ -72,8 +71,7 @@ - - + @@ -116,8 +114,7 @@ - - + @@ -139,8 +136,7 @@ - - +