From b1621db2c1487e33243ea08d7d0b116220217a57 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Fri, 23 Oct 2026 16:41:53 +0200
Subject: [PATCH 01/11] Introduce tests for the CorotationalFEMForceField
 tangent stifness

---
 .../FEM/Elastic/tests/CMakeLists.txt          |   1 +
 .../tests/CorotationalFEMForceField_test.cpp  | 202 ++++++++++++++++++
 2 files changed, 203 insertions(+)
 create mode 100644 Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt
index c6d6e9e1b00..7f4149bb8c4 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt
@@ -5,6 +5,7 @@ project(Sofa.Component.SolidMechanics.FEM.Elastic_test)
 set(SOURCE_FILES
     BaseTetrahedronFEMForceField_test.h
     BeamFEMForceField_test.cpp
+    CorotationalFEMForceField_test.cpp
     FastTetrahedralCorotationalForceField_test.cpp
     FEMSourceTermIntegrator_test.cpp
     HexahedralFEMForceField_test.cpp
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
new file mode 100644
index 00000000000..5fd913628e7
--- /dev/null
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
@@ -0,0 +1,202 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* This program is free software; you can redistribute it and/or modify it     *
+* under the terms of the GNU Lesser General Public License as published by    *
+* the Free Software Foundation; either version 2.1 of the License, or (at     *
+* your option) any later version.                                             *
+*                                                                             *
+* This program is distributed in the hope that it will be useful, but WITHOUT *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
+* for more details.                                                           *
+*                                                                             *
+* You should have received a copy of the GNU Lesser General Public License    *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#include <sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h>
+#include <sofa/component/solidmechanics/testing/ForceFieldTestCreation.h>
+#include <sofa/component/topology/container/constant/MeshTopology.h>
+
+#include <cmath>
+#include <numbers>
+#include <string>
+#include <type_traits>
+#include <vector>
+
+namespace sofa
+{
+
+using namespace sofa::component::solidmechanics::fem::elastic;
+
+namespace
+{
+
+using sofa::defaulttype::Vec2Types;
+using sofa::defaulttype::Vec3Types;
+
+template <class DataTypes>
+using Matrix = sofa::type::Mat<DataTypes::spatial_dimensions, DataTypes::spatial_dimensions, SReal>;
+
+constexpr SReal degreeToRadian = std::numbers::pi_v<SReal> / 180;
+
+/// Helpers for computing rotations in 2D and 3D
+Matrix<Vec2Types> rotation(SReal angle)
+{
+    return {{std::cos(angle), -std::sin(angle)}, {std::sin(angle), std::cos(angle)}};
+}
+
+Matrix<Vec3Types> rotation(const sofa::type::Vec3& axis, SReal angle)
+{
+    const auto K = sofa::type::crossProductMatrix(axis.normalized());
+    return Matrix<Vec3Types>::Identity() + K * std::sin(angle) + K * K * (1 - std::cos(angle));
+}
+
+}
+
+/// Checks the force derivatives against a finite difference of the force, on a single element.
+template <class ForceField>
+struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<ForceField>
+{
+    using DataTypes = typename ForceField::DataTypes;
+    using VecCoord = sofa::VecCoord_t<DataTypes>;
+    using VecDeriv = sofa::VecDeriv_t<DataTypes>;
+    using Coord = sofa::Coord_t<DataTypes>;
+    using Deriv = sofa::Deriv_t<DataTypes>;
+
+    /// Rest positions of the element.
+    VecCoord restPositions;
+
+    /// Expected force produced when the element is stretched by 50% along x
+    VecDeriv expectedStretchForces;
+
+    /// The rotation methods available for the element, except identity, which does not rotate.
+    std::vector<std::string> rotationMethods {"polar", "polar_affine", "stable_polar", "stable_polar_affine"};
+
+    CorotationalFEMForceFieldTangent_test()
+    {
+        const auto topology = sofa::core::objectmodel::New<sofa::component::topology::container::constant::MeshTopology>();
+        this->node->addObject(topology);
+
+        // E = 1 and nu = 0: the stretch stress is 0.5 along x only.
+        if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Quad>>)
+        {
+            topology->addQuad(0, 1, 2, 3);
+            restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0), Coord(1.0, 1.0), Coord(0.0, 1.0)};
+            expectedStretchForces = VecDeriv{Deriv(0.25, 0.0), Deriv(-0.25, 0.0), Deriv(-0.25, 0.0), Deriv(0.25, 0.0)};
+        }
+        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Triangle>>)
+        {
+            topology->addTriangle(0, 1, 2);
+            restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0), Coord(0.0, 1.0)};
+            expectedStretchForces = VecDeriv{Deriv(0.25, 0.0), Deriv(-0.25, 0.0), Deriv(0.0, 0.0)};
+        }
+        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>>)
+        {
+            topology->addHexa(0, 1, 2, 3, 4, 5, 6, 7);
+            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(1.0, 1.0, 0.0), Coord(0.0, 1.0, 0.0),
+                                     Coord(0.0, 0.0, 1.0), Coord(1.0, 0.0, 1.0), Coord(1.0, 1.0, 1.0), Coord(0.0, 1.0, 1.0)};
+            expectedStretchForces = VecDeriv{Deriv(0.125, 0.0, 0.0), Deriv(-0.125, 0.0, 0.0), Deriv(-0.125, 0.0, 0.0), Deriv(0.125, 0.0, 0.0),
+                                             Deriv(0.125, 0.0, 0.0), Deriv(-0.125, 0.0, 0.0), Deriv(-0.125, 0.0, 0.0), Deriv(0.125, 0.0, 0.0)};
+            rotationMethods.push_back("hexahedron");
+        }
+        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>>)
+        {
+            topology->addTetra(0, 1, 2, 3);
+            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0), Coord(0.0, 0.0, 1.0)};
+            expectedStretchForces = VecDeriv{Deriv(1.0 / 12, 0.0, 0.0), Deriv(-1.0 / 12, 0.0, 0.0), Deriv(0.0, 0.0, 0.0), Deriv(0.0, 0.0, 0.0)};
+            rotationMethods.push_back("triangle");
+        }
+        else
+        {
+            static_assert(sizeof(ForceField) == 0, "No element defined for this force field");
+        }
+
+        this->dof->writeOnlyRestPositions().wref() = restPositions;
+
+        this->force->d_youngModulus.setValue({1});
+        this->force->d_poissonRatio.setValue({0});
+
+        this->deltaRange = {5e5, 5e6};
+        this->flags &= ~sofa::ForceField_test<ForceField>::TEST_ADD_K_TO_MATRIX;
+    }
+
+    /// Rotation by 30 degrees. In 3D, the rotation is taken about the axis (1, 2, 3).
+    static Matrix<DataTypes> rotation()
+    {
+        if constexpr (DataTypes::spatial_dimensions == 2)
+        {
+            return sofa::rotation(30 * degreeToRadian);
+        }
+        else
+        {
+            return sofa::rotation(sofa::type::Vec3(1.0, 2.0, 3.0), 30 * degreeToRadian);
+        }
+    }
+
+    /// Stretch by 50% along x.
+    static Matrix<DataTypes> stretch()
+    {
+        if constexpr (DataTypes::spatial_dimensions == 2)
+        {
+            return {{1.5, 0.0}, {0.0, 1.0}};
+        }
+        else
+        {
+            return {{1.5, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
+        }
+    }
+
+    /// Moves the nodes by the deformation gradient F = R U and expects the forces f rotated by R, for
+    /// each rotation method.
+    void runTest(const Matrix<DataTypes>& R, const Matrix<DataTypes>& U, const VecDeriv& f)
+    {
+        VecCoord x(restPositions.size());
+        VecDeriv v(restPositions.size());
+        VecDeriv expectedForces(restPositions.size());
+        for (std::size_t j = 0; j < restPositions.size(); ++j)
+        {
+            x[j] = R * U * restPositions[j];
+            expectedForces[j] = R * f[j];
+        }
+
+        for (const auto& rotationMethod : rotationMethods)
+        {
+            SCOPED_TRACE("rotationMethod = " + rotationMethod);
+            this->force->findData("rotationMethod")->read(rotationMethod);
+            this->run_test(x, v, expectedForces);
+        }
+    }
+};
+
+using TestTypes = ::testing::Types<
+    CorotationalFEMForceField<Vec2Types, sofa::geometry::Quad>,
+    CorotationalFEMForceField<Vec2Types, sofa::geometry::Triangle>,
+    CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>,
+    CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>>;
+
+TYPED_TEST_SUITE(CorotationalFEMForceFieldTangent_test, TestTypes);
+
+TYPED_TEST(CorotationalFEMForceFieldTangent_test, rotated)
+{
+    using Matrix = Matrix<typename TestFixture::DataTypes>;
+    this->runTest(this->rotation(), Matrix::Identity(), typename TestFixture::VecDeriv(this->restPositions.size()));
+}
+
+TYPED_TEST(CorotationalFEMForceFieldTangent_test, stretched)
+{
+    using Matrix = Matrix<typename TestFixture::DataTypes>;
+    this->runTest(Matrix::Identity(), this->stretch(), this->expectedStretchForces);
+}
+
+TYPED_TEST(CorotationalFEMForceFieldTangent_test, stretchedRotated)
+{
+    this->runTest(this->rotation(), this->stretch(), this->expectedStretchForces);
+}
+
+}

