[FEM.Elastic] Introduce generic element-agnostic elastic FEM force field by alxbilger · Pull Request #5882 · sofa-framework/sofa

alxbilger · 2026-01-22T15:43:11Z

Coming from https://github.com/alxbilger/Elasticity

TODO:

add tests
add scenes

By submitting this pull request, I acknowledge that
I have read, understand, and agree SOFA Developer Certificate of Origin (DCO).

Reviewers will merge this pull-request only if

it builds with SUCCESS for all platforms on the CI.
it does not generate new warnings.
it does not generate new unit test failures.
it does not generate new scene test failures.
it does not break API compatibility.
it is more than 1 week old (or has fast-merge label).

fredroy · 2026-01-28T23:11:54Z

[ci-build][with-all-tests]

th-skam

I only touched some easy to read files. Will try to take a better look in the next days.
But I really like the overall work.

...nent/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VecView.h

hugtalbot · 2026-02-06T13:11:00Z

This PR deserves the attention of all @sofa-framework/reviewers

th-skam · 2026-02-09T10:12:41Z

...anics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/FullySymmetric4Tensor.h

+                    for (sofa::Size l = 0; l < spatial_dimensions; ++l)
+                    {
+                        const auto ijkl = callable(i, j, k, l);
+                        const auto jikl = callable(i, j, k, l);


Suggested change

const auto jikl = callable(i, j, k, l);

const auto jikl = callable(j, i, k, l);

th-skam · 2026-02-09T11:15:47Z

...stic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Hexahedron].h

+        return gradient;
+    }
+
+    static constexpr std::array<QuadraturePointAndWeight, 4> quadraturePoints()


Shouldn't we be using 8-point rule for hexas? And similarly 4-point rule for quads instead of 3 in the respective file?
I guess it's a performance/approximative trade-off but if the work coming from Elasticity is to be used for didacticism, perhaps we should have some comment here and mention that it is approximative.

...EM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h

th-skam · 2026-02-09T12:24:05Z

...anics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/FullySymmetric4Tensor.h

+ * @param youngModulus The Young's modulus of the material, representing its stiffness.
+ * @param poissonRatio The Poisson's ratio of the material, describing its deformation behavior.
+ * @return The isotropic elasticity tensor represented as an object of
+ * `ElasticityTensor<DataTypes>`.


Suggested change

* `ElasticityTensor<DataTypes>`.

* `FullySymmetric4Tensor<DataTypes>`.

...lidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/KroneckerDelta.h

...lidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/LameParameters.h

...anics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MajorSymmetric4Tensor.h

...echanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/StrainDisplacement.h

th-skam · 2026-02-09T16:51:55Z

...echanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/StrainDisplacement.h

+        auto row = spatial_dimensions;
+        for (sofa::Size i = 0; i < spatial_dimensions; ++i)
+        {
+            for (sofa::Size j = i + 1; j < spatial_dimensions; ++j)
+            {
+                B(row, ne * spatial_dimensions + i) = gradientShapeFunctions[ne][j];
+                B(row, ne * spatial_dimensions + j) = gradientShapeFunctions[ne][i];
+                ++row;
+            }
+        }


Suggested change

auto row = spatial_dimensions;

for (sofa::Size i = 0; i < spatial_dimensions; ++i)

{

for (sofa::Size j = i + 1; j < spatial_dimensions; ++j)

{

B(row, ne * spatial_dimensions + i) = gradientShapeFunctions[ne][j];

B(row, ne * spatial_dimensions + j) = gradientShapeFunctions[ne][i];

++row;

}

}

for (sofa::Size i = 0; i < spatial_dimensions; ++i)

{

for (sofa::Size j = i + 1; j < spatial_dimensions; ++j)

{

const auto row = tensorToVoigtIndex<DataTypes>(i, j);

B(row, ne * spatial_dimensions + i) = gradientShapeFunctions[ne][j];

B(row, ne * spatial_dimensions + j) = gradientShapeFunctions[ne][i];

}

}

