diff --git a/Sofa/framework/Simulation/Core/CMakeLists.txt b/Sofa/framework/Simulation/Core/CMakeLists.txt index cf32c217c79..b18c01d2992 100644 --- a/Sofa/framework/Simulation/Core/CMakeLists.txt +++ b/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 diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/DifferentialOperations.cpp b/Sofa/framework/Simulation/Core/src/sofa/simulation/DifferentialOperations.cpp new file mode 100644 index 00000000000..38ea8847b19 --- /dev/null +++ b/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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include +#include + +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 diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/DifferentialOperations.h b/Sofa/framework/Simulation/Core/src/sofa/simulation/DifferentialOperations.h new file mode 100644 index 00000000000..775a1517bca --- /dev/null +++ b/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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +#include +#include + +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); +}; + +} diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.cpp b/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.cpp index 51ca30266b2..c167fc7404d 100644 --- a/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.cpp +++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.cpp @@ -20,6 +20,7 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #include +#include 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); diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.h b/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.h index 80a12a60945..8ebccdb2367 100644 --- a/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraphMechanicalOperations.h +++ b/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 $