From 78a6beee710762a3681352c6f79a742c7071adc5 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Fri, 23 Oct 2026 16:50:03 +0200
Subject: [PATCH 02/11] Add embedded elements in the corotational test

---
 .../tests/CorotationalFEMForceField_test.cpp  | 24 ++++++++++++++++++-
 1 file changed, 23 insertions(+), 1 deletion(-)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
index 5fd913628e7..30cdbe82456 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
@@ -96,6 +96,25 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
             restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0), Coord(0.0, 1.0)};
             expectedStretchForces = VecDeriv{Deriv(0.25, 0.0), Deriv(-0.25, 0.0), Deriv(0.0, 0.0)};
         }
+        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Edge>>)
+        {
+            topology->addEdge(0, 1);
+            restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0)};
+            expectedStretchForces = VecDeriv{Deriv(0.5, 0.0), Deriv(-0.5, 0.0)};
+        }
+        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Quad>>)
+        {
+            topology->addQuad(0, 1, 2, 3);
+            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(1.0, 1.0, 0.0), Coord(0.0, 1.0, 0.0)};
+            expectedStretchForces = VecDeriv{Deriv(0.25, 0.0, 0.0), Deriv(-0.25, 0.0, 0.0), Deriv(-0.25, 0.0, 0.0), Deriv(0.25, 0.0, 0.0)};
+        }
+        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Triangle>>)
+        {
+            topology->addTriangle(0, 1, 2);
+            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0)};
+            expectedStretchForces = VecDeriv{Deriv(0.25, 0.0, 0.0), Deriv(-0.25, 0.0, 0.0), Deriv(0.0, 0.0, 0.0)};
+            rotationMethods.push_back("triangle");
+        }
         else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>>)
         {
             topology->addHexa(0, 1, 2, 3, 4, 5, 6, 7);
@@ -178,7 +197,10 @@ using TestTypes = ::testing::Types<
     CorotationalFEMForceField<Vec2Types, sofa::geometry::Quad>,
     CorotationalFEMForceField<Vec2Types, sofa::geometry::Triangle>,
     CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>,
-    CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>>;
+    CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>,
+    CorotationalFEMForceField<Vec2Types, sofa::geometry::Edge>,
+    CorotationalFEMForceField<Vec3Types, sofa::geometry::Quad>,
+    CorotationalFEMForceField<Vec3Types, sofa::geometry::Triangle>>;
 
 TYPED_TEST_SUITE(CorotationalFEMForceFieldTangent_test, TestTypes);
 

From f544e2a562b4eb1e0211d04052365e29e7d0adbc Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Fri, 23 Oct 2026 17:25:36 +0200
Subject: [PATCH 03/11] Use irregular element shapes and refactor tests

---
 .../tests/CorotationalFEMForceField_test.cpp  | 146 +++++++++++-------
 1 file changed, 93 insertions(+), 53 deletions(-)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
index 30cdbe82456..9ad0913e03b 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
@@ -45,18 +45,42 @@ using Matrix = sofa::type::Mat<DataTypes::spatial_dimensions, DataTypes::spatial
 
 constexpr SReal degreeToRadian = std::numbers::pi_v<SReal> / 180;
 
-/// Helpers for computing rotations in 2D and 3D
-Matrix<Vec2Types> rotation(SReal angle)
+/// Helpers for computing rotation matrices in 2D and 3D
+Matrix<Vec2Types> rotationMatrix2D(SReal angle)
 {
     return {{std::cos(angle), -std::sin(angle)}, {std::sin(angle), std::cos(angle)}};
 }
 
-Matrix<Vec3Types> rotation(const sofa::type::Vec3& axis, SReal angle)
+Matrix<Vec3Types> rotationMatrix3D(const sofa::type::Vec3& axis, SReal angle)
 {
     const auto K = sofa::type::crossProductMatrix(axis.normalized());
     return Matrix<Vec3Types>::Identity() + K * std::sin(angle) + K * K * (1 - std::cos(angle));
 }
 
+// The matrix for the applied rotation in 2D and 3D
+const Matrix<Vec2Types> prescribedRotationMatrix2D = rotationMatrix2D(30 * degreeToRadian);
+const Matrix<Vec3Types> prescribedRotationMatrix3D = rotationMatrix3D(sofa::type::Vec3(1.0, 2.0, 3.0), 30 * degreeToRadian);
+
+// Stretches by 50% along the axes 
+// (0.866, 0.5) in 2D and
+// (0.876, 0.420, -0.239) in 3D
+// This axis was chosen so that it does not align with any edge of the test elements. If placed
+// along an edge, the "triangle" and "hexahedron" rotation methods can extract the exact rotation
+// by coincidence.
+const Matrix<Vec2Types> prescribedStretchMatrix2D =
+    prescribedRotationMatrix2D * Matrix<Vec2Types>{{1.5, 0.0}, {0.0, 1.0}} * prescribedRotationMatrix2D.transposed();
+const Matrix<Vec3Types> prescribedStretchMatrix3D =
+    prescribedRotationMatrix3D * Matrix<Vec3Types>{{1.5, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}} * prescribedRotationMatrix3D.transposed();
+
+/// prescribedStretchMatrix2D within the plane z = 0, for the triangles and quads in 3D.
+const Matrix<Vec3Types> prescribedFlatStretchMatrix3D {
+    {prescribedStretchMatrix2D(0, 0), prescribedStretchMatrix2D(0, 1), 0.0},
+    {prescribedStretchMatrix2D(1, 0), prescribedStretchMatrix2D(1, 1), 0.0},
+    {0.0, 0.0, 1.0}};
+
+// An edge is stretched along its own line. any other direction involves rotation.
+const Matrix<Vec2Types> prescribedEdgeStretchMatrix2D {{1.5, 0.0}, {0.0, 1.0}};
+
 }
 
 /// Checks the force derivatives against a finite difference of the force, on a single element.
@@ -67,13 +91,12 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
     using VecCoord = sofa::VecCoord_t<DataTypes>;
     using VecDeriv = sofa::VecDeriv_t<DataTypes>;
     using Coord = sofa::Coord_t<DataTypes>;
-    using Deriv = sofa::Deriv_t<DataTypes>;
 
-    /// Rest positions of the element.
+    /// Rest positions of the element. 
     VecCoord restPositions;
 
-    /// Expected force produced when the element is stretched by 50% along x
-    VecDeriv expectedStretchForces;
+    Matrix<DataTypes> rotation;
+    Matrix<DataTypes> stretch;
 
     /// The rotation methods available for the element, except identity, which does not rotate.
     std::vector<std::string> rotationMethods {"polar", "polar_affine", "stable_polar", "stable_polar_affine"};
@@ -83,52 +106,60 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
         const auto topology = sofa::core::objectmodel::New<sofa::component::topology::container::constant::MeshTopology>();
         this->node->addObject(topology);
 
-        // E = 1 and nu = 0: the stretch stress is 0.5 along x only.
+        /// An irregular shape is chosen for all elements that follow. This skews the rigid fit of
+        /// the polar decomposition methods. Rigid and affine fits extract different rotations and
+        /// the tangent of each must follow how its own rotation changes.
         if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Quad>>)
         {
             topology->addQuad(0, 1, 2, 3);
-            restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0), Coord(1.0, 1.0), Coord(0.0, 1.0)};
-            expectedStretchForces = VecDeriv{Deriv(0.25, 0.0), Deriv(-0.25, 0.0), Deriv(-0.25, 0.0), Deriv(0.25, 0.0)};
+            restPositions = VecCoord{Coord(0.0, 0.0), Coord(3.0, 0.0), Coord(3.0, 1.0), Coord(0.0, 1.0)};
+            rotation = prescribedRotationMatrix2D;
+            stretch = prescribedStretchMatrix2D;
         }
         else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Triangle>>)
         {
             topology->addTriangle(0, 1, 2);
             restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0), Coord(0.0, 1.0)};
