From d9e46a56bc792b8aa39914b2cd2b84774a24a7ab Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 16 Oct 2026 11:19:02 +0900 Subject: [PATCH] [Mass] Make MeshMatrixMass positive definite on hexahedra and quads The non-lumped mass matrix only stores vertex and edge masses. On hexahedra it used rho*V * (1/20 per vertex, 1/40 per edge), whose element eigenvalues are rho*V/40 * {5, 3, 1, -1}: the checkerboard mode has negative mass, so any hexahedral MeshMatrixMass was indefinite (e.g. 42 negative eigenvalues out of 576 in TetrahedronHyperelasticityFEMForceField_invertedTets.scn). On quads, rho*A * (1/8, 1/16) was singular (eigenvalues rho*A/16 * {4, 2, 2, 0}). Start from the consistent mass matrices (trilinear hexa: rho*V/216 * (8, 4, 2, 1), bilinear quad: rho*A/36 * (4, 2, 1)) and add the couplings that cannot be stored (face and body diagonals, quad diagonals) to the vertex mass: - hexa: rho*V/216 * (15 vertex, 4 edge), eigenvalues rho*V/216 * {27, 19, 11, 3} - quad: rho*A/36 * (5 vertex, 2 edge), eigenvalues rho*A/36 * {9, 5, 5, 1} Total mass and row sums (the lumped masses, also used for the weight) are unchanged; the lumping coefficients are updated accordingly (27/15 and 9/5). Re-enable MeshMatrixMass_test.cpp, which was dropped from the CMake list in sums to totalMass and has the lumped masses as row sums, including after topological changes. --- .../sofa/component/mass/MeshMatrixMass.inl | 33 ++- Sofa/Component/Mass/tests/CMakeLists.txt | 1 + .../Mass/tests/MeshMatrixMass_test.cpp | 206 +++++++++++++++--- 3 files changed, 202 insertions(+), 38 deletions(-) diff --git a/Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl b/Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl index eeb7a34a52a..035abb77d32 100644 --- a/Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl +++ b/Sofa/Component/Mass/src/sofa/component/mass/MeshMatrixMass.inl @@ -323,6 +323,12 @@ void MeshMatrixMass::applyEdgeMassTriangleDestructi // --------------------------------------------------- //{ +// The consistent mass matrix of a bilinear quad is rho*A/36 * (4 on the vertices, 2 on the edges, +// 1 between opposite vertices). Only vertex and edge masses are stored, so the coupling between +// opposite vertices is added to the vertex mass: rho*A/36 * (5 on the vertices, 2 on the edges). +// This keeps the total mass and the row sums (rho*A/4 per vertex, which is the lumped mass) and +// gives a positive definite element matrix (eigenvalues rho*A/36 * {9, 5, 5, 1}). + /// Creation function for mass stored on vertices template template = 2, int > > @@ -363,7 +369,7 @@ void MeshMatrixMass::applyVertexMassQuadCreation(co const auto& pos3 = GeometricalTypes::getCPos(positions[q[3]]); const auto quadArea = sofa::geometry::Quad::area(pos0, pos1, pos2, pos3); - mass = (densityM[quadAdded[i]] * quadArea) / (typename DataTypes::Real(8.0)); + mass = (densityM[quadAdded[i]] * quadArea) * (typename DataTypes::Real(5.0) / typename DataTypes::Real(36.0)); // Adding mass for (unsigned int j = 0; j < 4; ++j) @@ -423,7 +429,7 @@ void MeshMatrixMass::applyEdgeMassQuadCreation(cons const auto& pos3 = GeometricalTypes::getCPos(positions[q[3]]); const auto quadArea = sofa::geometry::Quad::area(pos0, pos1, pos2, pos3); - mass = (densityM[quadAdded[i]] * quadArea) / (typename DataTypes::Real(16.0)); + mass = (densityM[quadAdded[i]] * quadArea) * (typename DataTypes::Real(2.0) / typename DataTypes::Real(36.0)); // Adding mass edges of concerne quad for (unsigned int j = 0; j < 4; ++j) @@ -467,7 +473,7 @@ void MeshMatrixMass::applyVertexMassQuadDestruction const auto& pos3 = GeometricalTypes::getCPos(positions[q[3]]); const auto quadArea = sofa::geometry::Quad::area(pos0, pos1, pos2, pos3); - mass = (densityM[quadRemoved[i]] * quadArea) / (typename DataTypes::Real(8.0)); + mass = (densityM[quadRemoved[i]] * quadArea) * (typename DataTypes::Real(5.0) / typename DataTypes::Real(36.0)); // Removing mass for (unsigned int j = 0; j < 4; ++j) @@ -516,7 +522,7 @@ void MeshMatrixMass::applyEdgeMassQuadDestruction(c const auto& pos3 = GeometricalTypes::getCPos(positions[q[3]]); const auto quadArea = sofa::geometry::Quad::area(pos0, pos1, pos2, pos3); - mass = (densityM[quadRemoved[i]] * quadArea) / (typename DataTypes::Real(16.0)); + mass = (densityM[quadRemoved[i]] * quadArea) * (typename DataTypes::Real(2.0) / typename DataTypes::Real(36.0)); // Removing mass edges of concerne quad for (unsigned int j = 0; j < 4; ++j) @@ -757,6 +763,13 @@ void MeshMatrixMass::applyEdgeMassTetrahedronDestru // --------------------------------------------------------- //{ +// The consistent mass matrix of a trilinear hexahedron is rho*V/216 * (8 on the vertices, 4 on the +// edges, 2 on the face diagonals, 1 on the body diagonals). Only vertex and edge masses are stored, +// so the face and body diagonal couplings are added to the vertex mass: +// rho*V/216 * (15 on the vertices, 4 on the edges). This keeps the total mass and the row sums +// (rho*V/8 per vertex, which is the lumped mass) and gives a positive definite element matrix +// (eigenvalues rho*V/216 * {27, 19, 11, 3}). + /// Creation function for mass stored on vertices template template = 3, int > > @@ -801,7 +814,7 @@ void MeshMatrixMass::applyVertexMassHexahedronCreat const auto& rpos7 = GeometricalTypes::getCPos(positions[h[7]]); const auto hexaVolume = sofa::geometry::Hexahedron::volume(rpos0, rpos1, rpos2, rpos3, rpos4, rpos5, rpos6, rpos7); - mass = (densityM[hexahedronAdded[i]] * hexaVolume) / (typename DataTypes::Real(20.0)); + mass = (densityM[hexahedronAdded[i]] * hexaVolume) * (typename DataTypes::Real(15.0) / typename DataTypes::Real(216.0)); // Adding mass for (unsigned int j = 0; j < 8; ++j) @@ -865,7 +878,7 @@ void MeshMatrixMass::applyEdgeMassHexahedronCreatio const auto& rpos7 = GeometricalTypes::getCPos(positions[h[7]]); const auto hexaVolume = sofa::geometry::Hexahedron::volume(rpos0, rpos1, rpos2, rpos3, rpos4, rpos5, rpos6, rpos7); - mass = (densityM[hexahedronAdded[i]] * hexaVolume) / (typename DataTypes::Real(40.0)); + mass = (densityM[hexahedronAdded[i]] * hexaVolume) * (typename DataTypes::Real(4.0) / typename DataTypes::Real(216.0)); // Adding mass edges of concerne triangle for (unsigned int j = 0; j < 12; ++j) @@ -913,7 +926,7 @@ void MeshMatrixMass::applyVertexMassHexahedronDestr const auto& rpos7 = GeometricalTypes::getCPos(positions[h[7]]); const auto hexaVolume = sofa::geometry::Hexahedron::volume(rpos0, rpos1, rpos2, rpos3, rpos4, rpos5, rpos6, rpos7); - mass = (densityM[hexahedronRemoved[i]] * hexaVolume) / (typename DataTypes::Real(20.0)); + mass = (densityM[hexahedronRemoved[i]] * hexaVolume) * (typename DataTypes::Real(15.0) / typename DataTypes::Real(216.0)); // Removing mass for (unsigned int j = 0; j < 8; ++j) @@ -966,7 +979,7 @@ void MeshMatrixMass::applyEdgeMassHexahedronDestruc const auto& rpos7 = GeometricalTypes::getCPos(positions[h[7]]); const auto hexaVolume = sofa::geometry::Hexahedron::volume(rpos0, rpos1, rpos2, rpos3, rpos4, rpos5, rpos6, rpos7); - mass = (densityM[hexahedronRemoved[i]] * hexaVolume) / (typename DataTypes::Real(40.0)); + mass = (densityM[hexahedronRemoved[i]] * hexaVolume) * (typename DataTypes::Real(4.0) / typename DataTypes::Real(216.0)); // Removing mass edges of concerne triangle for (unsigned int j = 0; j < 12; ++j) @@ -1425,7 +1438,7 @@ void MeshMatrixMass::computeMass() for (Index i = 0; i < n; ++i) quadsAdded.push_back(i); - m_massLumpingCoeff = 2.0; + m_massLumpingCoeff = 9.0 / 5.0; if (!isLumped()) { applyEdgeMassQuadCreation(quadsAdded, l_topology->getQuads(), emptyAncestors, emptyCoefficients); @@ -1463,7 +1476,7 @@ void MeshMatrixMass::computeMass() for (Index i = 0; i < n; ++i) hexahedraAdded.push_back(i); - m_massLumpingCoeff = 2.5; + m_massLumpingCoeff = 27.0 / 15.0; if (!isLumped()) { applyEdgeMassHexahedronCreation(hexahedraAdded, l_topology->getHexahedra(), emptyAncestors, emptyCoefficients); diff --git a/Sofa/Component/Mass/tests/CMakeLists.txt b/Sofa/Component/Mass/tests/CMakeLists.txt index 31b2ab3ed26..9d03a6a8127 100644 --- a/Sofa/Component/Mass/tests/CMakeLists.txt +++ b/Sofa/Component/Mass/tests/CMakeLists.txt @@ -5,6 +5,7 @@ project(Sofa.Component.Mass_test) set(SOURCE_FILES DiagonalMass_test.cpp FEMMass_test.cpp + MeshMatrixMass_test.cpp MassTestCreation[FEMMass].cpp MassTestCreation[MeshMatrixMass].cpp MassTestCreation[UniformMass].cpp diff --git a/Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp b/Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp index 01548d0f40d..7f8566009eb 100644 --- a/Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp +++ b/Sofa/Component/Mass/tests/MeshMatrixMass_test.cpp @@ -58,6 +58,8 @@ using sofa::simulation::SceneLoaderXML ; #include using std::string ; +#include + #include using sofa::testing::BaseTest; @@ -90,6 +92,8 @@ class MeshMatrixMass_test : public BaseTest typedef typename DataTypes::Coord Coord; typedef typename DataTypes::VecCoord VecCoord; typedef typename DataTypes::Real Real; + typedef typename DataTypes::Deriv Deriv; + typedef typename DataTypes::VecDeriv VecDeriv; typedef typename type::vector VecMass; typedef MeshMatrixMass TheMeshMatrixMass ; @@ -204,8 +208,8 @@ class MeshMatrixMass_test : public BaseTest static const MassType volume = 8.0; static const MassType expectedTotalMass = 1.0f; static const MassType expectedDensity = expectedTotalMass / volume; - static const VecMass expectedVMass(8, (MassType)(expectedDensity * volume * 1 / 20)); - static const VecMass expectedEMass(18, (MassType)(expectedDensity * volume * 1 / 40)); + static const VecMass expectedVMass(8, (MassType)(expectedDensity * volume * 15 / 216)); + static const VecMass expectedEMass(18, (MassType)(expectedDensity * volume * 4 / 216)); check(expectedTotalMass, expectedDensity, expectedVMass, expectedEMass); } @@ -245,10 +249,10 @@ class MeshMatrixMass_test : public BaseTest EXPECT_FLOATINGPOINT_EQ( mass->getTotalMass(), expectedTotalMass); EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], expectedDensity); - EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[0], (MassType)(expectedDensity * volumeElem * 1 / 20)); - EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[1], (MassType)(expectedDensity * volumeElem * 1 / 20 * 2)); // vertex shared by 2 hexa - EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[0], (MassType)(expectedDensity * volumeElem * 1 / 40 * 2)); // Edge shared by 2 hexa - EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[2], (MassType)(expectedDensity * volumeElem * 1 / 40)); + EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[0], (MassType)(expectedDensity * volumeElem * 15 / 216)); + EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[1], (MassType)(expectedDensity * volumeElem * 15 / 216 * 2)); // vertex shared by 2 hexa + EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[0], (MassType)(expectedDensity * volumeElem * 4 / 216 * 2)); // Edge shared by 2 hexa + EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[2], (MassType)(expectedDensity * volumeElem * 4 / 216)); } return ; @@ -287,10 +291,10 @@ class MeshMatrixMass_test : public BaseTest EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), expectedTotalMass); EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], expectedDensity); - EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[0], (MassType)(expectedDensity * volumeElem * 1 / 20)); - EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[1], (MassType)(expectedDensity * volumeElem * 1 / 20 * 2)); // vertex shared by 2 hexa - EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[0], (MassType)(expectedDensity * volumeElem * 1 / 40 * 2)); // Edge shared by 2 hexa - EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[2], (MassType)(expectedDensity * volumeElem * 1 / 40)); + EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[0], (MassType)(expectedDensity * volumeElem * 15 / 216)); + EXPECT_FLOATINGPOINT_EQ(mass->d_vertexMass.getValue()[1], (MassType)(expectedDensity * volumeElem * 15 / 216 * 2)); // vertex shared by 2 hexa + EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[0], (MassType)(expectedDensity * volumeElem * 4 / 216 * 2)); // Edge shared by 2 hexa + EXPECT_FLOATINGPOINT_EQ(mass->d_edgeMass.getValue()[2], (MassType)(expectedDensity * volumeElem * 4 / 216)); } return ; @@ -371,7 +375,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ( mass->getTotalMass(), 2.0_sreal); EXPECT_FLOATINGPOINT_EQ( mass->getMassDensity()[0], SReal(mass->getTotalMass() / 8.0)); // 8 hexa - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(mass->getMassDensity()[0] / 20)); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(mass->getMassDensity()[0] * 15 / 216)); } return ; @@ -404,7 +408,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ( mass->getTotalMass(), 2.0_sreal); EXPECT_FLOATINGPOINT_EQ( mass->getMassDensity()[0], SReal(mass->getTotalMass() / 8.0)); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(mass->getMassDensity()[0] / 20)); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(mass->getMassDensity()[0] * 15 / 216)); } return ; @@ -437,7 +441,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ( mass->getTotalMass(), 8.0_sreal); EXPECT_FLOATINGPOINT_EQ( mass->getMassDensity()[0], 1.0_sreal); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 0.05_sreal); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(1.0 * 15 / 216)); } return ; @@ -471,7 +475,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ( mass->getTotalMass(), 1.0_sreal); EXPECT_FLOATINGPOINT_EQ( mass->getMassDensity()[0], 0.125_sreal); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 0.00625_sreal); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(0.125 * 15 / 216)); } return ; @@ -502,7 +506,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), 1.0_sreal); EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], 0.125_sreal); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 0.00625_sreal); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(0.125 * 15 / 216)); } return ; @@ -533,7 +537,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), 1.0_sreal); EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], 0.125_sreal); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 0.00625_sreal); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(0.125 * 15 / 216)); } return ; @@ -565,7 +569,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), 1.0_sreal); EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], 0.125_sreal); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 0.00625_sreal); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(0.125 * 15 / 216)); } return ; @@ -596,7 +600,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), 1.0_sreal); EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], 0.125_sreal); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 0.00625_sreal); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(0.125 * 15 / 216)); } return ; @@ -631,7 +635,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), 1.0_sreal); EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], 0.125_sreal); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 0.00625_sreal); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(0.125 * 15 / 216)); } return ; @@ -664,7 +668,7 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ( mass->getMassCount(), 27 ); EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), 2.0_sreal); EXPECT_FLOATINGPOINT_EQ(mass->getMassDensity()[0], 0.25_sreal); - EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], 0.0125_sreal); + EXPECT_FLOATINGPOINT_EQ( mass->getVertexMass()[0], SReal(0.25 * 15 / 216)); } return ; @@ -1273,8 +1277,8 @@ class MeshMatrixMass_test : public BaseTest const VecMass& vMasses = mass->d_vertexMass.getValue(); const VecMass& eMasses = mass->d_edgeMass.getValue(); - static const MassType refValueV = (MassType)(expectedDensity * volumeElem * 1 / 20); - static const MassType refValueE = (MassType)(expectedDensity * volumeElem * 1 / 40); + static const MassType refValueV = (MassType)(expectedDensity * volumeElem * 15 / 216); + static const MassType refValueE = (MassType)(expectedDensity * volumeElem * 4 / 216); // check value at init EXPECT_EQ(vMasses.size(), 27); @@ -1302,6 +1306,10 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ(vMasses.size(), 26); EXPECT_EQ(eMasses.size(), 87); EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), 7.0_sreal); + if (!lumped) + { + checkPositiveDefinite(mass); + } // check vertex mass EXPECT_FLOATINGPOINT_EQ(vMasses[0], refValueV); // check update of Mass when removing tetra @@ -1584,8 +1592,8 @@ class MeshMatrixMass_test : public BaseTest const VecMass& vMasses = mass->d_vertexMass.getValue(); const VecMass& eMasses = mass->d_edgeMass.getValue(); - static const MassType refValueV = (MassType)(expectedDensity * volumeElem * 1 / 8); - static const MassType refValueE = (MassType)(expectedDensity * volumeElem * 1 / 16); + static const MassType refValueV = (MassType)(expectedDensity * volumeElem * 5 / 36); + static const MassType refValueE = (MassType)(expectedDensity * volumeElem * 2 / 36); // check value at init EXPECT_EQ(vMasses.size(), 9); @@ -1613,6 +1621,10 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ(vMasses.size(), 8); EXPECT_EQ(eMasses.size(), 14); EXPECT_FLOATINGPOINT_EQ(mass->getTotalMass(), expectedTotalMass - massElem); + if (!lumped) + { + checkPositiveDefinite(mass); + } // check vertex mass EXPECT_FLOATINGPOINT_EQ(vMasses[0], refValueV); // check update of Mass when removing tetra @@ -1922,6 +1934,118 @@ class MeshMatrixMass_test : public BaseTest EXPECT_EQ(eMasses.size(), 0); EXPECT_NEAR(mass->getTotalMass(), 0, 1e-4); } + + /// Assemble the scalar mass matrix (one row per vertex) from vertexMass and edgeMass, + /// as addMDx does when the mass is not lumped. + static Eigen::MatrixXd assembleScalarMassMatrix(TheMeshMatrixMass* mass) + { + const VecMass& vMasses = mass->d_vertexMass.getValue(); + const VecMass& eMasses = mass->d_edgeMass.getValue(); + const auto& edges = mass->l_topology->getEdges(); + + Eigen::MatrixXd M = Eigen::MatrixXd::Zero(vMasses.size(), vMasses.size()); + for (std::size_t i = 0; i < vMasses.size(); ++i) + { + M(i, i) += vMasses[i]; + } + for (std::size_t i = 0; i < edges.size(); ++i) + { + M(edges[i][0], edges[i][1]) += eMasses[i]; + M(edges[i][1], edges[i][0]) += eMasses[i]; + } + return M; + } + + /// Check that the assembled non-lumped mass matrix is positive definite and sums to totalMass + static void checkPositiveDefinite(TheMeshMatrixMass* mass) + { + const Eigen::MatrixXd M = assembleScalarMassMatrix(mass); + EXPECT_NEAR(M.sum(), mass->getTotalMass(), 1e-10 * mass->getTotalMass()); + + const Eigen::SelfAdjointEigenSolver eigenSolver(M, Eigen::EigenvaluesOnly); + const Eigen::VectorXd& eigenvalues = eigenSolver.eigenvalues(); + const auto nbNonPositive = (eigenvalues.array() <= 1e-12 * eigenvalues.maxCoeff()).count(); + EXPECT_EQ(nbNonPositive, 0) << "The mass matrix has " << nbNonPositive << " non-positive eigenvalues out of " + << eigenvalues.size() << " (min eigenvalue: " << eigenvalues.minCoeff() << ")"; + } + + /// Returns M * 1 (a column of ones per dof direction), i.e. the row sums of the mass matrix + static VecDeriv computeMOnes(TheMeshMatrixMass* mass) + { + const auto nbDofs = mass->d_vertexMass.getValue().size(); + Data ones(VecDeriv(nbDofs, Deriv(1, 1, 1))); + Data res(VecDeriv(nbDofs, Deriv())); + mass->addMDx(core::mechanicalparams::defaultInstance(), res, ones, 1.0); + return res.getValue(); + } + + /// Check that the non-lumped mass matrix is symmetric positive definite, sums to totalMass, + /// and that its row sums are the masses of the lumped matrix built on the same mesh. + void checkMassMatrix(const string& topologyNodes) + { + const auto createScene = [&topologyNodes](bool lumped) + { + const string scene = + "" + "" + " " + + topologyNodes + + " " + " " + ""; + const Node::SPtr root = SceneLoaderXML::loadFromMemory("massMatrix", scene.c_str()); + sofa::simulation::node::initRoot(root.get()); + return root; + }; + + const Node::SPtr root = createScene(false); + ASSERT_NE(root.get(), nullptr); + TheMeshMatrixMass* mass = root->getTreeObject(); + ASSERT_NE(mass, nullptr); + ASSERT_TRUE(mass->isComponentStateValid()); + ASSERT_FALSE(mass->isLumped()); + + EXPECT_NEAR(mass->getTotalMass(), 2.0, 1e-10); + checkPositiveDefinite(mass); + const Eigen::MatrixXd M = assembleScalarMassMatrix(mass); + + // M is applied with the same matrix as the one assembled above + const VecDeriv MOnes = computeMOnes(mass); + ASSERT_EQ(MOnes.size(), M.rows()); + for (std::size_t i = 0; i < MOnes.size(); ++i) + { + EXPECT_NEAR(MOnes[i][0], M.row(i).sum(), 1e-12); + } + + // the row sums are the lumped masses... + const Node::SPtr rootLumped = createScene(true); + ASSERT_NE(rootLumped.get(), nullptr); + TheMeshMatrixMass* lumpedMass = rootLumped->getTreeObject(); + ASSERT_NE(lumpedMass, nullptr); + ASSERT_TRUE(lumpedMass->isComponentStateValid()); + ASSERT_TRUE(lumpedMass->isLumped()); + EXPECT_NEAR(lumpedMass->getTotalMass(), 2.0, 1e-10); + const VecDeriv lumpedMasses = computeMOnes(lumpedMass); + ASSERT_EQ(lumpedMasses.size(), MOnes.size()); + + // ... and the masses used for the weight, whether lumped or not + const Data x(VecCoord(MOnes.size(), Coord())); + const Data v(VecDeriv(MOnes.size(), Deriv())); + Data weight(VecDeriv(MOnes.size(), Deriv())); + Data weightLumped(VecDeriv(MOnes.size(), Deriv())); + mass->addForce(core::mechanicalparams::defaultInstance(), weight, x, v); + lumpedMass->addForce(core::mechanicalparams::defaultInstance(), weightLumped, x, v); + + for (std::size_t i = 0; i < MOnes.size(); ++i) + { + EXPECT_NEAR(MOnes[i][0], lumpedMasses[i][0], 1e-12) << "vertex " << i; + EXPECT_NEAR(-weight.getValue()[i][1], lumpedMasses[i][0], 1e-12) << "vertex " << i; + EXPECT_NEAR(-weightLumped.getValue()[i][1], lumpedMasses[i][0], 1e-12) << "vertex " << i; + } + + sofa::simulation::node::unload(root); + sofa::simulation::node::unload(rootLumped); + } }; @@ -1967,8 +2091,8 @@ TEST_F(MeshMatrixMass3_test, singleQuad) static const MassType surface = 1.0; static const MassType expectedTotalMass = 1.0f; static const MassType density = expectedTotalMass / surface; - static const VecMass expectedVMass(4, (MassType)(density * surface * 1 / 8)); - static const VecMass expectedEMass(4, (MassType)(density * surface * 1 / 16)); + static const VecMass expectedVMass(4, (MassType)(density * surface * 5 / 36)); + static const VecMass expectedEMass(4, (MassType)(density * surface * 2 / 36)); runTest(positions, topologyContainer, geometryAlgorithms, expectedTotalMass, density, expectedVMass, expectedEMass); } @@ -2017,8 +2141,8 @@ TEST_F(MeshMatrixMass3_test, singleHexahedron) static const MassType volume = 1.0; static const MassType expectedTotalMass = 1.0f; static const MassType density = expectedTotalMass / volume; - static const VecMass expectedVMass(8, (MassType)(density * volume * 1 / 20)); - static const VecMass expectedEMass(12, (MassType)(density * volume * 1 / 40)); + static const VecMass expectedVMass(8, (MassType)(density * volume * 15 / 216)); + static const VecMass expectedEMass(12, (MassType)(density * volume * 4 / 216)); runTest(positions, topologyContainer, geometryAlgorithms, expectedTotalMass, density, expectedVMass, expectedEMass); } @@ -2213,4 +2337,30 @@ TEST_F(MeshMatrixMass3_test, checkTopologicalChanges_Edge_lumped) { } +TEST_F(MeshMatrixMass3_test, checkMassMatrix_Hexa) +{ + EXPECT_MSG_NOEMIT(Error); + checkMassMatrix(""); +} + +TEST_F(MeshMatrixMass3_test, checkMassMatrix_Quad) +{ + EXPECT_MSG_NOEMIT(Error); + checkMassMatrix(""); +} + +TEST_F(MeshMatrixMass3_test, checkMassMatrix_Tetra) +{ + EXPECT_MSG_NOEMIT(Error); + checkMassMatrix("" + ""); +} + +TEST_F(MeshMatrixMass3_test, checkMassMatrix_Triangle) +{ + EXPECT_MSG_NOEMIT(Error); + checkMassMatrix("" + ""); +} + } // namespace sofa