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

Open
lililii124 wants to merge 1 commit into
base: master
Choose a base branch
from
Open

[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
1 commit
Select commit
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
6 changes: 6 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,11 @@
# SOFA Changelog

## Unreleased

### Bug fixes

- [HyperElastic] Fix the Stable Neo-Hookean example by isolating the inverted configuration's topology, computing mass on the tetrahedral rest mesh, and reducing its time step. Add scene regression coverage and material derivative tests for inverted deformations ([#6220](https://github.com/sofa-framework/sofa/issues/6220)).



## [v26.06.00]( https://github.com/sofa-framework/sofa/tree/v26.06.00 )
Expand Down
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
15 changes: 9 additions & 6 deletions ...les/Component/SolidMechanics/FEM/TetrahedronHyperelasticityFEMForceField_invertedTets.scn
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
<?xml version="1.0" ?>
<Node name="root" dt="0.00005" showBoundingTree="0" gravity="0 0 0">
<Node name="root" dt="0.00001" 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] -->
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"/>

<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"/>
<MeshMatrixMass totalMass="1.0" topology="@Container" geometryState="@../mechObj" computeMassOnRest="true"/>

<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
1 change: 1 addition & 0 deletions examples/RegressionStateScenes.regression-tests
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,7 @@ Component/SolidMechanics/FEM/TetrahedronFEMForceField.scn 300 1e-4 0 1
Component/SolidMechanics/FEM/TetrahedronFEMForceField_assemble.scn 300 1e-4 0 1
Component/SolidMechanics/FEM/TetrahedronFEMForceField_plasticity.scn 300 1e-4 0 1
Component/SolidMechanics/FEM/TetrahedronHyperelasticityFEMForceField.scn 100 1e-4 0 1
Component/SolidMechanics/FEM/TetrahedronHyperelasticityFEMForceField_invertedTets.scn 100 1e-4 0 0
Component/SolidMechanics/FEM/TriangularFEMForceFieldOptim.scn 100 1e-4 0 1
Component/SolidMechanics/Spring/FastTriangularBendingSprings.scn 100 1e-4 0 1
Component/SolidMechanics/Spring/MeshSpringForceField.scn 300 1e-4 0 1
Expand Down
Loading