[Mass] Fix MeshMatrixMass totalMass when initialized from vertexMass/edgemass

Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl

@@ -1597,6 +1597,7 @@ void MeshMatrixMass<DataTypes, GeometricalTypes>::initFromVertexMass()
 {
     msg_info() << "vertexMass information is used";
 
+    // The given vertexMass are the lumped masses, i.e. the diagonal of the applied mass matrix
     const auto vertexMass = d_vertexMass.getValue();
     Real totalMassSave = 0.0;
     for(size_t i=0; i<vertexMass.size(); i++)
@@ -1614,7 +1615,8 @@ void MeshMatrixMass<DataTypes, GeometricalTypes>::initFromVertexMass()
     for(size_t i=0; i<vertexMassInfo.size(); i++)
     {
         volume += vertexMassInfo[i]*m_massLumpingCoeff;
-        vertexMassInfo[i] = vertexMass[i];
+        // the stored vertexMass is multiplied by m_massLumpingCoeff wherever the lumped mass is used
+        vertexMassInfo[i] = m_massLumpingCoeff > 0 ? vertexMass[i] / m_massLumpingCoeff : vertexMass[i];
     }
 
     //Update all computed values
@@ -1656,20 +1658,24 @@ void MeshMatrixMass<DataTypes, GeometricalTypes>::initFromVertexAndEdgeMass()
 
     const sofa::type::vector<MassType> vertexMass = d_vertexMass.getValue();
     const sofa::type::vector<MassType> edgeMass = d_edgeMass.getValue();
-    Real totalMassSave = 0.0;
-    for(size_t i=0; i<vertexMass.size(); i++)
-    {
-        totalMassSave += vertexMass[i];
-    }
-    for(size_t i=0; i<edgeMass.size(); i++)
-    {
-        totalMassSave += edgeMass[i];
-    }
+
     //Compute the volume
     setMassDensity(1.0);
 
     computeMass();
 
+    // Sum of the entries of the applied mass matrix
+    Real totalMassSave = std::accumulate(vertexMass.begin(), vertexMass.end(), Real(0));
+    if (!isLumped())
+    {
+        // Add mass per edges if not lumped, *2 as it is added to both edge vertices
+        totalMassSave += std::accumulate(edgeMass.begin(), edgeMass.end(), Real(0)) * 2.0;
+    }
+    else
+    {
+        totalMassSave *= m_massLumpingCoeff;
+    }
+
     helper::WriteAccessor<Data<MassVector> > vertexMassInfo = d_vertexMass;
     helper::WriteAccessor<Data<MassVector> > edgeMassInfo = d_edgeMass;
     //Compute volume = mass since massDensity = 1.0

Sofa/Component/Mass/tests/CMakeLists.txt

@@ -8,6 +8,7 @@ set(SOURCE_FILES
     MassTestCreation[FEMMass].cpp
     MassTestCreation[MeshMatrixMass].cpp
     MassTestCreation[UniformMass].cpp
+    MeshMatrixMass_test.cpp
     UniformMass_test.cpp
 )
 

Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp

@@ -55,6 +55,10 @@ using sofa::simulation::Node ;
 #include <sofa/simulation/common/SceneLoaderXML.h>
 using sofa::simulation::SceneLoaderXML ;
 
+#include <sofa/core/MechanicalParams.h>
+
+#include <numeric>
+#include <sstream>
 #include <string>
 using std::string ;
 
@@ -170,6 +174,162 @@ class MeshMatrixMass_test : public BaseTest
     }
 
 