-            expectedStretchForces = VecDeriv{Deriv(0.25, 0.0), Deriv(-0.25, 0.0), Deriv(0.0, 0.0)};
+            rotation = prescribedRotationMatrix2D;
+            stretch = prescribedStretchMatrix2D;
         }
         else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Edge>>)
         {
             topology->addEdge(0, 1);
             restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0)};
-            expectedStretchForces = VecDeriv{Deriv(0.5, 0.0), Deriv(-0.5, 0.0)};
+            rotation = prescribedRotationMatrix2D;
+            stretch = prescribedEdgeStretchMatrix2D;
         }
         else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Quad>>)
         {
             topology->addQuad(0, 1, 2, 3);
-            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(1.0, 1.0, 0.0), Coord(0.0, 1.0, 0.0)};
-            expectedStretchForces = VecDeriv{Deriv(0.25, 0.0, 0.0), Deriv(-0.25, 0.0, 0.0), Deriv(-0.25, 0.0, 0.0), Deriv(0.25, 0.0, 0.0)};
+            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(3.0, 0.0, 0.0), Coord(3.0, 1.0, 0.0), Coord(0.0, 1.0, 0.0)};
+            rotation = prescribedRotationMatrix3D;
+            stretch = prescribedFlatStretchMatrix3D;
         }
         else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Triangle>>)
         {
             topology->addTriangle(0, 1, 2);
             restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0)};
-            expectedStretchForces = VecDeriv{Deriv(0.25, 0.0, 0.0), Deriv(-0.25, 0.0, 0.0), Deriv(0.0, 0.0, 0.0)};
+            rotation = prescribedRotationMatrix3D;
+            stretch = prescribedFlatStretchMatrix3D;
             rotationMethods.push_back("triangle");
         }
         else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>>)
         {
             topology->addHexa(0, 1, 2, 3, 4, 5, 6, 7);
-            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(1.0, 1.0, 0.0), Coord(0.0, 1.0, 0.0),
-                                     Coord(0.0, 0.0, 1.0), Coord(1.0, 0.0, 1.0), Coord(1.0, 1.0, 1.0), Coord(0.0, 1.0, 1.0)};
-            expectedStretchForces = VecDeriv{Deriv(0.125, 0.0, 0.0), Deriv(-0.125, 0.0, 0.0), Deriv(-0.125, 0.0, 0.0), Deriv(0.125, 0.0, 0.0),
-                                             Deriv(0.125, 0.0, 0.0), Deriv(-0.125, 0.0, 0.0), Deriv(-0.125, 0.0, 0.0), Deriv(0.125, 0.0, 0.0)};
+            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(3.0, 0.0, 0.0), Coord(3.0, 2.0, 0.0), Coord(0.0, 2.0, 0.0),
+                                     Coord(0.0, 0.0, 1.0), Coord(3.0, 0.0, 1.0), Coord(3.0, 2.0, 1.0), Coord(0.0, 2.0, 1.0)};
+            rotation = prescribedRotationMatrix3D;
+            stretch = prescribedStretchMatrix3D;
             rotationMethods.push_back("hexahedron");
         }
         else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>>)
         {
             topology->addTetra(0, 1, 2, 3);
             restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0), Coord(0.0, 0.0, 1.0)};
-            expectedStretchForces = VecDeriv{Deriv(1.0 / 12, 0.0, 0.0), Deriv(-1.0 / 12, 0.0, 0.0), Deriv(0.0, 0.0, 0.0), Deriv(0.0, 0.0, 0.0)};
+            rotation = prescribedRotationMatrix3D;
+            stretch = prescribedStretchMatrix3D;
             rotationMethods.push_back("triangle");
         }
         else
@@ -145,50 +176,59 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
         this->flags &= ~sofa::ForceField_test<ForceField>::TEST_ADD_K_TO_MATRIX;
     }
 
-    /// Rotation by 30 degrees. In 3D, the rotation is taken about the axis (1, 2, 3).
-    static Matrix<DataTypes> rotation()
+    /// The forces of the force field at the positions x.
+    VecDeriv forcesAt(const VecCoord& x, const VecDeriv& v)
     {
-        if constexpr (DataTypes::spatial_dimensions == 2)
-        {
-            return sofa::rotation(30 * degreeToRadian);
-        }
-        else
-        {
-            return sofa::rotation(sofa::type::Vec3(1.0, 2.0, 3.0), 30 * degreeToRadian);
-        }
+        this->setupState(x, v);
+        sofa::simulation::node::initRoot(this->node.get());
+        core::MechanicalParams mparams;
+        this->computeForce(&mparams);
+        VecDeriv forces;
+        sofa::testing::copyFromData(forces, this->dof->readForces());
+        return forces;
     }
 
-    /// Stretch by 50% along x.
-    static Matrix<DataTypes> stretch()
-    {
-        if constexpr (DataTypes::spatial_dimensions == 2)
-        {
-            return {{1.5, 0.0}, {0.0, 1.0}};
-        }
-        else
-        {
-            return {{1.5, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
-        }
-    }
-
-    /// Moves the nodes by the deformation gradient F = R U and expects the forces f rotated by R, for
-    /// each rotation method.
-    void runTest(const Matrix<DataTypes>& R, const Matrix<DataTypes>& U, const VecDeriv& f)
+    /// Moves the nodes by the deformation gradient F = R U and checks each rotation method.
+    void runTest(const Matrix<DataTypes>& R, const Matrix<DataTypes>& U)
     {
+        // The element at rest, stretched by U, then rotated by R.
+        VecCoord stretchedPositions(restPositions.size());
         VecCoord x(restPositions.size());
         VecDeriv v(restPositions.size());
-        VecDeriv expectedForces(restPositions.size());
         for (std::size_t j = 0; j < restPositions.size(); ++j)
         {
-            x[j] = R * U * restPositions[j];
-            expectedForces[j] = R * f[j];
+            stretchedPositions[j] = U * restPositions[j];
+            x[j] = R * stretchedPositions[j];
         }
 
+        // The expected forces: the linear elastic forces of the stretch, computed by the identity
+        // method, rotated by R.
+        this->force->findData("rotationMethod")->read("identity");
+        VecDeriv expectedForces = forcesAt(stretchedPositions, v);
+        for (auto& force : expectedForces)
+        {
+            force = R * force;
+        }
+
+        // This test checks the consistency of the tangent. The elements have irregular shapes, on which
+        // only the polar decomposition of the affine fit extracts the rotation that produces the correct
+        // forces. For polar_affine and stable_polar_affine, the forces are therefore checked against the
+        // expected ones. The other methods can produce incorrect forces on irregular shapes, so their
+        // forces are checked against the ones they compute themselves, which bypasses the force check
+        // and leaves only the tangent checked.
         for (const auto& rotationMethod : rotationMethods)
         {
             SCOPED_TRACE("rotationMethod = " + rotationMethod);
             this->force->findData("rotationMethod")->read(rotationMethod);
-            this->run_test(x, v, expectedForces);
+
+            if (rotationMethod == "polar_affine" || rotationMethod == "stable_polar_affine")
+            {
+                this->run_test(x, v, expectedForces);
+            }
+            else
+            {
+                this->run_test(x, v, forcesAt(x, v));
+            }
         }
     }
 };
@@ -207,18 +247,18 @@ TYPED_TEST_SUITE(CorotationalFEMForceFieldTangent_test, TestTypes);
 TYPED_TEST(CorotationalFEMForceFieldTangent_test, rotated)
 {
     using Matrix = Matrix<typename TestFixture::DataTypes>;
-    this->runTest(this->rotation(), Matrix::Identity(), typename TestFixture::VecDeriv(this->restPositions.size()));
+    this->runTest(this->rotation, Matrix::Identity());
 }
 
 TYPED_TEST(CorotationalFEMForceFieldTangent_test, stretched)
 {
     using Matrix = Matrix<typename TestFixture::DataTypes>;
-    this->runTest(Matrix::Identity(), this->stretch(), this->expectedStretchForces);
+    this->runTest(Matrix::Identity(), this->stretch);
 }
 
 TYPED_TEST(CorotationalFEMForceFieldTangent_test, stretchedRotated)
 {
-    this->runTest(this->rotation(), this->stretch(), this->expectedStretchForces);
+    this->runTest(this->rotation, this->stretch);
 }
 
 }

