sofa-framework · alxbilger · Mar 31, 2026
diff --git a/...lastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp b/...lastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp
@@ -49,6 +49,9 @@ void registerElementCorotationalFEMForceField(sofa::core::ObjectFactory* factory
 
     factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear hexahedra using the corotational approach")
         .add< ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron> >(true));
+
+    factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear prisms using the corotational approach")
+        .add< ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism> >(true));
 }
 
 // template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
@@ -60,5 +63,6 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotational
 template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
 template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
 template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
 
 }
diff --git a/.../Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h b/.../Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h
@@ -198,6 +198,7 @@ extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorot
 extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
 extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
 extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
 #endif
 
 }  // namespace sofa::component::solidmechanics::fem::elastic
diff --git a/...c/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp b/...c/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp
@@ -49,6 +49,9 @@ void registerElementLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* fa
 
     factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear hexahedra assuming small strain")
         .add< ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron> >(true));
+
+    factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear prisms assuming small strain")
+        .add< ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism> >(true));
 }
 
 template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
@@ -60,5 +63,6 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallS
 template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
 template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
 template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
 
 }
diff --git a/...tic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h b/...tic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h
@@ -105,6 +105,7 @@ extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinea
 extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
 extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
 extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
 #endif
 
 }  // namespace sofa::component::solidmechanics::fem::elastic
diff --git a/Sofa/framework/FEM/CMakeLists.txt b/Sofa/framework/FEM/CMakeLists.txt
@@ -11,6 +11,7 @@ set(HEADER_FILES
     ${SOFAFEMSRC_ROOT}/FiniteElement[all].h
     ${SOFAFEMSRC_ROOT}/FiniteElement[Edge].h
     ${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].h
+    ${SOFAFEMSRC_ROOT}/FiniteElement[Prism].h
     ${SOFAFEMSRC_ROOT}/FiniteElement[Quad].h
     ${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].h
     ${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].h
@@ -21,6 +22,7 @@ set(SOURCE_FILES
 
     ${SOFAFEMSRC_ROOT}/FiniteElement[Edge].cpp
     ${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].cpp
+    ${SOFAFEMSRC_ROOT}/FiniteElement[Prism].cpp
     ${SOFAFEMSRC_ROOT}/FiniteElement[Quad].cpp
     ${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].cpp
     ${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].cpp

diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].cpp
@@ -0,0 +1,31 @@
+/******************************************************************************
+*                 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                             *
+******************************************************************************/
+#define SOFA_FEM_FINITE_ELEMENT_PRISM_CPP
+#include <sofa/fem/FiniteElement[Prism].h>
+#include <sofa/defaulttype/VecTypes.h>
+
+namespace sofa::fem
+{
+
+template struct SOFA_FEM_API FiniteElement<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>;
+
+}
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h
@@ -0,0 +1,97 @@
+/******************************************************************************
+*                 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                             *
+******************************************************************************/
+#pragma once
+#include <sofa/fem/FiniteElement.h>
+
+#if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP)
+#include <sofa/defaulttype/VecTypes.h>
+#endif
+
+namespace sofa::fem
+{
+
+template <class DataTypes>
+struct FiniteElement<sofa::geometry::Prism, DataTypes>
+{
+    FINITEELEMENT_HEADER(sofa::geometry::Prism, DataTypes, 3);
+    static_assert(spatial_dimensions == 3, "Prisms are only defined in 3D");
+
+    constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
+        {0, 0, 0},
+        {1, 0, 0},
+        {0, 1, 0},
+        {0, 0, 1},
+        {1, 0, 1},
+        {0, 1, 1},
+    }};
+
+    static const sofa::type::vector<TopologyElement>& getElementSequence(sofa::core::topology::BaseMeshTopology& topology)
+    {
+        return topology.getPrisms();
+    }
+
+    static constexpr sofa::type::Vec<NumberOfNodesInElement, Real> shapeFunctions(const sofa::type::Vec<TopologicalDimension, Real>& q)
+    {
+        return {
+            q[0] * q[2] - q[0] + q[1] * q[2] - q[1] - q[2] + 1,
+            (1 - q[2]) * q[0],
+            (1 - q[2]) * q[1],
+            (-q[0] - q[1] + 1) * q[2],
+            q[0] * q[2],
+            q[1] * q[2]
+        };
+    }
+
+    static constexpr sofa::type::Mat<NumberOfNodesInElement, TopologicalDimension, Real> gradientShapeFunctions(const sofa::type::Vec<TopologicalDimension, Real>& q)
+    {
+        SOFA_UNUSED(q);
+        return {
+            {q[2] - 1, q[2] - 1, q[0] + q[1] - 1},
+            {1 - q[2], 0, -q[0]},
+            {0, 1 - q[2], -q[1]},
+            {-q[2], -q[2], -q[0] - q[1] + 1},
+            {q[2], 0, q[0]},
+            {0, q[2], q[1]},
+        };
+    }
+
+    static constexpr std::array<QuadraturePointAndWeight, 2> quadraturePoints()
+    {
+        constexpr auto third = static_cast<Real>(1) / static_cast<Real>(3);
+        constexpr auto sqrt_3 = static_cast<Real>(0.57735026919); // 1/sqrt(3)
+        constexpr auto one = static_cast<Real>(1);
+        constexpr QuadraturePoint q0 {third, third, static_cast<Real>(0.5) * (one - sqrt_3)};
+        constexpr QuadraturePoint q1 {third, third, static_cast<Real>(0.5) * (one + sqrt_3)};
+
+        constexpr std::array<QuadraturePointAndWeight, 2> q {
+            std::make_pair(q0, 1./4.),
+            std::make_pair(q1, 1./4.),
+        };
+        return q;
+    }
+};
+
+#if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP)
+extern template struct SOFA_FEM_API FiniteElement<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>;
+#endif
+
+}
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h
@@ -23,6 +23,7 @@
 
 #include <sofa/fem/FiniteElement[Edge].h>
 #include <sofa/fem/FiniteElement[Hexahedron].h>
