From 315943ad77238e78d4d9c7d830af2385ca38f9d0 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 5 Mar 2026 09:38:40 +0100 Subject: [PATCH 1/2] move MappingGraph class from Sofa.Component.LinearSystem to Sofa.Simulation.Core --- Sofa/Component/LinearSystem/CMakeLists.txt | 12 +++++-- .../component/linearsystem/MappingGraph.h | 31 +++++++++++++++++++ .../linearsystem/MatrixLinearSystem.inl | 2 +- .../linearsystem/MatrixProjectionMethod.h | 8 ++--- .../linearsystem/MatrixProjectionMethod.inl | 16 +++++----- ...sembleGlobalVectorFromLocalVectorVisitor.h | 1 + ...atchFromGlobalVectorToLocalVectorVisitor.h | 1 + .../LinearSystem/tests/CMakeLists.txt | 1 - Sofa/framework/Simulation/Core/CMakeLists.txt | 2 ++ .../src/sofa/simulation}/MappingGraph.cpp | 4 +-- .../Core/src/sofa/simulation}/MappingGraph.h | 6 ++-- .../Simulation/Core/test/CMakeLists.txt | 6 +++- .../Core/test}/MappingGraph_test.cpp | 16 +++++----- 13 files changed, 76 insertions(+), 30 deletions(-) create mode 100644 Sofa/Component/LinearSystem/compat/sofa/component/linearsystem/MappingGraph.h rename Sofa/{Component/LinearSystem/src/sofa/component/linearsystem => framework/Simulation/Core/src/sofa/simulation}/MappingGraph.cpp (99%) rename Sofa/{Component/LinearSystem/src/sofa/component/linearsystem => framework/Simulation/Core/src/sofa/simulation}/MappingGraph.h (98%) rename Sofa/{Component/LinearSystem/tests => framework/Simulation/Core/test}/MappingGraph_test.cpp (96%) diff --git a/Sofa/Component/LinearSystem/CMakeLists.txt b/Sofa/Component/LinearSystem/CMakeLists.txt index 3dd8a6f69f7..7047c0c120b 100644 --- a/Sofa/Component/LinearSystem/CMakeLists.txt +++ b/Sofa/Component/LinearSystem/CMakeLists.txt @@ -17,7 +17,6 @@ set(HEADER_FILES ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/LinearSystemData.h ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MappedMassMatrixObserver.h ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MappedMassMatrixObserver.inl - ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MappingGraph.h ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MatrixFreeSystem.h ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MatrixLinearSystem.h ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MatrixLinearSystem.inl @@ -46,7 +45,6 @@ set(SOURCE_FILES ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/ConstantSparsityPatternSystem.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/ConstantSparsityProjectionMethod.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MappedMassMatrixObserver.cpp - ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MappingGraph.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MatrixLinearSystem.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MatrixProjectionMethod.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/TypedMatrixLinearSystem.cpp @@ -55,11 +53,21 @@ set(SOURCE_FILES ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/visitors/DispatchFromGlobalVectorToLocalVectorVisitor.cpp ) +set(DEPRECATED_DIR "compat/sofa/component/linearsystem") +set(DEPRECATED_HEADER_FILES + ${DEPRECATED_DIR}/MappingGraph.h +) + sofa_find_package(Sofa.Simulation.Core REQUIRED) add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES}) target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core Sofa.LinearAlgebra) +target_include_directories(${PROJECT_NAME} PUBLIC + $ + $ +) + sofa_create_package_with_targets( PACKAGE_NAME ${PROJECT_NAME} PACKAGE_VERSION ${Sofa_VERSION} diff --git a/Sofa/Component/LinearSystem/compat/sofa/component/linearsystem/MappingGraph.h b/Sofa/Component/LinearSystem/compat/sofa/component/linearsystem/MappingGraph.h new file mode 100644 index 00000000000..6aa664670b2 --- /dev/null +++ b/Sofa/Component/LinearSystem/compat/sofa/component/linearsystem/MappingGraph.h @@ -0,0 +1,31 @@ +/****************************************************************************** +* 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 +SOFA_HEADER_DEPRECATED("v26.06", "v26.12", "sofa/simulation/MappingGraph.h") + +namespace sofa::component::linearsystem +{ +using MappingGraph = sofa::simulation::MappingGraph; +} diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixLinearSystem.inl b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixLinearSystem.inl index fd728161787..ca14f33bac1 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixLinearSystem.inl +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixLinearSystem.inl @@ -239,7 +239,7 @@ inline sofa::type::vector retrieveAssociat inline sofa::type::vector retrieveAssociatedMechanicalState(BaseMapping* component) { - type::vector mstates = component->getMechFrom(); + type::vector mstates = component->getMechFrom(); //remove duplicates: it may happen for MultiMappings std::sort( mstates.begin(), mstates.end() ); diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixProjectionMethod.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixProjectionMethod.h index bf91160bbad..77854b28e79 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixProjectionMethod.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixProjectionMethod.h @@ -61,13 +61,13 @@ class MatrixProjectionMethod : public BaseMatrixProjectionMethod linearalgebra::BaseMatrix* globalMatrix) override; /// Given a Mechanical State and its matrix, identifies the nodes affected by the matrix - std::vector identifyAffectedDoFs(BaseMechanicalState* mstate, TMatrix* crs); + std::vector identifyAffectedDoFs(core::behavior::BaseMechanicalState* mstate, TMatrix* crs); /** * Build the jacobian matrices of mappings from a mapped state to its top most parents (in the * sense of mappings) */ - MappingJacobians computeJacobiansFrom(BaseMechanicalState* mstate, const core::MechanicalParams* mparams, const MappingGraph& mappingGraph, TMatrix* crs); + simulation::MappingJacobians computeJacobiansFrom(core::behavior::BaseMechanicalState* mstate, const core::MechanicalParams* mparams, const MappingGraph& mappingGraph, TMatrix* crs); core::objectmodel::BaseContext* getSolveContext(); @@ -87,7 +87,7 @@ class MatrixProjectionMethod : public BaseMatrixProjectionMethod void addMappedMatrixToGlobalMatrixEigen( sofa::type::fixed_array mstatePair, TMatrix* mappedMatrix, - sofa::type::fixed_array, 2> jacobians, + sofa::type::fixed_array, 2> jacobians, const MappingGraph& mappingGraph, linearalgebra::BaseMatrix* globalMatrix); @@ -100,7 +100,7 @@ class MatrixProjectionMethod : public BaseMatrixProjectionMethod Data d_areJacobiansConstant; ///< True if mapping jacobians are considered constant over time. They are computed only the first time. - std::optional, 2>> m_mappingJacobians; + std::optional, 2>> m_mappingJacobians; }; #if !defined(SOFA_COMPONENT_LINEARSYSTEM_EIGENMATRIXMAPPING_CPP) diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixProjectionMethod.inl b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixProjectionMethod.inl index 075d3e48477..c42828c59f2 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixProjectionMethod.inl +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MatrixProjectionMethod.inl @@ -70,7 +70,7 @@ template void MatrixProjectionMethod::addMappedMatrixToGlobalMatrixEigen( sofa::type::fixed_array mstatePair, TMatrix* mappedMatrix, - sofa::type::fixed_array, 2> jacobians, + sofa::type::fixed_array, 2> jacobians, const MappingGraph& mappingGraph, linearalgebra::BaseMatrix* globalMatrix) { if (!mappedMatrix) @@ -206,14 +206,14 @@ void MatrixProjectionMethod::computeMatrixJacobians(const core::Mechani { if (!m_mappingJacobians.has_value() || !d_areJacobiansConstant.getValue()) { - const MappingJacobians J0 = computeJacobiansFrom( + const simulation::MappingJacobians J0 = computeJacobiansFrom( this->l_mechanicalStates[0], mparams, mappingGraph, matrixToProject); - const MappingJacobians J1 = + const simulation::MappingJacobians J1 = (this->l_mechanicalStates[0] == this->l_mechanicalStates[1]) ? J0 : computeJacobiansFrom(this->l_mechanicalStates[1], mparams, mappingGraph, matrixToProject); - m_mappingJacobians.emplace(sofa::type::fixed_array, 2>({J0, J1})); + m_mappingJacobians.emplace(sofa::type::fixed_array, 2>({J0, J1})); } } @@ -237,7 +237,7 @@ void MatrixProjectionMethod::projectMatrixToGlobalMatrix(const core::Me template std::vector MatrixProjectionMethod::identifyAffectedDoFs( - BaseMechanicalState* mstate, TMatrix* crs) + core::behavior::BaseMechanicalState* mstate, TMatrix* crs) { const auto blockSize = mstate->getMatrixBlockSize(); std::set setAffectedDoFs; @@ -262,13 +262,13 @@ std::vector MatrixProjectionMethod::identifyAffectedDoFs( } template -MappingJacobians MatrixProjectionMethod::computeJacobiansFrom( - BaseMechanicalState* mstate, const core::MechanicalParams* mparams, +simulation::MappingJacobians MatrixProjectionMethod::computeJacobiansFrom( + core::behavior::BaseMechanicalState* mstate, const core::MechanicalParams* mparams, const MappingGraph& mappingGraph, TMatrix* crs) { core::ConstraintParams cparams(*mparams); - MappingJacobians jacobians(*mstate); + simulation::MappingJacobians jacobians(*mstate); if (!mappingGraph.hasAnyMappingInput(mstate)) { diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/visitors/AssembleGlobalVectorFromLocalVectorVisitor.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/visitors/AssembleGlobalVectorFromLocalVectorVisitor.h index 7cd7de402c1..a57341bad2d 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/visitors/AssembleGlobalVectorFromLocalVectorVisitor.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/visitors/AssembleGlobalVectorFromLocalVectorVisitor.h @@ -21,6 +21,7 @@ ******************************************************************************/ #pragma once +#include #include #include diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/visitors/DispatchFromGlobalVectorToLocalVectorVisitor.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/visitors/DispatchFromGlobalVectorToLocalVectorVisitor.h index 80dbb4d9c6f..45f7650c37e 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/visitors/DispatchFromGlobalVectorToLocalVectorVisitor.h +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/visitors/DispatchFromGlobalVectorToLocalVectorVisitor.h @@ -21,6 +21,7 @@ ******************************************************************************/ #pragma once +#include #include #include diff --git a/Sofa/Component/LinearSystem/tests/CMakeLists.txt b/Sofa/Component/LinearSystem/tests/CMakeLists.txt index 3d2591e119c..96c8a4e9b5e 100644 --- a/Sofa/Component/LinearSystem/tests/CMakeLists.txt +++ b/Sofa/Component/LinearSystem/tests/CMakeLists.txt @@ -4,7 +4,6 @@ project(Sofa.Component.LinearSystem_test) set(SOURCE_FILES MatrixLinearSystem_test.cpp - MappingGraph_test.cpp ) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) diff --git a/Sofa/framework/Simulation/Core/CMakeLists.txt b/Sofa/framework/Simulation/Core/CMakeLists.txt index 8a07c0a14ce..e16319d1d4e 100644 --- a/Sofa/framework/Simulation/Core/CMakeLists.txt +++ b/Sofa/framework/Simulation/Core/CMakeLists.txt @@ -66,6 +66,7 @@ set(HEADER_FILES ${SRC_ROOT}/SceneCheck.h ${SRC_ROOT}/SceneCheckRegistry.h ${SRC_ROOT}/SceneCheckMainRegistry.h + ${SRC_ROOT}/MappingGraph.h ${SRC_ROOT}/events/BuildConstraintSystemEndEvent.h ${SRC_ROOT}/events/SimulationInitDoneEvent.h @@ -168,6 +169,7 @@ set(SOURCE_FILES ${SRC_ROOT}/InitVisitor.cpp ${SRC_ROOT}/IntegrateBeginEvent.cpp ${SRC_ROOT}/IntegrateEndEvent.cpp + ${SRC_ROOT}/MappingGraph.cpp ${SRC_ROOT}/MechanicalOperations.cpp ${SRC_ROOT}/MechanicalVPrintVisitor.cpp ${SRC_ROOT}/MechanicalVisitor.cpp diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MappingGraph.cpp b/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraph.cpp similarity index 99% rename from Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MappingGraph.cpp rename to Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraph.cpp index bfdf653c2db..b9735b9e37f 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MappingGraph.cpp +++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraph.cpp @@ -19,7 +19,7 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#include +#include #include #include @@ -30,7 +30,7 @@ #include #include -namespace sofa::component::linearsystem +namespace sofa::simulation { core::objectmodel::BaseContext* MappingGraph::getRootNode() const diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MappingGraph.h b/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraph.h similarity index 98% rename from Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MappingGraph.h rename to Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraph.h index 482b716e3ae..75d0aac75c0 100644 --- a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/MappingGraph.h +++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/MappingGraph.h @@ -20,12 +20,12 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once -#include +#include #include #include -namespace sofa::component::linearsystem +namespace sofa::simulation { using core::behavior::BaseMechanicalState; @@ -35,7 +35,7 @@ using core::behavior::BaseMechanicalState; * * Graph must be built with the build() function. */ -class SOFA_COMPONENT_LINEARSYSTEM_API MappingGraph +class SOFA_SIMULATION_CORE_API MappingGraph { public: diff --git a/Sofa/framework/Simulation/Core/test/CMakeLists.txt b/Sofa/framework/Simulation/Core/test/CMakeLists.txt index e440df28395..fa559b2d6d1 100644 --- a/Sofa/framework/Simulation/Core/test/CMakeLists.txt +++ b/Sofa/framework/Simulation/Core/test/CMakeLists.txt @@ -4,6 +4,7 @@ project(Sofa.Simulation.Core_test) set(SOURCE_FILES Colors_test.cpp + MappingGraph_test.cpp ParallelForEach_test.cpp RequiredPlugin_test.cpp SceneCheckRegistry_test.cpp @@ -16,6 +17,9 @@ set(SOURCE_FILES ) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) -target_link_libraries(${PROJECT_NAME} Sofa.Testing Sofa.Simulation.Core Sofa.LinearAlgebra.Testing) +target_link_libraries(${PROJECT_NAME} Sofa.Testing Sofa.Simulation.Core Sofa.LinearAlgebra.Testing + Sofa.Component.StateContainer + Sofa.Component.Mapping.Linear +) add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME}) diff --git a/Sofa/Component/LinearSystem/tests/MappingGraph_test.cpp b/Sofa/framework/Simulation/Core/test/MappingGraph_test.cpp similarity index 96% rename from Sofa/Component/LinearSystem/tests/MappingGraph_test.cpp rename to Sofa/framework/Simulation/Core/test/MappingGraph_test.cpp index 7de0406d2b5..0d6dc401b5b 100644 --- a/Sofa/Component/LinearSystem/tests/MappingGraph_test.cpp +++ b/Sofa/framework/Simulation/Core/test/MappingGraph_test.cpp @@ -20,7 +20,7 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #include -#include +#include #include #include #include @@ -30,7 +30,7 @@ TEST(MappingGraph, noBuild) { - const sofa::component::linearsystem::MappingGraph graph; + const sofa::simulation::MappingGraph graph; EXPECT_FALSE(graph.isBuilt()); EXPECT_EQ(graph.getRootNode(), nullptr); @@ -47,7 +47,7 @@ TEST(MappingGraph, noBuild) TEST(MappingGraph, nullRootNode) { - sofa::component::linearsystem::MappingGraph graph; + sofa::simulation::MappingGraph graph; graph.build(sofa::core::MechanicalParams::defaultInstance(), nullptr); EXPECT_FALSE(graph.isBuilt()); @@ -67,7 +67,7 @@ TEST(MappingGraph, emptyRootNode) { const sofa::simulation::Node::SPtr root = sofa::core::objectmodel::New(); - sofa::component::linearsystem::MappingGraph graph; + sofa::simulation::MappingGraph graph; graph.build(sofa::core::MechanicalParams::defaultInstance(), root.get()); EXPECT_TRUE(graph.isBuilt()); @@ -91,7 +91,7 @@ TEST(MappingGraph, oneMechanicalObject) root->addObject(mstate); mstate->resize(10); - sofa::component::linearsystem::MappingGraph graph; + sofa::simulation::MappingGraph graph; graph.build(sofa::core::MechanicalParams::defaultInstance(), root.get()); EXPECT_TRUE(graph.isBuilt()); @@ -117,7 +117,7 @@ TEST(MappingGraph, twoMechanicalObject) root->addObject(mstate2); mstate2->resize(2); - sofa::component::linearsystem::MappingGraph graph; + sofa::simulation::MappingGraph graph; graph.build(sofa::core::MechanicalParams::defaultInstance(), root.get()); EXPECT_TRUE(graph.isBuilt()); @@ -151,7 +151,7 @@ TEST(MappingGraph, oneMapping) mapping->setFrom(mstate1.get()); mapping->setTo(mstate2.get()); - sofa::component::linearsystem::MappingGraph graph; + sofa::simulation::MappingGraph graph; graph.build(sofa::core::MechanicalParams::defaultInstance(), root.get()); EXPECT_TRUE(graph.isBuilt()); @@ -213,7 +213,7 @@ TEST(MappingGraph, diamondMapping) {"input", "@/left/left @/right/right"}, {"output", "@bottom"} }); - sofa::component::linearsystem::MappingGraph graph; + sofa::simulation::MappingGraph graph; graph.build(sofa::core::MechanicalParams::defaultInstance(), root.get()); EXPECT_TRUE(graph.isBuilt()); From fdd47c01e36bfb6ac00b8d080c224c1e81b6446e Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 5 Mar 2026 14:36:10 +0100 Subject: [PATCH 2/2] Detect mapped masses --- .../odesolver/forward/EulerExplicitSolver.cpp | 83 +++++++++++++++++-- .../odesolver/forward/EulerExplicitSolver.h | 2 + 2 files changed, 76 insertions(+), 9 deletions(-) diff --git a/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.cpp b/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.cpp index 9f33d9df4f0..10ed828a093 100644 --- a/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.cpp +++ b/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.cpp @@ -20,14 +20,15 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #include -#include -#include -#include #include +#include #include +#include #include #include - +#include +#include +#include #include using sofa::simulation::mechanicalvisitor::MechanicalGetNonDiagonalMassesCountVisitor; @@ -84,11 +85,9 @@ void EulerExplicitSolver::solve(const core::ExecParams* params, addSeparateGravity(&mop, dt, vResult); computeForce(&mop, f); - sofa::Size nbNonDiagonalMasses = 0; - MechanicalGetNonDiagonalMassesCountVisitor(&mop.mparams, &nbNonDiagonalMasses).execute(this->getContext()); - - // Mass matrix is diagonal, solution can thus be found by computing acc = f/m - if(nbNonDiagonalMasses == 0.) + // The inverse of the mass matrix is trivial to compute. Otherwise, it requires the + // assembly of a linear system and a linear solver to compute its solution. + if(isMassMatrixTriviallyInvertible(params)) { // acc = M^-1 * f computeAcceleration(&mop, acc, f); @@ -343,4 +342,70 @@ void EulerExplicitSolver::solveSystem(core::MultiVecDerivId solution, core::Mult l_linearSolver->getLinearSystem()->dispatchSystemSolution(solution); } +class MappedMassVisitor final : public simulation::BaseMechanicalVisitor +{ +public: + MappedMassVisitor(const sofa::core::ExecParams* params, sofa::simulation::MappingGraph* mappingGraph) + : BaseMechanicalVisitor(params), m_mappingGraph(mappingGraph) + {} + + Result fwdMass(simulation::Node*, sofa::core::behavior::BaseMass* mass) override + { + if (mass && m_mappingGraph) + { + m_hasMappedMass |= m_mappingGraph->hasAnyMappingInput(mass); + } + return Result::RESULT_CONTINUE; + } + + [[nodiscard]] bool hasMappedMass() const { return m_hasMappedMass; } + +private: + sofa::simulation::MappingGraph* m_mappingGraph { nullptr }; + bool m_hasMappedMass { false }; +}; + +class AllOfMassesAreDiagonalVisitor final : public simulation::BaseMechanicalVisitor +{ +public: + using simulation::BaseMechanicalVisitor::BaseMechanicalVisitor; + Result fwdMass(simulation::Node*, sofa::core::behavior::BaseMass* mass) override + { + if (mass) + { + m_allMassesAreDiagonal &= mass->isDiagonal(); + } + return Result::RESULT_CONTINUE; + } + + [[nodiscard]] bool areAllMassesAreDiagonal() const { return m_allMassesAreDiagonal; } + +private: + bool m_allMassesAreDiagonal { true }; +}; + +bool EulerExplicitSolver::isMassMatrixTriviallyInvertible(const core::ExecParams* params) +{ + sofa::simulation::MappingGraph mappingGraph; + mappingGraph.build(params, this->getContext()); + + // To achieve a diagonal global mass matrix in this system: + // 1) Each individual mass matrix must itself be diagonal. + // 2) The mapped masses must also be mapped through a diagonal Jacobian matrix. + + // Diagonal Jacobian matrices are uncommon and challenging to detect, which means that even if + // we identify a mapped mass, we cannot guarantee the global mass matrix will remain diagonal. + // Moreover, computing the inverse of a mapped mass would require a complex API. Therefore, this + // case is not supported without assembling the global mass matrix. + MappedMassVisitor mappedMassVisitor(params, &mappingGraph); + mappedMassVisitor.execute(this->getContext()); + if (mappedMassVisitor.hasMappedMass()) + return false; + + // At this stage, we now that we don't have any mapped mass. We can check if they are all diagonal. + AllOfMassesAreDiagonalVisitor allOfMassesAreDiagonalVisitor(params); + allOfMassesAreDiagonalVisitor.execute(this->getContext()); + return allOfMassesAreDiagonalVisitor.areAllMassesAreDiagonal(); +} + } // namespace sofa::component::odesolver::forward diff --git a/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.h b/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.h index d5c9d1b83f2..72173e477a9 100644 --- a/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.h +++ b/Sofa/Component/ODESolver/Forward/src/sofa/component/odesolver/forward/EulerExplicitSolver.h @@ -123,6 +123,8 @@ class SOFA_COMPONENT_ODESOLVER_FORWARD_API EulerExplicitSolver : public sofa::co void assembleSystemMatrix(sofa::simulation::common::MechanicalOperations* mop) const; void solveSystem(core::MultiVecDerivId solution, core::MultiVecDerivId rhs) const; + + bool isMassMatrixTriviallyInvertible(const core::ExecParams* params); }; } // namespace sofa::component::odesolver::forward