From 580e44757ddbe256a6714e7ee059095c24005ac9 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Fri, 23 Oct 2026 17:58:59 +0200
Subject: [PATCH 04/11] Circumvent the potential energy test

---
 .../FEM/Elastic/tests/CorotationalFEMForceField_test.cpp       | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
index 9ad0913e03b..a65a93fa0b6 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
@@ -227,7 +227,10 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
             }
             else
             {
+                // These methods are also expected to miscompute the potential energy on irregular shapes.
+                this->flags &= ~sofa::ForceField_test<ForceField>::TEST_POTENTIAL_ENERGY;
                 this->run_test(x, v, forcesAt(x, v));
+                this->flags |= sofa::ForceField_test<ForceField>::TEST_POTENTIAL_ENERGY;
             }
         }
     }

From 0f6fcedcd3dfab5a28e04504a9c63d5e00851d62 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Sat, 24 Oct 2026 01:27:16 +0200
Subject: [PATCH 05/11] Instantiate test based on rotatio method

---
 .../tests/CorotationalFEMForceField_test.cpp  | 131 +++++++++++-------
 1 file changed, 83 insertions(+), 48 deletions(-)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
index a65a93fa0b6..938fcb04d21 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
@@ -25,9 +25,7 @@
 
 #include <cmath>
 #include <numbers>
-#include <string>
 #include <type_traits>
-#include <vector>
 
 namespace sofa
 {
@@ -81,12 +79,35 @@ const Matrix<Vec3Types> prescribedFlatStretchMatrix3D {
 // An edge is stretched along its own line. any other direction involves rotation.
 const Matrix<Vec2Types> prescribedEdgeStretchMatrix2D {{1.5, 0.0}, {0.0, 1.0}};
 
+struct PolarMethod { static constexpr const char* key = "polar"; };
+struct PolarAffineMethod { static constexpr const char* key = "polar_affine"; };
+struct StablePolarMethod { static constexpr const char* key = "stable_polar"; };
+struct StablePolarAffineMethod { static constexpr const char* key = "stable_polar_affine"; };
+struct TriangleMethod { static constexpr const char* key = "triangle"; };
+struct HexahedronMethod { static constexpr const char* key = "hexahedron"; };
+
+using Quad2D = CorotationalFEMForceField<Vec2Types, sofa::geometry::Quad>;
+using Triangle2D = CorotationalFEMForceField<Vec2Types, sofa::geometry::Triangle>;
+using Edge2D = CorotationalFEMForceField<Vec2Types, sofa::geometry::Edge>;
+using Quad3D = CorotationalFEMForceField<Vec3Types, sofa::geometry::Quad>;
+using Triangle3D = CorotationalFEMForceField<Vec3Types, sofa::geometry::Triangle>;
+using Hexahedron3D = CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>;
+using Tetrahedron3D = CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>;
+
+template <class ForceFieldType, class MethodType>
+struct TestCase
+{
+    using ForceField = ForceFieldType;
+    using Method = MethodType;
+};
+
 }
 
 /// Checks the force derivatives against a finite difference of the force, on a single element.
-template <class ForceField>
-struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<ForceField>
+template <class TestCase>
+struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<typename TestCase::ForceField>
 {
+    using ForceField = typename TestCase::ForceField;
     using DataTypes = typename ForceField::DataTypes;
     using VecCoord = sofa::VecCoord_t<DataTypes>;
     using VecDeriv = sofa::VecDeriv_t<DataTypes>;
@@ -98,9 +119,6 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
     Matrix<DataTypes> rotation;
     Matrix<DataTypes> stretch;
 
-    /// The rotation methods available for the element, except identity, which does not rotate.
-    std::vector<std::string> rotationMethods {"polar", "polar_affine", "stable_polar", "stable_polar_affine"};
-
     CorotationalFEMForceFieldTangent_test()
     {
         const auto topology = sofa::core::objectmodel::New<sofa::component::topology::container::constant::MeshTopology>();
@@ -109,58 +127,55 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
         /// An irregular shape is chosen for all elements that follow. This skews the rigid fit of
         /// the polar decomposition methods. Rigid and affine fits extract different rotations and
         /// the tangent of each must follow how its own rotation changes.
-        if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Quad>>)
+        if constexpr (std::is_same_v<ForceField, Quad2D>)
         {
             topology->addQuad(0, 1, 2, 3);
             restPositions = VecCoord{Coord(0.0, 0.0), Coord(3.0, 0.0), Coord(3.0, 1.0), Coord(0.0, 1.0)};
             rotation = prescribedRotationMatrix2D;
             stretch = prescribedStretchMatrix2D;
         }
-        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Triangle>>)
+        else if constexpr (std::is_same_v<ForceField, Triangle2D>)
         {
             topology->addTriangle(0, 1, 2);
             restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0), Coord(0.0, 1.0)};
             rotation = prescribedRotationMatrix2D;
             stretch = prescribedStretchMatrix2D;
         }
-        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec2Types, sofa::geometry::Edge>>)
+        else if constexpr (std::is_same_v<ForceField, Edge2D>)
         {
             topology->addEdge(0, 1);
             restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0)};
             rotation = prescribedRotationMatrix2D;
             stretch = prescribedEdgeStretchMatrix2D;
         }
-        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Quad>>)
+        else if constexpr (std::is_same_v<ForceField, Quad3D>)
         {
             topology->addQuad(0, 1, 2, 3);
             restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(3.0, 0.0, 0.0), Coord(3.0, 1.0, 0.0), Coord(0.0, 1.0, 0.0)};
             rotation = prescribedRotationMatrix3D;
             stretch = prescribedFlatStretchMatrix3D;
         }
-        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Triangle>>)
+        else if constexpr (std::is_same_v<ForceField, Triangle3D>)
         {
             topology->addTriangle(0, 1, 2);
             restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0)};
             rotation = prescribedRotationMatrix3D;
             stretch = prescribedFlatStretchMatrix3D;
-            rotationMethods.push_back("triangle");
         }
-        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>>)
+        else if constexpr (std::is_same_v<ForceField, Hexahedron3D>)
         {
             topology->addHexa(0, 1, 2, 3, 4, 5, 6, 7);
             restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(3.0, 0.0, 0.0), Coord(3.0, 2.0, 0.0), Coord(0.0, 2.0, 0.0),
                                      Coord(0.0, 0.0, 1.0), Coord(3.0, 0.0, 1.0), Coord(3.0, 2.0, 1.0), Coord(0.0, 2.0, 1.0)};
             rotation = prescribedRotationMatrix3D;
             stretch = prescribedStretchMatrix3D;
-            rotationMethods.push_back("hexahedron");
         }
-        else if constexpr (std::is_same_v<ForceField, CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>>)
+        else if constexpr (std::is_same_v<ForceField, Tetrahedron3D>)
         {
             topology->addTetra(0, 1, 2, 3);
             restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(1.0, 0.0, 0.0), Coord(0.0, 1.0, 0.0), Coord(0.0, 0.0, 1.0)};
             rotation = prescribedRotationMatrix3D;
             stretch = prescribedStretchMatrix3D;
-            rotationMethods.push_back("triangle");
         }
         else
         {
@@ -177,7 +192,7 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
     }
 
     /// The forces of the force field at the positions x.
-    VecDeriv forcesAt(const VecCoord& x, const VecDeriv& v)
+    VecDeriv computeForcesAt(const VecCoord& x, const VecDeriv& v)
     {
         this->setupState(x, v);
         sofa::simulation::node::initRoot(this->node.get());
@@ -188,7 +203,8 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
         return forces;
     }
 
