[Simulation.Core] Differential operations

Sofa/framework/Simulation/Core/CMakeLists.txt

@@ -17,6 +17,7 @@ set(HEADER_FILES
     ${SRC_ROOT}/DefaultAnimationLoop.h
     ${SRC_ROOT}/DefaultVisualManagerLoop.h
     ${SRC_ROOT}/DeleteVisitor.h
+    ${SRC_ROOT}/DifferentialOperations.h
     ${SRC_ROOT}/ExportDotVisitor.h
     ${SRC_ROOT}/ExportGnuplotVisitor.h
     ${SRC_ROOT}/ExportVisualModelOBJVisitor.h
@@ -169,6 +170,7 @@ set(SOURCE_FILES
     ${SRC_ROOT}/DefaultAnimationLoop.cpp
     ${SRC_ROOT}/DefaultVisualManagerLoop.cpp
     ${SRC_ROOT}/DeleteVisitor.cpp
+    ${SRC_ROOT}/DifferentialOperations.cpp
     ${SRC_ROOT}/ExportDotVisitor.cpp
     ${SRC_ROOT}/ExportGnuplotVisitor.cpp
     ${SRC_ROOT}/ExportVisualModelOBJVisitor.cpp

Sofa/framework/Simulation/Core/src/sofa/simulation/DifferentialOperations.cpp

@@ -0,0 +1,77 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* This program is free software; you can redistribute it and/or modify it     *
+* under the terms of the GNU Lesser General Public License as published by    *
+* the Free Software Foundation; either version 2.1 of the License, or (at     *
+* your option) any later version.                                             *
+*                                                                             *
+* This program is distributed in the hope that it will be useful, but WITHOUT *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
+* for more details.                                                           *
+*                                                                             *
+* You should have received a copy of the GNU Lesser General Public License    *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#include <sofa/simulation/DifferentialOperations.h>
+#include <sofa/core/MechanicalParams.h>
+
+namespace sofa::simulation::common
+{
+
+void DifferentialOperations::pushforwardCoord(
+    const MappingGraph& mappingGraph, const core::MechanicalParams& mparams,
+    core::MultiVecCoordId coordVectorId)
+{
+    mappingGraph.algorithms.traverseTopDown_([&](core::BaseMapping& mapping)
+    {
+        mapping.apply(&mparams, coordVectorId, coordVectorId);
+    });
+}
+
+void DifferentialOperations::pushforwardTangent(
+    const MappingGraph& mappingGraph, const core::MechanicalParams& mparams,
+    core::MultiVecDerivId tangentVectorId)
+{
+    mappingGraph.algorithms.traverseTopDown_([&](core::BaseMapping& mapping)
+    {
+        mapping.applyJ(&mparams, tangentVectorId, tangentVectorId);
+    });
+}
+
+void DifferentialOperations::pullbackCotangent(
+    const MappingGraph& mappingGraph, const core::MechanicalParams& mparams,
+    core::MultiVecDerivId cotangentVectorId, bool ignoreMappingFlag)
+{
+    mappingGraph.algorithms.traverseBottomUp_([&](core::BaseMapping& mapping)
+    {
+        if (mapping.areForcesMapped() || ignoreMappingFlag)
+        {
+            mapping.applyJT(&mparams, cotangentVectorId, cotangentVectorId);
+        }
+    });
+}
+
+void DifferentialOperations::pullbackCotangentTangent(
+    const MappingGraph& mappingGraph, const core::MechanicalParams& mparams,
+    core::MultiVecDerivId cotangentTangentVectorId)
+{
+    mappingGraph.algorithms.traverseBottomUp_(
+        [&](core::BaseMapping& mapping)
+        {
+            mapping.applyJT(&mparams, cotangentTangentVectorId, cotangentTangentVectorId);
+            if (mparams.kFactor() != 0)
+            {
+                // ideally, the child cotangent vector must be provided here, instead of implicitly getting it in the mapping
+                mapping.applyDJT(&mparams, cotangentTangentVectorId, cotangentTangentVectorId);
+            }
+        });
+}
+
+}  // namespace sofa::simulation::common

Sofa/framework/Simulation/Core/src/sofa/simulation/DifferentialOperations.h

@@ -0,0 +1,121 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* This program is free software; you can redistribute it and/or modify it     *
+* under the terms of the GNU Lesser General Public License as published by    *
+* the Free Software Foundation; either version 2.1 of the License, or (at     *
+* your option) any later version.                                             *
+*                                                                             *
+* This program is distributed in the hope that it will be useful, but WITHOUT *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
+* for more details.                                                           *
+*                                                                             *
+* You should have received a copy of the GNU Lesser General Public License    *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+
+#include <sofa/simulation/config.h>
+
+#include <sofa/core/MultiVecId.h>
+#include <sofa/simulation/MappingGraph.h>
+
+namespace sofa::simulation::common
+{
+
+/**
+ * @brief Differential and geometric operations over a MappingGraph.
+ *
+ * Provides static operations for pushing forward and pulling back multi-vectors
+ * through the hierarchy of mappings defined within a MappingGraph.
+ *
+ * In physical simulations, primal quantities (coordinates, tangent vectors such as velocities
+ * or displacements) are propagated top-to-bottom through pushforward operations (applying
+ * the mapping function or its Jacobian $\mathbf{J}$). Dual quantities (cotangent vectors such
+ * as forces) are transformed bottom-to-top through pullback operations (using the transpose
+ * Jacobian $\mathbf{J}^T$).
+ */
+class SOFA_SIMULATION_CORE_API DifferentialOperations
+{
+public:
+
+    /**
+     * @brief Push forward primal coordinates top-down through the mapping graph.
+     *
+     * Traverses the mapping graph from top to bottom and applies each mapping function
+     * to transform primal configuration coordinates (e.g. positions) across coordinate frames.
+     *
+     * @param mappingGraph The mapping graph containing the topology and mappings to traverse.
+     * @param mparams Mechanical parameters associated with the current operation.
+     * @param coordVectorId Identifier of the coordinate multi-vector to push forward.
+     */
+    static void pushforwardCoord(
+        const MappingGraph& mappingGraph, const core::MechanicalParams& mparams,
+        core::MultiVecCoordId coordVectorId);
+
+    /**
+     * @brief Push forward tangent vectors top-down via the mapping Jacobian.
+     *
+     * Traverses the mapping graph from top to bottom and applies the Jacobian matrix $\mathbf{J}$
+     * of each mapping to propagate tangent quantities (e.g. velocities, displacements $dx$)
+     * via Jacobian-vector products ($\mathbf{J}v$).
+     *
+     * @param mappingGraph The mapping graph containing the topology and mappings to traverse.
+     * @param mparams Mechanical parameters associated with the current operation.
+     * @param tangentVectorId Identifier of the derivative multi-vector representing tangent quantities to push forward.
+     */
+    static void pushforwardTangent(
+        const MappingGraph& mappingGraph, const core::MechanicalParams& mparams,
+        core::MultiVecDerivId tangentVectorId);
+
+    /**
+     * @brief Pull back cotangent vectors bottom-up via the transpose of the mapping Jacobian.
+     *
+     * Traverses the mapping graph from bottom to top and applies the transposed Jacobian $\mathbf{J}^T$
+     * of each mapping to pull back dual quantities (e.g. forces, momentum) via
+     * vector-Jacobian products ($v^T\mathbf{J}$).
+     *
+     * @param mappingGraph The mapping graph containing the topology and mappings to traverse.
+     * @param mparams Mechanical parameters associated with the current operation.
+     * @param cotangentVectorId Identifier of the derivative multi-vector representing cotangent quantities to pull back.
+     * @param ignoreMappingFlag If true, forces are pulled back through all mappings regardless of
+     *                          whether their internal force-mapping flag (`areForcesMapped()`) is enabled.
+     */
+    static void pullbackCotangent(
+        const MappingGraph& mappingGraph, const core::MechanicalParams& mparams,
+        core::MultiVecDerivId cotangentVectorId,
+        bool ignoreMappingFlag = true);
+
+    /**
+     * @brief Pull back the differential (tangent variation) of a cotangent vector bottom-up.
+     *
+     * Computes the total differential of a vector-Jacobian product (VJP), corresponding to
+     * the variation of a pulled-back cotangent quantity (such as force variation @f$ df_x = d(\mathbf{J}^T f_y) @f$):
+     * @f[
+     *     df_x = \mathbf{J}^T df_y + d(\mathbf{J}^T) f_y
+     * @f]
+     *
+     * Traverses the mapping graph from bottom to top and accumulates:
+     * 1. The first-order pullback of the cotangent differential via the transpose Jacobian:
+     *    @f$ \mathbf{J}^T df_y @f$ (`applyJT()`).
+     * 2. When the stiffness factor @c kFactor is non-zero, the geometric stiffness contribution
+     *    arising from the non-linear variation of the mapping's Jacobian (a Hessian-vector contraction):
+     *    @f$ d(\mathbf{J}^T) f_y @f$ (`applyDJT()`).
+     *
+     * @param mappingGraph The mapping graph containing the topology and mappings to traverse.
+     * @param mparams Mechanical parameters associated with the current operation (including @c kFactor).
+     * @param cotangentTangentVectorId Identifier of the derivative multi-vector storing the cotangent
+     *                                 tangent variations (e.g. @f$ df @f$) to pull back.
+     */
+    static void pullbackCotangentTangent(
+        const MappingGraph& mappingGraph, const core::MechanicalParams& mparams,
+        core::MultiVecDerivId cotangentTangentVectorId);
+};
+
+}

Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.cpp

@@ -20,6 +20,7 @@
 * Contact information: contact@sofa-framework.org                             *
 ******************************************************************************/
 #include <sofa/simulation/MappingGraphMechanicalOperations.h>
+#include <sofa/simulation/DifferentialOperations.h>
 
 namespace sofa::simulation::common
 {
@@ -41,7 +42,7 @@ void MappingGraphMechanicalOperations::projectResponse(const MappingGraph& mappi
 void MappingGraphMechanicalOperations::computeForce(const MappingGraph& mappingGraph,
                                                     core::MultiVecDerivId result,
                                                     bool clearForceBefore,
-                                                    bool accumulateForcesFromMappedStates,
+                                                    bool pullbackForces,
                                                     TaskScheduler* taskScheduler)
 {
     //assumes the mapping graph is valid and properly initialized
@@ -79,16 +80,14 @@ void MappingGraphMechanicalOperations::computeForce(const MappingGraph& mappingG
         forceField.addForce(&mparams, result);
     }, sofa::simulation::VisitorApplication::ALL_NODES, taskScheduler);
 
-    if (accumulateForcesFromMappedStates)
+    if (pullbackForces)
     {
         /**
-         * Compute f_in += J^T * f_out using mappings in the mapping graph. This operation must be
-         * performed in a bottom-up order to ensure correct force accumulation.
+         * Pull back force f_in += J^T * f_out using the mappings in the mapping graph. Force is a
+         * dual quantity (a cotangent vector). That is why this operation must be performed in a
+         * bottom-up order.
          */
-        mappingGraph.algorithms.traverseBottomUp_([&](core::BaseMapping& mapping)
-        {
-            mapping.applyJT(&mparams, result, result);
-        });
+        DifferentialOperations::pullbackCotangent(mappingGraph, mparams, result);
     }
 }
 void MappingGraphMechanicalOperations::addMBKv(const MappingGraph& mappingGraph,
@@ -124,14 +123,7 @@ void MappingGraphMechanicalOperations::addMBKv(const MappingGraph& mappingGraph,
 
     if (accumulate)
     {
-        mappingGraph.algorithms.traverseBottomUp_([&](core::BaseMapping& mapping)
-        {
-            mapping.applyJT(&mparams, df, df);
-            if( mparams.kFactor() != 0 )
-            {
-                mapping.applyDJT(&mparams, df, df);
-            }
-        });
+        DifferentialOperations::pullbackCotangentTangent(mappingGraph, mparams, df);
     }
 
     mparams.setDx(dx);

Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.h

@@ -53,7 +53,7 @@ class SOFA_SIMULATION_CORE_API MappingGraphMechanicalOperations : public Mechani
     using MechanicalOperations::projectResponse;
 
     /// Compute the current force (given the latest propagated position and velocity)
-    void computeForce(const MappingGraph& mappingGraph, core::MultiVecDerivId result, bool clearForceBefore, bool accumulateForcesFromMappedStates, TaskScheduler* taskScheduler);
+    void computeForce(const MappingGraph& mappingGraph, core::MultiVecDerivId result, bool clearForceBefore, bool pullbackForces, TaskScheduler* taskScheduler);
     using MechanicalOperations::computeForce;
 
     /// accumulate $ df += (m M + b B + k K) velocity $