[SofaExplicitOdeSolver] Introduce visitor to know the number of non-diagonal mass matrices

SofaKernel/modules/SofaExplicitOdeSolver/src/SofaExplicitOdeSolver/EulerSolver.cpp

@@ -27,6 +27,9 @@
 #include <sofa/core/behavior/MultiMatrix.h>
 #include <sofa/helper/AdvancedTimer.h>
 
+#include <sofa/simulation/mechanicalvisitor/MechanicalGetNonDiagonalMassesCountVisitor.h>
+using sofa::simulation::mechanicalvisitor::MechanicalGetNonDiagonalMassesCountVisitor;
+
 //#define SOFA_NO_VMULTIOP
 
 namespace sofa::component::odesolver
@@ -42,7 +45,6 @@ int EulerExplicitSolverClass = core::RegisterObject("A simple explicit time inte
 
 EulerExplicitSolver::EulerExplicitSolver()
     : d_symplectic( initData( &d_symplectic, true, "symplectic", "If true, the velocities are updated before the positions and the method is symplectic (more robust). If false, the positions are updated before the velocities (standard Euler, less robust).") )
-    , d_optimizedForDiagonalMatrix(initData(&d_optimizedForDiagonalMatrix, true, "optimizedForDiagonalMatrix", "If true, solution to the system Ax=b can be directly found by computing x = f/m. Must be set to false if M is sparse."))
     , d_threadSafeVisitor(initData(&d_threadSafeVisitor, false, "threadSafeVisitor", "If true, do not use realloc and free visitors in fwdInteractionForceField."))
 {
 }
@@ -73,8 +75,11 @@ void EulerExplicitSolver::solve(const core::ExecParams* params,
     addSeparateGravity(&mop, dt, vResult);
     computeForce(&mop, f);
 
+    SReal nbNonDiagonalMasses = 0;
+    MechanicalGetNonDiagonalMassesCountVisitor(&mop.mparams, &nbNonDiagonalMasses).execute(this->getContext());
+
     // Mass matrix is diagonal, solution can thus be found by computing acc = f/m
-    if(d_optimizedForDiagonalMatrix.getValue())
+    if(nbNonDiagonalMasses == 0.)
     {
         // acc = M^-1 * f
         computeAcceleration(&mop, acc, f);
@@ -240,6 +245,20 @@ void EulerExplicitSolver::init()
     reinit();
 }
 
+void EulerExplicitSolver::parse(sofa::core::objectmodel::BaseObjectDescription* arg)
+{
+    Inherit1::parse(arg);
+
+    const char* val = arg->getAttribute("optimizedForDiagonalMatrix",nullptr) ;
+    if(val)
+    {
+        msg_deprecated() << "The attribute 'optimizedForDiagonalMatrix' is deprecated since SOFA 21.06." << msgendl
+                         << "This data was previously used to solve the system more efficiently in case the "
+                            "global mass matrix were diagonal." << msgendl
+                         << "Now, this property is detected automatically.";
+    }
+}
+
 void EulerExplicitSolver::addSeparateGravity(sofa::simulation::common::MechanicalOperations* mop, SReal dt, core::MultiVecDerivId v)
 {
     sofa::helper::ScopedAdvancedTimer timer("addSeparateGravity");

SofaKernel/modules/SofaExplicitOdeSolver/src/SofaExplicitOdeSolver/EulerSolver.h

@@ -81,7 +81,6 @@ class SOFA_SOFAEXPLICITODESOLVER_API EulerExplicitSolver : public sofa::core::be
     void solve(const core::ExecParams* params, SReal dt, sofa::core::MultiVecCoordId xResult, sofa::core::MultiVecDerivId vResult) override;
 
     Data<bool> d_symplectic; ///< If true, the velocities are updated before the positions and the method is symplectic (more robust). If false, the positions are updated before the velocities (standard Euler, less robust).
-    Data<bool> d_optimizedForDiagonalMatrix; ///< If M matrix is sparse (MeshMatrixMass), must be set to false (function addMDx() will compute the mass). Else, if true, solution to the system Ax=b can be directly found by computing x = f/m. The function accFromF() in the mass API will be used.
     Data<bool> d_threadSafeVisitor;
 
     /// Given an input derivative order (0 for position, 1 for velocity, 2 for acceleration),
@@ -95,6 +94,8 @@ class SOFA_SOFAEXPLICITODESOLVER_API EulerExplicitSolver : public sofa::core::be
 
     void init() override ;
 
+    void parse(sofa::core::objectmodel::BaseObjectDescription* arg) override;
+
 protected:
 
     /// Update state variable (new position and velocity) based on the computed acceleration

SofaKernel/modules/SofaSimulationCore/CMakeLists.txt

@@ -99,6 +99,7 @@ set(HEADER_FILES
     ${SRC_ROOT}/mechanicalvisitor/MechanicalGetDimensionVisitor.h
     ${SRC_ROOT}/mechanicalvisitor/MechanicalGetMatrixDimensionVisitor.h
     ${SRC_ROOT}/mechanicalvisitor/MechanicalGetMomentumVisitor.h
+    ${SRC_ROOT}/mechanicalvisitor/MechanicalGetNonDiagonalMassesCountVisitor.h
     ${SRC_ROOT}/mechanicalvisitor/MechanicalIntegrateConstraintVisitor.h
     ${SRC_ROOT}/mechanicalvisitor/MechanicalIntegrationVisitor.h
     ${SRC_ROOT}/mechanicalvisitor/MechanicalMultiVectorFromBaseVectorVisitor.h
@@ -220,6 +221,7 @@ set(SOURCE_FILES
     ${SRC_ROOT}/mechanicalvisitor/MechanicalGetDimensionVisitor.cpp
     ${SRC_ROOT}/mechanicalvisitor/MechanicalGetMatrixDimensionVisitor.cpp
     ${SRC_ROOT}/mechanicalvisitor/MechanicalGetMomentumVisitor.cpp
+    ${SRC_ROOT}/mechanicalvisitor/MechanicalGetNonDiagonalMassesCountVisitor.cpp
     ${SRC_ROOT}/mechanicalvisitor/MechanicalIntegrateConstraintVisitor.cpp
     ${SRC_ROOT}/mechanicalvisitor/MechanicalIntegrationVisitor.cpp
     ${SRC_ROOT}/mechanicalvisitor/MechanicalMultiVectorFromBaseVectorVisitor.cpp

SofaKernel/modules/SofaSimulationCore/src/sofa/simulation/mechanicalvisitor/MechanicalGetNonDiagonalMassesCountVisitor.cpp

@@ -0,0 +1,35 @@
+/******************************************************************************
+*                 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/mechanicalvisitor/MechanicalGetNonDiagonalMassesCountVisitor.h>
+#include <sofa/core/behavior/BaseMass.h>
+
+namespace sofa::simulation::mechanicalvisitor
+{
+
+Visitor::Result MechanicalGetNonDiagonalMassesCountVisitor::fwdMass(VisitorContext* ctx, core::behavior::BaseMass* mass)
+{
+    *ctx->nodeData += !mass->isDiagonal();
+    return RESULT_CONTINUE;
+}
+
+}

SofaKernel/modules/SofaSimulationCore/src/sofa/simulation/mechanicalvisitor/MechanicalGetNonDiagonalMassesCountVisitor.h

@@ -0,0 +1,52 @@
+/******************************************************************************
+*                 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/MechanicalVisitor.h>
+
+namespace sofa::simulation::mechanicalvisitor
+{
+
+/** Count the number of masses which are not diagonal */
+class SOFA_SIMULATION_CORE_API MechanicalGetNonDiagonalMassesCountVisitor : public MechanicalVisitor
+{
+public:
+    MechanicalGetNonDiagonalMassesCountVisitor(const sofa::core::MechanicalParams* mparams, SReal* result)
+            : MechanicalVisitor(mparams)
+    {
+        rootData = result;
+    }
+
+    Result fwdMass(VisitorContext* ctx, sofa::core::behavior::BaseMass* mass) override;
+
+    /// Return a class name for this visitor
+    /// Only used for debugging / profiling purposes
+    const char* getClassName() const override { return "MechanicalGetNonDiagonalMassesCountVisitor";}
+
+    bool writeNodeData() const override
+    {
+        return true;
+    }
+
+};
+
+}