-    /// Moves the nodes by the deformation gradient F = R U and checks each rotation method.
+    /// Moves the nodes by the deformation gradient F = R U and runs the ForceField test there.
+    /// The consistency of the tangent stiffness is checked therein.
     void runTest(const Matrix<DataTypes>& R, const Matrix<DataTypes>& U)
     {
         // The element at rest, stretched by U, then rotated by R.
@@ -204,46 +220,65 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<Forc
         // The expected forces: the linear elastic forces of the stretch, computed by the identity
         // method, rotated by R.
         this->force->findData("rotationMethod")->read("identity");
-        VecDeriv expectedForces = forcesAt(stretchedPositions, v);
+        VecDeriv expectedForces = computeForcesAt(stretchedPositions, v);
         for (auto& force : expectedForces)
         {
             force = R * force;
         }
 
-        // This test checks the consistency of the tangent. The elements have irregular shapes, on which
-        // only the polar decomposition of the affine fit extracts the rotation that produces the correct
-        // forces. For polar_affine and stable_polar_affine, the forces are therefore checked against the
-        // expected ones. The other methods can produce incorrect forces on irregular shapes, so their
-        // forces are checked against the ones they compute themselves, which bypasses the force check
-        // and leaves only the tangent checked.
-        for (const auto& rotationMethod : rotationMethods)
+        using Method = typename TestCase::Method;
+        this->force->findData("rotationMethod")->read(Method::key);
+
+        // For irregular elements, only the polar decomposition of the affine fit extracts the
+        // rotation that produces the correct forces. 
+        if constexpr (std::is_same_v<Method, PolarAffineMethod> || std::is_same_v<Method, StablePolarAffineMethod>)
+        {
+            // Therefore the expected force is used as reference.
+            this->run_test(x, v, expectedForces);
+        }
+        else
         {
-            SCOPED_TRACE("rotationMethod = " + rotationMethod);
-            this->force->findData("rotationMethod")->read(rotationMethod);
-
-            if (rotationMethod == "polar_affine" || rotationMethod == "stable_polar_affine")
-            {
-                this->run_test(x, v, expectedForces);
-            }
-            else
-            {
-                // These methods are also expected to miscompute the potential energy on irregular shapes.
-                this->flags &= ~sofa::ForceField_test<ForceField>::TEST_POTENTIAL_ENERGY;
-                this->run_test(x, v, forcesAt(x, v));
-                this->flags |= sofa::ForceField_test<ForceField>::TEST_POTENTIAL_ENERGY;
-            }
+            // The rest of the methods produce incorrect forces on irregular shapes. So the check is
+            // implicitly bypassed by using the force they compute themselves as reference. These
+            // methods are also expected to miscompute the potential energy on irregular shapes.
+            this->flags &= ~sofa::ForceField_test<ForceField>::TEST_POTENTIAL_ENERGY;
+            this->run_test(x, v, computeForcesAt(x, v));
         }
     }
 };
 
 using TestTypes = ::testing::Types<
-    CorotationalFEMForceField<Vec2Types, sofa::geometry::Quad>,
-    CorotationalFEMForceField<Vec2Types, sofa::geometry::Triangle>,
-    CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>,
-    CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>,
-    CorotationalFEMForceField<Vec2Types, sofa::geometry::Edge>,
-    CorotationalFEMForceField<Vec3Types, sofa::geometry::Quad>,
-    CorotationalFEMForceField<Vec3Types, sofa::geometry::Triangle>>;
+    TestCase<Quad2D, PolarMethod>,
+    TestCase<Quad2D, PolarAffineMethod>,
+    TestCase<Quad2D, StablePolarMethod>,
+    TestCase<Quad2D, StablePolarAffineMethod>,
+    TestCase<Triangle2D, PolarMethod>,
+    TestCase<Triangle2D, PolarAffineMethod>,
+    TestCase<Triangle2D, StablePolarMethod>,
+    TestCase<Triangle2D, StablePolarAffineMethod>,
+    TestCase<Edge2D, PolarMethod>,
+    TestCase<Edge2D, PolarAffineMethod>,
+    TestCase<Edge2D, StablePolarMethod>,
+    TestCase<Edge2D, StablePolarAffineMethod>,
+    TestCase<Quad3D, PolarMethod>,
+    TestCase<Quad3D, PolarAffineMethod>,
+    TestCase<Quad3D, StablePolarMethod>,
+    TestCase<Quad3D, StablePolarAffineMethod>,
+    TestCase<Triangle3D, PolarMethod>,
+    TestCase<Triangle3D, PolarAffineMethod>,
+    TestCase<Triangle3D, StablePolarMethod>,
+    TestCase<Triangle3D, StablePolarAffineMethod>,
+    TestCase<Triangle3D, TriangleMethod>,
+    TestCase<Hexahedron3D, PolarMethod>,
+    TestCase<Hexahedron3D, PolarAffineMethod>,
+    TestCase<Hexahedron3D, StablePolarMethod>,
+    TestCase<Hexahedron3D, StablePolarAffineMethod>,
+    TestCase<Hexahedron3D, HexahedronMethod>,
+    TestCase<Tetrahedron3D, PolarMethod>,
+    TestCase<Tetrahedron3D, PolarAffineMethod>,
+    TestCase<Tetrahedron3D, StablePolarMethod>,
+    TestCase<Tetrahedron3D, StablePolarAffineMethod>,
+    TestCase<Tetrahedron3D, TriangleMethod>>;
 
 TYPED_TEST_SUITE(CorotationalFEMForceFieldTangent_test, TestTypes);
 

From 0f0f593444b4ce75ee2f907bcb7821555444bec9 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Sat, 24 Oct 2026 01:29:54 +0200
Subject: [PATCH 06/11] Add back the tangent computation in Corotational force
 field and helpers

---
 .../fem/elastic/CorotationalFEMForceField.h   |  6 +++
 .../fem/elastic/CorotationalFEMForceField.inl | 47 +++++++++++++++++--
 .../impl/rotations/RotationMethodsContainer.h | 30 ++++++++++++
 Sofa/framework/Type/src/sofa/type/Mat.h       |  7 +++
 4 files changed, 87 insertions(+), 3 deletions(-)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h
index 4f610a84d83..21142175371 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h
@@ -171,6 +171,8 @@ class CorotationalFEMForceField :
 
     StressVoigtVector computeStress(const DeformationGradient& F, sofa::Size elementId) override;
 