+#include <sofa/fem/FiniteElement[Prism].h>
 #include <sofa/fem/FiniteElement[Quad].h>
 #include <sofa/fem/FiniteElement[Tetrahedron].h>
 #include <sofa/fem/FiniteElement[Triangle].h>
diff --git a/Sofa/framework/FEM/test/FiniteElement_test.cpp b/Sofa/framework/FEM/test/FiniteElement_test.cpp
@@ -83,6 +83,11 @@ TEST(FiniteElement, hexa3dWeights)
     testSumWeights<sofa::geometry::Hexahedron, sofa::defaulttype::Vec3Types>(8);
 }
 
+TEST(FiniteElement, prism3dWeights)
+{
+    testSumWeights<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>(1 / 2.);
+}
+
 /**
  * Checks that the sum of the gradients of shape functions is zero at the evaluation point.
  */
@@ -158,6 +163,10 @@ TEST(FiniteElement, hexa3dGradientShapeFunctions)
     testGradientShapeFunctions<sofa::geometry::Hexahedron, sofa::defaulttype::Vec3Types>(sofa::type::Vec3(1., 1., 1.));
 }
 
-
+TEST(FiniteElement, prism3dGradientShapeFunctions)
+{
+    testGradientShapeFunctions<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>(sofa::type::Vec3(0., 0., 0.));
+    testGradientShapeFunctions<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>(sofa::type::Vec3(1., 1., 1.));
+}
 
 }
diff --git a/examples/Component/Topology/Mapping/Hexa2PrismTopologicalMapping.scn b/examples/Component/Topology/Mapping/Hexa2PrismTopologicalMapping.scn
@@ -1,24 +1,49 @@
 <?xml version="1.0"?>
