diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt
index fddd6033a0d..eeda2d5edeb 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt
@@ -13,8 +13,12 @@ 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}/NonConstantSourceTerm.h
${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.h
${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.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
@@ -35,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
@@ -68,7 +73,10 @@ 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}/NonConstantSourceTerm.cpp
${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.cpp
+ ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.cpp
${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.cpp
${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.cpp
${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMForceField.cpp
@@ -79,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/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/FEMSourceTermIntegrator.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.cpp
new file mode 100644
index 00000000000..f7a908ed703
--- /dev/null
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.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_INTEGRATOR_CPP
+
+#include
+
+#include
+#include
+#include
+
+namespace sofa::component::solidmechanics::fem::elastic
+{
+
+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 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 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/FEMSourceTermIntegrator.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h
new file mode 100644
index 00000000000..2cf86b8e808
--- /dev/null
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h
@@ -0,0 +1,237 @@
+/******************************************************************************
+* 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
+#include
+#include
+#include
+
+#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP)
+#include
+#endif
+
+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. 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).
+ */
+template
+class FEMSourceTermIntegrator :
+ public sofa::core::behavior::ForceField,
+ public virtual sofa::core::behavior::TopologyAccessor
+{
+public:
+ using DataTypes = TDataTypes;
+ using ElementType = TElementType;
+ 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;
+ 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:
+
+ /**
+ * @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 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.
+ *
+ * This method performs several initialization steps:
+ * 1. Initializes the base force field.
+ * 2. Initializes the topology accessor.
+ * 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 nodal source term to the RHS vector.
+ *
+ * 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, used only by displacement-dependent terms.
+ * @param v The current velocities (unused in this implementation).
+ */
+ void addForce(
+ const sofa::core::MechanicalParams* mparams,
+ DataVecDeriv& f,
+ const DataVecCoord& x,
+ const DataVecDeriv& v) override;
+
+ /**
+ * @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,
+ DataVecDeriv& df,
+ const DataVecDeriv& dx) override;
+
+ /**
+ * @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;
+
+ using sofa::core::behavior::ForceField::getPotentialEnergy;
+ /**
+ * @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 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:
+
+ /**
+ * @brief Default constructor.
+ */
+ 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();
+
+ /**
+ * @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 VecDeriv& nodalSourceTerm,
+ VecDeriv& result) const;
+
+ /**
+ * @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.
+ *
+ * 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.
+ */
+ VecDeriv m_constantForce;
+};
+
+#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/FEMSourceTermIntegrator.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl
new file mode 100644
index 00000000000..d15a730aa02
--- /dev/null
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl
@@ -0,0 +1,412 @@
+/******************************************************************************
+* 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
+
+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_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},
+ [this](const sofa::core::DataTracker&)
+ {
+ if (!this->isComponentStateInvalid() && this->l_topology && this->mstate && !l_constantSources.empty())
+ {
+ assembleGlobalMatrix();
+ assembleConstantForce();
+ }
+
+ return this->getComponentState();
+ }, {});
+}
+
+template
+void FEMSourceTermIntegrator::init()
+{
+ sofa::core::behavior::ForceField::init();
+
+ if (!this->isComponentStateInvalid())
+ {
+ sofa::core::behavior::TopologyAccessor::init();
+ }
+
+ if (!this->isComponentStateInvalid())
+ {
+ this->validateSources();
+ }
+
+ if (!this->isComponentStateInvalid() && this->l_topology && this->mstate && !l_constantSources.empty())
+ {
+ this->assembleGlobalMatrix();
+ this->assembleConstantForce();
+ }
+
+ if (!this->isComponentStateInvalid())
+ {
+ this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid);
+ }
+}
+
+template
+void FEMSourceTermIntegrator::validateSources()
+{
+ // Gather all matching source components in Context if a link is left empty
+ auto fallbackToContext = [this](auto& link)
+ {
+ using SourceType = typename std::remove_reference_t::DestType;
+
+ if (link.empty())
+ {
+ const auto sourcesInContext = this->getContext()->template getObjects(
+ sofa::core::objectmodel::BaseContext::Local);
+
+ 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() && l_nonConstantSources.empty(), this)
+ << "No source term linked, and none found in the current context '"
+ << this->getContext()->getName() << "'. This component has zero force contribution.";
+}
+
+template
+void FEMSourceTermIntegrator::assembleGlobalMatrix()
+{
+ const auto& elements = FiniteElement::getElementSequence(*this->l_topology);
+ sofa::type::vector elementMatrices;
+
+ // 1. compute the geometry-only matrix of each element
+ calculateElementMatrix(elements, elementMatrices);
+
+ // 2. scatter the element matrices into the global matrix
+ initializeGlobalMatrix(elements, elementMatrices);
+}
+
+template
+void FEMSourceTermIntegrator::calculateElementMatrix(
+ const auto& elements, sofa::type::vector& elementMatrices)
+{
+ const auto restPositionsAccessor = this->mstate->readRestPositions();
+ elementMatrices.resize(elements.size());
+
+ const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue());
+
+ for (sofa::Index elementId = 0; elementId < elements.size(); ++elementId)
+ {
+ const auto& element = elements[elementId];
+ auto& elementMatrix = elementMatrices[elementId];
+
+ const std::array elementNodesRestCoordinates =
+ 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] : quadratureRule)
+ {
+ 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);
+
+ const auto NT_N = sofa::type::dyad(N, N);
+
+ elementMatrix += (weight * detJ) * NT_N;
+ }
+ }
+}
+
+template
+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 (sofa::Index elementId = 0; elementId < elements.size(); ++elementId)
+ {
+ const auto& element = elements[elementId];
+ const auto& elementMatrix = elementMatrices[elementId];
+
+ 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 FEMSourceTermIntegrator::applyGlobalMatrix(
+ 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)
+ {
+ 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
+ VecDeriv sourceTerms(size, Deriv{});
+
+ for (const auto& source : l_constantSources)
+ {
+ for (sofa::Index i = 0; i < size; ++i)
+ sourceTerms[i] += source->getNodeProperty(i);
+ }
+
+ m_constantForce.assign(size, Deriv{});
+ applyGlobalMatrix(sourceTerms, m_constantForce);
+}
+
+template
+void FEMSourceTermIntegrator::addForce(const sofa::core::MechanicalParams* mparams,
+ DataVecDeriv& f,
+ const DataVecCoord& x,
+ const DataVecDeriv& v)
+{
+ SOFA_UNUSED(mparams);
+ SOFA_UNUSED(v);
+
+ if (this->isComponentStateInvalid())
+ {
+ return;
+ }
+
+ auto forceAccessor = sofa::helper::getWriteAccessor(f);
+
+ for (sofa::Index i = 0; i < m_constantForce.size(); ++i)
+ {
+ 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,
+ DataVecDeriv& df,
+ const DataVecDeriv& 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)
+{
+ 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 DataVecCoord& x) const
+{
+ SOFA_UNUSED(mparams);
+
+ 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/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 2c37e45f155..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
@@ -27,16 +27,20 @@ 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 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);
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);
@@ -67,16 +71,20 @@ const char* getModuleVersion()
void registerObjects(sofa::core::ObjectFactory* factory)
{
registerBeamFEMForceField(factory);
+ registerConstantSourceTerm(factory);
registerCorotationalFEMForceField(factory);
+ registerFEMSourceTermIntegrator(factory);
registerFastTetrahedralCorotationalForceField(factory);
registerHexahedralFEMForceField(factory);
registerHexahedralFEMForceFieldAndMass(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/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..b038c7bb24d
--- /dev/null
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp
@@ -0,0 +1,312 @@
+/******************************************************************************
+* 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
+#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;
+
+// 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:
+ 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();
+ }
+
+ 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.
+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);
+}
+
+// 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