+    sofa::Data<bool> d_consistentTangent;
+
 protected:
 
     void beforeElementForce(const sofa::core::MechanicalParams* mparams,
@@ -192,6 +194,10 @@ class CorotationalFEMForceField :
     sofa::type::vector<RotationMatrix> m_rotations;
     sofa::type::vector<RotationMatrix> m_initialRotationsTransposed;
 
+    /// Element stiffness plus the stiffness due to the element rotation, in the element frame,
+    /// accumulated in computeElementsForces. Empty if consistentTangent is off.
+    sofa::type::vector<typename trait::ElementHessian> m_tangents;
+
     /// Local displacement of the element nodes in the element's own rotated frame. Flat vector.
     ElementDisplacement computeElementLocalDisplacement(
         const std::array<sofa::Coord_t<DataTypes>, trait::NumberOfNodesInElement>& nodes,
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl
index fc5fc9fd809..7e9bbdffd72 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl
@@ -32,7 +32,8 @@ namespace sofa::component::solidmechanics::fem::elastic
 
 template <class DataTypes, class ElementType>
 CorotationalFEMForceField<DataTypes, ElementType>::CorotationalFEMForceField()
-    : m_rotationMethods(this)
+    : d_consistentTangent(initData(&d_consistentTangent, true, "consistentTangent", "Include the change of the element rotations in the force derivatives"))
+    , m_rotationMethods(this)
 {
     this->addUpdateCallback("selectRotationMethod", {&this->m_rotationMethods.d_rotationMethod},
         [this](const sofa::core::DataTracker&)
@@ -74,6 +75,7 @@ void CorotationalFEMForceField<DataTypes, ElementType>::beforeElementForce(
     SOFA_UNUSED(x);
     const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
     m_rotations.resize(elements.size(), RotationMatrix::Identity());
+    m_tangents.resize(d_consistentTangent.getValue() ? elements.size() : 0);
 }
 
 template <class DataTypes, class ElementType>
@@ -129,6 +131,43 @@ void CorotationalFEMForceField<DataTypes, ElementType>::computeElementsForces(
         auto& elementForce = elementForces[elementId];
         elementForce = stiffnessMatrix * displacement;
 
+        if (!m_tangents.empty())
+        {
+            typename RotationMethods<DataTypes, ElementType>::RotationJacobian jacobian;
+            if (!m_rotationMethods.computeRotationJacobian(jacobian, elementRotation, elementNodesCoordinates, restElementNodesCoordinates))
+            {
+                msg_warning_once() << "The rotation jacobian of an element could not be computed. Its tangent ignores the change of rotation.";
+            }
+
+            const auto t0 = translation(restElementNodesCoordinates);
+            std::array<decltype(sofa::type::crossProductMatrix(sofa::Coord_t<DataTypes>()).transposed()), trait::NumberOfNodesInElement> localPositions;
+            sofa::Coord_t<DataTypes> nodeValue { sofa::type::NOINIT };
+            for (sofa::Size c = 0; c < trait::NumberOfNodesInElement; ++c)
+            {
+                displacement.getsub(c * DIM, nodeValue);
+                localPositions[c] = sofa::type::crossProductMatrix(nodeValue + restElementNodesCoordinates[c] - t0).transposed();
+            }
+
+            auto& tangent = m_tangents[elementId];
+            tangent = stiffnessMatrix;
+            RotationMatrix block { sofa::type::NOINIT };
+            for (sofa::Size a = 0; a < trait::NumberOfNodesInElement; ++a)
+            {
+                elementForce.getsub(a * DIM, nodeValue);
+                auto W = sofa::type::crossProductMatrix(nodeValue).transposed();
+                for (sofa::Size c = 0; c < trait::NumberOfNodesInElement; ++c)
+                {
+                    stiffnessMatrix.getsub(a * DIM, c * DIM, block);
+                    W -= block * localPositions[c];
+                }
+                for (sofa::Size b = 0; b < trait::NumberOfNodesInElement; ++b)
+                {
+                    tangent.getsub(a * DIM, b * DIM, block);
+                    tangent.setsub(a * DIM, b * DIM, block + W * jacobian[b]);
+                }
+            }
+        }
+
         sofa::type::Vec<DIM, sofa::Real_t<DataTypes>> nodeForce { sofa::type::NOINIT };
         for (sofa::Size n = 0; n < trait::NumberOfNodesInElement; ++n)
         {
@@ -149,6 +188,7 @@ void CorotationalFEMForceField<DataTypes, ElementType>::computeElementsForcesDer
     SOFA_UNUSED(mparams);
     const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
     auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
+    const auto& elementTangentStiffness = m_tangents.empty() ? elementStiffness.ref() : m_tangents;
 
     for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
     {
@@ -163,7 +203,7 @@ void CorotationalFEMForceField<DataTypes, ElementType>::computeElementsForcesDer
             element_dx.setsub(n * trait::spatial_dimensions, elementRotation.multTranspose(nodeDx[element[n]]));
         }
 
-        const auto& stiffnessMatrix = elementStiffness[elementId];
+        const auto& stiffnessMatrix = elementTangentStiffness[elementId];
 
         auto& df = elementForcesDeriv[elementId];
         df = stiffnessMatrix * element_dx;
@@ -194,7 +234,8 @@ void CorotationalFEMForceField<DataTypes, ElementType>::buildStiffnessMatrix(
         return;
     }
 
-    auto elementStiffnessIt = elementStiffness.begin();
+    const auto& elementTangentStiffness = m_tangents.empty() ? elementStiffness.ref() : m_tangents;
+    auto elementStiffnessIt = elementTangentStiffness.begin();
     auto rotationMatrixIt = m_rotations.begin();
     for (const auto& element : elements)
     {
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/RotationMethodsContainer.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/RotationMethodsContainer.h
index f4f84998014..322a46c3d8b 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/RotationMethodsContainer.h
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/RotationMethodsContainer.h
@@ -71,6 +71,7 @@ struct RotationMethodsContainer
 
     using RotationMatrix = sofa::type::Mat<DataTypes::spatial_dimensions, DataTypes::spatial_dimensions, sofa::Real_t<DataTypes>>;
     static constexpr std::size_t NumberOfNodesInElement = ElementType::NumberOfNodes;
+    using RotationJacobian = std::array<decltype(sofa::type::crossProductMatrix(sofa::Coord_t<DataTypes>())), NumberOfNodesInElement>;
 
     /**
      * Computes the current rotation state relative to the initial configuration.
@@ -92,6 +93,35 @@ struct RotationMethodsContainer
             m_rotationComputer);
     }
 
+    /**
+     * Computes the jacobian of the rotation matrix. Zero for methods that do not provide one.
+     *
+     * @param jacobian Output: Derivative of the rotation with respect to each node position, in the element frame
+     * @param rotationMatrix Current rotation matrix, as computed by computeRotation
+     * @param nodesPosition Current positions of element nodes
+     * @param nodesRestPosition Rest positions of element nodes (initial configuration)
+     * @return False if the method failed to compute the jacobian
+     */
+    bool computeRotationJacobian(RotationJacobian& jacobian, const RotationMatrix& rotationMatrix,
+        const std::array<sofa::Coord_t<DataTypes>, NumberOfNodesInElement>& nodesPosition,
+        const std::array<sofa::Coord_t<DataTypes>, NumberOfNodesInElement>& nodesRestPosition)
+    {
+        return std::visit(
+            [&](auto& rotationComputer)
+            {
+                if constexpr (requires { rotationComputer.computeRotationJacobian(jacobian, rotationMatrix, nodesPosition, nodesRestPosition); })
+                {
+                    return rotationComputer.computeRotationJacobian(jacobian, rotationMatrix, nodesPosition, nodesRestPosition);
+                }
+                else
+                {
+                    jacobian = {};
+                    return true;
+                }
+            },
+            m_rotationComputer);
+    }
+
     struct RotationMethodsItems final : sofa::helper::SelectableItem<RotationMethodsItems>
     {
         using sofa::helper::SelectableItem<RotationMethodsItems>::SelectableItem;
diff --git a/Sofa/framework/Type/src/sofa/type/Mat.h b/Sofa/framework/Type/src/sofa/type/Mat.h
index 32db3108afc..897eb090654 100644
--- a/Sofa/framework/Type/src/sofa/type/Mat.h
+++ b/Sofa/framework/Type/src/sofa/type/Mat.h
@@ -1463,6 +1463,13 @@ constexpr Mat<3, 3, Real> crossProductMatrix(const Vec<3, Real>& v) noexcept
     return res;
 }
 
+/// crossProductMatrix(v) * x = v.cross(x), a scalar in 2D
+template<class Real>
+constexpr Mat<1, 2, Real> crossProductMatrix(const Vec<2, Real>& v) noexcept
+{
+    return {{-v[1], v[0]}};
+}
+
 
 /// return a * b^T
 template<sofa::Size L,class Real>

From 79b2a4854b8577b14e1a6dd445dd904134db67f5 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Sat, 24 Oct 2026 01:30:28 +0200
Subject: [PATCH 07/11] Apply the tangent fix to the rigid fit

---
 .../impl/rotations/PolarDecomposition.h       | 66 +++++++++++++++++++
 1 file changed, 66 insertions(+)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h
index f9b413d911a..c66090e912d 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h
@@ -41,6 +41,7 @@ template <class DataTypes, class ElementType, Fit FitType, Algorithm AlgorithmTy
 struct PolarDecomposition
 {
     using RotationMatrix = sofa::type::Mat<DataTypes::spatial_dimensions, DataTypes::spatial_dimensions, sofa::Real_t<DataTypes>>;
+    using StretchMatrix = RotationMatrix;
 
     /// The remainder of the dimensions of an embedded element within a higher-dimensional space.
     // TODO: embedded elements (Codimension > 0) are not handled.
@@ -108,6 +109,71 @@ struct PolarDecomposition
         }
     }
 
+    /// Computes the jacobian of the rotation matrix
+    template<std::size_t NumberOfNodesInElement, class JacobianBlock>
+    bool computeRotationJacobian(std::array<JacobianBlock, NumberOfNodesInElement>& jacobian,
+        const RotationMatrix& rotationMatrix,
+        const std::array<sofa::Coord_t<DataTypes>, NumberOfNodesInElement>& nodesPosition,
+        const std::array<sofa::Coord_t<DataTypes>, NumberOfNodesInElement>& nodesRestPosition)
+    {
+        if constexpr (DataTypes::spatial_dimensions == 1)
+        {
+            // In 1D, there is no rotation.
+            return true;
+        }
+        else if constexpr (Codimension > 0 || FitType == Fit::Affine)
+        {
+            // TODO: the affine fit and embedded elements. Without a jacobian, the tangent treats the
+            // rotation as constant.
+            jacobian = {};
+            return false;
+        }
+        else
+        {
+            const auto t = computeCentroid(nodesPosition);
+            const auto t0 = computeCentroid(nodesRestPosition);
+
+            sofa::type::Mat<NumberOfNodesInElement, DataTypes::spatial_dimensions, sofa::Real_t<DataTypes>> P(sofa::type::NOINIT);
+            sofa::type::Mat<NumberOfNodesInElement, DataTypes::spatial_dimensions, sofa::Real_t<DataTypes>> Q(sofa::type::NOINIT);
+
+            for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j)
+            {
+                P[j] = nodesPosition[j] - t;
+                Q[j] = nodesRestPosition[j] - t0;
+            }
+
+            const auto H = P.multTranspose(Q);
+
+            // S and the offsets q_k are divided by the trace of M so that small elements do not test A as
+            // singular. The jacobian is unchanged.
+            const auto inverseTraceM = static_cast<sofa::Real_t<DataTypes>>(1) / sofa::type::trace(Q.multTranspose(Q));
+            const StretchMatrix S = rotationMatrix.multTranspose(H) * inverseTraceM;
+
+            sofa::type::Mat<JacobianBlock::nbLines, JacobianBlock::nbLines, sofa::Real_t<DataTypes>> A;
+            if constexpr (DataTypes::spatial_dimensions == 3)
+            {
+                A = StretchMatrix::Identity() * sofa::type::trace(S) - S;
+            }
+            else
+            {
+                A(0, 0) = sofa::type::trace(S);
+            }
+
+            decltype(A) inverseA(sofa::type::NOINIT);
+            if (!inverseA.invert(A))
+            {
+                jacobian = {};
+                return false;
+            }
+
+            for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j)
+            {
+                jacobian[j] = inverseA * sofa::type::crossProductMatrix(Q[j] * inverseTraceM);
+            }
+            return true;
+        }
+    }
+
     static constexpr sofa::helper::Item getItem()
     {
         if constexpr (AlgorithmType == Algorithm::Standard)

From ed9986c4afc7deac9ddde6711798da320163cda7 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Tue, 27 Oct 2026 09:15:17 +0100
Subject: [PATCH 08/11] Apply the affine fit tangent

---
 .../impl/rotations/PolarDecomposition.h       | 30 ++++++++++++-------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h
index c66090e912d..14d0560affd 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h
@@ -121,10 +121,9 @@ struct PolarDecomposition
             // In 1D, there is no rotation.
             return true;
         }
-        else if constexpr (Codimension > 0 || FitType == Fit::Affine)
+        else if constexpr (Codimension > 0)
         {
-            // TODO: the affine fit and embedded elements. Without a jacobian, the tangent treats the
-            // rotation as constant.
+            // TODO: embedded elements.
             jacobian = {};
             return false;
         }
@@ -142,12 +141,23 @@ struct PolarDecomposition
                 Q[j] = nodesRestPosition[j] - t0;
             }
 
-            const auto H = P.multTranspose(Q);
+            const auto M = Q.multTranspose(Q);
+            const auto traceM = sofa::type::trace(M);
 
-            // S and the offsets q_k are divided by the trace of M so that small elements do not test A as
-            // singular. The jacobian is unchanged.
-            const auto inverseTraceM = static_cast<sofa::Real_t<DataTypes>>(1) / sofa::type::trace(Q.multTranspose(Q));
-            const StretchMatrix S = rotationMatrix.multTranspose(H) * inverseTraceM;
+            const auto w = [&Q, &M, traceM]()
+            {
+                if constexpr (FitType == Fit::Affine)
+                {
+                    return Q * (M / traceM).inverted() / traceM;
+                }
+                else
+                {
+                    return Q / traceM;
+                }
+            }();
+
+            const auto F = P.multTranspose(w);
+            const StretchMatrix S = rotationMatrix.multTranspose(F);
 
             sofa::type::Mat<JacobianBlock::nbLines, JacobianBlock::nbLines, sofa::Real_t<DataTypes>> A;
             if constexpr (DataTypes::spatial_dimensions == 3)
@@ -166,9 +176,9 @@ struct PolarDecomposition
                 return false;
             }
 
-            for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j)
+            for (sofa::Size k = 0; k < NumberOfNodesInElement; ++k)
             {
-                jacobian[j] = inverseA * sofa::type::crossProductMatrix(Q[j] * inverseTraceM);
+                jacobian[k] = inverseA * sofa::type::crossProductMatrix(w[k]);
             }
             return true;
         }

