Skip to content

[HyperElastic] Fix inverted Stable Neo-Hookean scene and add regression coverage #6377

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
hugtalbot merged 5 commits into from
Nov 5, 2026
Merged

[HyperElastic] Fix inverted Stable Neo-Hookean scene and add regression coverage #6377

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d2eaff8
[HyperElastic] FIX inverted Stable Neo-Hookean scene and add regressi…
lililii124 Oct 10, 2026
1b6e093
Restore original Stable Neo-Hookean mass and timestep after review
lililii124 Oct 26, 2026
3406ea4
Update CHANGELOG.md
lililii124 Oct 30, 2026
64d2fbc
Merge branch 'master' into fix/6220-stable-neohookean-regression
lililii124 Oct 30, 2026
aeff0c0
Remove inverted example from state regression suite
lililii124 Nov 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
108 changes: 95 additions & 13 deletions Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/Material_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,12 +72,36 @@ generatePositiveDefiniteMatrix(sofa::testing::LinearCongruentialRandomGenerator&
return W;
}

void computeStrainInformation(StrainInformation<defaulttype::Vec3Types>& strain)
StrainInformation<defaulttype::Vec3Types>::MatrixSym
generateTestDeformationTensor(sofa::testing::LinearCongruentialRandomGenerator& lcg, bool inverted)
{
if (!inverted)
{
return generatePositiveDefiniteMatrix(lcg);
}

// C = F^T F for F = {{1.2, 0.2, 0.1}, {0, 0.9, 0.15}, {0, 0, -0.8}}.
// Keep the inverted sample away from singularity for finite differences.
StrainInformation<defaulttype::Vec3Types>::MatrixSym C;
C(0, 0) = 1.44;
C(0, 1) = 0.24;
C(0, 2) = 0.12;
C(1, 1) = 0.85;
C(1, 2) = 0.155;
C(2, 2) = 0.6725;
return C;
}

void computeStrainInformation(StrainInformation<defaulttype::Vec3Types>& strain, bool inverted = false)
{
strain.trC = sofa::type::trace(strain.deformationTensor);

//det(C)=J^2
// det(C) = J^2 does not retain the orientation of the deformation.
strain.J = std::sqrt(sofa::type::determinant(strain.deformationTensor));
if (inverted)
{
strain.J = -strain.J;
}
}

SReal perturbedStrainEnergy(
Expand All @@ -97,7 +121,7 @@ SReal perturbedStrainEnergy(
}

strain.deformationTensor[i] += h;
computeStrainInformation(strain);
computeStrainInformation(strain, strain_0.J < 0);

return material.getStrainEnergy(&strain, materialParameters);
}
Expand All @@ -119,7 +143,7 @@ StrainInformation<defaulttype::Vec3Types>::MatrixSym perturbedPK2(
}

strain.deformationTensor[i] += h;
computeStrainInformation(strain);
computeStrainInformation(strain, strain_0.J < 0);