+    /// Load a scene made of a 2x2x2 box (hexahedra or tetrahedra) with a MeshMatrixMass using the given attributes.
+    /// The scene replaces the root of the test, so that it is unloaded at tear down.
+    TheMeshMatrixMass* loadBoxScene(const bool hexa, const string& massAttributes)
+    {
+        string topology;
+        if (hexa)
+        {
+            topology =
+                "    <RegularGridTopology name='grid' n='3 3 3' min='0 0 0' max='2 2 2' p0='0 0 0' />               "
+                "    <Node name='Hexa' >                                                                            "
+                "        <MechanicalObject position='@../grid.position' />                                          "
+                "        <HexahedronSetTopologyContainer name='Container' src='@../grid' />                         "
+                "        <HexahedronSetGeometryAlgorithms template='Vec3d' name='GeomAlgo' />                       ";
+        }
+        else
+        {
+            topology =
+                "    <RequiredPlugin pluginName='Sofa.Component.Topology.Mapping'/>                                 "
+                "    <MechanicalObject />                                                                           "
+                "    <RegularGridTopology name='grid' n='2 2 2' min='0 0 0' max='2 2 2' p0='0 0 0' />               "
+                "    <Node name='Tetra' >                                                                           "
+                "        <MechanicalObject src='@../grid'/>                                                         "
+                "        <TetrahedronSetTopologyContainer name='Container' />                                       "
+                "        <TetrahedronSetTopologyModifier name='Modifier' />                                         "
+                "        <TetrahedronSetGeometryAlgorithms template='Vec3d' name='GeomAlgo' />                      "
+                "        <Hexa2TetraTopologicalMapping input='@../grid' output='@Container' />                      ";
+        }
+
+        const string scene =
+                "<?xml version='1.0'?>                                                                              "
+                "<Node  name='Root' gravity='0 0 0' time='0' animate='0'   >                                        "
+                "    <DefaultAnimationLoop />                                                                       "
+                + topology +
+                "        <MeshMatrixMass name='m_mass' " + massAttributes + " />                                   "
+                "    </Node>                                                                                        "
+                "</Node>                                                                                            ";
+
+        if (root != nullptr)
+            sofa::simulation::node::unload(root);
+
+        root = SceneLoaderXML::loadFromMemory("loadBoxScene", scene.c_str());
+        if (root == nullptr)
+            return nullptr;
+
+        sofa::simulation::node::initRoot(root.get());
+        return root->getTreeObject<TheMeshMatrixMass>();
+    }
+
+    static string toString(const VecMass& values)
+    {
+        std::ostringstream oss;
+        oss.precision(17);
+        oss << values;
+        return oss.str();
+    }
+
+    /// Row sums of the mass matrix actually applied by the component: M * 1, computed with addMDx
+    static VecMass computeAppliedRowSums(TheMeshMatrixMass* mass, const std::size_t nbVertices)
+    {
+        typename DataTypes::Deriv one;
+        one[0] = 1;
+
+        core::objectmodel::Data<typename DataTypes::VecDeriv> dx;
+        dx.setValue(typename DataTypes::VecDeriv(nbVertices, one));
+        core::objectmodel::Data<typename DataTypes::VecDeriv> res;
+        res.setValue(typename DataTypes::VecDeriv(nbVertices));
+
+        mass->addMDx(core::mechanicalparams::defaultInstance(), res, dx, 1.0);
+
+        VecMass rowSums(nbVertices);
+        for (std::size_t i = 0; i < nbVertices; ++i)
+            rowSums[i] = res.getValue()[i][0];
+        return rowSums;
+    }
+
+    /// Check that the total mass is the sum of all the entries of the mass matrix, both the matrix
+    /// assembled from vertexMass (diagonal) and edgeMass (both off-diagonal entries), and the matrix
+    /// actually applied by addMDx (which uses the lumped masses if lumping is enabled).
+    /// The box has a volume of 8, so the density must be totalMass / 8.
+    void checkTotalMassMatchesMassMatrix(TheMeshMatrixMass* mass)
+    {
+        const auto& vertexMass = mass->d_vertexMass.getValue();
+        const auto& edgeMass = mass->d_edgeMass.getValue();
+        const MassType totalMass = mass->getTotalMass();
+        const MassType tolerance = 1e-10 * totalMass;
+
+        if (!mass->isLumped())
+        {
+            const MassType assembledSum = std::accumulate(vertexMass.begin(), vertexMass.end(), MassType(0))
+                + 2 * std::accumulate(edgeMass.begin(), edgeMass.end(), MassType(0));
+            EXPECT_NEAR(totalMass, assembledSum, tolerance);
+        }
+
+        const VecMass rowSums = computeAppliedRowSums(mass, vertexMass.size());
+        EXPECT_NEAR(totalMass, std::accumulate(rowSums.begin(), rowSums.end(), MassType(0)), tolerance);
+
+        EXPECT_NEAR(mass->getMassDensity()[0], totalMass / 8, tolerance);
+    }
+
+    /// With lumping and vertexMass only, the given vertexMass are the lumped masses:
+    /// they are the diagonal of the applied mass matrix and sum to the total mass.
+    void check_VertexMass_Lumping_TotalMassMatchesMassMatrix(const bool hexa)
+    {
+        TheMeshMatrixMass* mass = loadBoxScene(hexa, "massDensity='1.0' lumping='1'");
+        ASSERT_NE(mass, nullptr);
+        const std::size_t nbVertices = mass->d_vertexMass.getValue().size();
+
+        VecMass vertexMass(nbVertices);
+        for (std::size_t i = 0; i < nbVertices; ++i)
+            vertexMass[i] = MassType(1.0 + 0.1 * i);
+
+        mass = loadBoxScene(hexa, "lumping='1' vertexMass='" + toString(vertexMass) + "'");
+        ASSERT_NE(mass, nullptr);
+        ASSERT_TRUE(mass->isLumped());
+
+        const MassType expectedTotalMass = std::accumulate(vertexMass.begin(), vertexMass.end(), MassType(0));
+        EXPECT_NEAR(mass->getTotalMass(), expectedTotalMass, 1e-10 * expectedTotalMass);
+
+        const VecMass rowSums = computeAppliedRowSums(mass, nbVertices);
+        for (std::size_t i = 0; i < nbVertices; ++i)
+            EXPECT_NEAR(rowSums[i], vertexMass[i], 1e-10 * vertexMass[i]);
+
+        checkTotalMassMatchesMassMatrix(mass);
+    }
+
+    /// With vertexMass and edgeMass, each edge mass is applied twice: M(i,j) and M(j,i)
+    void check_VertexAndEdgeMass_TotalMassMatchesMassMatrix(const bool hexa, const bool lumped)
+    {
+        TheMeshMatrixMass* mass = loadBoxScene(hexa, "massDensity='1.0'");
+        ASSERT_NE(mass, nullptr);
+        const std::size_t nbVertices = mass->d_vertexMass.getValue().size();
+        const std::size_t nbEdges = mass->d_edgeMass.getValue().size();
+
+        VecMass vertexMass(nbVertices);
+        for (std::size_t i = 0; i < nbVertices; ++i)
+            vertexMass[i] = MassType(1.0 + 0.1 * i);
+        VecMass edgeMass(nbEdges);
+        for (std::size_t i = 0; i < nbEdges; ++i)
+            edgeMass[i] = MassType(0.05 + 0.01 * i);
+
+        mass = loadBoxScene(hexa, string("lumping='") + (lumped ? "1" : "0") + "'"
+            + " vertexMass='" + toString(vertexMass) + "' edgeMass='" + toString(edgeMass) + "'");
+        ASSERT_NE(mass, nullptr);
+        ASSERT_EQ(mass->isLumped(), lumped);
+
+        if (!lumped)
+        {
+            const MassType expectedTotalMass = std::accumulate(vertexMass.begin(), vertexMass.end(), MassType(0))
+                + 2 * std::accumulate(edgeMass.begin(), edgeMass.end(), MassType(0));
+            EXPECT_NEAR(mass->getTotalMass(), expectedTotalMass, 1e-10 * expectedTotalMass);
+        }
+
+        checkTotalMassMatchesMassMatrix(mass);
+    }
+
+
     //---------------------------------------------------------------
     // HEXA topology
     //---------------------------------------------------------------
