[SofaCUDA] ElementFEMForceField: Generic CUDA implementation

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h

@@ -55,6 +55,10 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f
     using StrainDisplacement = typename trait::StrainDisplacement;
     using Real = typename trait::Real;
 
+public:
+    using AssembledStiffnessMatrix = sofa::type::Mat<
+        trait::NumberOfDofsInElement, trait::NumberOfDofsInElement, Real>;
+
 protected:
 
     BaseElementLinearFEMForceField();
@@ -70,6 +74,14 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f
      * List of precomputed element stiffness matrices
      */
     sofa::Data<sofa::type::vector<ElementStiffness> > d_elementStiffness;
+
+    /**
+     * Contiguous buffer of assembled stiffness matrices (one per element).
+     * Extracted from d_elementStiffness for cache-friendly access in the hot path.
+     * This avoids loading the full FactorizedElementStiffness struct (~15 KB per element)
+     * when only the assembled matrix (~4.6 KB) is needed.
+     */
+    sofa::type::vector<AssembledStiffnessMatrix> m_assembledStiffnessMatrices;
 };
 
 #if !defined(ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP)

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl

@@ -100,6 +100,13 @@ void BaseElementLinearFEMForceField<DataTypes, ElementType>::precomputeElementSt
             const std::array<sofa::Coord_t<DataTypes>, trait::NumberOfNodesInElement> nodesCoordinates = extractNodesVectorFromGlobalVector(element, restPositionAccessor.ref());
             elementStiffness[elementId] = integrate<DataTypes, ElementType, trait::matrixVectorProductType>(nodesCoordinates, elasticityTensor);
         });
+
+    // Extract assembled matrices into a contiguous buffer for cache-friendly access
+    m_assembledStiffnessMatrices.resize(elements.size());
+    for (std::size_t i = 0; i < elements.size(); ++i)
+    {
+        m_assembledStiffnessMatrices[i] = elementStiffness[i].getAssembledMatrix();
+    }
 }
 
 }

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl

@@ -86,7 +86,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsFo
 {
     const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
     auto restPositionAccessor = this->mstate->readRestPositions();
-    auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
+    const auto& assembledMatrices = this->m_assembledStiffnessMatrices;
 
     for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
     {
@@ -112,7 +112,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsFo
             transformedDisplacement = elementRotation.multTranspose(elementNodesCoordinates[j] - t) - (restElementNodesCoordinates[j] - t0);
         }
 
-        const auto& stiffnessMatrix = elementStiffness[elementId];
+        const auto& stiffnessMatrix = assembledMatrices[elementId];
 
         auto& elementForce = elementForces[elementId];
         elementForce = stiffnessMatrix * displacement;
@@ -134,7 +134,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsFo
     const sofa::VecDeriv_t<DataTypes>& nodeDx)
 {
     const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
-    auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
+    const auto& assembledMatrices = this->m_assembledStiffnessMatrices;
 
     for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
     {
@@ -150,7 +150,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsFo
             rotated_dx = elementRotation.multTranspose(nodeDx[element[n]]);
         }
 
-        const auto& stiffnessMatrix = elementStiffness[elementId];
+        const auto& stiffnessMatrix = assembledMatrices[elementId];
 
         auto& df = elementForcesDeriv[elementId];
         df = stiffnessMatrix * element_dx;

Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl

@@ -49,12 +49,12 @@ void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::computeEleme
 {
     const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
     auto restPositionAccessor = this->mstate->readRestPositions();
-    auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
+    const auto& assembledMatrices = this->m_assembledStiffnessMatrices;
 
     for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
     {
         const auto& element = elements[elementId];
-        const auto& stiffnessMatrix = elementStiffness[elementId];
+        const auto& stiffnessMatrix = assembledMatrices[elementId];
 
         typename trait::ElementDisplacement displacement{ sofa::type::NOINIT };
 
@@ -79,12 +79,12 @@ void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::computeEleme
     const sofa::VecDeriv_t<DataTypes>& nodeDx)
 {
     const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
-    auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
+    const auto& assembledMatrices = this->m_assembledStiffnessMatrices;
 
     for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
     {
         const auto& element = elements[elementId];
-        const auto& stiffnessMatrix = elementStiffness[elementId];
+        const auto& stiffnessMatrix = assembledMatrices[elementId];
 
         const std::array<sofa::Coord_t<DataTypes>, trait::NumberOfNodesInElement> elementNodesDx =
             extractNodesVectorFromGlobalVector(element, nodeDx);

applications/plugins/SofaCUDA/Component/CMakeLists.txt

@@ -39,6 +39,11 @@ set(HEADER_FILES
 
 
     ### solidmechanics
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementFEMKernelUtils.cuh
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementCorotationalFEMForceField.h
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementCorotationalFEMForceField.inl
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementLinearSmallStrainFEMForceField.h
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementLinearSmallStrainFEMForceField.inl
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaHexahedronFEMForceField.h
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaHexahedronFEMForceField.inl
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/hyperelastic/CudaStandardTetrahedralFEMForceField.h
@@ -111,6 +116,8 @@ set(SOURCE_FILES
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/mass/CudaUniformMass.cpp
 
     ### Solidmechanics
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementCorotationalFEMForceField.cpp
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementLinearSmallStrainFEMForceField.cpp
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaHexahedronFEMForceField.cpp
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/hyperelastic/CudaStandardTetrahedralFEMForceField.cpp
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/tensormass/CudaTetrahedralTensorMassForceField.cpp
@@ -181,6 +188,8 @@ set(CUDA_SOURCES
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/mass/CudaUniformMass.cu
 
     ### solidmechanics
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementCorotationalFEMForceField.cu
+    ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementLinearSmallStrainFEMForceField.cu
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/elastic/CudaHexahedronFEMForceField.cu
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/fem/hyperelastic/CudaStandardTetrahedralFEMForceField.cu
     ${SOFACUDA_COMPONENT_SOURCE_DIR}/SofaCUDA/component/solidmechanics/tensormass/CudaTetrahedralTensorMassForceField.cu

applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/init.cpp

@@ -90,6 +90,8 @@ extern void registerPlaneForceField(sofa::core::ObjectFactory* factory);
 extern void registerSphereForceField(sofa::core::ObjectFactory* factory);
 
 // component::solidmechanics::fem::elastic
+extern void registerElementCorotationalFEMForceField(sofa::core::ObjectFactory* factory);
+extern void registerElementLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* factory);
 extern void registerHexahedronFEMForceField(sofa::core::ObjectFactory* factory);
 extern void registerTetrahedronFEMForceField(sofa::core::ObjectFactory* factory);
 extern void registerTriangularFEMForceFieldOptim(sofa::core::ObjectFactory* factory);
@@ -224,6 +226,8 @@ void registerObjects(sofa::core::ObjectFactory* factory)
     registerLinearForceField(factory);
     registerPlaneForceField(factory);
     registerSphereForceField(factory);
+    registerElementCorotationalFEMForceField(factory);
+    registerElementLinearSmallStrainFEMForceField(factory);
     registerHexahedronFEMForceField(factory);
     registerTetrahedronFEMForceField(factory);
     registerTriangularFEMForceFieldOptim(factory);

applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaElementCorotationalFEMForceField.cpp

@@ -0,0 +1,101 @@
+/******************************************************************************
+*                 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 <SofaCUDA/component/config.h>
+
+#include <sofa/gpu/cuda/CudaTypes.h>
+#include <SofaCUDA/component/solidmechanics/fem/elastic/CudaElementCorotationalFEMForceField.inl>
+#include <sofa/core/ObjectFactory.h>
+
+namespace sofa::component::solidmechanics::fem::elastic
+{
+
+using namespace sofa::gpu::cuda;
+
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Edge>;
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Triangle>;
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Quad>;
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Tetrahedron>;
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Hexahedron>;
+
+#ifdef SOFA_GPU_CUDA_DOUBLE
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Edge>;
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Triangle>;
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Quad>;
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Tetrahedron>;
+template class SOFACUDA_COMPONENT_API CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Hexahedron>;
+#endif
+
+} // namespace sofa::component::solidmechanics::fem::elastic
+
+namespace sofa::gpu::cuda
+{
+
+void registerElementCorotationalFEMForceField(sofa::core::ObjectFactory* factory)
+{
+    using namespace sofa::component::solidmechanics::fem::elastic;
+
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA for EdgeCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Edge> >()
+    );
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA for TriangleCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Triangle> >()
+    );
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA for QuadCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Quad> >()
+    );
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA for TetrahedronCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Tetrahedron> >()
+    );
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA for HexahedronCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3fTypes, sofa::geometry::Hexahedron> >()
+    );
+
+#ifdef SOFA_GPU_CUDA_DOUBLE
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA (double) for EdgeCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Edge> >()
+    );
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA (double) for TriangleCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Triangle> >()
+    );
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA (double) for QuadCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Quad> >()
+    );
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA (double) for TetrahedronCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Tetrahedron> >()
+    );
+    factory->registerObjects(sofa::core::ObjectRegistrationData(
+        "Supports GPU-side computations using CUDA (double) for HexahedronCorotationalFEMForceField")
+        .add< CudaElementCorotationalFEMForceField<CudaVec3dTypes, sofa::geometry::Hexahedron> >()
+    );
+#endif
+}
+
+} // namespace sofa::gpu::cuda