StrainInformation<defaulttype::Vec3Types>::MatrixSym PK2;
material.deriveSPKTensor(&strain, materialParameters, PK2);
Expand All @@ -128,17 +152,19 @@ StrainInformation<defaulttype::Vec3Types>::MatrixSym perturbedPK2(

void testSecondPiolaKirchhoffFromStrainEnergyDensityFunction(
HyperelasticMaterial<defaulttype::Vec3Types>& material,
const MaterialParameters<defaulttype::Vec3Types>& materialParameters)
const MaterialParameters<defaulttype::Vec3Types>& materialParameters,
bool inverted = false)
{
using Real = defaulttype::Vec3Types::Coord::value_type;

StrainInformation<defaulttype::Vec3Types> strain;

//random right Cauchy-Green deformation tensor
sofa::testing::LinearCongruentialRandomGenerator lcg(96547);
strain.deformationTensor = generatePositiveDefiniteMatrix(lcg);
strain.deformationTensor = generateTestDeformationTensor(lcg, inverted);

computeStrainInformation(strain);
computeStrainInformation(strain, inverted);
ASSERT_EQ(strain.J < 0, inverted);

StrainInformation<defaulttype::Vec3Types>::MatrixSym PK2;
material.deriveSPKTensor(&strain, materialParameters, PK2);
Expand Down Expand Up @@ -167,7 +193,8 @@ void testSecondPiolaKirchhoffFromStrainEnergyDensityFunction(

void testElasticityTensorFromSecondPiolaKirchhoff(
HyperelasticMaterial<defaulttype::Vec3Types>& material,
const MaterialParameters<defaulttype::Vec3Types>& materialParameters)
const MaterialParameters<defaulttype::Vec3Types>& materialParameters,
bool inverted = false)
{
using Real = defaulttype::Vec3Types::Coord::value_type;
using MatrixSym = StrainInformation<defaulttype::Vec3Types>::MatrixSym;
Expand All @@ -177,9 +204,10 @@ void testElasticityTensorFromSecondPiolaKirchhoff(

//random right Cauchy-Green deformation tensor
sofa::testing::LinearCongruentialRandomGenerator lcg(96547);
strain.deformationTensor = generatePositiveDefiniteMatrix(lcg);
strain.deformationTensor = generateTestDeformationTensor(lcg, inverted);

computeStrainInformation(strain);
computeStrainInformation(strain, inverted);
ASSERT_EQ(strain.J < 0, inverted);

Matrix6 elasticityTensor;
material.ElasticityTensor(&strain, materialParameters, elasticityTensor);
Expand Down Expand Up @@ -212,7 +240,8 @@ void testElasticityTensorFromSecondPiolaKirchhoff(

void testApplyElasticityTensor(
HyperelasticMaterial<defaulttype::Vec3Types>& material,
const MaterialParameters<defaulttype::Vec3Types>& materialParameters)
const MaterialParameters<defaulttype::Vec3Types>& materialParameters,
bool inverted = false)
{
using MatrixSym = StrainInformation<defaulttype::Vec3Types>::MatrixSym;
using Matrix6 = HyperelasticMaterial<defaulttype::Vec3Types>::Matrix6;
Expand All @@ -221,7 +250,7 @@ void testApplyElasticityTensor(

//random right Cauchy-Green deformation tensor
sofa::testing::LinearCongruentialRandomGenerator lcg(96547);
strain.deformationTensor = generatePositiveDefiniteMatrix(lcg);
strain.deformationTensor = generateTestDeformationTensor(lcg, inverted);

//another random-generated symmetric second-order tensor
const StrainInformation<defaulttype::Vec3Types>::MatrixSym e =
Expand All @@ -232,7 +261,8 @@ void testApplyElasticityTensor(
eVec[i] = e[i];
}

computeStrainInformation(strain);
computeStrainInformation(strain, inverted);
ASSERT_EQ(strain.J < 0, inverted);

//first method to obtain the application of the elasticity tensor: matrix-vector product
Matrix6 elasticityTensor;
Expand Down Expand Up @@ -313,6 +343,40 @@ TEST(HyperelasticMaterial, PK2_StableNeoHookean)
testSecondPiolaKirchhoffFromStrainEnergyDensityFunction(material, materialParameters);
}

TEST(HyperelasticMaterial, PK2_StableNeoHookean_Inverted)
{
MaterialParameters<defaulttype::Vec3Types> materialParameters;
materialParameters.parameterArray = { 1., 1.};

StableNeoHookean<defaulttype::Vec3Types> material{};
testSecondPiolaKirchhoffFromStrainEnergyDensityFunction(material, materialParameters, true);
}

TEST(HyperelasticMaterial, StableNeoHookean_Reflection)
{
MaterialParameters<defaulttype::Vec3Types> materialParameters;
materialParameters.parameterArray = { 1., 1.};

StableNeoHookean<defaulttype::Vec3Types> material{};
StrainInformation<defaulttype::Vec3Types> strain;
strain.deformationTensor.identity();
computeStrainInformation(strain);
const auto restEnergy = material.getStrainEnergy(&strain, materialParameters);

// F = diag(1, 1, -1) has the same C as the rest state, but J = -1.
strain.J = -1;
StrainInformation<defaulttype::Vec3Types>::MatrixSym stress;
material.deriveSPKTensor(&strain, materialParameters, stress);

EXPECT_DOUBLE_EQ(material.getStrainEnergy(&strain, materialParameters) - restEnergy, 6.);
for (sofa::Size i = 0; i < 3; ++i)
{
for (sofa::Size j = 0; j <= i; ++j)
{
EXPECT_DOUBLE_EQ(stress(i, j), i == j ? 6. : 0.);
}
}
}


TEST(HyperelasticMaterial, ElasticityTensor_StVenantKirchhoff)
Expand Down Expand Up @@ -351,6 +415,15 @@ TEST(HyperelasticMaterial, ElasticityTensor_StableNeoHookean)
testElasticityTensorFromSecondPiolaKirchhoff(material, materialParameters);
}

TEST(HyperelasticMaterial, ElasticityTensor_StableNeoHookean_Inverted)
{
MaterialParameters<defaulttype::Vec3Types> materialParameters;
materialParameters.parameterArray = { 1., 1.};

StableNeoHookean<defaulttype::Vec3Types> material{};
testElasticityTensorFromSecondPiolaKirchhoff(material, materialParameters, true);
}

TEST(HyperelasticMaterial, ElasticityTensor_ArrudaBoyce)
{
MaterialParameters<defaulttype::Vec3Types> materialParameters;
Expand Down Expand Up @@ -409,6 +482,15 @@ TEST(HyperelasticMaterial, applyElasticityTensor_StableNeoHookean)
testApplyElasticityTensor(material, materialParameters);
}

TEST(HyperelasticMaterial, applyElasticityTensor_StableNeoHookean_Inverted)
{
MaterialParameters<defaulttype::Vec3Types> materialParameters;
materialParameters.parameterArray = { 1., 1.};

StableNeoHookean<defaulttype::Vec3Types> material{};
testApplyElasticityTensor(material, materialParameters, true);
}

TEST(HyperelasticMaterial, applyElasticityTensor_ArrudaBoyce)
{
MaterialParameters<defaulttype::Vec3Types> materialParameters;
Expand Down
17 changes: 10 additions & 7 deletions ...les/Component/SolidMechanics/FEM/TetrahedronHyperelasticityFEMForceField_invertedTets.scn
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
<?xml version="1.0" ?>
<?xml version="1.0" ?>
<Node name="root" dt="0.00005" showBoundingTree="0" gravity="0 0 0">
<Node name="plugins">
<RequiredPlugin pluginName="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedProjectiveConstraint] -->
<RequiredPlugin pluginName="Sofa.Component.Engine.Select"/> <!-- Needed to use components [BoxROI] -->
<RequiredPlugin pluginName="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [SparseLDLSolver] -->
<RequiredPlugin pluginName="Sofa.Component.LinearSystem"/> <!-- Needed to use components [ConstantSparsityPatternSystem] -->
<RequiredPlugin pluginName="Sofa.Component.Mass"/> <!-- Needed to use components [MeshMatrixMass] -->
<RequiredPlugin pluginName="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
<RequiredPlugin pluginName="Sofa.Component.IntegrationScheme.Backward"/> <!-- Needed to use components [EulerImplicitIntegrationScheme] -->
<RequiredPlugin pluginName="Sofa.Component.SolidMechanics.FEM.HyperElastic"/> <!-- Needed to use components [TetrahedronHyperelasticityFEMForceField] -->
<RequiredPlugin pluginName="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
Expand All @@ -25,21 +25,24 @@
<SparseLDLSolver template="CompressedRowSparseMatrixd"/>

<RegularGridTopology name="hexaGrid" min="0 0 0" max="1 1 2.7" n="6 6 16" p0="0 0 0"/>
<RegularGridTopology name="hexaGridRest" min="0 0 0" max="1 1 -2.7" n="6 6 16" p0="0 0 0"/>
<Node name="invertedConfiguration">
<RegularGridTopology name="hexaGrid" min="0 0 0" max="1 1 -2.7" n="6 6 16" p0="0 0 0"/>
</Node>

<MechanicalObject name="mechObj" rest_position="@hexaGrid.position" position="@hexaGridRest.position"/>
<MeshMatrixMass totalMass="1.0"/>
<MechanicalObject name="mechObj" rest_position="@hexaGrid.position" position="@invertedConfiguration/hexaGrid.position"/>
<!-- Preserve the uniform nodal mass used by the original inverted-element example. -->
<UniformMass totalMass="1.0"/>

<Node name="tetras">
<TetrahedronSetTopologyContainer name="Container"/>
<TetrahedronSetTopologyModifier name="Modifier" />
<TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />
<Hexa2TetraTopologicalMapping name="default28" input="@../hexaGrid" output="@Container" swapping="true"/>

<TetrahedronHyperelasticityFEMForceField name="FEM" ParameterSet="1644295.30201342 33557.0469798658" materialName="StableNeoHookean"/>
<TetrahedronHyperelasticityFEMForceField name="FEM" topology="@Container" ParameterSet="1644295.30201342 33557.0469798658" materialName="StableNeoHookean"/>
</Node>

<BoxROI drawBoxes="1" box="0 0 0 1 1 0.05" name="box"/>
<BoxROI drawBoxes="1" box="0 0 0 1 1 0.05" name="box" position="@mechObj.rest_position"/>
<FixedProjectiveConstraint indices="@box.indices"/>
</Node>

Expand Down
Loading