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(); }