From 5a2282ccb6d6062e2241f823edc0732d2dc56ae0 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 2 Feb 2026 15:35:27 +0100 Subject: [PATCH 1/8] Refactor `ForceFieldTestCreation`: modularize functions for clarity --- .../testing/ForceFieldTestCreation.h | 199 ++++++++++-------- 1 file changed, 117 insertions(+), 82 deletions(-) diff --git a/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h b/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h index 34b22f048da..7750fa18cd5 100644 --- a/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h +++ b/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h @@ -146,6 +146,103 @@ struct ForceField_test : public BaseSimulationTest, NumericTest(node); } + void setupState( const VecCoord& x, const VecDeriv& v) + { + std::size_t n = x.size(); + this->dof->resize(static_cast(n)); + typename DOF::WriteVecCoord xdof = this->dof->writePositions(); + sofa::testing::copyToData( xdof, x ); + typename DOF::WriteVecDeriv vdof = this->dof->writeVelocities(); + sofa::testing::copyToData( vdof, v ); + } + + void computeForce(core::MechanicalParams* mparams) const + { + MechanicalResetForceVisitor resetForce(mparams, core::vec_id::write_access::force); + node->execute(resetForce); + MechanicalComputeForceVisitor computeForce( mparams, core::vec_id::write_access::force ); + this->node->execute(computeForce); + } + + void checkForce(const VecDeriv& ef, const std::string& message = "") + { + typename DOF::ReadVecDeriv f= this->dof->readForces(); + EXPECT_LT( this->vectorMaxDiff(f,ef), errorMax*this->epsilon() ) << message; + } + + void checkPotentialEnergy(const core::MechanicalParams* mparams, SReal potentialEnergyBefore, const VecDeriv& curF, const VecDeriv& dX) + { + if( !(flags & TEST_POTENTIAL_ENERGY) ) return; + + // Get potential energy after displacement of dofs + SReal potentialEnergyAfterDisplacement = force->getPotentialEnergy(mparams); + + // Check getPotentialEnergy() we should have dE = -dX.F + + // Compute dE = E(x+dx)-E(x) + SReal differencePotentialEnergy = potentialEnergyAfterDisplacement-potentialEnergyBefore; + + // Compute the expected difference of potential energy: -dX.F (dot product between applied displacement and Force) + SReal expectedDifferencePotentialEnergy = 0; + for( unsigned i=0; i errorFactorPotentialEnergy*errorMax*this->epsilon() ){ + ADD_FAILURE()<<"dPotentialEnergy differs from -dX.F (threshold=" << errorFactorPotentialEnergy*errorMax*this->epsilon() << ")" << std::endl + << "dPotentialEnergy is " << differencePotentialEnergy << std::endl + << "-dX.F is " << expectedDifferencePotentialEnergy << std::endl + << "Failed seed number = " << this->seed << std::endl; + } + } + + void checkComputeDf(core::MechanicalParams* mparams, const VecDeriv& dX, const VecDeriv& changeOfForce) + { + MechanicalResetForceVisitor resetForce(mparams, core::vec_id::write_access::force); + node->execute(resetForce); + dof->vRealloc( mparams, core::vec_id::write_access::dx); // dx is not allocated by default + typename DOF::WriteVecDeriv wdx = dof->writeDx(); + sofa::testing::copyToData ( wdx, dX ); + MechanicalComputeDfVisitor computeDf( mparams, core::vec_id::write_access::force ); + node->execute(computeDf); + VecDeriv dF; + sofa::testing::copyFromData( dF, dof->readForces() ); + + EXPECT_LE(this->vectorMaxDiff(changeOfForce, dF), errorMax * this->epsilon()) << + "dF differs from change of force\n" + "Failed seed number = " << this->seed; + } + + void checkStiffnessMatrix(core::MechanicalParams* mparams, const VecDeriv& dX, const VecDeriv& changeOfForce) + { + typedef sofa::linearalgebra::EigenBaseSparseMatrix Sqmat; + const std::size_t n = dX.size(); + const sofa::SignedIndex matrixSize = static_cast(n * DataTypes::deriv_total_size); + Sqmat K( matrixSize, matrixSize); + sofa::core::behavior::SingleMatrixAccessor accessor( &K ); + mparams->setKFactor(1.0); + force->addKToMatrix( mparams, &accessor); + K.compress(); + + modeling::Vector dx; + sofa::testing::data_traits::VecDeriv_to_Vector( dx, dX ); + + modeling::Vector Kdx = K * dx; + modeling::Vector df; + sofa::testing::data_traits::VecDeriv_to_Vector( df, changeOfForce ); + + if( debug ){ + std::cout << " dX = " << dX << std::endl; + std::cout << " change of force = " << changeOfForce << std::endl; + std::cout << " Kdx = " << Kdx.transpose() << std::endl; + } + + EXPECT_LE( this->vectorMaxDiff(Kdx, df), errorMax * this->epsilon() ) << + "Kdx differs from change of force" + "\nFailed seed number = " << this->seed; + } + /** * @brief Given positions and velocities, checks that the expected forces are obtained, and that a small change of positions generates the corresponding change of forces. * @param x positions @@ -169,33 +266,29 @@ struct ForceField_test : public BaseSimulationTest, NumericTestdof->resize(static_cast(n)); - typename DOF::WriteVecCoord xdof = this->dof->writePositions(); - sofa::testing::copyToData( xdof, x ); - typename DOF::WriteVecDeriv vdof = this->dof->writeVelocities(); - sofa::testing::copyToData( vdof, v ); + setupState(x, v); // init scene and compute force if (initScene) { sofa::simulation::node::initRoot(this->node.get()); } + core::MechanicalParams mparams; mparams.setKFactor(1.0); - MechanicalResetForceVisitor resetForce(&mparams, core::vec_id::write_access::force); - node->execute(resetForce); - MechanicalComputeForceVisitor computeForce( &mparams, core::vec_id::write_access::force ); - this->node->execute(computeForce); + + computeForce(&mparams); // check force - typename DOF::ReadVecDeriv f= this->dof->readForces(); - if(debug){ + if(debug) + { + typename DOF::ReadVecDeriv f= this->dof->readForces(); std::cout << "run_test, x = " << x << std::endl; std::cout << " v = " << v << std::endl; std::cout << " expected f = " << ef << std::endl; std::cout << " actual f = " << f.ref() << std::endl; } - ASSERT_LT( this->vectorMaxDiff(f,ef), errorMax*this->epsilon() ); + checkForce(ef); if( !checkStiffness ) return; @@ -206,18 +299,21 @@ struct ForceField_test : public BaseSimulationTest, NumericTestreadForces() ); // Get potential Energy before applying a displacement to dofs - SReal potentialEnergyBeforeDisplacement = (flags & TEST_POTENTIAL_ENERGY) ? ((const core::behavior::BaseForceField*)force.get())->getPotentialEnergy(&mparams) : 0; + SReal potentialEnergyBeforeDisplacement = (flags & TEST_POTENTIAL_ENERGY) ? force->getPotentialEnergy(&mparams) : 0; // change position VecDeriv dX(n); - for( unsigned i=0; iepsilon(), deltaRange.second * this->epsilon() ); // todo: better random, with negative values - xdof[i] += dX[i]; + { + typename DOF::WriteVecCoord xdof = this->dof->writePositions(); + for( unsigned i=0; iepsilon(), deltaRange.second * this->epsilon() ); // todo: better random, with negative values + xdof[i] += dX[i]; + } } // compute new force and difference between previous force - node->execute(resetForce); - node->execute(computeForce); + computeForce(&mparams); + VecDeriv newF; sofa::testing::copyFromData( newF, dof->readForces() ); VecDeriv changeOfForce(curF); @@ -225,72 +321,11 @@ struct ForceField_test : public BaseSimulationTest, NumericTestgetPotentialEnergy(&mparams); - - // Check getPotentialEnergy() we should have dE = -dX.F - - // Compute dE = E(x+dx)-E(x) - SReal differencePotentialEnergy = potentialEnergyAfterDisplacement-potentialEnergyBeforeDisplacement; - - // Compute the expected difference of potential energy: -dX.F (dot product between applied displacement and Force) - SReal expectedDifferencePotentialEnergy = 0; - for( unsigned i=0; i errorFactorPotentialEnergy*errorMax*this->epsilon() ){ - ADD_FAILURE()<<"dPotentialEnergy differs from -dX.F (threshold=" << errorFactorPotentialEnergy*errorMax*this->epsilon() << ")" << std::endl - << "dPotentialEnergy is " << differencePotentialEnergy << std::endl - << "-dX.F is " << expectedDifferencePotentialEnergy << std::endl - << "Failed seed number = " << this->seed << std::endl; - } - } - - - // check computeDf: compare its result to actual change - node->execute(resetForce); - dof->vRealloc( &mparams, core::vec_id::write_access::dx); // dx is not allocated by default - typename DOF::WriteVecDeriv wdx = dof->writeDx(); - sofa::testing::copyToData ( wdx, dX ); - MechanicalComputeDfVisitor computeDf( &mparams, core::vec_id::write_access::force ); - node->execute(computeDf); - VecDeriv dF; - sofa::testing::copyFromData( dF, dof->readForces() ); - - EXPECT_LE(this->vectorMaxDiff(changeOfForce, dF), errorMax * this->epsilon()) << - "dF differs from change of force\n" - "Failed seed number = " << this->seed; - - // check stiffness matrix: compare its product with dx to actual force change - typedef sofa::linearalgebra::EigenBaseSparseMatrix Sqmat; - const sofa::SignedIndex matrixSize = static_cast(n * DataTypes::deriv_total_size); - Sqmat K( matrixSize, matrixSize); - sofa::core::behavior::SingleMatrixAccessor accessor( &K ); - mparams.setKFactor(1.0); - force->addKToMatrix( &mparams, &accessor); - K.compress(); - // cout << "stiffness: " << K << endl; - modeling::Vector dx; - sofa::testing::data_traits::VecDeriv_to_Vector( dx, dX ); + checkPotentialEnergy(&mparams, potentialEnergyBeforeDisplacement, curF, dX); - modeling::Vector Kdx = K * dx; - if( debug ){ - std::cout << " dX = " << dX << std::endl; - std::cout << " newF = " << newF << std::endl; - std::cout << " change of force = " << changeOfForce << std::endl; - std::cout << " addDforce = " << dF << std::endl; - std::cout << " Kdx = " << Kdx.transpose() << std::endl; - } + checkComputeDf(&mparams, dX, changeOfForce); - modeling::Vector df; - sofa::testing::data_traits::VecDeriv_to_Vector( df, changeOfForce ); - EXPECT_LE( this->vectorMaxDiff(Kdx, df), errorMax * this->epsilon() ) << - "Kdx differs from change of force" - "\nFailed seed number = " << this->seed; + checkStiffnessMatrix(&mparams, dX, changeOfForce); } From 327be3f6facd286de343d7e449376b63e70f9a41 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 2 Feb 2026 15:52:44 +0100 Subject: [PATCH 2/8] Enhance `ForceFieldTestCreation`: add `checkBuildStiffnessMatrix` and refactor stiffness matrix checks --- .../testing/ForceFieldTestCreation.h | 59 +++++++++++++++++-- 1 file changed, 53 insertions(+), 6 deletions(-) diff --git a/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h b/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h index 7750fa18cd5..d3504e1da9b 100644 --- a/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h +++ b/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h @@ -36,6 +36,7 @@ using sofa::testing::NumericTest; #include #include #include +#include #include using sofa::simulation::mechanicalvisitor::MechanicalComputeDfVisitor; @@ -189,7 +190,8 @@ struct ForceField_test : public BaseSimulationTest, NumericTest errorFactorPotentialEnergy*errorMax*this->epsilon() ){ + if( absoluteErrorPotentialEnergy> errorFactorPotentialEnergy*errorMax*this->epsilon() ) + { ADD_FAILURE()<<"dPotentialEnergy differs from -dX.F (threshold=" << errorFactorPotentialEnergy*errorMax*this->epsilon() << ")" << std::endl << "dPotentialEnergy is " << differencePotentialEnergy << std::endl << "-dX.F is " << expectedDifferencePotentialEnergy << std::endl @@ -214,7 +216,7 @@ struct ForceField_test : public BaseSimulationTest, NumericTest Date: Mon, 2 Feb 2026 18:02:36 +0100 Subject: [PATCH 4/8] use concept for ForceField_test template parameter --- .../solidmechanics/testing/ForceFieldTestCreation.h | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h b/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h index d3504e1da9b..fedceca6fa5 100644 --- a/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h +++ b/Sofa/Component/SolidMechanics/Testing/src/sofa/component/solidmechanics/testing/ForceFieldTestCreation.h @@ -22,12 +22,7 @@ #pragma once #include -using sofa::testing::BaseSimulationTest; - #include -using sofa::testing::NumericTest; - -#include #include #include #include @@ -58,8 +53,8 @@ namespace sofa { * @author François Faure, 2014 * */ -template -struct ForceField_test : public BaseSimulationTest, NumericTest +template requires std::is_base_of_v +struct ForceField_test : public sofa::testing::BaseSimulationTest, public sofa::testing::NumericTest { typedef _ForceFieldType ForceField; typedef typename ForceField::DataTypes DataTypes; From ae01a2ea5349a827f3823eee40c19de21d1c2e9e Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 2 Feb 2026 18:10:03 +0100 Subject: [PATCH 5/8] [Component] Expose `getPotentialEnergy` in `QuadPressureForceField` --- .../component/mechanicalload/QuadPressureForceField.h | 3 ++- .../mechanicalload/QuadPressureForceField.inl | 10 +++++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/QuadPressureForceField.h b/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/QuadPressureForceField.h index 5880a219c41..3e049f67a05 100644 --- a/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/QuadPressureForceField.h +++ b/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/QuadPressureForceField.h @@ -118,7 +118,8 @@ class QuadPressureForceField : public core::behavior::ForceField void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /*matrix*/) override; void buildDampingMatrix(core::behavior::DampingMatrix* /*matrix*/) final; - SReal getPotentialEnergy(const core::MechanicalParams* /*mparams*/, const DataVecCoord& /* x */) const override { msg_warning() << "Method getPotentialEnergy not implemented yet."; return 0.0; } + using Inherit1::getPotentialEnergy; + SReal getPotentialEnergy(const core::MechanicalParams* /*mparams*/, const DataVecCoord& /* x */) const override; void draw(const core::visual::VisualParams* vparams) override; diff --git a/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/QuadPressureForceField.inl b/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/QuadPressureForceField.inl index 0b869038c8c..45428a8e616 100644 --- a/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/QuadPressureForceField.inl +++ b/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/QuadPressureForceField.inl @@ -232,7 +232,15 @@ void QuadPressureForceField::buildDampingMatrix(core::behavior::Dampi // No damping in this ForceField } -template +template +SReal QuadPressureForceField::getPotentialEnergy(const core::MechanicalParams*, + const DataVecCoord&) const +{ + msg_warning() << "Method getPotentialEnergy not implemented yet."; + return 0.0; +} + +template void QuadPressureForceField::draw(const core::visual::VisualParams* vparams) { const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); From 8f04f7df9ca0b1c233f23837e11cbda4a5f28572 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 2 Feb 2026 18:10:50 +0100 Subject: [PATCH 6/8] [Component] Expose `getPotentialEnergy` in `PenalityContactForceField` --- .../collision/response/contact/PenalityContactForceField.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Sofa/Component/Collision/Response/Contact/src/sofa/component/collision/response/contact/PenalityContactForceField.h b/Sofa/Component/Collision/Response/Contact/src/sofa/component/collision/response/contact/PenalityContactForceField.h index 1ad75539831..e7f485fa4f2 100644 --- a/Sofa/Component/Collision/Response/Contact/src/sofa/component/collision/response/contact/PenalityContactForceField.h +++ b/Sofa/Component/Collision/Response/Contact/src/sofa/component/collision/response/contact/PenalityContactForceField.h @@ -123,6 +123,7 @@ class PenalityContactForceField : public core::behavior::PairInteractionForceFie void buildDampingMatrix(core::behavior::DampingMatrix* /*matrix*/) final; + using Inherit::getPotentialEnergy; SReal getPotentialEnergy(const sofa::core::MechanicalParams*, const DataVecCoord&, const DataVecCoord& ) const override; const type::vector< Contact >& getContact() const { return contacts.getValue();} From 158e6eb3acdf304f7b2ea20b4ca8a43b53fdf3dd Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 2 Feb 2026 18:10:50 +0100 Subject: [PATCH 7/8] [Component] Expose `getPotentialEnergy` in `TrianglePressureForceField` --- .../sofa/component/mechanicalload/TrianglePressureForceField.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/TrianglePressureForceField.h b/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/TrianglePressureForceField.h index 8aeffb3ad10..039cb74fa87 100644 --- a/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/TrianglePressureForceField.h +++ b/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/TrianglePressureForceField.h @@ -117,6 +117,7 @@ class TrianglePressureForceField : public core::behavior::ForceField void buildDampingMatrix(core::behavior::DampingMatrix* /*matrix*/) final; + using Inherit1::getPotentialEnergy; SReal getPotentialEnergy(const core::MechanicalParams* /*mparams*/, const DataVecCoord& /* x */) const override; void draw(const core::visual::VisualParams* vparams) override; From 6caa59cd643caf919f1c86e69c2bc6999e3a8d0d Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 2 Feb 2026 18:14:35 +0100 Subject: [PATCH 8/8] [Component] Expose `getPotentialEnergy` in `SurfacePressureForceField` --- .../sofa/component/mechanicalload/SurfacePressureForceField.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/SurfacePressureForceField.h b/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/SurfacePressureForceField.h index e866cd70f32..d5846d28abd 100644 --- a/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/SurfacePressureForceField.h +++ b/Sofa/Component/MechanicalLoad/src/sofa/component/mechanicalload/SurfacePressureForceField.h @@ -103,6 +103,7 @@ class SurfacePressureForceField : public core::behavior::ForceField void addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override; void buildStiffnessMatrix(core::behavior::StiffnessMatrix* matrix) override; void buildDampingMatrix(core::behavior::DampingMatrix* /*matrix*/) final; + using Inherit1::getPotentialEnergy; SReal getPotentialEnergy(const core::MechanicalParams* /*mparams*/, const DataVecCoord& /* x */) const override; void draw(const core::visual::VisualParams* vparams) override;