-<Node name="root" dt="0.01" gravity="0 -9 0">
+<Node name="root" dt="0.01" gravity="0 -9.81 0">
+
     <Node name="plugins">
+        <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedProjectiveConstraint] -->
+        <RequiredPlugin name="Sofa.Component.Engine.Select"/> <!-- Needed to use components [BoxROI] -->
+        <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [SparseLDLSolver] -->
+        <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <!-- Needed to use components [CGLinearSolver] -->
+        <RequiredPlugin name="Sofa.Component.LinearSolver.Ordering"/> <!-- Needed to use components [NaturalOrderingMethod] -->
+        <RequiredPlugin name="Sofa.Component.LinearSystem"/> <!-- Needed to use components [ConstantSparsityPatternSystem] -->
+        <RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [MeshMatrixMass] -->
+        <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
+        <RequiredPlugin name="Sofa.Component.SolidMechanics.FEM.Elastic"/>
         <RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
-        <RequiredPlugin name="Sofa.Component.Topology.Container.Constant"/> <!-- Needed to use components [MeshTopology] -->
+        <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <!-- Needed to use components [TetrahedronSetGeometryAlgorithms,TetrahedronSetTopologyContainer,TetrahedronSetTopologyModifier] -->
         <RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
-        <RequiredPlugin name="Sofa.Component.Topology.Mapping"/> <!-- Needed to use components [Hexa2PrismTopologicalMapping] -->
-        <RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualMesh] -->
+        <RequiredPlugin name="Sofa.Component.Topology.Mapping"/> <!-- Needed to use components [Hexa2TetraTopologicalMapping] -->
+        <RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [LineAxis,VisualGrid,VisualStyle] -->
+        <RequiredPlugin name="Sofa.GL.Component.Rendering3D"/> <!-- Needed to use components [OglSceneFrame] -->
     </Node>
 
     <DefaultAnimationLoop/>
+    <VisualStyle displayFlags="showBehaviorModels showForceFields" />
+
+    <VisualGrid size="0.1"/>
+    <LineAxis size="0.1"/>
+    <OglSceneFrame/>
+
+    <EulerImplicitSolver name="backward Euler" rayleighStiffness="0.1" rayleighMass="0.1" />
+    <SparseLDLSolver/>
 
-    <Node name="grid">
-        <RegularGridTopology name="grid" min="-5 -5 0" max="5 5 40" n="5 5 20"/>
-        <MechanicalObject template="Vec3" name="state" position="@grid.position"/>
+    <RegularGridTopology name="grid" min="-0.01 -0.01 0" max="0.01 0.01 0.2" n="5 5 30"/>
+    <MechanicalObject template="Vec3" name="state" showObject="true"/>
 
-        <Node name="prisms">
-            <MeshTopology name="prism_topology"/>
-            <Hexa2PrismTopologicalMapping input="@grid" output="@prism_topology" />
-            <VisualMesh position="@../state.position" topology="@prism_topology" enable="true"/>
-        </Node>
+    <MeshMatrixMass massDensity="1100"/>
 
+    <Node name="prisms">
+        <MeshTopology name="prism_topology"/>
+        <Hexa2PrismTopologicalMapping input="@grid" output="@prism_topology" />
+
+        <PrismCorotationalFEMForceField name="FEM" youngModulus="2e6" poissonRatio="0.45" topology="@prism_topology"
+                                        rotationMethod="polar" computeForceStrategy="sequenced" computeForceDerivStrategy="sequenced"/>
+
+        <VisualMesh position="@../state.position" topology="@prism_topology" enable="true"/>
     </Node>
+
+    <BoxROI template="Vec3" name="box_roi" box="-0.011 -0.011 -0.0001   0.011 0.011 0.0001" drawBoxes="1" />
+    <FixedProjectiveConstraint template="Vec3" indices="@box_roi.indices" />
 </Node>