From a4c29967462b525d1622f292c08c03533f365068 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 10 Nov 2026 10:15:16 +0100 Subject: [PATCH 1/7] Turn mu and alpha into vectors --- .../hyperelasticmaterials/OgdenMaterial.h | 6 ++- .../hyperelasticmaterials/OgdenMaterial.inl | 39 +++++++++++++++---- 2 files changed, 35 insertions(+), 10 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.h index 14a0d2ae033..13e5ea61f34 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.h +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.h @@ -50,10 +50,12 @@ class OgdenMaterial : public: // Ogden material constants - sofa::Data m_mu; - sofa::Data m_alpha; + sofa::Data> m_mu; + sofa::Data> m_alpha; sofa::Data m_kappa; + void init() override; + StressTensor secondPiolaKirchhoffStress(Strain& strain) override; ElasticityTensor elasticityTensor(Strain& strain) override; diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl index dc8f8e5aba4..5e55b02eb65 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl @@ -24,6 +24,7 @@ #include #include #include +#include #include namespace sofa::component::solidmechanics::fem::hyperelastic @@ -31,14 +32,36 @@ namespace sofa::component::solidmechanics::fem::hyperelastic template OgdenMaterial::OgdenMaterial() - : m_mu(initData(&m_mu, static_cast(1e2), "mu", + : m_mu(initData(&m_mu, sofa::type::vector{static_cast(1e2)}, "mu", "Material constant relevant to the shear modulus")) - , m_alpha(initData(&m_alpha, static_cast(1.5), "alpha", + , m_alpha(initData(&m_alpha, sofa::type::vector{static_cast(1.5)}, "alpha", "Material constant exponent related to strain-stiffening")) , m_kappa(initData(&m_kappa, static_cast(1e3), "kappa", "Material constant related to the bulk modulus")) {} +template +void OgdenMaterial::init() +{ + const auto& mu = m_mu.getValue(); + const auto& alpha = m_alpha.getValue(); + + if (mu.empty() || mu.size() != alpha.size()) + { + msg_error() << "mu and alpha must be non-empty and of the same size (got " << mu.size() + << " and " << alpha.size() << ")"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } + + if (std::any_of(alpha.begin(), alpha.end(), [](Real a) { return a == static_cast(0); })) + { + msg_error() << "alpha values must be non-zero"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } + + PK2HyperelasticMaterial::init(); +} + template auto OgdenMaterial::secondPiolaKirchhoffStress(Strain& strain) -> StressTensor { @@ -52,8 +75,8 @@ auto OgdenMaterial::secondPiolaKirchhoffStress(Strain& str const auto S_isochoric = [this, J, &strain, &C_1, &C] () { - const Real mu = m_mu.getValue(); - const Real alpha = m_alpha.getValue(); + const Real mu = m_mu.getValue()[0]; + const Real alpha = m_alpha.getValue()[0]; const Real FJ = pow(J, -alpha / static_cast(spatial_dimensions)); @@ -112,8 +135,8 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El { using EigenMatrix = Eigen::Matrix; - const Real mu = m_mu.getValue(); - const Real alpha = m_alpha.getValue(); + const Real mu = m_mu.getValue()[0]; + const Real alpha = m_alpha.getValue()[0]; const Real kappa = m_kappa.getValue(); const auto& C = strain.getRightCauchyGreenTensor(); @@ -223,8 +246,8 @@ auto OgdenMaterial::strainEnergyDensity(Strain& strain) -> { using EigenMatrix = Eigen::Matrix; - const Real mu = m_mu.getValue(); - const Real alpha = m_alpha.getValue(); + const Real mu = m_mu.getValue()[0]; + const Real alpha = m_alpha.getValue()[0]; const Real kappa = m_kappa.getValue(); const auto& C = strain.getRightCauchyGreenTensor(); From 2c18674b3f650475194b63eb8ece4e3dd8d56170 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 10 Nov 2026 10:23:08 +0100 Subject: [PATCH 2/7] accumulate strain density in a loop --- .../hyperelasticmaterials/OgdenMaterial.inl | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl index 5e55b02eb65..ceff5aec8de 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl @@ -246,8 +246,8 @@ auto OgdenMaterial::strainEnergyDensity(Strain& strain) -> { using EigenMatrix = Eigen::Matrix; - const Real mu = m_mu.getValue()[0]; - const Real alpha = m_alpha.getValue()[0]; + const auto& mu = m_mu.getValue(); + const auto& alpha = m_alpha.getValue(); const Real kappa = m_kappa.getValue(); const auto& C = strain.getRightCauchyGreenTensor(); @@ -265,13 +265,17 @@ auto OgdenMaterial::strainEnergyDensity(Strain& strain) -> dmsg_warning("OgdenMaterial") << "EigenSolver iterations failed to converge"; const auto eigenvalues = EigenProblemSolver.eigenvalues(); - const Real aBy2 = static_cast(0.5) * alpha; - Real trCaBy2{static_cast(0)}; - for (sofa::Index n = 0; n < spatial_dimensions; ++n) - trCaBy2 += pow(eigenvalues(n), aBy2); + Real W_isochoric{static_cast(0)}; + for (sofa::Index p = 0; p < mu.size(); ++p) + { + const Real aBy2 = static_cast(0.5) * alpha[p]; + Real trCaBy2{static_cast(0)}; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + trCaBy2 += pow(eigenvalues(n), aBy2); - const Real FJ = pow(J, -alpha / static_cast(spatial_dimensions)); - const Real W_isochoric = mu / (alpha * alpha) * (FJ * trCaBy2 - static_cast(spatial_dimensions)); + const Real FJ = pow(J, -alpha[p] / static_cast(spatial_dimensions)); + W_isochoric += mu[p] / (alpha[p] * alpha[p]) * (FJ * trCaBy2 - static_cast(spatial_dimensions)); + } const Real logJ = log(J); const Real W_volumetric = static_cast(0.5) * kappa * logJ * logJ; From d60e6d9cf5f4fa5168fe1b241ac81474228410f6 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 10 Nov 2026 10:26:03 +0100 Subject: [PATCH 3/7] turn SPK into a loop --- .../hyperelasticmaterials/OgdenMaterial.inl | 49 ++++++++++--------- 1 file changed, 26 insertions(+), 23 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl index ceff5aec8de..1ee3f972494 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl @@ -75,10 +75,8 @@ auto OgdenMaterial::secondPiolaKirchhoffStress(Strain& str const auto S_isochoric = [this, J, &strain, &C_1, &C] () { - const Real mu = m_mu.getValue()[0]; - const Real alpha = m_alpha.getValue()[0]; - - const Real FJ = pow(J, -alpha / static_cast(spatial_dimensions)); + const auto& mu = m_mu.getValue(); + const auto& alpha = m_alpha.getValue(); EigenMatrix CEigen; for (sofa::Index m = 0; m < spatial_dimensions; ++m) @@ -92,29 +90,34 @@ auto OgdenMaterial::secondPiolaKirchhoffStress(Strain& str const auto eigenvectors = EigenProblemSolver.eigenvectors(); const auto eigenvalues = EigenProblemSolver.eigenvalues(); - // Precompute trace of C^(alpha/2) and eigenvalue powers: lambda^(alpha/2 - 1) - const Real aBy2 = static_cast(0.5) * alpha; - const Real aBy2Minus1 = aBy2 - static_cast(1); - Real trCaBy2{static_cast(0)}; - typename DataTypes::Coord eigenvaluePowers{}; - for (sofa::Index n = 0; n < spatial_dimensions; ++n) - { - trCaBy2 += pow(eigenvalues(n), aBy2); - eigenvaluePowers(n) = pow(eigenvalues(n), aBy2Minus1); - } - - const Real coefS1 = FJ * mu / (alpha * static_cast(2)); - const Real coefS2 = -FJ * mu / (alpha * static_cast(2 * spatial_dimensions)) * trCaBy2; - StressTensor S{}; - for (sofa::Index i = 0; i < spatial_dimensions; ++i) + for (sofa::Index p = 0; p < mu.size(); ++p) { - for (sofa::Index j = 0; j < spatial_dimensions; ++j) + const Real FJ = pow(J, -alpha[p] / static_cast(spatial_dimensions)); + + // Precompute trace of C^(alpha/2) and eigenvalue powers: lambda^(alpha/2 - 1) + const Real aBy2 = static_cast(0.5) * alpha[p]; + const Real aBy2Minus1 = aBy2 - static_cast(1); + Real trCaBy2{static_cast(0)}; + typename DataTypes::Coord eigenvaluePowers{}; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) { - for (sofa::Index n = 0; n < spatial_dimensions; ++n) - S(i,j) += coefS1 * eigenvaluePowers(n) * eigenvectors(i,n) * eigenvectors(j,n); - S(i,j) += coefS2 * C_1(j, i); + trCaBy2 += pow(eigenvalues(n), aBy2); + eigenvaluePowers(n) = pow(eigenvalues(n), aBy2Minus1); + } + + const Real coefS1 = FJ * mu[p] / (alpha[p] * static_cast(2)); + const Real coefS2 = -FJ * mu[p] / (alpha[p] * static_cast(2 * spatial_dimensions)) * trCaBy2; + + for (sofa::Index i = 0; i < spatial_dimensions; ++i) + { + for (sofa::Index j = 0; j < spatial_dimensions; ++j) + { + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + S(i,j) += coefS1 * eigenvaluePowers(n) * eigenvectors(i,n) * eigenvectors(j,n); + S(i,j) += coefS2 * C_1(j, i); + } } } From ba60da46e8c1060c8fa21b651aba167a1f35ba33 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 10 Nov 2026 11:00:06 +0100 Subject: [PATCH 4/7] Loop over the elasticity tensor --- .../hyperelasticmaterials/OgdenMaterial.inl | 126 +++++++++--------- 1 file changed, 65 insertions(+), 61 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl index 1ee3f972494..5aa4e2c46d0 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl @@ -138,8 +138,8 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El { using EigenMatrix = Eigen::Matrix; - const Real mu = m_mu.getValue()[0]; - const Real alpha = m_alpha.getValue()[0]; + const auto& mu = m_mu.getValue(); + const auto& alpha = m_alpha.getValue(); const Real kappa = m_kappa.getValue(); const auto& C = strain.getRightCauchyGreenTensor(); @@ -147,7 +147,6 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El const Real J = strain.getDeterminantDeformationGradient(); assert(J > 0); - const Real FJ = pow(J, -alpha / static_cast(spatial_dimensions)); EigenMatrix CEigen; for (sofa::Index m = 0; m < spatial_dimensions; ++m) @@ -160,31 +159,6 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El const auto eigenvectors = EigenProblemSolver.eigenvectors(); const auto eigenvalues = EigenProblemSolver.eigenvalues(); - // Precompute trace of C^(alpha/2) and eigenvalue powers: lambda^(a/2 - 1) and lambda^(a/2 - 2) - const Real aBy2 = static_cast(0.5) * alpha; - const Real aBy2Minus1 = aBy2 - static_cast(1); - Real trCaBy2{static_cast(0)}; - typename DataTypes::Coord eigenvaluePowers1{}, eigenvaluePowers2{}; - for (sofa::Index n = 0; n < spatial_dimensions; ++n) - { - trCaBy2 += pow(eigenvalues(n), aBy2); - eigenvaluePowers1(n) = pow(eigenvalues(n), aBy2Minus1); - eigenvaluePowers2(n) = pow(eigenvalues(n), aBy2Minus1 - static_cast(1)); - } - - const Real coefSpectral = FJ * mu / alpha * aBy2Minus1; - const Real coefRotational = FJ * mu / alpha; - - // Precompute C^(alpha1/2 - 1) from eigenbasis: V * D * V^T; D_i = lambda_i^(alpha1/2 - 1) - RightCauchyGreenTensor D, EigenBasis; - for (sofa::Index n = 0; n < spatial_dimensions; ++n) - { - for (sofa::Index m = 0; m < spatial_dimensions; ++m) - EigenBasis(n,m) = eigenvectors(n,m); - D(n,n) = eigenvaluePowers1(n); - } - const RightCauchyGreenTensor CaBy2Minus1(EigenBasis*D*(EigenBasis.transposed())); - return ElasticityTensor( [&](sofa::Index i, sofa::Index j, sofa::Index k, sofa::Index l) { @@ -193,44 +167,74 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El { Real sum{static_cast(0)}; - // Spectral decomposition contributions; derivatives of eigenvalues and eigenvectors - for (sofa::Index n = 0; n < spatial_dimensions; ++n) + for (sofa::Index p = 0; p < mu.size(); ++p) { - sum += coefSpectral * eigenvaluePowers2(n) * eigenvectors(i,n) * eigenvectors(j,n) - * eigenvectors(k,n) * eigenvectors(l,n); + const Real FJ = pow(J, -alpha[p] / static_cast(spatial_dimensions)); - for (sofa::Index m = 0; m < spatial_dimensions; ++m) + // Precompute trace of C^(alpha/2) and eigenvalue powers: lambda^(a/2 - 1) and lambda^(a/2 - 2) + const Real aBy2 = static_cast(0.5) * alpha[p]; + const Real aBy2Minus1 = aBy2 - static_cast(1); + Real trCaBy2{static_cast(0)}; + typename DataTypes::Coord eigenvaluePowers1{}, eigenvaluePowers2{}; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) { - // This sum is taken over only for m != n - if (m == n) continue; - - // Eigenvalue multiplicity (val1 == val2) causes an indeterminate form; - // In that case, the limit of the expression is used - constexpr Real tol = static_cast(1e-8); - const Real maxEval = std::max(std::fabs(eigenvalues(n)), std::fabs(eigenvalues(m))); - const bool isDegenerate = std::fabs(eigenvalues(n) - eigenvalues(m)) < maxEval * tol; - const Real eigenvalueMultiplier = isDegenerate - ? aBy2Minus1 * eigenvaluePowers2(n) - : (eigenvaluePowers1(n) - eigenvaluePowers1(m)) / (eigenvalues(n) - eigenvalues(m)); - - sum += static_cast(0.5) * coefRotational * eigenvalueMultiplier * - (eigenvectors(i,n) * eigenvectors(j,m) * eigenvectors(k,m) * eigenvectors(l,n) - + eigenvectors(i,n) * eigenvectors(j,m) * eigenvectors(k,n) * eigenvectors(l,m)); + trCaBy2 += pow(eigenvalues(n), aBy2); + eigenvaluePowers1(n) = pow(eigenvalues(n), aBy2Minus1); + eigenvaluePowers2(n) = pow(eigenvalues(n), aBy2Minus1 - static_cast(1)); } - } - // Remaining terms grouped to highlight respective contribution from - // differentiating coefS1, coefS2 and C^(-1) from S_isochoric - sum += FJ * mu * ( - // Contribution from derivative of F(J) - trCaBy2 / static_cast(2 * spatial_dimensions * spatial_dimensions) * C_1(j,i) * C_1(l,k) - - // Contribution from derivatives of C^(alpha/2 - 1) and tr(C^(alpha/2)) - (CaBy2Minus1(j,i) * C_1(l,k) + C_1(j,i) * CaBy2Minus1(k,l)) - / static_cast(2 * spatial_dimensions) + - // Contribution from derivative of C^(-1) - trCaBy2 / (static_cast(2 * spatial_dimensions) * alpha) - * (C_1(j,k) * C_1(l,i) + C_1(j,l) * C_1(k,i)) - ); + const Real coefSpectral = FJ * mu[p] / alpha[p] * aBy2Minus1; + const Real coefRotational = FJ * mu[p] / alpha[p]; + + // Precompute C^(alpha1/2 - 1) from eigenbasis: V * D * V^T; D_i = lambda_i^(alpha1/2 - 1) + RightCauchyGreenTensor D, EigenBasis; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + { + for (sofa::Index m = 0; m < spatial_dimensions; ++m) + EigenBasis(n,m) = eigenvectors(n,m); + D(n,n) = eigenvaluePowers1(n); + } + const RightCauchyGreenTensor CaBy2Minus1(EigenBasis*D*(EigenBasis.transposed())); + + // Spectral decomposition contributions; derivatives of eigenvalues and eigenvectors + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + { + sum += coefSpectral * eigenvaluePowers2(n) * eigenvectors(i,n) * eigenvectors(j,n) + * eigenvectors(k,n) * eigenvectors(l,n); + + for (sofa::Index m = 0; m < spatial_dimensions; ++m) + { + // This sum is taken over only for m != n + if (m == n) continue; + + // Eigenvalue multiplicity (val1 == val2) causes an indeterminate form; + // In that case, the limit of the expression is used + constexpr Real tol = static_cast(1e-8); + const Real maxEval = std::max(std::fabs(eigenvalues(n)), std::fabs(eigenvalues(m))); + const bool isDegenerate = std::fabs(eigenvalues(n) - eigenvalues(m)) < maxEval * tol; + const Real eigenvalueMultiplier = isDegenerate + ? aBy2Minus1 * eigenvaluePowers2(n) + : (eigenvaluePowers1(n) - eigenvaluePowers1(m)) / (eigenvalues(n) - eigenvalues(m)); + + sum += static_cast(0.5) * coefRotational * eigenvalueMultiplier * + (eigenvectors(i,n) * eigenvectors(j,m) * eigenvectors(k,m) * eigenvectors(l,n) + + eigenvectors(i,n) * eigenvectors(j,m) * eigenvectors(k,n) * eigenvectors(l,m)); + } + } + + // Remaining terms grouped to highlight respective contribution from + // differentiating coefS1, coefS2 and C^(-1) from S_isochoric + sum += FJ * mu[p] * ( + // Contribution from derivative of F(J) + trCaBy2 / static_cast(2 * spatial_dimensions * spatial_dimensions) * C_1(j,i) * C_1(l,k) - + // Contribution from derivatives of C^(alpha/2 - 1) and tr(C^(alpha/2)) + (CaBy2Minus1(j,i) * C_1(l,k) + C_1(j,i) * CaBy2Minus1(k,l)) + / static_cast(2 * spatial_dimensions) + + // Contribution from derivative of C^(-1) + trCaBy2 / (static_cast(2 * spatial_dimensions) * alpha[p]) + * (C_1(j,k) * C_1(l,i) + C_1(j,l) * C_1(k,i)) + ); + } return sum; }(); From 53fc40ec1fbe75d5ecde197ac077c070f8ebaa6a Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 10 Nov 2026 11:07:15 +0100 Subject: [PATCH 5/7] Apply 3-term test in ogden --- .../FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp index 3e414acdeca..1138f10582f 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp @@ -85,6 +85,12 @@ class BaseMaterialTest : public sofa::testing::BaseTest linearMaterial->d_youngModulus.setValue(10.); } + if (auto* ogdenMaterial = dynamic_cast*>(material.get())) + { + ogdenMaterial->m_mu.setValue({80., 25., 40.}); + ogdenMaterial->m_alpha.setValue({1.3, 3.7, -1.8}); + } + material->init(); } From 6626d05df890dddd332548b3cdb05ca36e8bded2 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 10 Nov 2026 11:13:11 +0100 Subject: [PATCH 6/7] fix warnings in HyperelasticMaterials_test --- .../FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp index 1138f10582f..22459373ff7 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/tests/HyperelasticMaterials_test.cpp @@ -101,9 +101,9 @@ class BaseMaterialTest : public sofa::testing::BaseTest Real J; do { - for (int j = 0; j < spatial_dimensions; j++) + for (sofa::Size j = 0; j < spatial_dimensions; j++) { - for (int i = 0; i < spatial_dimensions; i++) + for (sofa::Size i = 0; i < spatial_dimensions; i++) { F(i, j) = lcg.generateInRange(-1., 1.); } @@ -307,7 +307,6 @@ class PK2MaterialTest : public BaseMaterialTest> const auto& C = strain.getRightCauchyGreenTensor(); - const auto S = this->material->secondPiolaKirchhoffStress(strain); const auto elasticityTensor = this->material->elasticityTensor(strain); // Small perturbation for finite difference From fbef6b20625b335cc9df3091882f0b2db0694640 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 10 Nov 2026 11:19:43 +0100 Subject: [PATCH 7/7] vectorial precomputation of shared terms --- .../hyperelasticmaterials/OgdenMaterial.inl | 76 ++++++++++--------- 1 file changed, 42 insertions(+), 34 deletions(-) diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl index 5aa4e2c46d0..4b33645c743 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl @@ -159,6 +159,38 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El const auto eigenvectors = EigenProblemSolver.eigenvectors(); const auto eigenvalues = EigenProblemSolver.eigenvalues(); + const auto nbTerms = mu.size(); + sofa::type::vector FJ(nbTerms), aBy2Minus1(nbTerms), trCaBy2(nbTerms); + sofa::type::vector eigenvaluePowers1(nbTerms), eigenvaluePowers2(nbTerms); + sofa::type::vector CaBy2Minus1(nbTerms); + + RightCauchyGreenTensor EigenBasis; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + for (sofa::Index m = 0; m < spatial_dimensions; ++m) + EigenBasis(n,m) = eigenvectors(n,m); + + // Precompute terms that are shared over the upcoming ijkl loop + for (sofa::Index p = 0; p < nbTerms; ++p) + { + FJ[p] = pow(J, -alpha[p] / static_cast(spatial_dimensions)); + + // Trace of C^(alpha/2) and eigenvalue powers: lambda^(a/2 - 1) and lambda^(a/2 - 2) + const Real aBy2 = static_cast(0.5) * alpha[p]; + aBy2Minus1[p] = aBy2 - static_cast(1); + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + { + trCaBy2[p] += pow(eigenvalues(n), aBy2); + eigenvaluePowers1[p](n) = pow(eigenvalues(n), aBy2Minus1[p]); + eigenvaluePowers2[p](n) = pow(eigenvalues(n), aBy2Minus1[p] - static_cast(1)); + } + + // C^(alpha1/2 - 1) from eigenbasis: V * D * V^T; D_i = lambda_i^(alpha1/2 - 1) + RightCauchyGreenTensor D; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + D(n,n) = eigenvaluePowers1[p](n); + CaBy2Minus1[p] = EigenBasis*D*(EigenBasis.transposed()); + } + return ElasticityTensor( [&](sofa::Index i, sofa::Index j, sofa::Index k, sofa::Index l) { @@ -167,39 +199,15 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El { Real sum{static_cast(0)}; - for (sofa::Index p = 0; p < mu.size(); ++p) + for (sofa::Index p = 0; p < nbTerms; ++p) { - const Real FJ = pow(J, -alpha[p] / static_cast(spatial_dimensions)); - - // Precompute trace of C^(alpha/2) and eigenvalue powers: lambda^(a/2 - 1) and lambda^(a/2 - 2) - const Real aBy2 = static_cast(0.5) * alpha[p]; - const Real aBy2Minus1 = aBy2 - static_cast(1); - Real trCaBy2{static_cast(0)}; - typename DataTypes::Coord eigenvaluePowers1{}, eigenvaluePowers2{}; - for (sofa::Index n = 0; n < spatial_dimensions; ++n) - { - trCaBy2 += pow(eigenvalues(n), aBy2); - eigenvaluePowers1(n) = pow(eigenvalues(n), aBy2Minus1); - eigenvaluePowers2(n) = pow(eigenvalues(n), aBy2Minus1 - static_cast(1)); - } - - const Real coefSpectral = FJ * mu[p] / alpha[p] * aBy2Minus1; - const Real coefRotational = FJ * mu[p] / alpha[p]; - - // Precompute C^(alpha1/2 - 1) from eigenbasis: V * D * V^T; D_i = lambda_i^(alpha1/2 - 1) - RightCauchyGreenTensor D, EigenBasis; - for (sofa::Index n = 0; n < spatial_dimensions; ++n) - { - for (sofa::Index m = 0; m < spatial_dimensions; ++m) - EigenBasis(n,m) = eigenvectors(n,m); - D(n,n) = eigenvaluePowers1(n); - } - const RightCauchyGreenTensor CaBy2Minus1(EigenBasis*D*(EigenBasis.transposed())); + const Real coefSpectral = FJ[p] * mu[p] / alpha[p] * aBy2Minus1[p]; + const Real coefRotational = FJ[p] * mu[p] / alpha[p]; // Spectral decomposition contributions; derivatives of eigenvalues and eigenvectors for (sofa::Index n = 0; n < spatial_dimensions; ++n) { - sum += coefSpectral * eigenvaluePowers2(n) * eigenvectors(i,n) * eigenvectors(j,n) + sum += coefSpectral * eigenvaluePowers2[p](n) * eigenvectors(i,n) * eigenvectors(j,n) * eigenvectors(k,n) * eigenvectors(l,n); for (sofa::Index m = 0; m < spatial_dimensions; ++m) @@ -213,8 +221,8 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El const Real maxEval = std::max(std::fabs(eigenvalues(n)), std::fabs(eigenvalues(m))); const bool isDegenerate = std::fabs(eigenvalues(n) - eigenvalues(m)) < maxEval * tol; const Real eigenvalueMultiplier = isDegenerate - ? aBy2Minus1 * eigenvaluePowers2(n) - : (eigenvaluePowers1(n) - eigenvaluePowers1(m)) / (eigenvalues(n) - eigenvalues(m)); + ? aBy2Minus1[p] * eigenvaluePowers2[p](n) + : (eigenvaluePowers1[p](n) - eigenvaluePowers1[p](m)) / (eigenvalues(n) - eigenvalues(m)); sum += static_cast(0.5) * coefRotational * eigenvalueMultiplier * (eigenvectors(i,n) * eigenvectors(j,m) * eigenvectors(k,m) * eigenvectors(l,n) @@ -224,14 +232,14 @@ auto OgdenMaterial::elasticityTensor(Strain& strain) -> El // Remaining terms grouped to highlight respective contribution from // differentiating coefS1, coefS2 and C^(-1) from S_isochoric - sum += FJ * mu[p] * ( + sum += FJ[p] * mu[p] * ( // Contribution from derivative of F(J) - trCaBy2 / static_cast(2 * spatial_dimensions * spatial_dimensions) * C_1(j,i) * C_1(l,k) - + trCaBy2[p] / static_cast(2 * spatial_dimensions * spatial_dimensions) * C_1(j,i) * C_1(l,k) - // Contribution from derivatives of C^(alpha/2 - 1) and tr(C^(alpha/2)) - (CaBy2Minus1(j,i) * C_1(l,k) + C_1(j,i) * CaBy2Minus1(k,l)) + (CaBy2Minus1[p](j,i) * C_1(l,k) + C_1(j,i) * CaBy2Minus1[p](k,l)) / static_cast(2 * spatial_dimensions) + // Contribution from derivative of C^(-1) - trCaBy2 / (static_cast(2 * spatial_dimensions) * alpha[p]) + trCaBy2[p] / (static_cast(2 * spatial_dimensions) * alpha[p]) * (C_1(j,k) * C_1(l,i) + C_1(j,l) * C_1(k,i)) ); }