@@ -330,10 +490,11 @@ class MeshMatrixMass_test : public BaseTest
             EXPECT_EQ(mass->isLumped(), true);
             EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), expectedTotalMass);
             EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], expectedDensity);
-            EXPECT_EQ(mass->getVertexMass()[0], 1.0);
 
-            EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[0], 1.0_sreal);
-            EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[1], 1.0_sreal);
+            // the given vertexMass are the lumped masses applied on each vertex
+            const VecMass lumpedMasses = computeAppliedRowSums(mass, 27);
+            EXPECT_FLOATINGPOINT_EQ(lumpedMasses[0], 1.0_sreal);
+            EXPECT_FLOATINGPOINT_EQ(lumpedMasses[1], 1.0_sreal);
             EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[0], 0.0_sreal);
             EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[2], 0.0_sreal);
         }
@@ -827,7 +988,8 @@ class MeshMatrixMass_test : public BaseTest
             EXPECT_EQ( mass->getMassCount(), 8 );
             EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), 8.0_sreal);
             EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], 1.0_sreal);
-            EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 1.0_sreal);
+            // the given vertexMass are the lumped masses applied on each vertex
+            EXPECT_FLOATINGPOINT_EQ(computeAppliedRowSums(mass, 8)[0], 1.0_sreal);
         }
 
         return ;
@@ -2113,6 +2275,31 @@ TEST_F(MeshMatrixMass3_test, check_VertexMass_Lumping_Initialization_Tetra){
 }
 
 
+TEST_F(MeshMatrixMass3_test, check_VertexMass_Lumping_TotalMassMatchesMassMatrix_Hexa) {
+    check_VertexMass_Lumping_TotalMassMatchesMassMatrix(true);
+}
+
+TEST_F(MeshMatrixMass3_test, check_VertexMass_Lumping_TotalMassMatchesMassMatrix_Tetra) {
+    check_VertexMass_Lumping_TotalMassMatchesMassMatrix(false);
+}
+
+TEST_F(MeshMatrixMass3_test, check_VertexAndEdgeMass_TotalMassMatchesMassMatrix_Hexa_sparse) {
+    check_VertexAndEdgeMass_TotalMassMatchesMassMatrix(true, false);
+}
+
+TEST_F(MeshMatrixMass3_test, check_VertexAndEdgeMass_TotalMassMatchesMassMatrix_Tetra_sparse) {
+    check_VertexAndEdgeMass_TotalMassMatchesMassMatrix(false, false);
+}
+
+TEST_F(MeshMatrixMass3_test, check_VertexAndEdgeMass_TotalMassMatchesMassMatrix_Hexa_lumped) {
+    check_VertexAndEdgeMass_TotalMassMatchesMassMatrix(true, true);
+}
+
+TEST_F(MeshMatrixMass3_test, check_VertexAndEdgeMass_TotalMassMatchesMassMatrix_Tetra_lumped) {
+    check_VertexAndEdgeMass_TotalMassMatchesMassMatrix(false, true);
+}
+
+
 TEST_F(MeshMatrixMass3_test, check_DoubleDeclaration_TotalMassAndMassDensity_Tetra){
     check_DoubleDeclaration_TotalMassAndMassDensity_Tetra();
 }