From ed45a6dc14dfb8575a5f671bfc34267bdc7eef2e Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 16 Oct 2026 11:36:30 +0900 Subject: [PATCH] [Mass] Fix MeshMatrixMass totalMass when initialized from vertexMass/edgeMass - initFromVertexAndEdgeMass: each edge mass is applied twice (M(i,j) and M(j,i)), so totalMass = sum(vertexMass) + 2 * sum(edgeMass). When lumped, totalMass = coeff * sum(vertexMass), which is what addMDx applies. - initFromVertexMass (lumping): the given vertexMass are the lumped masses (behavior before #2305, and what totalMass and massDensity already assume). Store vertexMass / massLumpingCoeff, since the lumping coefficient is applied wherever the lumped mass is used. - Re-enable MeshMatrixMass_test.cpp and add tests checking that totalMass is the sum of the mass matrix entries (assembled and applied through addMDx) for both init paths, on tetra and hexa meshes. Behavior changes: - lumping=true with vertexMass only: the applied mass is now exactly the given vertexMass, instead of massLumpingCoeff (2 for triangles/quads, 2.5 for tetrahedra/hexahedra) times it. Free fall under gravity is unchanged, but objects are lighter relative to stiffness, constraints and contacts. d_vertexMass now holds vertexMass / massLumpingCoeff after init. totalMass and massDensity are unchanged. - vertexMass + edgeMass: the dynamics are unchanged, only the reported totalMass and massDensity change. The corrected massDensity also gives the right mass to elements added by topological changes. - No scene in the repository sets vertexMass or edgeMass on a MeshMatrixMass, so no regression reference is affected. --- .../sofa/component/mass/MeshMatrixMass.inl | 26 ++- Sofa/Component/Mass/tests/CMakeLists.txt | 1 + .../Mass/tests/MeshMatrixMass_test.cpp | 195 +++++++++++++++++- 3 files changed, 208 insertions(+), 14 deletions(-) diff --git a/Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl b/Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl index 79af714bbff..8bf6b6851e9 100644 --- a/Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl +++ b/Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl @@ -1597,6 +1597,7 @@ void MeshMatrixMass::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::initFromVertexMass() for(size_t i=0; i 0 ? vertexMass[i] / m_massLumpingCoeff : vertexMass[i]; } //Update all computed values @@ -1656,20 +1658,24 @@ void MeshMatrixMass::initFromVertexAndEdgeMass() const sofa::type::vector vertexMass = d_vertexMass.getValue(); const sofa::type::vector edgeMass = d_edgeMass.getValue(); - Real totalMassSave = 0.0; - for(size_t i=0; i > vertexMassInfo = d_vertexMass; helper::WriteAccessor > edgeMassInfo = d_edgeMass; //Compute volume = mass since massDensity = 1.0 diff --git a/Sofa/Component/Mass/tests/CMakeLists.txt b/Sofa/Component/Mass/tests/CMakeLists.txt index 31b2ab3ed26..2876b34cec8 100644 --- a/Sofa/Component/Mass/tests/CMakeLists.txt +++ b/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 ) diff --git a/Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp b/Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp index 01548d0f40d..3e2a51c6061 100644 --- a/Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp +++ b/Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp @@ -55,6 +55,10 @@ using sofa::simulation::Node ; #include using sofa::simulation::SceneLoaderXML ; +#include + +#include +#include #include 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 = + " " + " " + " " + " " + " "; + } + else + { + topology = + " " + " " + " " + " " + " " + " " + " " + " " + " "; + } + + const string scene = + " " + " " + " " + + topology + + " " + " " + " "; + + 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(); + } + + 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 dx; + dx.setValue(typename DataTypes::VecDeriv(nbVertices, one)); + core::objectmodel::Data 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(); }