From ee0819173c227ea290e600e03a2bdb6305034af7 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Tue, 27 Oct 2026 09:50:01 +0100
Subject: [PATCH 09/11] Rotation jacobian in triangle and hexahedron rotation
 methods

---
 .../impl/rotations/HexahedronRotation.h       | 39 +++++++++++++++++++
 .../elastic/impl/rotations/TriangleRotation.h | 28 +++++++++++++
 2 files changed, 67 insertions(+)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/HexahedronRotation.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/HexahedronRotation.h
index 996ba550949..134890a6ee5 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/HexahedronRotation.h
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/HexahedronRotation.h
@@ -64,6 +64,45 @@ struct HexahedronRotation
         rotationMatrix = currentRotation.multTranspose(initialRotationMatrix);
     }
 
+    template<std::size_t NumberOfNodesInElement, class JacobianBlock>
+    bool computeRotationJacobian(std::array<JacobianBlock, NumberOfNodesInElement>& jacobian,
+        const RotationMatrix& rotationMatrix,
+        const std::array<sofa::Coord_t<DataTypes>, NumberOfNodesInElement>& nodesPosition,
+        const std::array<sofa::Coord_t<DataTypes>, NumberOfNodesInElement>& nodesRestPosition)
+    {
+        using Real = sofa::Real_t<DataTypes>;
+        SOFA_UNUSED(nodesRestPosition);
+
+        RotationMatrix currentFrame(sofa::type::NOINIT);
+        computeOrientationFromHexahedron(currentFrame, nodesPosition);
+        const RotationMatrix restFrame = currentFrame * rotationMatrix;
+
+        const auto a = (nodesPosition[1] - nodesPosition[0] +
+                        nodesPosition[2] - nodesPosition[3] +
+                        nodesPosition[5] - nodesPosition[4] +
+                        nodesPosition[6] - nodesPosition[7]) * 0.25;
+        const auto b = (nodesPosition[3] - nodesPosition[0] +
+                        nodesPosition[2] - nodesPosition[1] +
+                        nodesPosition[7] - nodesPosition[4] +
+                        nodesPosition[6] - nodesPosition[5]) * 0.25;
+        const auto lengthA = a.norm();
+        const auto b1 = b * currentFrame[0];
+        const auto b2 = b * currentFrame[1];
+
+        const RotationMatrix Ga {{0, 0, -b1 / (lengthA * b2)}, {0, 0, -1 / lengthA}, {0, 1 / lengthA, 0}};
+        const RotationMatrix Gb {{0, 0, 1 / b2}, {0, 0, 0}, {0, 0, 0}};
+
+        constexpr std::array<Real, 8> alpha {-0.25, 0.25, 0.25, -0.25, -0.25, 0.25, 0.25, -0.25};
+        constexpr std::array<Real, 8> beta {-0.25, -0.25, 0.25, 0.25, -0.25, -0.25, 0.25, 0.25};
+
+        jacobian = {};
+        for (std::size_t k = 0; k < 8; ++k)
+        {
+            jacobian[k] = restFrame.multTranspose((Ga * alpha[k] + Gb * beta[k]) * restFrame);
+        }
+        return true;
+    }
+
     static constexpr sofa::helper::Item getItem()
     {
         return {"hexahedron", "Compute the rotation based on two average edges in the hexahedron"};
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/TriangleRotation.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/TriangleRotation.h
index 4247b921dba..c86cbe2e4f6 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/TriangleRotation.h
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/TriangleRotation.h
@@ -45,6 +45,34 @@ struct TriangleRotation
         rotationMatrix = currentRotation.multTranspose(initialRotationMatrix);
     }
 
