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

Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt

@@ -7,12 +7,20 @@ set(HEADER_FILES
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/config.h.in
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/init.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/fwd.h
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseElementLinearFEMForceField.h
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseElementLinearFEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.inl
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ElementCorotationalFEMForceField.h
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ElementCorotationalFEMForceField.inl
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ElementLinearSmallStrainFEMForceField.h
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ElementLinearSmallStrainFEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.inl
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMForceField.h
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/HexahedralFEMForceField.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/HexahedralFEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/HexahedralFEMForceFieldAndMass.h
@@ -37,13 +45,46 @@ set(HEADER_FILES
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceField.inl
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceFieldOptim.h
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceFieldOptim.inl
+
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/finiteelement/FiniteElement.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/finiteelement/FiniteElement[all].h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/finiteelement/FiniteElement[Edge].h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/finiteelement/FiniteElement[Hexahedron].h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/finiteelement/FiniteElement[Quad].h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/finiteelement/FiniteElement[Tetrahedron].h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/finiteelement/FiniteElement[Triangle].h
+
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/ComputeStrategy.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/ElementStiffnessMatrix.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/FullySymmetric4Tensor.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/IsotropicElasticityTensor.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/KroneckerDelta.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/LameParameters.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/MatrixTools.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/StrainDisplacement.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/SymmetricTensor.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/VectorTools.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/VecView.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/VoigtNotation.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/trait.h
+
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/HexahedronRotation.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/IdentityRotation.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/PolarDecomposition.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/RotationMethodsContainer.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/StablePolarDecomposition.h
+   ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/TriangleRotation.h
 )
 
 set(SOURCE_FILES
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/init.cpp
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseElementLinearFEMForceField.cpp
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.cpp
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.cpp
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ElementCorotationalFEMForceField.cpp
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/ElementLinearSmallStrainFEMForceField.cpp
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.cpp
+    ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMForceField.cpp
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/HexahedralFEMForceField.cpp
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/HexahedralFEMForceFieldAndMass.cpp
     ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/HexahedronFEMForceField.cpp

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

@@ -0,0 +1,37 @@
+/******************************************************************************
+*                 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: [email protected]                             *
+******************************************************************************/
+#define ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP
+#include <sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl>
+#include <sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[all].h>
+
+namespace sofa::component::solidmechanics::fem::elastic
+{
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec2Types, sofa::geometry::Edge>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Edge>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec2Types, sofa::geometry::Triangle>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Triangle>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec2Types, sofa::geometry::Quad>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
+template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
+}

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

@@ -0,0 +1,87 @@
+/******************************************************************************
+*                 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: [email protected]                             *
+******************************************************************************/
+#pragma once
+
+#include <sofa/component/solidmechanics/fem/elastic/config.h>
+#include <sofa/component/solidmechanics/fem/elastic/impl/trait.h>
+#include <sofa/component/solidmechanics/fem/elastic/BaseLinearElasticityFEMForceField.h>
+
+#if !defined(ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP)
+#include <sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[all].h>
+#endif
+
+namespace sofa::component::solidmechanics::fem::elastic
+{
+
+/**
+ * A base class for all element-based linear elastic force fields.
+ *
+ * It stores precomputed stiffness matrices (one per element) that are derived from:
+ *   - The initial configuration of the mechanical model
+ *   - Material properties (Young's modulus, Poisson's ratio)
+ */
+template <class DataTypes, class ElementType>
+class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::fem::elastic::BaseLinearElasticityFEMForceField<DataTypes>
+{
+public:
+    SOFA_ABSTRACT_CLASS(
+        SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType),
+        sofa::component::solidmechanics::fem::elastic::BaseLinearElasticityFEMForceField<DataTypes>);
+
+    void init() override;
+
+private:
+    using trait = sofa::component::solidmechanics::fem::elastic::trait<DataTypes, ElementType>;
+    using ElementStiffness = typename trait::ElementStiffness;
+    using ElasticityTensor = typename trait::ElasticityTensor;
+    using StrainDisplacement = typename trait::StrainDisplacement;
+
+protected:
+
+    BaseElementLinearFEMForceField();
+
+    /**
+     * With linear small strain, the element stiffness matrix is constant, so it can be precomputed.
+     */
+    void precomputeElementStiffness();
+
+public:
+
+    /**
+     * List of precomputed element stiffness matrices
+     */
+    sofa::Data<sofa::type::vector<ElementStiffness> > d_elementStiffness;
+};
+
+#if !defined(ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP)
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec2Types, sofa::geometry::Edge>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Edge>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec2Types, sofa::geometry::Triangle>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Triangle>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec2Types, sofa::geometry::Quad>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
+extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
+#endif
+
+}  // namespace sofa::component::solidmechanics::fem::elastic

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

@@ -0,0 +1,100 @@
+/******************************************************************************
+*                 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: [email protected]                             *
+******************************************************************************/
+#pragma once
+
+#include <sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h>
+#include <sofa/component/solidmechanics/fem/elastic/BaseLinearElasticityFEMForceField.inl>
+#include <sofa/component/solidmechanics/fem/elastic/impl/LameParameters.h>
+#include <sofa/component/solidmechanics/fem/elastic/impl/VectorTools.h>
+
+#include <ranges>
+
+#include <sofa/component/solidmechanics/fem/elastic/FEMForceField.h>
+
+namespace sofa::component::solidmechanics::fem::elastic
+{
+
+template <class DataTypes, class ElementType>
+BaseElementLinearFEMForceField<DataTypes, ElementType>::BaseElementLinearFEMForceField()
+    : d_elementStiffness(initData(&d_elementStiffness, "elementStiffness", "List of stiffness matrices per element"))
+{
+    this->addUpdateCallback("precomputeStiffness", {&this->d_youngModulus, &this->d_poissonRatio},
+    [this](const sofa::core::DataTracker& )
+    {
+        precomputeElementStiffness();
+        return this->getComponentState();
+    }, {});
+}
+
+template <class DataTypes, class ElementType>
+void BaseElementLinearFEMForceField<DataTypes, ElementType>::init()
+{
+    sofa::component::solidmechanics::fem::elastic::BaseLinearElasticityFEMForceField<DataTypes>::init();
+
+    if (!this->isComponentStateInvalid())
+    {
+        this->precomputeElementStiffness();
+    }
+}
+
+template <class DataTypes, class ElementType>
+void BaseElementLinearFEMForceField<DataTypes, ElementType>::precomputeElementStiffness()
+{
+    if (!this->l_topology)
+        return;
+
+    if (this->isComponentStateInvalid())
+        return;
+
+    if (!this->mstate)
+        return;
+
+    const auto youngModulusAccessor = sofa::helper::ReadAccessor(this->d_youngModulus);
+    const auto poissonRatioAccessor = sofa::helper::ReadAccessor(this->d_poissonRatio);
+
+    auto restPositionAccessor = this->mstate->readRestPositions();
+
+    const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
+
+    auto elementStiffness = sofa::helper::getWriteOnlyAccessor(d_elementStiffness);
+    elementStiffness.resize(elements.size());
+
+    SCOPED_TIMER("precomputeStiffness");
+    sofa::helper::IotaView indices {static_cast<decltype(elements.size())>(0ul), elements.size()};
+    std::for_each(indices.begin(), indices.end(),
+        [&](const auto elementId)
+        {
+            const auto& element = elements[elementId];
+
+            const auto youngModulus = this->getYoungModulusInElement(elementId);
+            const auto poissonRatio = this->getPoissonRatioInElement(elementId);
+
+            const auto [mu, lambda] = sofa::component::solidmechanics::fem::elastic::toLameParameters<DataTypes>(youngModulus, poissonRatio);
+
+            const auto elasticityTensor = makeIsotropicElasticityTensor<DataTypes>(mu, lambda);
+
+            const std::array<sofa::Coord_t<DataTypes>, trait::NumberOfNodesInElement> nodesCoordinates = extractNodesVectorFromGlobalVector(element, restPositionAccessor.ref());
+            elementStiffness[elementId] = integrate<DataTypes, ElementType, trait::matrixVectorProductType>(nodesCoordinates, elasticityTensor);
+        });
+}
+
+}