From 4af45091cb7d79d1e516a98e24b52c58aa3770e8 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Wed, 26 Aug 2026 08:25:28 +0900 Subject: [PATCH 1/5] add scene checking the projection of mapped matrices --- .../MatrixProjectionEquivalence.scn | 90 +++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 examples/Validation/MatrixProjectionEquivalence.scn diff --git a/examples/Validation/MatrixProjectionEquivalence.scn b/examples/Validation/MatrixProjectionEquivalence.scn new file mode 100644 index 00000000000..f51698ce92d --- /dev/null +++ b/examples/Validation/MatrixProjectionEquivalence.scn @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 1dafe96d2436eeaefdd7aa924f1d57bac05a9178 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Thu, 27 Aug 2026 07:10:22 +0900 Subject: [PATCH 2/5] convert test scene in unit test --- .../LinearSystem/tests/CMakeLists.txt | 2 + .../tests/MatrixProjectionMethod_test.cpp | 198 ++++++++++++++++++ .../MatrixProjectionEquivalence.scn | 90 -------- 3 files changed, 200 insertions(+), 90 deletions(-) create mode 100644 Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp delete mode 100644 examples/Validation/MatrixProjectionEquivalence.scn diff --git a/Sofa/Component/LinearSystem/tests/CMakeLists.txt b/Sofa/Component/LinearSystem/tests/CMakeLists.txt index 96c8a4e9b5e..e225b77aa1b 100644 --- a/Sofa/Component/LinearSystem/tests/CMakeLists.txt +++ b/Sofa/Component/LinearSystem/tests/CMakeLists.txt @@ -4,6 +4,7 @@ project(Sofa.Component.LinearSystem_test) set(SOURCE_FILES MatrixLinearSystem_test.cpp + MatrixProjectionMethod_test.cpp ) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) @@ -13,4 +14,5 @@ target_link_libraries(${PROJECT_NAME} Sofa.Testing Sofa.Component.StateContainer Sofa.Component.Mapping.Linear Sofa.Component.SolidMechanics.Spring + Sofa.Simulation.Graph ) diff --git a/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp b/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp new file mode 100644 index 00000000000..a850427e256 --- /dev/null +++ b/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp @@ -0,0 +1,198 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace sofa +{ + +using DataTypes = defaulttype::Vec3Types; +using MatrixType = linearalgebra::FullMatrix; +using VectorType = linearalgebra::FullVector; +using MatrixSystem = component::linearsystem::MatrixLinearSystem; +using MechanicalObject3 = component::statecontainer::MechanicalObject; +using Spring3 = component::solidmechanics::spring::SpringForceField; +using IdentityMapping3 = component::mapping::linear::IdentityMapping; + +constexpr SReal springStiffness = 100_sreal; +constexpr SReal springRestLength = 1_sreal; + +/// A particle, alone in its node, optionally carrying an identity-mapped copy of itself +struct Particle +{ + simulation::Node::SPtr node; + MechanicalObject3::SPtr state; ///< the degrees of freedom of the system + MechanicalObject3::SPtr mapped; ///< the identity-mapped copy, if any + + /// The state the force field is applied on + MechanicalObject3* target() const + { + return mapped ? mapped.get() : state.get(); + } +}; + +static Particle createParticle(simulation::Node* parent, const std::string& name, + const type::Vec3& position, bool withMapping) +{ + Particle particle; + particle.node = parent->createChild(name); + + particle.state = core::objectmodel::New(); + particle.state->setName("dofs"); + particle.node->addObject(particle.state); + particle.state->resize(1); + particle.state->writePositions()[0] = position; + + if (withMapping) + { + const auto mappedNode = particle.node->createChild(name + "_mapped"); + + particle.mapped = core::objectmodel::New(); + particle.mapped->setName("dofs"); + mappedNode->addObject(particle.mapped); + particle.mapped->resize(1); + particle.mapped->writePositions()[0] = position; + + const auto mapping = core::objectmodel::New(); + mapping->setModels(particle.state.get(), particle.mapped.get()); + mappedNode->addObject(mapping); + } + + return particle; +} + +/// Assembles the matrix of a system made of two independent particles coupled by a +/// spring. When `throughMappings` is true, the spring is applied on identity-mapped +/// copies of the particles rather than on the particles themselves. In both cases the +/// degrees of freedom of the system are the two particles, so the matrix is 6x6. +static void assembleMatrix(bool throughMappings, MatrixType& result) +{ + const simulation::Node::SPtr root = simulation::getSimulation()->createNewGraph("root"); + + const MatrixSystem::SPtr linearSystem = core::objectmodel::New(); + root->addObject(linearSystem); + + const Particle particle0 = createParticle(root.get(), "p0", {0_sreal, 0_sreal, 0_sreal}, throughMappings); + const Particle particle1 = createParticle(root.get(), "p1", {2_sreal, 0_sreal, 0_sreal}, throughMappings); + + const auto spring = core::objectmodel::New(particle0.target(), particle1.target()); + spring->setName("spring"); + root->addObject(spring); + spring->addSpring(0, 0, springStiffness, 0_sreal, springRestLength); + + simulation::node::initRoot(root.get()); + + auto mparams = *core::MechanicalParams::defaultInstance(); + mparams.setKFactor(1_sreal); + + // force fields usually pre-compute elements required by the assembly in addForce + core::MultiVecDerivId forceId = core::vec_id::write_access::externalForce; + static_cast(spring.get())->addForce(&mparams, forceId); + + linearSystem->buildSystemMatrix(&mparams); + + const MatrixType* matrix = linearSystem->getSystemMatrix(); + ASSERT_NE(matrix, nullptr); + + result.resize(matrix->rowSize(), matrix->colSize()); + for (MatrixType::Index i = 0; i < matrix->rowSize(); ++i) + { + for (MatrixType::Index j = 0; j < matrix->colSize(); ++j) + { + result.set(i, j, matrix->element(i, j)); + } + } + + simulation::node::unload(root); +} + +/// A force field applied on mapped states must be projected into the global matrix as +/// J^T K J. Here the mappings are identities, so J = I and the projected matrix must be +/// exactly the matrix obtained by applying the same force field directly on the degrees +/// of freedom of the system. +/// +/// The comparison is sensitive in both directions: a missing coupling term and a +/// spurious contribution both break the equality. +TEST(MatrixProjectionMethod, mappedForceFieldMatchesNonMapped) +{ + MatrixType reference, projected; + assembleMatrix(false, reference); + assembleMatrix(true, projected); + + // two particles of 3 degrees of freedom each + ASSERT_EQ(reference.rowSize(), 6); + ASSERT_EQ(reference.colSize(), 6); + ASSERT_EQ(projected.rowSize(), reference.rowSize()); + ASSERT_EQ(projected.colSize(), reference.colSize()); + + static constexpr SReal tolerance = 1e-12_sreal; + + for (MatrixType::Index i = 0; i < reference.rowSize(); ++i) + { + for (MatrixType::Index j = 0; j < reference.colSize(); ++j) + { + EXPECT_NEAR(projected.element(i, j), reference.element(i, j), tolerance) + << "at (" << i << ", " << j << ")"; + } + } +} + +/// Guards the test above: it would still pass if both matrices were empty. The spring +/// couples the two particles, so both off-diagonal blocks must be non-zero, i.e. the +/// projection must produce the coupling terms and not only the diagonal ones. +TEST(MatrixProjectionMethod, mappedForceFieldProducesCouplingTerms) +{ + MatrixType projected; + assembleMatrix(true, projected); + + ASSERT_EQ(projected.rowSize(), 6); + ASSERT_EQ(projected.colSize(), 6); + + SReal diagonalBlocks = 0_sreal; + SReal offDiagonalBlocks = 0_sreal; + + for (MatrixType::Index i = 0; i < 6; ++i) + { + for (MatrixType::Index j = 0; j < 6; ++j) + { + const bool sameParticle = (i < 3) == (j < 3); + (sameParticle ? diagonalBlocks : offDiagonalBlocks) += std::abs(projected.element(i, j)); + } + } + + EXPECT_GT(diagonalBlocks, 0_sreal); + EXPECT_GT(offDiagonalBlocks, 0_sreal); +} + +} diff --git a/examples/Validation/MatrixProjectionEquivalence.scn b/examples/Validation/MatrixProjectionEquivalence.scn deleted file mode 100644 index f51698ce92d..00000000000 --- a/examples/Validation/MatrixProjectionEquivalence.scn +++ /dev/null @@ -1,90 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - From add5c43069562b935b2f74c359c47ec1c2cb3a99 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Thu, 27 Aug 2026 07:15:26 +0900 Subject: [PATCH 3/5] use sofa simpleapi --- .../tests/MatrixProjectionMethod_test.cpp | 153 ++++++++---------- 1 file changed, 69 insertions(+), 84 deletions(-) diff --git a/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp b/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp index a850427e256..51c01061557 100644 --- a/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp +++ b/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp @@ -21,95 +21,76 @@ ******************************************************************************/ #include -#include -#include - -#include -#include -#include - +#include #include -#include +#include +#include #include +#include #include #include namespace sofa { -using DataTypes = defaulttype::Vec3Types; -using MatrixType = linearalgebra::FullMatrix; -using VectorType = linearalgebra::FullVector; -using MatrixSystem = component::linearsystem::MatrixLinearSystem; -using MechanicalObject3 = component::statecontainer::MechanicalObject; -using Spring3 = component::solidmechanics::spring::SpringForceField; -using IdentityMapping3 = component::mapping::linear::IdentityMapping; - -constexpr SReal springStiffness = 100_sreal; -constexpr SReal springRestLength = 1_sreal; +using Matrix = linearalgebra::FullMatrix; +using MatrixIndex = linearalgebra::BaseMatrix::Index; -/// A particle, alone in its node, optionally carrying an identity-mapped copy of itself -struct Particle +/// Assembles the matrix of a scene made of two independent particles coupled by a +/// spring. The spring is applied either directly on the degrees of freedom, or on +/// identity-mapped copies of them: +/// +/// root +/// |- MatrixLinearSystem +/// |- p0 - dofs (+ p0/mapped/dofs and an IdentityMapping) +/// |- p1 - dofs (+ p1/mapped/dofs and an IdentityMapping) +/// |- SpringForceField (on the dofs, or on the mapped ones) +/// +/// In both cases the degrees of freedom of the system are the two particles, so the +/// assembled matrix is 6x6. +static void assembleMatrix(bool throughMappings, Matrix& result) { - simulation::Node::SPtr node; - MechanicalObject3::SPtr state; ///< the degrees of freedom of the system - MechanicalObject3::SPtr mapped; ///< the identity-mapped copy, if any + const auto plugins = testing::makeScopedPlugin({ + Sofa.Component.LinearSystem, + Sofa.Component.Mapping.Linear, + Sofa.Component.SolidMechanics.Spring, + Sofa.Component.StateContainer}); - /// The state the force field is applied on - MechanicalObject3* target() const - { - return mapped ? mapped.get() : state.get(); - } -}; + const simulation::Node::SPtr root = simulation::getSimulation()->createNewGraph("root"); -static Particle createParticle(simulation::Node* parent, const std::string& name, - const type::Vec3& position, bool withMapping) -{ - Particle particle; - particle.node = parent->createChild(name); + const auto linearSystem = simpleapi::createObject(root, "MatrixLinearSystem", + {{"template", "FullMatrix"}}); - particle.state = core::objectmodel::New(); - particle.state->setName("dofs"); - particle.node->addObject(particle.state); - particle.state->resize(1); - particle.state->writePositions()[0] = position; + static const std::array positions { "0 0 0", "2 0 0" }; + std::array springTargets; - if (withMapping) + for (std::size_t i = 0; i < 2; ++i) { - const auto mappedNode = particle.node->createChild(name + "_mapped"); + const std::string name = "p" + std::to_string(i); - particle.mapped = core::objectmodel::New(); - particle.mapped->setName("dofs"); - mappedNode->addObject(particle.mapped); - particle.mapped->resize(1); - particle.mapped->writePositions()[0] = position; + const auto particle = simpleapi::createChild(root, name); + simpleapi::createObject(particle, "MechanicalObject", + {{"name", "dofs"}, {"template", "Vec3"}, {"position", positions[i]}}); - const auto mapping = core::objectmodel::New(); - mapping->setModels(particle.state.get(), particle.mapped.get()); - mappedNode->addObject(mapping); - } - - return particle; -} + springTargets[i] = "@/" + name + "/dofs"; -/// Assembles the matrix of a system made of two independent particles coupled by a -/// spring. When `throughMappings` is true, the spring is applied on identity-mapped -/// copies of the particles rather than on the particles themselves. In both cases the -/// degrees of freedom of the system are the two particles, so the matrix is 6x6. -static void assembleMatrix(bool throughMappings, MatrixType& result) -{ - const simulation::Node::SPtr root = simulation::getSimulation()->createNewGraph("root"); - - const MatrixSystem::SPtr linearSystem = core::objectmodel::New(); - root->addObject(linearSystem); + if (throughMappings) + { + const auto mapped = simpleapi::createChild(particle, "mapped"); + simpleapi::createObject(mapped, "MechanicalObject", + {{"name", "dofs"}, {"template", "Vec3"}, {"position", positions[i]}}); + simpleapi::createObject(mapped, "IdentityMapping", {}); - const Particle particle0 = createParticle(root.get(), "p0", {0_sreal, 0_sreal, 0_sreal}, throughMappings); - const Particle particle1 = createParticle(root.get(), "p1", {2_sreal, 0_sreal, 0_sreal}, throughMappings); + springTargets[i] = "@/" + name + "/mapped/dofs"; + } + } - const auto spring = core::objectmodel::New(particle0.target(), particle1.target()); - spring->setName("spring"); - root->addObject(spring); - spring->addSpring(0, 0, springStiffness, 0_sreal, springRestLength); + const auto spring = simpleapi::createObject(root, "SpringForceField", + {{"name", "spring"}, + {"object1", springTargets[0]}, + {"object2", springTargets[1]}, + // index1 index2 stiffness damping restLength + {"spring", "0 0 100 0 1"}}); simulation::node::initRoot(root.get()); @@ -117,18 +98,22 @@ static void assembleMatrix(bool throughMappings, MatrixType& result) mparams.setKFactor(1_sreal); // force fields usually pre-compute elements required by the assembly in addForce + auto* forceField = dynamic_cast(spring.get()); + ASSERT_NE(forceField, nullptr); core::MultiVecDerivId forceId = core::vec_id::write_access::externalForce; - static_cast(spring.get())->addForce(&mparams, forceId); + forceField->addForce(&mparams, forceId); - linearSystem->buildSystemMatrix(&mparams); + auto* system = dynamic_cast(linearSystem.get()); + ASSERT_NE(system, nullptr); + system->buildSystemMatrix(&mparams); - const MatrixType* matrix = linearSystem->getSystemMatrix(); + const linearalgebra::BaseMatrix* matrix = system->getSystemBaseMatrix(); ASSERT_NE(matrix, nullptr); result.resize(matrix->rowSize(), matrix->colSize()); - for (MatrixType::Index i = 0; i < matrix->rowSize(); ++i) + for (MatrixIndex i = 0; i < matrix->rowSize(); ++i) { - for (MatrixType::Index j = 0; j < matrix->colSize(); ++j) + for (MatrixIndex j = 0; j < matrix->colSize(); ++j) { result.set(i, j, matrix->element(i, j)); } @@ -137,16 +122,16 @@ static void assembleMatrix(bool throughMappings, MatrixType& result) simulation::node::unload(root); } -/// A force field applied on mapped states must be projected into the global matrix as -/// J^T K J. Here the mappings are identities, so J = I and the projected matrix must be -/// exactly the matrix obtained by applying the same force field directly on the degrees -/// of freedom of the system. +/// A force field acting on mapped states is projected into the global matrix as +/// J^T K J. The mappings here are identities, so J = I, and the projection must +/// reproduce exactly the matrix obtained by applying the same force field directly on +/// the degrees of freedom of the system. /// /// The comparison is sensitive in both directions: a missing coupling term and a /// spurious contribution both break the equality. TEST(MatrixProjectionMethod, mappedForceFieldMatchesNonMapped) { - MatrixType reference, projected; + Matrix reference, projected; assembleMatrix(false, reference); assembleMatrix(true, projected); @@ -158,9 +143,9 @@ TEST(MatrixProjectionMethod, mappedForceFieldMatchesNonMapped) static constexpr SReal tolerance = 1e-12_sreal; - for (MatrixType::Index i = 0; i < reference.rowSize(); ++i) + for (MatrixIndex i = 0; i < reference.rowSize(); ++i) { - for (MatrixType::Index j = 0; j < reference.colSize(); ++j) + for (MatrixIndex j = 0; j < reference.colSize(); ++j) { EXPECT_NEAR(projected.element(i, j), reference.element(i, j), tolerance) << "at (" << i << ", " << j << ")"; @@ -173,7 +158,7 @@ TEST(MatrixProjectionMethod, mappedForceFieldMatchesNonMapped) /// projection must produce the coupling terms and not only the diagonal ones. TEST(MatrixProjectionMethod, mappedForceFieldProducesCouplingTerms) { - MatrixType projected; + Matrix projected; assembleMatrix(true, projected); ASSERT_EQ(projected.rowSize(), 6); @@ -182,9 +167,9 @@ TEST(MatrixProjectionMethod, mappedForceFieldProducesCouplingTerms) SReal diagonalBlocks = 0_sreal; SReal offDiagonalBlocks = 0_sreal; - for (MatrixType::Index i = 0; i < 6; ++i) + for (MatrixIndex i = 0; i < 6; ++i) { - for (MatrixType::Index j = 0; j < 6; ++j) + for (MatrixIndex j = 0; j < 6; ++j) { const bool sameParticle = (i < 3) == (j < 3); (sameParticle ? diagonalBlocks : offDiagonalBlocks) += std::abs(projected.element(i, j)); From 34c5fb6e36da339d4c6fdf676dda81e8ae311429 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Thu, 27 Aug 2026 07:20:41 +0900 Subject: [PATCH 4/5] add ConstantSparsityPatternSystem and ConstantSparsityProjectionMethod tests --- .../tests/MatrixProjectionMethod_test.cpp | 149 ++++++++++++------ 1 file changed, 102 insertions(+), 47 deletions(-) diff --git a/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp b/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp index 51c01061557..3cb93d04f8d 100644 --- a/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp +++ b/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp @@ -36,19 +36,42 @@ namespace sofa using Matrix = linearalgebra::FullMatrix; using MatrixIndex = linearalgebra::BaseMatrix::Index; +/// The linear system and the projection method under test. Both have a variant taking +/// advantage of a constant sparsity pattern, and both variants take a different code +/// path once the pattern has been built, so every combination is tested. +struct Parameters +{ + std::string linearSystem; + std::string projectionMethod; +}; + +static std::string testName(const ::testing::TestParamInfo& info) +{ + return info.param.linearSystem + "_" + info.param.projectionMethod; +} + +/// The assembled matrix is 6x6: two independent particles of 3 degrees of freedom. +static constexpr MatrixIndex matrixSize = 6; + +/// The number of times the matrix is assembled. The systems and projection methods +/// relying on a constant sparsity pattern only use it from the second assembly on, so +/// more than one is required to cover them. +static constexpr unsigned int nbAssembly = 3; + /// Assembles the matrix of a scene made of two independent particles coupled by a /// spring. The spring is applied either directly on the degrees of freedom, or on /// identity-mapped copies of them: /// /// root -/// |- MatrixLinearSystem +/// |- /// |- p0 - dofs (+ p0/mapped/dofs and an IdentityMapping) /// |- p1 - dofs (+ p1/mapped/dofs and an IdentityMapping) /// |- SpringForceField (on the dofs, or on the mapped ones) +/// |- (only when the spring is applied on mapped states) /// -/// In both cases the degrees of freedom of the system are the two particles, so the -/// assembled matrix is 6x6. -static void assembleMatrix(bool throughMappings, Matrix& result) +/// The matrix is assembled `nbAssembly` times, and every assembly is returned. +static void assembleMatrix(const Parameters& parameters, bool throughMappings, + std::array& results) { const auto plugins = testing::makeScopedPlugin({ Sofa.Component.LinearSystem, @@ -58,8 +81,8 @@ static void assembleMatrix(bool throughMappings, Matrix& result) const simulation::Node::SPtr root = simulation::getSimulation()->createNewGraph("root"); - const auto linearSystem = simpleapi::createObject(root, "MatrixLinearSystem", - {{"template", "FullMatrix"}}); + const auto linearSystem = simpleapi::createObject(root, parameters.linearSystem, + {{"template", "CompressedRowSparseMatrixd"}}); static const std::array positions { "0 0 0", "2 0 0" }; std::array springTargets; @@ -85,6 +108,21 @@ static void assembleMatrix(bool throughMappings, Matrix& result) } } + if (throughMappings) + { + // one projection method per ordered pair of mapped states: the contribution of + // the spring is split into as many mapped matrices + for (const auto& first : springTargets) + { + for (const auto& second : springTargets) + { + simpleapi::createObject(root, parameters.projectionMethod, + {{"template", "CompressedRowSparseMatrixd"}, + {"mechanicalStates", first + " " + second}}); + } + } + } + const auto spring = simpleapi::createObject(root, "SpringForceField", {{"name", "spring"}, {"object1", springTargets[0]}, @@ -97,31 +135,41 @@ static void assembleMatrix(bool throughMappings, Matrix& result) auto mparams = *core::MechanicalParams::defaultInstance(); mparams.setKFactor(1_sreal); - // force fields usually pre-compute elements required by the assembly in addForce auto* forceField = dynamic_cast(spring.get()); ASSERT_NE(forceField, nullptr); - core::MultiVecDerivId forceId = core::vec_id::write_access::externalForce; - forceField->addForce(&mparams, forceId); auto* system = dynamic_cast(linearSystem.get()); ASSERT_NE(system, nullptr); - system->buildSystemMatrix(&mparams); - const linearalgebra::BaseMatrix* matrix = system->getSystemBaseMatrix(); - ASSERT_NE(matrix, nullptr); - - result.resize(matrix->rowSize(), matrix->colSize()); - for (MatrixIndex i = 0; i < matrix->rowSize(); ++i) + for (unsigned int assembly = 0; assembly < nbAssembly; ++assembly) { - for (MatrixIndex j = 0; j < matrix->colSize(); ++j) + // force fields usually pre-compute elements required by the assembly in addForce + core::MultiVecDerivId forceId = core::vec_id::write_access::externalForce; + forceField->addForce(&mparams, forceId); + + system->buildSystemMatrix(&mparams); + + const linearalgebra::BaseMatrix* matrix = system->getSystemBaseMatrix(); + ASSERT_NE(matrix, nullptr); + ASSERT_EQ(matrix->rowSize(), matrixSize); + ASSERT_EQ(matrix->colSize(), matrixSize); + + Matrix& result = results[assembly]; + result.resize(matrixSize, matrixSize); + for (MatrixIndex i = 0; i < matrixSize; ++i) { - result.set(i, j, matrix->element(i, j)); + for (MatrixIndex j = 0; j < matrixSize; ++j) + { + result.set(i, j, matrix->element(i, j)); + } } } simulation::node::unload(root); } +class MatrixProjectionMethodTest : public ::testing::TestWithParam {}; + /// A force field acting on mapped states is projected into the global matrix as /// J^T K J. The mappings here are identities, so J = I, and the projection must /// reproduce exactly the matrix obtained by applying the same force field directly on @@ -129,26 +177,23 @@ static void assembleMatrix(bool throughMappings, Matrix& result) /// /// The comparison is sensitive in both directions: a missing coupling term and a /// spurious contribution both break the equality. -TEST(MatrixProjectionMethod, mappedForceFieldMatchesNonMapped) +TEST_P(MatrixProjectionMethodTest, mappedForceFieldMatchesNonMapped) { - Matrix reference, projected; - assembleMatrix(false, reference); - assembleMatrix(true, projected); - - // two particles of 3 degrees of freedom each - ASSERT_EQ(reference.rowSize(), 6); - ASSERT_EQ(reference.colSize(), 6); - ASSERT_EQ(projected.rowSize(), reference.rowSize()); - ASSERT_EQ(projected.colSize(), reference.colSize()); + std::array reference, projected; + assembleMatrix(GetParam(), false, reference); + assembleMatrix(GetParam(), true, projected); static constexpr SReal tolerance = 1e-12_sreal; - for (MatrixIndex i = 0; i < reference.rowSize(); ++i) + for (unsigned int assembly = 0; assembly < nbAssembly; ++assembly) { - for (MatrixIndex j = 0; j < reference.colSize(); ++j) + for (MatrixIndex i = 0; i < matrixSize; ++i) { - EXPECT_NEAR(projected.element(i, j), reference.element(i, j), tolerance) - << "at (" << i << ", " << j << ")"; + for (MatrixIndex j = 0; j < matrixSize; ++j) + { + EXPECT_NEAR(projected[assembly].element(i, j), reference[assembly].element(i, j), tolerance) + << "at (" << i << ", " << j << ") of assembly " << assembly; + } } } } @@ -156,28 +201,38 @@ TEST(MatrixProjectionMethod, mappedForceFieldMatchesNonMapped) /// Guards the test above: it would still pass if both matrices were empty. The spring /// couples the two particles, so both off-diagonal blocks must be non-zero, i.e. the /// projection must produce the coupling terms and not only the diagonal ones. -TEST(MatrixProjectionMethod, mappedForceFieldProducesCouplingTerms) +TEST_P(MatrixProjectionMethodTest, mappedForceFieldProducesCouplingTerms) { - Matrix projected; - assembleMatrix(true, projected); + std::array projected; + assembleMatrix(GetParam(), true, projected); - ASSERT_EQ(projected.rowSize(), 6); - ASSERT_EQ(projected.colSize(), 6); - - SReal diagonalBlocks = 0_sreal; - SReal offDiagonalBlocks = 0_sreal; - - for (MatrixIndex i = 0; i < 6; ++i) + for (unsigned int assembly = 0; assembly < nbAssembly; ++assembly) { - for (MatrixIndex j = 0; j < 6; ++j) + SReal diagonalBlocks = 0_sreal; + SReal offDiagonalBlocks = 0_sreal; + + for (MatrixIndex i = 0; i < matrixSize; ++i) { - const bool sameParticle = (i < 3) == (j < 3); - (sameParticle ? diagonalBlocks : offDiagonalBlocks) += std::abs(projected.element(i, j)); + for (MatrixIndex j = 0; j < matrixSize; ++j) + { + const bool sameParticle = (i < 3) == (j < 3); + (sameParticle ? diagonalBlocks : offDiagonalBlocks) += + std::abs(projected[assembly].element(i, j)); + } } - } - EXPECT_GT(diagonalBlocks, 0_sreal); - EXPECT_GT(offDiagonalBlocks, 0_sreal); + EXPECT_GT(diagonalBlocks, 0_sreal) << "at assembly " << assembly; + EXPECT_GT(offDiagonalBlocks, 0_sreal) << "at assembly " << assembly; + } } +INSTANTIATE_TEST_SUITE_P(MatrixProjectionMethod, MatrixProjectionMethodTest, + ::testing::ValuesIn(std::vector{ + {"MatrixLinearSystem", "MatrixProjectionMethod"}, + {"MatrixLinearSystem", "ConstantSparsityProjectionMethod"}, + {"ConstantSparsityPatternSystem", "MatrixProjectionMethod"}, + {"ConstantSparsityPatternSystem", "ConstantSparsityProjectionMethod"}, + }), + testName); + } From 1ec6c3c7c77c40582be81923382f96c3c846dcfa Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Thu, 27 Aug 2026 07:24:20 +0900 Subject: [PATCH 5/5] add fullmatrix template --- .../tests/MatrixProjectionMethod_test.cpp | 27 ++++++++++++------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp b/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp index 3cb93d04f8d..3f76ed87474 100644 --- a/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp +++ b/Sofa/Component/LinearSystem/tests/MatrixProjectionMethod_test.cpp @@ -36,18 +36,21 @@ namespace sofa using Matrix = linearalgebra::FullMatrix; using MatrixIndex = linearalgebra::BaseMatrix::Index; -/// The linear system and the projection method under test. Both have a variant taking -/// advantage of a constant sparsity pattern, and both variants take a different code -/// path once the pattern has been built, so every combination is tested. +/// The container of the global matrix, the linear system and the projection method under +/// test. The linear system and the projection method both have a variant taking advantage +/// of a constant sparsity pattern, and both variants take a different code path once the +/// pattern has been built, so every combination is tested. struct Parameters { + std::string matrixTemplate; std::string linearSystem; std::string projectionMethod; }; static std::string testName(const ::testing::TestParamInfo& info) { - return info.param.linearSystem + "_" + info.param.projectionMethod; + return info.param.matrixTemplate + "_" + info.param.linearSystem + "_" + + info.param.projectionMethod; } /// The assembled matrix is 6x6: two independent particles of 3 degrees of freedom. @@ -82,7 +85,7 @@ static void assembleMatrix(const Parameters& parameters, bool throughMappings, const simulation::Node::SPtr root = simulation::getSimulation()->createNewGraph("root"); const auto linearSystem = simpleapi::createObject(root, parameters.linearSystem, - {{"template", "CompressedRowSparseMatrixd"}}); + {{"template", parameters.matrixTemplate}}); static const std::array positions { "0 0 0", "2 0 0" }; std::array springTargets; @@ -111,7 +114,8 @@ static void assembleMatrix(const Parameters& parameters, bool throughMappings, if (throughMappings) { // one projection method per ordered pair of mapped states: the contribution of - // the spring is split into as many mapped matrices + // the spring is split into as many mapped matrices. Contrary to the global + // matrix, the mapped matrices are always in the CRS format. for (const auto& first : springTargets) { for (const auto& second : springTargets) @@ -228,10 +232,13 @@ TEST_P(MatrixProjectionMethodTest, mappedForceFieldProducesCouplingTerms) INSTANTIATE_TEST_SUITE_P(MatrixProjectionMethod, MatrixProjectionMethodTest, ::testing::ValuesIn(std::vector{ - {"MatrixLinearSystem", "MatrixProjectionMethod"}, - {"MatrixLinearSystem", "ConstantSparsityProjectionMethod"}, - {"ConstantSparsityPatternSystem", "MatrixProjectionMethod"}, - {"ConstantSparsityPatternSystem", "ConstantSparsityProjectionMethod"}, + {"FullMatrix", "MatrixLinearSystem", "MatrixProjectionMethod"}, + {"FullMatrix", "MatrixLinearSystem", "ConstantSparsityProjectionMethod"}, + {"CompressedRowSparseMatrixd", "MatrixLinearSystem", "MatrixProjectionMethod"}, + {"CompressedRowSparseMatrixd", "MatrixLinearSystem", "ConstantSparsityProjectionMethod"}, + // ConstantSparsityPatternSystem only exists for the CRS format + {"CompressedRowSparseMatrixd", "ConstantSparsityPatternSystem", "MatrixProjectionMethod"}, + {"CompressedRowSparseMatrixd", "ConstantSparsityPatternSystem", "ConstantSparsityProjectionMethod"}, }), testName);