+    template<std::size_t NumberOfNodesInElement, class JacobianBlock>
+    bool computeRotationJacobian(std::array<JacobianBlock, NumberOfNodesInElement>& jacobian,
+        const RotationMatrix& rotationMatrix,
+        const std::array<sofa::Coord_t<DataTypes>, NumberOfNodesInElement>& nodesPosition,
+        const std::array<sofa::Coord_t<DataTypes>, NumberOfNodesInElement>& nodesRestPosition)
+    {
+        SOFA_UNUSED(nodesRestPosition);
+
+        RotationMatrix currentFrame(sofa::type::NOINIT);
+        computeRotationFrom3Points(currentFrame, {nodesPosition[0], nodesPosition[1], nodesPosition[2]});
+        const RotationMatrix restFrame = currentFrame * rotationMatrix;
+
+        const auto a = nodesPosition[1] - nodesPosition[0];
+        const auto b = nodesPosition[2] - nodesPosition[0];
+        const auto lengthA = a.norm();
+        const auto b1 = b * currentFrame[0];
+        const auto b2 = b * currentFrame[1];
+
+        const RotationMatrix Ga {{0, 0, -b1 / (lengthA * b2)}, {0, 0, -1 / lengthA}, {0, 1 / lengthA, 0}};
+        const RotationMatrix Gb {{0, 0, 1 / b2}, {0, 0, 0}, {0, 0, 0}};
+
+        jacobian = {};
+        jacobian[0] = restFrame.multTranspose((-Ga - Gb) * restFrame);
+        jacobian[1] = restFrame.multTranspose(Ga * restFrame);
+        jacobian[2] = restFrame.multTranspose(Gb * restFrame);
+        return true;
+    }
+
     static constexpr sofa::helper::Item getItem()
     {
         return {"triangle", "Compute the rotation based on the Gram-Schmidt frame alignment"};

From 1ad5d486ccb03793ed693c55f7ddf51d2682bc31 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Fri, 6 Nov 2026 11:14:18 +0100
Subject: [PATCH 10/11] Remove tests with embedded elements & polar
 decomposition

---
 .../tests/CorotationalFEMForceField_test.cpp  | 33 +------------------
 1 file changed, 1 insertion(+), 32 deletions(-)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
index 938fcb04d21..bde218de10f 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
@@ -70,15 +70,12 @@ const Matrix<Vec2Types> prescribedStretchMatrix2D =
 const Matrix<Vec3Types> prescribedStretchMatrix3D =
     prescribedRotationMatrix3D * Matrix<Vec3Types>{{1.5, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}} * prescribedRotationMatrix3D.transposed();
 
-/// prescribedStretchMatrix2D within the plane z = 0, for the triangles and quads in 3D.
+/// prescribedStretchMatrix2D within the plane z = 0, for the triangles in 3D.
 const Matrix<Vec3Types> prescribedFlatStretchMatrix3D {
     {prescribedStretchMatrix2D(0, 0), prescribedStretchMatrix2D(0, 1), 0.0},
     {prescribedStretchMatrix2D(1, 0), prescribedStretchMatrix2D(1, 1), 0.0},
     {0.0, 0.0, 1.0}};
 
-// An edge is stretched along its own line. any other direction involves rotation.
-const Matrix<Vec2Types> prescribedEdgeStretchMatrix2D {{1.5, 0.0}, {0.0, 1.0}};
-
 struct PolarMethod { static constexpr const char* key = "polar"; };
 struct PolarAffineMethod { static constexpr const char* key = "polar_affine"; };
 struct StablePolarMethod { static constexpr const char* key = "stable_polar"; };
@@ -88,8 +85,6 @@ struct HexahedronMethod { static constexpr const char* key = "hexahedron"; };
 
 using Quad2D = CorotationalFEMForceField<Vec2Types, sofa::geometry::Quad>;
 using Triangle2D = CorotationalFEMForceField<Vec2Types, sofa::geometry::Triangle>;
-using Edge2D = CorotationalFEMForceField<Vec2Types, sofa::geometry::Edge>;
-using Quad3D = CorotationalFEMForceField<Vec3Types, sofa::geometry::Quad>;
 using Triangle3D = CorotationalFEMForceField<Vec3Types, sofa::geometry::Triangle>;
 using Hexahedron3D = CorotationalFEMForceField<Vec3Types, sofa::geometry::Hexahedron>;
 using Tetrahedron3D = CorotationalFEMForceField<Vec3Types, sofa::geometry::Tetrahedron>;
@@ -141,20 +136,6 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<type
             rotation = prescribedRotationMatrix2D;
             stretch = prescribedStretchMatrix2D;
         }
-        else if constexpr (std::is_same_v<ForceField, Edge2D>)
-        {
-            topology->addEdge(0, 1);
-            restPositions = VecCoord{Coord(0.0, 0.0), Coord(1.0, 0.0)};
-            rotation = prescribedRotationMatrix2D;
-            stretch = prescribedEdgeStretchMatrix2D;
-        }
-        else if constexpr (std::is_same_v<ForceField, Quad3D>)
-        {
-            topology->addQuad(0, 1, 2, 3);
-            restPositions = VecCoord{Coord(0.0, 0.0, 0.0), Coord(3.0, 0.0, 0.0), Coord(3.0, 1.0, 0.0), Coord(0.0, 1.0, 0.0)};
-            rotation = prescribedRotationMatrix3D;
-            stretch = prescribedFlatStretchMatrix3D;
-        }
         else if constexpr (std::is_same_v<ForceField, Triangle3D>)
         {
             topology->addTriangle(0, 1, 2);
@@ -256,18 +237,6 @@ using TestTypes = ::testing::Types<
     TestCase<Triangle2D, PolarAffineMethod>,
     TestCase<Triangle2D, StablePolarMethod>,
     TestCase<Triangle2D, StablePolarAffineMethod>,
-    TestCase<Edge2D, PolarMethod>,
-    TestCase<Edge2D, PolarAffineMethod>,
-    TestCase<Edge2D, StablePolarMethod>,
-    TestCase<Edge2D, StablePolarAffineMethod>,
-    TestCase<Quad3D, PolarMethod>,
-    TestCase<Quad3D, PolarAffineMethod>,
-    TestCase<Quad3D, StablePolarMethod>,
-    TestCase<Quad3D, StablePolarAffineMethod>,
-    TestCase<Triangle3D, PolarMethod>,
-    TestCase<Triangle3D, PolarAffineMethod>,
-    TestCase<Triangle3D, StablePolarMethod>,
-    TestCase<Triangle3D, StablePolarAffineMethod>,
     TestCase<Triangle3D, TriangleMethod>,
     TestCase<Hexahedron3D, PolarMethod>,
     TestCase<Hexahedron3D, PolarAffineMethod>,

From 6ccfa65ce49aeb9d5769c0f4ae7856344dc67f9e Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 9 Nov 2026 11:34:26 +0100
Subject: [PATCH 11/11] add option to symmetrize tangent

---
 .../solidmechanics/fem/elastic/CorotationalFEMForceField.h  | 1 +
 .../fem/elastic/CorotationalFEMForceField.inl               | 6 ++++++
 .../FEM/Elastic/tests/CorotationalFEMForceField_test.cpp    | 1 +
 3 files changed, 8 insertions(+)

diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h
index 21142175371..f7b91c7e180 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h
@@ -172,6 +172,7 @@ class CorotationalFEMForceField :
     StressVoigtVector computeStress(const DeformationGradient& F, sofa::Size elementId) override;
 
     sofa::Data<bool> d_consistentTangent;
+    sofa::Data<bool> d_symmetrizeTangent;
 
 protected:
 
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl
index 7e9bbdffd72..0b2d6a4041d 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl
@@ -33,6 +33,7 @@ namespace sofa::component::solidmechanics::fem::elastic
 template <class DataTypes, class ElementType>
 CorotationalFEMForceField<DataTypes, ElementType>::CorotationalFEMForceField()
     : d_consistentTangent(initData(&d_consistentTangent, true, "consistentTangent", "Include the change of the element rotations in the force derivatives"))
+    , d_symmetrizeTangent(initData(&d_symmetrizeTangent, true, "symmetrizeTangent", "Use the symmetric part of the consistent tangent"))
     , m_rotationMethods(this)
 {
     this->addUpdateCallback("selectRotationMethod", {&this->m_rotationMethods.d_rotationMethod},
@@ -166,6 +167,11 @@ void CorotationalFEMForceField<DataTypes, ElementType>::computeElementsForces(
                     tangent.setsub(a * DIM, b * DIM, block + W * jacobian[b]);
                 }
             }
+
+            if (d_symmetrizeTangent.getValue())
+            {
+                tangent = (tangent + tangent.transposed()) / 2;
+            }
         }
 
         sofa::type::Vec<DIM, sofa::Real_t<DataTypes>> nodeForce { sofa::type::NOINIT };
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
index bde218de10f..14269e6098b 100644
--- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
+++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CorotationalFEMForceField_test.cpp
@@ -167,6 +167,7 @@ struct CorotationalFEMForceFieldTangent_test : public sofa::ForceField_test<type
 
         this->force->d_youngModulus.setValue({1});
         this->force->d_poissonRatio.setValue({0});
+        this->force->d_symmetrizeTangent.setValue(false);
 
         this->deltaRange = {5e5, 5e6};
         this->flags &= ~sofa::ForceField_test<ForceField>::TEST_ADD_K_TO_MATRIX;