There's a mismatch here between what's in VoigtNotation.h which states:
For 3D: (0,0)->0, (1,1)->1, (2,2)->2, (1,2)->3, (0,2)->4, (0,1)->5

But if you unroll the loop here, you would get
row = 3 -> (0,1), 4 -> (0,2), 5 -> (1,2) so 3 and 5 from voigt are swapped.

So the B^T C B in ElementStiffnessMatrix.h mixes them up.
And I assume we're getting away with it because we're dealing with isotropic problems but if the same code is to be used in the future we won't. Better to have a test to verify this now.

bakpaul · 2026-02-13T10:53:45Z

...stic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/StablePolarDecomposition.h

+
+        const auto H = P.multTranspose(Q);
+
+        sofa::helper::Decompose<sofa::Real_t<DataTypes>>::polarDecomposition_stable(H, rotationMatrix);


What's the diff with the regular one ? Why keep both ?

No idea. As usual, I want to keep compatibility with legacy SOFA. TetrahedronFEMForceField offers the option to select both polar decomposition. SVD uses the stable polar decomposition.

bakpaul · 2026-02-13T11:05:01Z

...anics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/FullySymmetric4Tensor.h

+                        const auto klij = callable(k, l, i, j);
+                        const auto ijlk = callable(i, j, l, k);


Isn't the second minor symetry missing ? That being said, testing the first one and the major one seems sufficient.

.../SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MatrixTools.h

bakpaul · 2026-02-13T11:08:44Z

.../SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MatrixTools.h

This should be purely added to the Core Matrix API

hugtalbot · 2026-02-24T08:10:30Z

I insist on the fact that this PR deserves the attention of all @sofa-framework/reviewers

alxbilger added pr: status to review To notify reviewers to review this pull-request pr: new feature Implement a new feature pr: highlighted in next release Highlight this contribution in the notes of the upcoming release labels Jan 22, 2026

hugtalbot added the topic for next dev-meeting PR to be discussed in sofa-dev meeting label Jan 23, 2026

th-skam reviewed Feb 2, 2026

View reviewed changes

...nent/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VecView.h Show resolved Hide resolved

alxbilger removed the topic for next dev-meeting PR to be discussed in sofa-dev meeting label Feb 5, 2026

th-skam requested changes Feb 9, 2026

View reviewed changes

bakpaul reviewed Feb 13, 2026

View reviewed changes

alxbilger added 6 commits March 2, 2026 15:18

[FEM.Elastic] Introduce generic element-agnostic elastic FEM force field

43e25bd

add header

f0dec93

register in the factory

ae89822

replace not supported iota_view with the SOFA alternative

7d3dc69

missing typename

e4ec49d

missing typename

ca8639b

alxbilger force-pushed the elasticity branch from fa404e6 to ca8639b Compare March 2, 2026 14:20

alxbilger added 10 commits March 2, 2026 15:56

const ref

4ba00b0

const ref

cf1d54c

replace include for cstddef

f9a3f77

missing listing of MajorSymmetric4Tensor in CMakeLists.txt

f9b2244

fix wrong namespace

3366077

make MajorSymmetric4Tensor constructor explicit

bd2baa0

MajorSymmetric4Tensor is not used

25c345e

remove determinantSquareMatrix as it is now integrated into SOFA

6ec7552

function can accept reference

00ba53c

result must be initialized before being used

32d7a96

	const auto jikl = callable(i, j, k, l);
	const auto jikl = callable(j, i, k, l);

	* `ElasticityTensor<DataTypes>`.
	* `FullySymmetric4Tensor<DataTypes>`.


		const auto H = P.multTranspose(Q);

		sofa::helper::Decompose<sofa::Real_t<DataTypes>>::polarDecomposition_stable(H, rotationMatrix);

		const auto klij = callable(k, l, i, j);
		const auto ijlk = callable(i, j, l, k);

Conversation

alxbilger commented Jan 22, 2026 • edited by th-skam Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

fredroy commented Jan 28, 2026

Uh oh!

th-skam left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

hugtalbot commented Feb 6, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hugtalbot commented Feb 24, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

alxbilger commented Jan 22, 2026 •

edited by th-skam

Loading