diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index 2602f434a89..28716799cb0 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/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 diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.cpp new file mode 100644 index 00000000000..2030ac9cb0d --- /dev/null +++ b/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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h new file mode 100644 index 00000000000..5163da659da --- /dev/null +++ b/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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +#if !defined(ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP) +#include +#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 BaseElementLinearFEMForceField : public sofa::component::solidmechanics::fem::elastic::BaseLinearElasticityFEMForceField +{ +public: + SOFA_ABSTRACT_CLASS( + SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType), + sofa::component::solidmechanics::fem::elastic::BaseLinearElasticityFEMForceField); + + void init() override; + +private: + using trait = sofa::component::solidmechanics::fem::elastic::trait; + 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 > d_elementStiffness; +}; + +#if !defined(ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseElementLinearFEMForceField; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl new file mode 100644 index 00000000000..0b90dc96768 --- /dev/null +++ b/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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include + +#include + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +BaseElementLinearFEMForceField::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 +void BaseElementLinearFEMForceField::init() +{ + sofa::component::solidmechanics::fem::elastic::BaseLinearElasticityFEMForceField::init(); + + if (!this->isComponentStateInvalid()) + { + this->precomputeElementStiffness(); + } +} + +template +void BaseElementLinearFEMForceField::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(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(youngModulus, poissonRatio); + + const auto elasticityTensor = makeIsotropicElasticityTensor(mu, lambda); + + const std::array, trait::NumberOfNodesInElement> nodesCoordinates = extractNodesVectorFromGlobalVector(element, restPositionAccessor.ref()); + elementStiffness[elementId] = integrate(nodesCoordinates, elasticityTensor); + }); +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp new file mode 100644 index 00000000000..36b65ff9d12 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp @@ -0,0 +1,64 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define ELASTICITY_COMPONENT_ELEMENT_COROTATIONAL_FEM_FORCE_FIELD_CPP + +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerElementCorotationalFEMForceField(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear beams using the corotational approach") + // .add< ElementCorotationalFEMForceField >() + .add< ElementCorotationalFEMForceField >() + .add< ElementCorotationalFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear triangles using the corotational approach") + .add< ElementCorotationalFEMForceField >() + .add< ElementCorotationalFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear quads using the corotational approach") + .add< ElementCorotationalFEMForceField >() + .add< ElementCorotationalFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear tetrahedra using the corotational approach") + .add< ElementCorotationalFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear hexahedra using the corotational approach") + .add< ElementCorotationalFEMForceField >(true)); +} + +// template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h new file mode 100644 index 00000000000..7bfc3c303f0 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h @@ -0,0 +1,203 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#if !defined(ELASTICITY_COMPONENT_ELEMENT_COROTATIONAL_FEM_FORCE_FIELD_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @brief A container for rotation computation methods in corotational formulations + * + * This class provides element-specific rotation computation strategies for common element types. + * + * @tparam DataTypes The data type used throughout the simulation (e.g., sofa::defaulttype::Vec3Types for 3D) + * @tparam ElementType The element type (e.g., sofa::geometry::Triangle, sofa::geometry::Tetrahedron) + * + * Inherits from RotationMethodsContainer with pre-defined rotation strategies. + * + * The class is specialized for some elements because some rotation strategies can be used only + * for specific elements. + */ +template +struct RotationMethods : RotationMethodsContainer, PolarDecomposition, IdentityRotation +> +{ + using Inherit = RotationMethodsContainer, PolarDecomposition, IdentityRotation>; + explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) + {} +}; + +template +using Vec3Real = sofa::defaulttype::StdVectorTypes, sofa::type::Vec<3, Real>, Real>; + +//partial specialization for linear triangle in 3D +template +struct RotationMethods, sofa::geometry::Triangle> : RotationMethodsContainer, sofa::geometry::Triangle, + StablePolarDecomposition>, PolarDecomposition>, IdentityRotation, TriangleRotation> +> +{ + using Inherit = RotationMethodsContainer, sofa::geometry::Triangle, + StablePolarDecomposition>, PolarDecomposition>, IdentityRotation, TriangleRotation> >; + + explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) + { + this->d_rotationMethod.setValue(TriangleRotation>::getItem().key); + } +}; + +//partial specialization for linear tetrahedron +template +struct RotationMethods : RotationMethodsContainer, PolarDecomposition, IdentityRotation, TriangleRotation +> +{ + using Inherit = RotationMethodsContainer, PolarDecomposition, IdentityRotation, TriangleRotation >; + + explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) + { + this->d_rotationMethod.setValue(TriangleRotation::getItem().key); + } +}; + +//partial specialization for linear hexahedron +template +struct RotationMethods : RotationMethodsContainer, PolarDecomposition, IdentityRotation, HexahedronRotation +> +{ + using Inherit = RotationMethodsContainer, PolarDecomposition, IdentityRotation, HexahedronRotation >; + + explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) + { + this->d_rotationMethod.setValue(HexahedronRotation::getItem().key); + } +}; + + +template +class ElementCorotationalFEMForceField : + public BaseElementLinearFEMForceField, + public FEMForceField +{ +public: + SOFA_CLASS2( + SOFA_TEMPLATE2(ElementCorotationalFEMForceField, DataTypes, ElementType), + SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType), + SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType)); + + /** + * The purpose of this function is to register the name of this class according to the provided + * pattern. + * + * Example: ElementCorotationalFEMForceField will produce + * the class name "EdgeCorotationalFEMForceField". + */ + static const std::string GetCustomClassName() + { + return std::string(sofa::geometry::elementTypeToString(ElementType::Element_type)) + + "CorotationalFEMForceField"; + } + + static const std::string GetCustomTemplateName() { return DataTypes::Name(); } + +private: + using trait = sofa::component::solidmechanics::fem::elastic::trait; + using ElementForce = typename trait::ElementForce; + using RotationMatrix = sofa::type::Mat>; + + +public: + + ElementCorotationalFEMForceField(); + + void init() override; + + void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; + + SReal getPotentialEnergy(const sofa::core::MechanicalParams*, const sofa::DataVecCoord_t& x) const override; + + const sofa::type::vector& getElementRotations() const { return m_rotations; } + +protected: + + void beforeElementForce(const sofa::core::MechanicalParams* mparams, + sofa::type::vector& f, + const sofa::VecCoord_t& x) override; + + void computeElementsForces( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& f, + const sofa::VecCoord_t& x) override; + + void computeElementsForcesDeriv( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForcesDeriv, + const sofa::VecDeriv_t& nodeDx) override; + + sofa::type::vector m_rotations; + sofa::type::vector m_initialRotationsTransposed; + + sofa::Coord_t translation(const std::array, trait::NumberOfNodesInElement>& nodes) const; + static sofa::Coord_t computeCentroid(const std::array, trait::NumberOfNodesInElement>& nodes); + + RotationMethods m_rotationMethods; + + void computeRotations(sofa::type::vector& rotations, + const sofa::VecCoord_t& nodePositions, + const sofa::VecCoord_t& nodeRestPositions); + void computeInitialRotations(); +}; + + + +#if !defined(ELASTICITY_COMPONENT_ELEMENT_COROTATIONAL_FEM_FORCE_FIELD_CPP) +// extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl new file mode 100644 index 00000000000..d618649030d --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl @@ -0,0 +1,277 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +ElementCorotationalFEMForceField::ElementCorotationalFEMForceField() + : m_rotationMethods(this) +{ + this->addUpdateCallback("selectRotationMethod", {&this->m_rotationMethods.d_rotationMethod}, + [this](const sofa::core::DataTracker&) + { + m_rotationMethods.selectRotationMethod(); + computeInitialRotations(); + return this->getComponentState(); + }, + {}); + + m_rotationMethods.selectRotationMethod(); +} + +template +void ElementCorotationalFEMForceField::init() +{ + BaseElementLinearFEMForceField::init(); + FEMForceField::init(); + + if (!this->isComponentStateInvalid()) + { + m_rotationMethods.selectRotationMethod(); + computeInitialRotations(); + } + + if (!this->isComponentStateInvalid()) + { + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + } + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + } +} + +template +void ElementCorotationalFEMForceField::beforeElementForce( + const sofa::core::MechanicalParams* mparams, sofa::type::vector& f, + const sofa::VecCoord_t& x) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + m_rotations.resize(elements.size(), RotationMatrix::Identity()); +} + +template +void ElementCorotationalFEMForceField::computeElementsForces( + const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForces, const sofa::VecCoord_t& nodePositions) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + auto restPositionAccessor = this->mstate->readRestPositions(); + auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness); + + for (std::size_t elementId = range.start; elementId < range.end; ++elementId) + { + const auto& element = elements[elementId]; + + const std::array, trait::NumberOfNodesInElement> elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, nodePositions); + const std::array, trait::NumberOfNodesInElement> restElementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, restPositionAccessor.ref()); + + auto& elementInitialRotationTransposed = this->m_initialRotationsTransposed[elementId]; + auto& elementRotation = this->m_rotations[elementId]; + + m_rotationMethods.computeRotation(elementRotation, elementInitialRotationTransposed, elementNodesCoordinates, restElementNodesCoordinates); + + const auto t = translation(elementNodesCoordinates); + const auto t0 = translation(restElementNodesCoordinates); + + typename trait::ElementDisplacement displacement(sofa::type::NOINIT); + for (sofa::Size j = 0; j < trait::NumberOfNodesInElement; ++j) + { + VecView> transformedDisplacement(displacement, j * trait::spatial_dimensions); + transformedDisplacement = elementRotation.multTranspose(elementNodesCoordinates[j] - t) - (restElementNodesCoordinates[j] - t0); + } + + const auto& stiffnessMatrix = elementStiffness[elementId]; + + auto& elementForce = elementForces[elementId]; + elementForce = stiffnessMatrix * displacement; + + for (sofa::Size i = 0; i < trait::NumberOfNodesInElement; ++i) + { + VecView> nodeForce(elementForce, i * trait::spatial_dimensions); + nodeForce = elementRotation * nodeForce; + } + } +} + + +template +void ElementCorotationalFEMForceField::computeElementsForcesDeriv( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForcesDeriv, + const sofa::VecDeriv_t& nodeDx) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness); + + for (std::size_t elementId = range.start; elementId < range.end; ++elementId) + { + const auto& element = elements[elementId]; + + const auto& elementRotation = m_rotations[elementId]; + + sofa::type::Vec> element_dx(sofa::type::NOINIT); + + for (sofa::Size n = 0; n < trait::NumberOfNodesInElement; ++n) + { + VecView> rotated_dx(element_dx, n * trait::spatial_dimensions); + rotated_dx = elementRotation.multTranspose(nodeDx[element[n]]); + } + + const auto& stiffnessMatrix = elementStiffness[elementId]; + + auto& df = elementForcesDeriv[elementId]; + df = stiffnessMatrix * element_dx; + + for (sofa::Size n = 0; n < trait::NumberOfNodesInElement; ++n) + { + VecView> nodedForce(df, n * trait::spatial_dimensions); + nodedForce = elementRotation * nodedForce; + } + } +} + +template +void ElementCorotationalFEMForceField::buildStiffnessMatrix( + sofa::core::behavior::StiffnessMatrix* matrix) +{ + auto dfdx = matrix->getForceDerivativeIn(this->sofa::core::behavior::ForceField::mstate) + .withRespectToPositionsIn(this->sofa::core::behavior::ForceField::mstate); + + sofa::type::Mat> localMatrix(sofa::type::NOINIT); + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness); + + if (m_rotations.size() < elements.size() || elementStiffness.size() < elements.size()) + { + return; + } + + auto elementStiffnessIt = elementStiffness.begin(); + auto rotationMatrixIt = m_rotations.begin(); + for (const auto& element : elements) + { + const auto& elementRotation = *rotationMatrixIt++; + const auto elementRotation_T = elementRotation.transposed(); + + const auto& stiffnessMatrix = *elementStiffnessIt++; + + for (sofa::Index n1 = 0; n1 < trait::NumberOfNodesInElement; ++n1) + { + for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2) + { + stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); // extract the submatrix corresponding to the coupling of nodes n1 and n2 + dfdx(element[n1] * trait::spatial_dimensions, element[n2] * trait::spatial_dimensions) += -elementRotation * localMatrix * elementRotation_T; + } + } + } +} + +template +SReal ElementCorotationalFEMForceField::getPotentialEnergy( + const sofa::core::MechanicalParams*, + const sofa::DataVecCoord_t& x) const +{ + return 0; +} + +template +auto ElementCorotationalFEMForceField::translation( + const std::array, trait::NumberOfNodesInElement>& nodes) const -> sofa::Coord_t +{ + // return nodes[0]; + return computeCentroid(nodes); +} + +template +auto ElementCorotationalFEMForceField::computeCentroid( + const std::array, trait::NumberOfNodesInElement>& nodes) -> sofa::Coord_t +{ + sofa::Coord_t centroid; + for (const auto node : nodes) + { + centroid += node; + } + centroid /= static_cast>(trait::NumberOfNodesInElement); + return centroid; +} + +template +void ElementCorotationalFEMForceField::computeRotations( + sofa::type::vector& rotations, + const sofa::VecCoord_t& nodePositions, + const sofa::VecCoord_t& nodeRestPositions) +{ + if (!this->l_topology) + return; + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + sofa::helper::IotaView indices {static_cast(0ul), elements.size()}; + + rotations.resize(elements.size(), RotationMatrix::Identity()); + if (m_initialRotationsTransposed.size() < elements.size()) + { + m_initialRotationsTransposed.resize(elements.size(), RotationMatrix::Identity()); + } + + std::for_each(indices.begin(), indices.end(), + [&](const auto elementId) + { + const auto& element = elements[elementId]; + + const std::array, trait::NumberOfNodesInElement> elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, nodePositions); + const std::array, trait::NumberOfNodesInElement> restElementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, nodeRestPositions); + + m_rotationMethods.computeRotation(rotations[elementId], m_initialRotationsTransposed[elementId], elementNodesCoordinates, restElementNodesCoordinates); + } + ); +} + +template +void ElementCorotationalFEMForceField::computeInitialRotations() +{ + auto restPositionAccessor = this->sofa::core::behavior::ForceField::mstate->readRestPositions(); + computeRotations(m_initialRotationsTransposed, restPositionAccessor.ref(), restPositionAccessor.ref()); + + for (auto& rotation : m_initialRotationsTransposed) + { + rotation.transpose(); + } +} + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp new file mode 100644 index 00000000000..227b95654d3 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp @@ -0,0 +1,64 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define ELASTICITY_COMPONENT_ELEMENT_LINEAR_SMALL_STRAIN_FEM_FORCE_FIELD_CPP + +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerElementLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear beams assuming small strain") + .add< ElementLinearSmallStrainFEMForceField >() + .add< ElementLinearSmallStrainFEMForceField >() + .add< ElementLinearSmallStrainFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear triangles assuming small strain") + .add< ElementLinearSmallStrainFEMForceField >() + .add< ElementLinearSmallStrainFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear quads assuming small strain") + .add< ElementLinearSmallStrainFEMForceField >() + .add< ElementLinearSmallStrainFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear tetrahedra assuming small strain") + .add< ElementLinearSmallStrainFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear hexahedra assuming small strain") + .add< ElementLinearSmallStrainFEMForceField >(true)); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h new file mode 100644 index 00000000000..f676e561cb3 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h @@ -0,0 +1,111 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +#include + +#if !defined(ELASTICITY_COMPONENT_ELEMENT_LINEAR_SMALL_STRAIN_FEM_FORCE_FIELD_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +class ElementLinearSmallStrainFEMForceField : + public BaseElementLinearFEMForceField, + public FEMForceField +{ +public: + SOFA_CLASS2( + SOFA_TEMPLATE2(ElementLinearSmallStrainFEMForceField, DataTypes, ElementType), + SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType), + SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType)); + + /** + * The purpose of this function is to register the name of this class according to the provided + * pattern. + * + * Example: ElementLinearSmallStrainFEMForceField will produce + * the class name "EdgeLinearSmallStrainFEMForceField". + */ + static const std::string GetCustomClassName() + { + return std::string(sofa::geometry::elementTypeToString(ElementType::Element_type)) + + "LinearSmallStrainFEMForceField"; + } + + static const std::string GetCustomTemplateName() { return DataTypes::Name(); } + +private: + using trait = sofa::component::solidmechanics::fem::elastic::trait; + using ElementStiffness = typename trait::ElementStiffness; + using ElasticityTensor = typename trait::ElasticityTensor; + using ElementDisplacement = typename trait::ElementDisplacement; + using StrainDisplacement = typename trait::StrainDisplacement; + using ElementForce = typename trait::ElementForce; + +public: + void init() override; + + void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; + + SReal getPotentialEnergy(const sofa::core::MechanicalParams*, const sofa::DataVecCoord_t& x) const override; + + using sofa::core::behavior::ForceField::addKToMatrix; + // almost deprecated, but here for compatibility with unit tests + void addKToMatrix(sofa::linearalgebra::BaseMatrix* matrix, SReal kFact, unsigned& offset) override; + +protected: + + void computeElementsForces( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForces, + const sofa::VecCoord_t& nodePositions) override; + + void computeElementsForcesDeriv( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForcesDeriv, + const sofa::VecDeriv_t& nodeDx) override; + +}; + + +#if !defined(ELASTICITY_COMPONENT_ELEMENT_LINEAR_SMALL_STRAIN_FEM_FORCE_FIELD_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl new file mode 100644 index 00000000000..26046d28f00 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl @@ -0,0 +1,183 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +void ElementLinearSmallStrainFEMForceField::init() +{ + BaseElementLinearFEMForceField::init(); + FEMForceField::init(); + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + } +} + + +template +void ElementLinearSmallStrainFEMForceField::computeElementsForces( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForces, + const sofa::VecCoord_t& nodePositions) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + auto restPositionAccessor = this->mstate->readRestPositions(); + auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness); + + for (std::size_t elementId = range.start; elementId < range.end; ++elementId) + { + const auto& element = elements[elementId]; + const auto& stiffnessMatrix = elementStiffness[elementId]; + + typename trait::ElementDisplacement displacement{ sofa::type::NOINIT }; + + for (sofa::Size j = 0; j < trait::NumberOfNodesInElement; ++j) + { + const auto nodeId = element[j]; + for (sofa::Size dim = 0; dim < trait::spatial_dimensions; ++dim) + { + displacement[j * trait::spatial_dimensions + dim] = nodePositions[nodeId][dim] - restPositionAccessor[nodeId][dim]; + } + } + + elementForces[elementId] = stiffnessMatrix * displacement; + } +} + +template +void ElementLinearSmallStrainFEMForceField::computeElementsForcesDeriv( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForcesDeriv, + const sofa::VecDeriv_t& nodeDx) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness); + + for (std::size_t elementId = range.start; elementId < range.end; ++elementId) + { + const auto& element = elements[elementId]; + const auto& stiffnessMatrix = elementStiffness[elementId]; + + const std::array, trait::NumberOfNodesInElement> elementNodesDx = + extractNodesVectorFromGlobalVector(element, nodeDx); + + sofa::type::Vec> element_dx(sofa::type::NOINIT); + for (sofa::Size nodeId = 0; nodeId < trait::NumberOfNodesInElement; ++nodeId) + { + const auto& dx = elementNodesDx[nodeId]; + for (sofa::Size dim = 0; dim < trait::spatial_dimensions; ++dim) + { + element_dx[nodeId * trait::spatial_dimensions + dim] = dx[dim]; + } + } + + elementForcesDeriv[elementId] = stiffnessMatrix * element_dx; + } +} + +template +void ElementLinearSmallStrainFEMForceField::buildStiffnessMatrix( + sofa::core::behavior::StiffnessMatrix* matrix) +{ + if (this->isComponentStateInvalid()) + return; + + auto dfdx = matrix->getForceDerivativeIn(this->mstate) + .withRespectToPositionsIn(this->mstate); + + sofa::type::Mat> localMatrix(sofa::type::NOINIT); + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + + const auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness); + + if (elementStiffness.size() < elements.size()) + { + return; + } + + auto elementStiffnessIt = elementStiffness.begin(); + for (const auto& element : elements) + { + const auto& stiffnessMatrix = *elementStiffnessIt++; + + for (sofa::Index n1 = 0; n1 < trait::NumberOfNodesInElement; ++n1) + { + for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2) + { + stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2 + dfdx(element[n1] * trait::spatial_dimensions, element[n2] * trait::spatial_dimensions) += -localMatrix; + } + } + } +} + +template +SReal ElementLinearSmallStrainFEMForceField::getPotentialEnergy( + const sofa::core::MechanicalParams*, + const sofa::DataVecCoord_t& x) const +{ + return 0; +} + +template +void ElementLinearSmallStrainFEMForceField::addKToMatrix( + sofa::linearalgebra::BaseMatrix* matrix, SReal kFact, unsigned& offset) +{ + if (this->isComponentStateInvalid()) + return; + + using LocalMatType = sofa::type::Mat>; + LocalMatType localMatrix{sofa::type::NOINIT}; + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness); + auto elementStiffnessIt = elementStiffness.begin(); + for (const auto& element : elements) + { + const auto& stiffnessMatrix = *elementStiffnessIt++; + + for (sofa::Index n1 = 0; n1 < trait::NumberOfNodesInElement; ++n1) + { + for (sofa::Index n2 = 0; n2 < trait::NumberOfNodesInElement; ++n2) + { + stiffnessMatrix.getAssembledMatrix().getsub(trait::spatial_dimensions * n1, trait::spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2 + + const auto value = (-static_cast>(kFact)) * static_cast>(localMatrix); + matrix->add( + offset + element[n1] * trait::spatial_dimensions, + offset + element[n2] * trait::spatial_dimensions, value); + } + } + } +} + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.cpp new file mode 100644 index 00000000000..fafab3a6a6b --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define ELASTICITY_COMPONENT_FEM_FORCEFIELD_CPP +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h new file mode 100644 index 00000000000..1e51da58349 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h @@ -0,0 +1,142 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include + +#if !defined(ELASTICITY_COMPONENT_FEM_FORCEFIELD_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +class FEMForceField : + public virtual sofa::core::behavior::ForceField, + public virtual sofa::core::behavior::TopologyAccessor, + public virtual sofa::simulation::TaskSchedulerUser +{ +public: + SOFA_CLASS3(SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType), + sofa::core::behavior::ForceField, + sofa::core::behavior::TopologyAccessor, + sofa::simulation::TaskSchedulerUser); + +private: + using trait = sofa::component::solidmechanics::fem::elastic::trait; + using ElementForce = typename trait::ElementForce; + +public: + void init() override; + + void addForce( + const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& f, + const sofa::DataVecCoord_t& x, + const sofa::DataVecDeriv_t& v) override; + + void addDForce(const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& df, + const sofa::DataVecDeriv_t& dx) override; + + void draw(const sofa::core::visual::VisualParams*) override; + + sofa::Data d_computeForceStrategy; + sofa::Data d_computeForceDerivStrategy; + + sofa::Data> d_elementSpace; + +protected: + + FEMForceField(); + + /// Methods related to addForce + /// @{ + void computeElementsForces(const sofa::core::MechanicalParams* mparams, + sofa::type::vector& f, + const sofa::VecCoord_t& x); + + virtual void beforeElementForce(const sofa::core::MechanicalParams* mparams, + sofa::type::vector& f, + const sofa::VecCoord_t& x) {} + + virtual void computeElementsForces( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& f, + const sofa::VecCoord_t& x) = 0; + + void dispatchElementForcesToNodes( + const sofa::type::vector& elements, + sofa::VecDeriv_t& nodeForces); + /// @} + + + /// Methods related to addDForce + /// @{ + void computeElementsForcesDeriv(const sofa::core::MechanicalParams* mparams, + sofa::type::vector& df, + const sofa::VecDeriv_t& dx); + + virtual void computeElementsForcesDeriv( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& df, + const sofa::VecDeriv_t& dx) = 0; + + /** + * Force derivatives were computed at the element level. This function dispatches the force + * derivatives from the elements to the nodes. + */ + void dispatchElementForcesDerivToNodes(const sofa::core::MechanicalParams* mparams, + const sofa::type::vector& elements, + sofa::VecDeriv_t& nodeForcesDeriv); + /// @} + + sofa::simulation::ForEachExecutionPolicy getExecutionPolicy(const sofa::Data& strategy) const; + + sofa::type::vector>> m_elementForce; + sofa::type::vector>> m_elementDForce; + + sofa::core::visual::DrawElementMesh m_drawMesh; +}; + +#if !defined(ELASTICITY_COMPONENT_FEM_FORCEFIELD_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl new file mode 100644 index 00000000000..41b436a3ca2 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl @@ -0,0 +1,239 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +FEMForceField::FEMForceField() + : d_computeForceStrategy(initData(&d_computeForceStrategy, "computeForceStrategy", std::string("The compute strategy used to compute the forces.\n" + ComputeStrategy::dataDescription()).c_str())) + , d_computeForceDerivStrategy(initData(&d_computeForceDerivStrategy, "computeForceDerivStrategy", std::string("The compute strategy used to compute the forces derivatives.\n" + ComputeStrategy::dataDescription()).c_str())) + , d_elementSpace(initData(&d_elementSpace, static_cast>(0.125), "elementSpace", "When rendering, the space between elements")) +{ + d_elementSpace.setGroup("Visualization"); + + d_computeForceStrategy.setGroup("Multithreading"); + d_computeForceDerivStrategy.setGroup("Multithreading"); +} + +template +void FEMForceField::init() +{ + sofa::core::behavior::ForceField::init(); + + if (!this->isComponentStateInvalid()) + { + TopologyAccessor::init(); + } + + if (!this->isComponentStateInvalid()) + { + this->initTaskScheduler(); + } +} + +template +void FEMForceField::addForce( + const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& f, + const sofa::DataVecCoord_t& x, + const sofa::DataVecDeriv_t& v) +{ + SOFA_UNUSED(mparams); + SOFA_UNUSED(v); + + auto positionAccessor = sofa::helper::getReadAccessor(x); + + if (this->l_topology == nullptr) return; + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + m_elementForce.resize(elements.size()); + + this->computeElementsForces(mparams, m_elementForce, positionAccessor.ref()); + + auto forceAccessor = sofa::helper::getWriteOnlyAccessor(f); + if (forceAccessor.size() < positionAccessor.size()) + { + forceAccessor.resize(positionAccessor.size()); + } + + // dispatch the element force to the degrees of freedom. + // this operation is done outside the compute strategy because it is not thread-safe. + dispatchElementForcesToNodes(elements, forceAccessor.wref()); +} + +template +void FEMForceField::computeElementsForces( + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& f, + const sofa::VecCoord_t& x) +{ + SCOPED_TIMER("ElementForces"); + this->beforeElementForce(mparams, f, x); + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + + sofa::simulation::forEachRange(getExecutionPolicy(d_computeForceStrategy), *this->m_taskScheduler, + static_cast(0), elements.size(), [this, mparams, &f, &x](const auto& range) + { + SCOPED_TIMER_TR("ElementForcesRange"); + this->computeElementsForces(range, mparams, f, x); + }); +} + +template +void FEMForceField::dispatchElementForcesToNodes( + const sofa::type::vector& elements, + sofa::VecDeriv_t& nodeForces) +{ + SCOPED_TIMER("DispatchElementForces"); + + for (sofa::Size i = 0; i < elements.size(); ++i) + { + const auto& element = elements[i]; + const auto& elementForce = m_elementForce[i]; + + for (sofa::Size j = 0; j < trait::NumberOfNodesInElement; ++j) + { + auto& nodeForce = nodeForces[element[j]]; + for (sofa::Size k = 0; k < trait::spatial_dimensions; ++k) + { + nodeForce[k] -= elementForce[j * trait::spatial_dimensions + k]; + } + } + } +} + +template +void FEMForceField::addDForce( + const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& df, + const sofa::DataVecDeriv_t& dx) +{ + if (this->isComponentStateInvalid()) + return; + + auto dfAccessor = sofa::helper::getWriteAccessor(df); + auto dxAccessor = sofa::helper::getReadAccessor(dx); + if (dxAccessor.size() != dfAccessor.size()) + { + dfAccessor.resize(dxAccessor.size()); + } + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + + if (m_elementDForce.size() != elements.size()) + { + m_elementDForce.resize(elements.size()); + } + + this->computeElementsForcesDeriv(mparams, m_elementDForce, dxAccessor.ref()); + + // dispatch the element dforce to the degrees of freedom. + // this operation is done outside the compute strategy because it is not thread-safe. + dispatchElementForcesDerivToNodes(mparams, elements, dfAccessor.wref()); +} + +template +void FEMForceField::computeElementsForcesDeriv( + const sofa::core::MechanicalParams* mparams, sofa::type::vector& df, + const sofa::VecDeriv_t& dx) +{ + SCOPED_TIMER("ElementForcesDeriv"); + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + + sofa::simulation::forEachRange(getExecutionPolicy(d_computeForceDerivStrategy), *this->m_taskScheduler, + static_cast(0), elements.size(), + [this, mparams, &df, &dx](const sofa::simulation::Range& range) + { + SCOPED_TIMER_TR("ElementForcesDerivRange"); + this->computeElementsForcesDeriv(range, mparams, df, dx); + }); +} + +template +void FEMForceField::dispatchElementForcesDerivToNodes(const sofa::core::MechanicalParams* mparams, + const sofa::type::vector& elements, + sofa::VecDeriv_t& nodeForcesDeriv) +{ + SCOPED_TIMER("DispatchElementForcesDeriv"); + + const auto kFactor = static_cast>(sofa::core::mechanicalparams::kFactorIncludingRayleighDamping( + mparams, this->rayleighStiffness.getValue())); + + for (std::size_t elementId = 0; elementId < elements.size(); ++elementId) + { + const auto& element = elements[elementId]; + const auto& elementDForce = m_elementDForce[elementId]; + + for (sofa::Size i = 0; i < trait::NumberOfNodesInElement; ++i) + { + const auto nodeId = element[i]; + auto& df = nodeForcesDeriv[nodeId]; + + for (sofa::Size dim = 0; dim < trait::spatial_dimensions; ++dim) + { + df[dim] -= kFactor * elementDForce[i * trait::spatial_dimensions + dim]; + } + } + } +} + + +template +sofa::simulation::ForEachExecutionPolicy FEMForceField::getExecutionPolicy( + const sofa::Data& strategy) const +{ + auto computeForceStrategyAccessor = sofa::helper::getReadAccessor(d_computeForceStrategy); + const auto& computeForceStrategy = computeForceStrategyAccessor->key(); + + return (computeForceStrategy == parallelComputeStrategy) + ? sofa::simulation::ForEachExecutionPolicy::PARALLEL + : sofa::simulation::ForEachExecutionPolicy::SEQUENTIAL; +} + +template +void FEMForceField::draw(const sofa::core::visual::VisualParams* vparams) +{ + if (!vparams->displayFlags().getShowForceFields()) + return; + + if (!this->l_topology) + return; + + const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); + + if (vparams->displayFlags().getShowWireFrame()) + vparams->drawTool()->setPolygonMode(0, true); + + const auto& x = this->mstate->read(sofa::core::vec_id::read_access::position)->getValue(); + + m_drawMesh.elementSpace = d_elementSpace.getValue(); + m_drawMesh.drawAllElements(vparams->drawTool(), x, this->l_topology.get()); +} + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement.h new file mode 100644 index 00000000000..4c04b76052b --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement.h @@ -0,0 +1,84 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct FiniteElement; + +#define FINITEELEMENT_HEADER(ElType, DataTypes, dimension) \ + using Coord = sofa::Coord_t;\ + using Real = sofa::Real_t;\ + using ElementType = ElType;\ + using TopologyElement = sofa::topology::Element;\ + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions;\ + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes;\ + static constexpr sofa::Size TopologicalDimension = dimension;\ + using ReferenceCoord = sofa::type::Vec;\ + using ShapeFunctionType = std::function;\ + using QuadraturePoint = ReferenceCoord; \ + using QuadraturePointAndWeight = std::pair + + +template +struct FiniteElementHelper +{ + using FiniteElement = sofa::component::solidmechanics::fem::elastic::FiniteElement; + using Coord = typename FiniteElement::Coord; + using Real = typename FiniteElement::Real; + + static constexpr sofa::Size spatial_dimensions = FiniteElement::spatial_dimensions; + static constexpr sofa::Size NumberOfNodesInElement = FiniteElement::NumberOfNodesInElement; + static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; + + // gradient of shape functions in the reference element evaluated at the quadrature point + static constexpr auto gradientShapeFunctionAtQuadraturePoints() + { + using Gradient = sofa::type::Mat; + constexpr auto quadraturePoints = FiniteElement::quadraturePoints(); + + std::array gradients; + std::transform(quadraturePoints.begin(), quadraturePoints.end(), gradients.begin(), + [](const auto& qp) { return FiniteElement::gradientShapeFunctions(qp.first); }); + return gradients; + } + + // jacobian of the mapping from the reference space to the physical space, evaluated where the + // gradient of the shape functions has been evaluated. + static constexpr auto jacobianFromReferenceToPhysical( + const std::array& elementNodesCoordinates, + const sofa::type::Mat& gradientShapeFunctionInReferenceElement) + { + sofa::type::Mat jacobian; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + jacobian += sofa::type::dyad(elementNodesCoordinates[i], gradientShapeFunctionInReferenceElement[i]); + } + return jacobian; + } + +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Edge].h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Edge].h new file mode 100644 index 00000000000..29cdbf9944d --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Edge].h @@ -0,0 +1,56 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::Edge, DataTypes, 1); + + constexpr static std::array referenceElementNodes {{ReferenceCoord{-1}, ReferenceCoord{1}}}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getEdges(); + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + SOFA_UNUSED(q); + return {{-static_cast(0.5)}, {static_cast(0.5)}}; + } + + static constexpr std::array quadraturePoints() + { + constexpr sofa::type::Vec q0(static_cast(0)); + return { + std::make_pair(q0, static_cast(2)) + }; + } + +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Hexahedron].h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Hexahedron].h new file mode 100644 index 00000000000..153c0003f65 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Hexahedron].h @@ -0,0 +1,88 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::Hexahedron, DataTypes, 3); + static_assert(spatial_dimensions == 3, "Hexahedrons are only defined in 3D"); + + constexpr static std::array referenceElementNodes {{ + {-1, -1, -1}, + {1, -1, -1}, + {1, 1, -1}, + {-1, 1, -1}, + {-1, -1, 1}, + {1, -1, 1}, + {1, 1, 1}, + {-1, 1, 1}, + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getHexahedra(); + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + const auto [x, y, z] = q; + sofa::type::Mat gradient(sofa::type::NOINIT); + using Line = typename sofa::type::Mat::Line; + + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + const auto& [xref, yref, zref] = referenceElementNodes[i]; + gradient[i] = 1./8. * Line( + xref * (1 + y * yref) * (1 + z * zref), + yref * (1 + x * xref) * (1 + z * zref), + zref * (1 + x * xref) * (1 + y * yref)); + } + + return gradient; + } + + static constexpr std::array quadraturePoints() + { + constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.) + constexpr Real sqrt3 = 1.73205080757; //sqrt(3.) + + constexpr sofa::type::Vec q0(0., sqrt2_3, -1./sqrt3); + constexpr sofa::type::Vec q1(0., -sqrt2_3, -1./sqrt3); + constexpr sofa::type::Vec q2(-sqrt2_3, 0., 1./sqrt3); + constexpr sofa::type::Vec q3(sqrt2_3, 0., 1./sqrt3); + + constexpr std::array q { + std::make_pair(q0, static_cast(2)), + std::make_pair(q1, static_cast(2)), + std::make_pair(q2, static_cast(2)), + std::make_pair(q3, static_cast(2)) + }; + + return q; + } +}; +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Quad].h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Quad].h new file mode 100644 index 00000000000..2d7fc812fa7 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Quad].h @@ -0,0 +1,74 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::Quad, DataTypes, 2); + static_assert(spatial_dimensions > 1, "Quads cannot be defined in 1D"); + + constexpr static std::array referenceElementNodes {{ + {-1, -1}, + {1, -1}, + {1, 1}, + {-1, 1} + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getQuads(); + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + return { + {1 / static_cast(4) * (-static_cast(1) + q[1]), 1 / static_cast(4) * (-static_cast(1) + q[0])}, + {1 / static_cast(4) * ( static_cast(1) - q[1]), 1 / static_cast(4) * (-static_cast(1) - q[0])}, + {1 / static_cast(4) * ( static_cast(1) + q[1]), 1 / static_cast(4) * ( static_cast(1) + q[0])}, + {1 / static_cast(4) * (-static_cast(1) - q[1]), 1 / static_cast(4) * ( static_cast(1) - q[0])} + }; + } + + static constexpr std::array quadraturePoints() + { + constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.) + constexpr Real sqrt6 = 2.44948974278; //sqrt(6.) + constexpr Real sqrt2 = 1.41421356237; //sqrt(2.) + + constexpr sofa::type::Vec q0(sqrt2_3, 0.); + constexpr sofa::type::Vec q1(-1/sqrt6, -1./sqrt2); + constexpr sofa::type::Vec q2(-1/sqrt6, 1./sqrt2); + + return { + std::make_pair(q0, 4./3.), + std::make_pair(q1, 4./3.), + std::make_pair(q2, 4./3.), + }; + } +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Tetrahedron].h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Tetrahedron].h new file mode 100644 index 00000000000..028689ce3d2 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Tetrahedron].h @@ -0,0 +1,65 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::Tetrahedron, DataTypes, 3); + static_assert(spatial_dimensions == 3, "Tetrahedrons are only defined in 3D"); + + constexpr static std::array referenceElementNodes {{ + {0, 0, 0}, + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getTetrahedra(); + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + SOFA_UNUSED(q); + return { + {-1, -1, -1}, + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }; + } + + static constexpr std::array quadraturePoints() + { + constexpr sofa::type::Vec q0(1./4., 1./4., 1./4.); + constexpr std::array q { std::make_pair(q0, 1./6.) }; + return q; + } +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Triangle].h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Triangle].h new file mode 100644 index 00000000000..5e6d61909cd --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[Triangle].h @@ -0,0 +1,62 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::Triangle, DataTypes, 2); + static_assert(spatial_dimensions > 1, "Triangles cannot be defined in 1D"); + + constexpr static std::array referenceElementNodes {{ + {0, 0}, + {1, 0}, + {0, 1}}}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getTriangles(); + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + SOFA_UNUSED(q); + return { + {-1, -1}, + {1, 0}, + {0, 1} + }; + } + + static constexpr std::array quadraturePoints() + { + return { + std::make_pair(sofa::type::Vec(1./3., 1./3.), 1./2.) + }; + } +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[all].h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[all].h new file mode 100644 index 00000000000..ebbf7d8fb91 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/finiteelement/FiniteElement[all].h @@ -0,0 +1,28 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include +#include diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ComputeStrategy.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ComputeStrategy.h new file mode 100644 index 00000000000..99970790c20 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ComputeStrategy.h @@ -0,0 +1,36 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +constexpr std::string_view parallelComputeStrategy = "parallel"; +constexpr std::string_view sequencedComputeStrategy = "sequenced"; + +MAKE_SELECTABLE_ITEMS(ComputeStrategy, + sofa::helper::Item{parallelComputeStrategy, "The algorithm is executed in parallel"}, + sofa::helper::Item{sequencedComputeStrategy, "The algorithm is executed sequentially"}, +); + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h new file mode 100644 index 00000000000..cf90443d52c --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h @@ -0,0 +1,205 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +using ElementStiffness = sofa::type::Mat< + ElementType::NumberOfNodes * DataTypes::spatial_dimensions, + ElementType::NumberOfNodes * DataTypes::spatial_dimensions, + sofa::Real_t +>; + +/** + * Specifies the type of matrix-vector product to be used with the stiffness matrix. + */ +enum class MatrixVectorProductType +{ + /** + * The matrix-vector product is computed using the factorization of the matrix + */ + Factorization, + + /** + * The matrix-vector product is computed using the dense matrix representation + */ + Dense +}; + +/** + * Represents an element stiffness matrix. It contains the dense matrix representation, but also + * its factorization. Both the factorization or the dense matrix can be used for the product of the + * stiffness matrix with a vector. Although it gives the same result, it has an impact on the + * performances. + */ +template +struct FactorizedElementStiffness +{ +private: + using Real = sofa::Real_t; + using FiniteElement = sofa::component::solidmechanics::fem::elastic::FiniteElement; + static constexpr auto NbQuadraturePoints = FiniteElement::quadraturePoints().size(); + static constexpr sofa::Size NumberOfIndependentElements = symmetric_tensor::NumberOfIndependentElements; + + /// The factors of the stiffness matrix + /// @{ + std::array, NbQuadraturePoints> B; + IsotropicElasticityTensor elasticityTensor; + std::array factors; + /// @} + + /** + * The assembled stiffness matrix + */ + sofa::type::Mat< + ElementType::NumberOfNodes * DataTypes::spatial_dimensions, + ElementType::NumberOfNodes * DataTypes::spatial_dimensions, + sofa::Real_t> stiffnessMatrix; + +public: + const StrainDisplacement& getB(std::size_t i) const { return B[i]; } + Real getFactor(std::size_t i) const { return factors[i]; } + const IsotropicElasticityTensor& getElasticityTensor() const { return elasticityTensor; } + + void setElasticityTensor(const FullySymmetric4Tensor& elasticityTensor_) + { + for (std::size_t i = 0; i < NumberOfIndependentElements; ++i) + { + for (std::size_t j = 0; j < NumberOfIndependentElements; ++j) + { + this->elasticityTensor(i, j) = elasticityTensor_.toVoigtMatSym()(i, j); + } + } + } + + void addFactor(std::size_t quadraturePointIndex, + const Real factor, + const StrainDisplacement& B_) + { + this->factors[quadraturePointIndex] = factor; + this->B[quadraturePointIndex] = B_; + + const auto K_i = factor * B_.multTranspose(elasticityTensor.toMat() * B_.B); + stiffnessMatrix += K_i; + } + + const auto& getAssembledMatrix() const { return stiffnessMatrix; } + + using Vec = sofa::type::Vec< + ElementType::NumberOfNodes * DataTypes::spatial_dimensions, + sofa::Real_t>; + + inline Vec operator*(const Vec& v) const + { + if constexpr (matrixVectorProductType == MatrixVectorProductType::Factorization) + { + if constexpr (NbQuadraturePoints > 1) + { + Vec result; + for (std::size_t i = 0; i < NbQuadraturePoints; ++i) + { + const auto& B = this->B[i]; + result += B.multTranspose(this->factors[i] * (this->elasticityTensor * (B * v))); + } + return result; + } + else + { + return this->B[0].multTranspose(this->factors[0] * (this->elasticityTensor * (this->B[0] * v))); + } + } + else + { + return this->stiffnessMatrix * v; + } + } + + friend std::ostream& operator<<(std::ostream& os, const FactorizedElementStiffness& obj) + { + return os << obj.getAssembledMatrix(); + } + + friend std::istream& operator>>(std::istream& in, FactorizedElementStiffness& obj) + { + return in; + } +}; + +template +FactorizedElementStiffness integrate( + const std::array, ElementType::NumberOfNodes>& nodesCoordinates, + const FullySymmetric4Tensor& elasticityTensor) +{ + using Real = sofa::Real_t; + using FiniteElement = FiniteElement; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; + + FactorizedElementStiffness K; + K.setElasticityTensor(elasticityTensor); + + std::size_t quadraturePointIndex = 0; + for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) + { + // gradient of shape functions in the reference element evaluated at the quadrature point + const sofa::type::Mat dN_dq_ref = + FiniteElement::gradientShapeFunctions(quadraturePoint); + + // jacobian of the mapping from the reference space to the physical space, evaluated at the + // quadrature point + sofa::type::Mat jacobian; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + jacobian += sofa::type::dyad(nodesCoordinates[i], dN_dq_ref[i]); + + const auto detJ = sofa::component::solidmechanics::fem::elastic::absGeneralizedDeterminant(jacobian); + const sofa::type::Mat J_inv = + sofa::component::solidmechanics::fem::elastic::inverse(jacobian); + + // gradient of the shape functions in the physical element evaluated at the quadrature point + sofa::type::Mat dN_dq(sofa::type::NOINIT); + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + dN_dq[i] = J_inv.transposed() * dN_dq_ref[i]; + + const auto B = makeStrainDisplacement(dN_dq); + K.addFactor(quadraturePointIndex, weight * detJ, B); + + ++quadraturePointIndex; + } + return K; +} + + + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/FullySymmetric4Tensor.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/FullySymmetric4Tensor.h new file mode 100644 index 00000000000..fd20bec8107 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/FullySymmetric4Tensor.h @@ -0,0 +1,156 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * A class to represent a fully symmetric 4th rank tensor. + * + * It's a tensor having both minor and major symmetries. Given the indices i,j,k,l, a fully symmetric + * tensor C has the following properties: + * C(i,j,k,l) = C(k,l,i,j) (major symmetry) + * C(i,j,k,l) = C(j,i,k,l) (minor symmetry) + * C(i,j,k,l) = C(i,j,l,k) (minor symmetry) + */ +template +class FullySymmetric4Tensor +{ +private: + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + static constexpr sofa::Size NumberOfIndependentElements = symmetric_tensor::NumberOfIndependentElements; + using Real = sofa::Real_t; + +public: + FullySymmetric4Tensor() = default; + + template + explicit FullySymmetric4Tensor(Callable callable) : m_matrix(sofa::type::NOINIT) + { + fill(callable); + } + + template + void fill(Callable callable) + { + SCOPED_TIMER_TR("fillFullySymmetric4Tensor"); + +#ifndef NDEBUG + checkSymmetry(callable); +#endif + + for (sofa::Size a = 0; a < NumberOfIndependentElements; ++a) + { + const auto [i, j] = toTensorIndices(a); + for (sofa::Size b = a; b < NumberOfIndependentElements; ++b) // the Voigt representation is symmetric, that is why b starts at a + { + const auto [k, l] = toTensorIndices(b); + m_matrix(a, b) = callable(i, j, k, l); + } + } + } + + Real& operator()(sofa::Size i, sofa::Size j, sofa::Size k, sofa::Size l) + { + const auto a = tensorToVoigtIndex(i, j); + const auto b = tensorToVoigtIndex(k, l); + return m_matrix(a, b); + } + + Real operator()(sofa::Size i, sofa::Size j, sofa::Size k, sofa::Size l) const + { + const auto a = tensorToVoigtIndex(i, j); + const auto b = tensorToVoigtIndex(k, l); + return m_matrix(a, b); + } + + const sofa::type::MatSym& toVoigtMatSym() const + { + return m_matrix; + } + +private: + sofa::type::MatSym m_matrix; + + template + static void checkSymmetry(Callable callable) + { + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + for (sofa::Size j = 0; j < spatial_dimensions; ++j) + { + for (sofa::Size k = 0; k < spatial_dimensions; ++k) + { + for (sofa::Size l = 0; l < spatial_dimensions; ++l) + { + const auto ijkl = callable(i, j, k, l); + const auto jikl = callable(i, j, k, l); + const auto klij = callable(k, l, i, j); + const auto ijlk = callable(i, j, l, k); + constexpr auto max_precision{std::numeric_limits::digits10 + 1}; + msg_error_when(std::abs(ijkl - klij) > 1e-6, "ElasticityTensor") << "No major symmetry (ij) <-> (kl) " << std::setprecision(max_precision) << ijkl << " != " << klij; + msg_error_when(std::abs(ijkl - jikl) > 1e-6, "ElasticityTensor") << "No minor symmetry i <-> j " << std::setprecision(max_precision) << ijkl << " != " << jikl; + msg_error_when(std::abs(ijkl - ijlk) > 1e-6, "ElasticityTensor") << "No minor symmetry k <-> l " << std::setprecision(max_precision) << ijkl << " != " << ijlk; + } + } + } + } + } +}; + + +/** + * @brief Creates an isotropic elasticity tensor for given material properties. + * + * This function constructs and returns an elasticity tensor for an isotropic material + * characterized by its Young's modulus and Poisson's ratio. It computes the tensor + * using the Lamé parameters, which are derived from the given material properties. + * + * @param youngModulus The Young's modulus of the material, representing its stiffness. + * @param poissonRatio The Poisson's ratio of the material, describing its deformation behavior. + * @return The isotropic elasticity tensor represented as an object of + * `ElasticityTensor`. + */ +template +FullySymmetric4Tensor makeIsotropicElasticityTensor(sofa::Real_t mu, sofa::Real_t lambda) +{ + using Real = sofa::Real_t; + + return FullySymmetric4Tensor{ + [mu, lambda](sofa::Index i, sofa::Index j, sofa::Index k, sofa::Index l) + { + return mu * (kroneckerDelta(i, k) * kroneckerDelta(j, l) + kroneckerDelta(i, l) * kroneckerDelta(j, k)) + + lambda * kroneckerDelta(i, j) * kroneckerDelta(k, l); + }}; +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/IsotropicElasticityTensor.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/IsotropicElasticityTensor.h new file mode 100644 index 00000000000..c7fcc6f649f --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/IsotropicElasticityTensor.h @@ -0,0 +1,94 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +constexpr sofa::type::Vec isotropicElasticityTensorProduct(const sofa::type::Mat& tensor, const sofa::type::Vec& v) +{ + return tensor * v; +} + +/** + * Specialization for 3D where the operations containing known zeros in the elasticity tensor are + * omitted. + */ +template +constexpr sofa::type::Vec<6, real> isotropicElasticityTensorProduct(const sofa::type::Mat<6, 6, real>& tensor, const sofa::type::Vec<6, real>& v) +{ + sofa::type::Vec<6, real> result { sofa::type::NOINIT }; + + result[0] = tensor(0, 0) * v[0] + tensor(0, 1) * v[1] + tensor(0, 2) * v[2]; + result[1] = tensor(1, 0) * v[0] + tensor(1, 1) * v[1] + tensor(1, 2) * v[2]; + result[2] = tensor(2, 0) * v[0] + tensor(2, 1) * v[1] + tensor(2, 2) * v[2]; + result[3] = tensor(3, 3) * v[3]; + result[4] = tensor(4, 4) * v[4]; + result[5] = tensor(5, 5) * v[5]; + + return result; +} + +template +struct IsotropicElasticityTensor +{ + static constexpr auto spatial_dimensions = DataType::spatial_dimensions; + static constexpr auto NbIndependentElements = symmetric_tensor::NumberOfIndependentElements; + + explicit IsotropicElasticityTensor(const sofa::type::Mat>& mat) : C(mat) {} + IsotropicElasticityTensor() = default; + + sofa::type::Vec> operator*(const sofa::type::Vec>& v) const + { + return isotropicElasticityTensorProduct(C, v); + } + + template + sofa::type::Mat> operator*(const sofa::type::Mat>& v) const noexcept + { + return C * v; + } + + sofa::Real_t operator()(sofa::Size i, sofa::Size j) const + { + return C(i, j); + } + + sofa::Real_t& operator()(sofa::Size i, sofa::Size j) + { + return C(i, j); + } + + const sofa::type::Mat>& toMat() const + { + return C; + } + +private: + sofa::type::Mat> C; + +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/KroneckerDelta.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/KroneckerDelta.h new file mode 100644 index 00000000000..9f7934ff07d --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/KroneckerDelta.h @@ -0,0 +1,34 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +constexpr static Real kroneckerDelta(std::size_t i, std::size_t j) +{ + return static_cast(i == j); +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/LameParameters.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/LameParameters.h new file mode 100644 index 00000000000..acc3a43acc0 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/LameParameters.h @@ -0,0 +1,53 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ +/** + * @brief Converts Young's modulus and Poisson's ratio to Lamé parameters. + * + * This function calculates and returns the two Lamé parameters, μ (shear modulus) + * and λ, derived from the given Young’s modulus and Poisson’s ratio of a material. + * These parameters are fundamental in describing isotropic elastic behavior. + * + * @param youngModulus The Young's modulus of the material, representing its stiffness. + * @param poissonRatio The Poisson's ratio of the material, describing its deformation behavior. + * @return A pair containing the calculated Lamé parameters: + * - First: μ (shear modulus), + * - Second: λ. + */ +template +std::pair, sofa::Real_t> +toLameParameters(sofa::Real_t youngModulus, sofa::Real_t poissonRatio) +{ + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + const auto mu = youngModulus / (2 * (1 + poissonRatio)); + const auto lambda = youngModulus * poissonRatio / ((1 + poissonRatio) * (1 - (spatial_dimensions - 1) * poissonRatio)); + return std::make_pair(mu, lambda); +} +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MajorSymmetric4Tensor.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MajorSymmetric4Tensor.h new file mode 100644 index 00000000000..cdc795ba5eb --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MajorSymmetric4Tensor.h @@ -0,0 +1,96 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +namespace elasticity +{ + +/** + * A class to represent a major symmetric 4th rank tensor. + * + * Given the indices i,j,k,l, a major symmetric tensor C has the following properties: + * C(i,j,k,l) = C(k,l,i,j) (major symmetry) + */ +template +class MajorSymmetric4Tensor +{ +private: + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + static constexpr sofa::Size spatial_dimension_square = spatial_dimensions * spatial_dimensions; + using Real = sofa::Real_t; + +public: + MajorSymmetric4Tensor() = delete; + + template + MajorSymmetric4Tensor(Callable callable) : m_matrix(sofa::type::NOINIT) + { + fill(callable); + } + + template + void fill(Callable callable) + { + SCOPED_TIMER_TR("fillMajorSymmetric4Tensor"); + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + for (sofa::Size j = 0; j < spatial_dimensions; ++j) + { + const auto ij = i * spatial_dimensions + j; + for (sofa::Size k = 0; k < spatial_dimensions; ++k) + { + for (sofa::Size l = 0; l < spatial_dimensions; ++l) + { + const auto kl = k * spatial_dimensions + l; + if (kl <= ij) + { + m_matrix(ij, kl) = callable(i, j, k, l); + } + } + } + } + } + } + + Real& operator()(sofa::Size i, sofa::Size j, sofa::Size k, sofa::Size l) + { + const auto a = i * spatial_dimensions + j; + const auto b = k * spatial_dimensions + l; + return m_matrix(a, b); + } + + Real operator()(sofa::Size i, sofa::Size j, sofa::Size k, sofa::Size l) const + { + const auto a = i * spatial_dimensions + j; + const auto b = k * spatial_dimensions + l; + return m_matrix(a, b); + } + +private: + sofa::type::MatSym m_matrix; +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MatrixTools.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MatrixTools.h new file mode 100644 index 00000000000..06dc755041c --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/MatrixTools.h @@ -0,0 +1,118 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * Alias for a sofa::type::Mat + * If T is sofa::type::Mat and L==1 && C==1, the alias is the scalar type. + * Otherwise, the alias is T itself. + * + * Example: static_cast>(matrix) + */ +template +using ScalarOrMatrix = std::conditional_t< + (T::nbLines==1 && T::nbCols==1), typename T::Real, T>; + +template +real determinantSquareMatrix(const sofa::type::Mat& mat) +{ + if constexpr (N == 1) + { + return mat(0, 0); + } + else if constexpr (N == 2) + { + return mat(0, 0) * mat(1, 1) - mat(0, 1) * mat(1, 0); + } + else + { + real det = 0; + for (size_t p = 0; p < N; ++p) + { + sofa::type::Mat submat; + for (size_t i = 1; i < N; ++i) + { + size_t colIndex = 0; + for (size_t j = 0; j < N; ++j) + { + if (j == p) continue; + submat(i - 1, colIndex++) = mat(i, j); + } + } + det += ((p % 2 == 0) ? 1 : -1) * mat(0, p) * determinantSquareMatrix(submat); + } + return det; + } +} + +template +real determinant(const sofa::type::Mat& mat) +{ + return determinantSquareMatrix(mat); +} + +template +real absGeneralizedDeterminant(const sofa::type::Mat& mat) +{ + if constexpr (L == C) + { + return std::abs(determinantSquareMatrix(mat)); + } + else + { + return std::sqrt(determinantSquareMatrix(mat.multTranspose(mat))); + } +} + +template +sofa::type::Mat leftPseudoInverse(const sofa::type::Mat& mat) +{ + return mat.multTranspose(mat).inverted() * mat.transposed(); +} + +/** + * Computes the inverse of a given matrix. + * For square matrices (L == C), the standard matrix inverse is computed. + * For non-square matrices, the left pseudo-inverse is returned. + * + * @param mat The input matrix of size LxC to be inverted or pseudo-inverted. + * @return A matrix of size CxL representing the inverse or left pseudo-inverse of the input matrix. + */ +template +sofa::type::Mat inverse(const sofa::type::Mat& mat) +{ + if constexpr (L == C) + { + return mat.inverted(); + } + else + { + return leftPseudoInverse(mat); + } +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/StrainDisplacement.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/StrainDisplacement.h new file mode 100644 index 00000000000..addc9c6263d --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/StrainDisplacement.h @@ -0,0 +1,183 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct StrainDisplacement; + + +template +sofa::type::Vec strainDisplacementVectorProduct(const sofa::type::Mat& B, const sofa::type::Vec& v) +{ + return B * v; +} + +template +sofa::type::Vec strainDisplacementTransposedVectorProduct( + const sofa::type::Mat& B, const sofa::type::Vec& v) +{ + return B.multTranspose(v); +} + +/** + * Specialization for a linear tetrahedron where the operations containing known zeros in the strain + * displacement tensor are omitted. + */ +template +sofa::type::Vec<6, real> strainDisplacementVectorProduct(const sofa::type::Mat<6, 12, real>& B, const sofa::type::Vec<12, real>& v) +{ + sofa::type::Vec<6, real> Bv{sofa::type::NOINIT}; + + Bv[0] = B(0, 0) * v[0] + B(0, 3) * v[3] + B(0, 6) * v[6] + B(0, 9) * v[9]; + Bv[1] = B(1, 1) * v[1] + B(1, 4) * v[4] + B(1, 7) * v[7] + B(1, 10) * v[10]; + Bv[2] = B(2, 2) * v[2] + B(2, 5) * v[5] + B(2, 8) * v[8] + B(2, 11) * v[11]; + Bv[3] = B(3, 0) * v[0] + B(3, 1) * v[1] + B(3, 3) * v[3] + B(3, 4) * v[4] + B(3, 6) * v[6] + B(3, 7) * v[7] + B(3, 9) * v[9] + B(3, 10) * v[10]; + Bv[4] = B(4, 0) * v[0] + B(4, 2) * v[2] + B(4, 3) * v[3] + B(4, 5) * v[5] + B(4, 6) * v[6] + B(4, 8) * v[8] + B(4, 9) * v[9] + B(4, 11) * v[11]; + Bv[5] = B(5, 1) * v[1] + B(5, 2) * v[2] + B(5, 4) * v[4] + B(5, 5) * v[5] + B(5, 7) * v[7] + B(5, 8) * v[8] + B(5, 10) * v[10] + B(5, 11) * v[11]; + + return Bv; +} + +/** + * Specialization for a linear tetrahedron where the operations containing known zeros in the strain + * displacement tensor are omitted. + */ +template +sofa::type::Vec<12, real> strainDisplacementTransposedVectorProduct(const sofa::type::Mat<6, 12, real>& B, const sofa::type::Vec<6, real>& v) +{ + sofa::type::Vec<12, real> B_Tv { sofa::type::NOINIT }; + + B_Tv[0] = B(0, 0) * v[0] + B(3, 0) * v[3] + B(4, 0) * v[4]; + B_Tv[1] = B(1, 1) * v[1] + B(3, 1) * v[3] + B(5, 1) * v[5]; + B_Tv[2] = B(2, 2) * v[2] + B(4, 2) * v[4] + B(5, 2) * v[5]; + B_Tv[3] = B(0, 3) * v[0] + B(3, 3) * v[3] + B(4, 3) * v[4]; + B_Tv[4] = B(1, 4) * v[1] + B(3, 4) * v[3] + B(5, 4) * v[5]; + B_Tv[5] = B(2, 5) * v[2] + B(4, 5) * v[4] + B(5, 5) * v[5]; + B_Tv[6] = B(0, 6) * v[0] + B(3, 6) * v[3] + B(4, 6) * v[4]; + B_Tv[7] = B(1, 7) * v[1] + B(3, 7) * v[3] + B(5, 7) * v[5]; + B_Tv[8] = B(2, 8) * v[2] + B(4, 8) * v[4] + B(5, 8) * v[5]; + B_Tv[9] = B(0, 9) * v[0] + B(3, 9) * v[3] + B(4, 9) * v[4]; + B_Tv[10] = B(1, 10) * v[1] + B(3, 10) * v[3] + B(5, 10) * v[5]; + B_Tv[11] = B(2, 11) * v[2] + B(4, 11) * v[4] + B(5, 11) * v[5]; + + return B_Tv; +} + + +template +struct StrainDisplacement +{ + using Real = sofa::Real_t; + static constexpr auto nbLines = symmetric_tensor::NumberOfIndependentElements; + static constexpr auto nbColumns = ElementType::NumberOfNodes * DataTypes::spatial_dimensions; + + constexpr Real& operator()(sofa::Size i, sofa::Size j) noexcept + { + return B(i, j); + } + + constexpr const Real& operator()(sofa::Size i, sofa::Size j) const noexcept + { + return B(i, j); + } + + constexpr sofa::type::Vec operator*(const sofa::type::Vec& v) const + { + return strainDisplacementVectorProduct(B, v); + } + + template + constexpr sofa::type::Mat multTranspose(const sofa::type::Mat& v) const noexcept + { + return B.multTranspose(v); + } + + constexpr sofa::type::Vec multTranspose(const sofa::type::Vec& v) const noexcept + { + return strainDisplacementTransposedVectorProduct(B, v); + } + + sofa::type::Mat> B; + + friend std::ostream& operator<<(std::ostream& out, const StrainDisplacement& B) + { + return operator<<(out, B.B); + } +}; + +template +sofa::type::Mat::nbColumns, real> +operator*(const sofa::type::Mat& A, const StrainDisplacement& B) +{ + return A * B.B; +} + +/** + * Creates a strain-displacement matrix (B-matrix) for finite element calculations. + * + * This function constructs a strain-displacement matrix based on the provided gradient of shape + * functions. The matrix is filled according to spatial dimensions and the number of nodes in the + * element. + * + * @param gradientShapeFunctions A matrix containing the gradient of the shape functions. + * + * @return A strain-displacement matrix of type StrainDisplacement, + * constructed for the finite element with the given gradient shape functions. + */ +template +StrainDisplacement makeStrainDisplacement( + const sofa::type::Mat > gradientShapeFunctions) +{ + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + + StrainDisplacement B; + for (sofa::Size ne = 0; ne < NumberOfNodesInElement; ++ne) + { + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + B(i, ne * spatial_dimensions + i) = gradientShapeFunctions[ne][i]; + } + + auto row = spatial_dimensions; + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + for (sofa::Size j = i + 1; j < spatial_dimensions; ++j) + { + B(row, ne * spatial_dimensions + i) = gradientShapeFunctions[ne][j]; + B(row, ne * spatial_dimensions + j) = gradientShapeFunctions[ne][i]; + ++row; + } + } + } + + return B; +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/SymmetricTensor.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/SymmetricTensor.h new file mode 100644 index 00000000000..30f613d7e35 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/SymmetricTensor.h @@ -0,0 +1,32 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +namespace sofa::component::solidmechanics::fem::elastic::symmetric_tensor +{ + +/// The number of independent elements in a symmetric 2nd-order tensor of size (N x N) +template +constexpr sofa::Size NumberOfIndependentElements = N * (N + 1) / 2; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VecView.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VecView.h new file mode 100644 index 00000000000..daa2cf05542 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VecView.h @@ -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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct VecView +{ + typedef sofa::Size Size; + typedef ValueType value_type; + typedef sofa::Size size_type; + typedef std::ptrdiff_t difference_type; + + template + VecView(sofa::type::Vec& vec) + : m_data(vec.data()) + {} + + template + VecView(sofa::type::Vec& vec, sofa::Size i) + : m_data(&vec.elems.data()[i]) + {} + + value_type operator[](sofa::Size i) const + { + return m_data[i]; + } + + value_type& operator[](sofa::Size i) + { + return m_data[i]; + } + + sofa::type::Vec operator-() const + { + sofa::type::Vec res(sofa::type::NOINIT); + for (sofa::Size i = 0; i < N; ++i) + res[i] = -m_data[i]; + return res; + } + + sofa::type::Vec toVec() const + { + sofa::type::Vec res(sofa::type::NOINIT); + for (sofa::Size i = 0; i < N; ++i) + res[i] = m_data[i]; + return res; + } + + VecView& operator=(const sofa::type::Vec& vec) + { + for (sofa::Size i = 0; i < N; ++i) + m_data[i] = vec[i]; + return *this; + } + +private: + ValueType* m_data { nullptr }; +}; + + +template +sofa::type::Vec operator*(const sofa::type::Mat& mat, const VecView& vec) +{ + sofa::type::Vec res(sofa::type::NOINIT); + + for (sofa::Size i = 0; i < L; ++i) + { + for (sofa::Size j = 0; j < C; ++j) + { + res[i] += mat(i, j) * vec[j]; + } + } + + return res; +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VectorTools.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VectorTools.h new file mode 100644 index 00000000000..3a6722febe4 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VectorTools.h @@ -0,0 +1,63 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +namespace detail +{ + +template +auto extractRefNodesVectorFromGlobalVector(const sofa::topology::Element& element, const VecCoord& vector, std::index_sequence) +{ + using Coord = typename VecCoord::value_type; + using CoordRef = std::reference_wrapper; + return std::array{ std::cref(vector[element[I]])... }; +} + +template +auto extractNodesVectorFromGlobalVector(const sofa::topology::Element& element, const VecCoord& vector, std::index_sequence) +{ + using Coord = typename VecCoord::value_type; + return std::array{ vector[element[I]]... }; +} + +} + +template +auto extractRefNodesVectorFromGlobalVector( + const sofa::topology::Element& element, const VecCoord& vector) +{ + return detail::extractRefNodesVectorFromGlobalVector(element, vector, std::make_index_sequence{}); +} + +template +std::array +extractNodesVectorFromGlobalVector(const sofa::topology::Element& element, const VecCoord& vector) +{ + return detail::extractNodesVectorFromGlobalVector(element, vector, std::make_index_sequence{}); +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VoigtNotation.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VoigtNotation.h new file mode 100644 index 00000000000..ec6830be8b6 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/VoigtNotation.h @@ -0,0 +1,117 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * Converts a Voigt index to the corresponding tensor indices (i, j) for a symmetric tensor. + * + * The Voigt index is a 0-indexed integer that uniquely identifies an independent component of a symmetric tensor. + * + * @tparam DataTypes The data type (e.g., sofa::defaulttype::Vec2Types, sofa::defaulttype::Vec3Types) + * @param voigtIndex The Voigt index (0 <= voigtIndex < N, where N is the number of independent components) + * @return A pair of tensor indices (i, j) such that the symmetric tensor component (i, j) corresponds to the Voigt index + * @pre voigtIndex must be less than symmetric_tensor::NumberOfIndependentElements + * @note For 3D: (0,0) -> 0, (1,1) -> 1, (2,2) -> 2, (1,2) -> 3, (0,2) -> 4, (0,1) -> 5 + * For 2D: (0,0) -> 0, (1,1) -> 1, (0,1) -> 2 + */ +template +constexpr auto toTensorIndices(std::size_t voigtIndex) +{ + assert(voigtIndex < symmetric_tensor::NumberOfIndependentElements); + if constexpr (DataTypes::spatial_dimensions == 3) + { + constexpr std::array voigt3d { + std::make_pair(0, 0), + std::make_pair(1, 1), + std::make_pair(2, 2), + std::make_pair(1, 2), + std::make_pair(0, 2), + std::make_pair(0, 1) + }; + assert(voigtIndex < voigt3d.size()); + return voigt3d[voigtIndex]; + } + else if constexpr (DataTypes::spatial_dimensions == 2) + { + constexpr std::array voigt2d { + std::make_pair(0, 0), + std::make_pair(1, 1), + std::make_pair(0, 1) + }; + assert(voigtIndex < voigt2d.size()); + return voigt2d[voigtIndex]; + } + else + { + return std::make_pair(0, 0); + } +} + +/** + * Converts tensor indices (i, j) to the corresponding Voigt index for a symmetric tensor. + * + * This function handles symmetric tensors by mapping (i, j) to an index in the Voigt convention. + * + * @tparam DataTypes The data type (e.g., sofa::defaulttype::Vec2Types, sofa::defaulttype::Vec3Types) + * @param i First tensor index + * @param j Second tensor index + * @return Voigt index (0 <= index < N, where N is the number of independent components) + * @pre i and j must be valid spatial dimensions (0 <= i, j < DataTypes::spatial_dimensions) + * @note For 3D: (0,0)->0, (1,1)->1, (2,2)->2, (1,2)->3, (0,2)->4, (0,1)->5 + * For 2D: (0,0)->0, (1,1)->1, (0,1)->2 + */ +template +constexpr std::size_t tensorToVoigtIndex(std::size_t i, std::size_t j) +{ + assert(i < DataTypes::spatial_dimensions); + assert(j < DataTypes::spatial_dimensions); + if (i == j) + return i; + return symmetric_tensor::NumberOfIndependentElements - i - j; +} + +static_assert(tensorToVoigtIndex(0,0) == 0); +static_assert(tensorToVoigtIndex(0,1) == 5); +static_assert(tensorToVoigtIndex(0,2) == 4); +static_assert(tensorToVoigtIndex(1,0) == 5); +static_assert(tensorToVoigtIndex(1,1) == 1); +static_assert(tensorToVoigtIndex(1,2) == 3); +static_assert(tensorToVoigtIndex(2,0) == 4); +static_assert(tensorToVoigtIndex(2,1) == 3); +static_assert(tensorToVoigtIndex(2,2) == 2); + +static_assert(tensorToVoigtIndex(0,0) == 0); +static_assert(tensorToVoigtIndex(0,1) == 2); +static_assert(tensorToVoigtIndex(1,0) == 2); +static_assert(tensorToVoigtIndex(1,1) == 1); + +static_assert(tensorToVoigtIndex(0,0) == 0); + + +} 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 new file mode 100644 index 00000000000..996ba550949 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/HexahedronRotation.h @@ -0,0 +1,113 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class HexahedronRotation + * @brief Computes the orientation matrix for 8-node hexahedral elements + * + * This rotation method calculates the element's orientation relative to its rest configuration + * by computing an orthonormal basis from two average edges of the hexahedron. The resulting + * rotation matrix is then used to transform the initial rotation matrix. + * + * @tparam DataTypes The data type used throughout the simulation (e.g., sofa::defaulttype::Vec3Types for 3D) + */ +template +struct HexahedronRotation +{ + using RotationMatrix = sofa::type::Mat>; + + /** + * Computes the rotation matrix for the hexahedral element + * + * @param rotationMatrix Output: Current rotation matrix (relative to initial configuration) + * @param initialRotationMatrix The initial rotation matrix (from rest configuration) + * @param nodesPosition Current positions of all nodes + * @param nodesRestPosition Rest positions of all nodes (initial configuration) + */ + template + void computeRotation(RotationMatrix& rotationMatrix, const RotationMatrix& initialRotationMatrix, + const std::array, NumberOfNodesInElement>& nodesPosition, + const std::array, NumberOfNodesInElement>& nodesRestPosition) + { + SOFA_UNUSED(nodesRestPosition); + + RotationMatrix currentRotation(sofa::type::NOINIT); + computeOrientationFromHexahedron(currentRotation, nodesPosition); + + rotationMatrix = currentRotation.multTranspose(initialRotationMatrix); + } + + static constexpr sofa::helper::Item getItem() + { + return {"hexahedron", "Compute the rotation based on two average edges in the hexahedron"}; + } + +private: + + /** + * Computes the orthonormal basis for the hexahedron + * + * @param rotationMatrix Output: Orthonormal basis matrix (x-axis, y-axis, z-axis) + * @param hexahedronNodes Current positions of all 8 nodes + */ + template + void computeOrientationFromHexahedron(RotationMatrix& rotationMatrix, + const std::array, NumberOfNodesInElement>& hexahedronNodes) + { + using Coord = sofa::Coord_t; + + // Compute average edge vectors + Coord xAxis = (hexahedronNodes[1] - hexahedronNodes[0] + + hexahedronNodes[2] - hexahedronNodes[3] + + hexahedronNodes[5] - hexahedronNodes[4] + + hexahedronNodes[6] - hexahedronNodes[7]) * 0.25; + + Coord yAxis = (hexahedronNodes[3] - hexahedronNodes[0] + + hexahedronNodes[2] - hexahedronNodes[1] + + hexahedronNodes[7] - hexahedronNodes[4] + + hexahedronNodes[6] - hexahedronNodes[5]) * 0.25; + + // Normalize the x-axis + xAxis.normalize(); + + // Compute z-axis as cross product of x and y + Coord zAxis = cross(xAxis, yAxis).normalized(); + + // Recompute y-axis to be orthogonal to x and z + yAxis = cross(zAxis, xAxis); + + // Set the orthonormal basis matrix + rotationMatrix[0] = xAxis; + rotationMatrix[1] = yAxis; + rotationMatrix[2] = zAxis; + } +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/IdentityRotation.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/IdentityRotation.h new file mode 100644 index 00000000000..1d462e91326 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/IdentityRotation.h @@ -0,0 +1,41 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +namespace sofa::component::solidmechanics::fem::elastic +{ + +struct IdentityRotation +{ + template + void computeRotation(RotationMatrix& rotationMatrix, const RotationMatrix& initialRotationMatrix, const NotUsed1&, const NotUsed2&) + { + rotationMatrix.identity(); + } + + static constexpr sofa::helper::Item getItem() + { + return {"identity", "Identity rotation. Equivalent to the linear small strain FEM."}; + } +}; + +} 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 new file mode 100644 index 00000000000..dd07960e9d1 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h @@ -0,0 +1,78 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct PolarDecomposition +{ + using RotationMatrix = sofa::type::Mat>; + + template + sofa::Coord_t computeCentroid(const std::array, NumberOfNodesInElement>& nodes) + { + sofa::Coord_t centroid; + for (const auto node : nodes) + { + centroid += node; + } + centroid /= static_cast>(NumberOfNodesInElement); + return centroid; + } + + template + void computeRotation(RotationMatrix& rotationMatrix, const RotationMatrix& initialRotationMatrix, + const std::array, NumberOfNodesInElement>& nodesPosition, + const std::array, NumberOfNodesInElement>& nodesRestPosition) + { + SOFA_UNUSED(initialRotationMatrix); + + const auto t = computeCentroid(nodesPosition); + const auto t0 = computeCentroid(nodesRestPosition); + + sofa::type::Mat> P(sofa::type::NOINIT); + sofa::type::Mat> 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); + + sofa::helper::Decompose>::polarDecomposition(H, rotationMatrix); + } + + static constexpr sofa::helper::Item getItem() + { + return {"polar", "Polar decomposition"}; + } +}; + +} 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 new file mode 100644 index 00000000000..ea708f32316 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/RotationMethodsContainer.h @@ -0,0 +1,145 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class RotationMethodsContainer + * @brief A container for rotation computation methods in corotational formulations + * + * This class provides a flexible mechanism to select and execute different rotation computation + * strategies for elements in corotational finite element simulations. The rotation methods + * calculate the element's rotation relative to its initial configuration, which is essential + * for corotational force fields. + * + * @tparam DataTypes The data type used throughout the simulation (e.g., sofa::defaulttype::Vec3Types for 3D) + * @tparam ElementType The element type (e.g., sofa::geometry::Triangle, sofa::geometry::Tetrahedron) + * @tparam Methods... A pack of rotation method classes implementing the `computeRotation` interface + * + * Key Features: + * - Supports multiple rotation computation strategies via std::variant + * - Provides GUI-driven method selection through a Data-driven mechanism + * + * Technical Details: + * - The active method is stored in `m_rotationComputer` (std::variant) + * - Method selection is driven by the `d_rotationMethod` data object (index-based) + * - The `selectRotationMethod` function updates the active method based on the selected index + * - The `MAKE_SELECTABLE_ITEMS` macro generates GUI items for method selection + * + * Why Multiple Methods? + * Corotational formulations often require different rotation strategies depending on computational efficiency requirements. + * + * This container enables seamless switching between methods without reconfiguring the entire simulation. + */ +template +struct RotationMethodsContainer +{ +private: + std::variant m_rotationComputer; + +public: + + static constexpr auto NumberOfMethods = std::variant_size_v; + + explicit RotationMethodsContainer(sofa::core::objectmodel::BaseObject* parent) + : d_rotationMethod(parent->initData(&d_rotationMethod, "rotationMethod", ("The method used to compute the element rotations.\n" + RotationMethodsItems::dataDescription()).c_str())) + {} + + using RotationMatrix = sofa::type::Mat>; + static constexpr std::size_t NumberOfNodesInElement = ElementType::NumberOfNodes; + + /** + * Computes the current rotation state relative to the initial configuration. + * + * @param rotationMatrix Output: Current rotation matrix (relative to initial config) + * @param initialRotationMatrix Initial rotation matrix (for reference) + * @param nodesPosition Current positions of element nodes + * @param nodesRestPosition Rest positions of element nodes (initial configuration) + */ + void computeRotation(RotationMatrix& rotationMatrix, const RotationMatrix& initialRotationMatrix, + const std::array, NumberOfNodesInElement>& nodesPosition, + const std::array, NumberOfNodesInElement>& nodesRestPosition) + { + std::visit( + [&](auto& rotationComputer) + { + rotationComputer.computeRotation(rotationMatrix, initialRotationMatrix, nodesPosition, nodesRestPosition); + }, + m_rotationComputer); + } + + struct RotationMethodsItems final : sofa::helper::SelectableItem + { + using sofa::helper::SelectableItem::SelectableItem; + using sofa::helper::SelectableItem::operator=; + using sofa::helper::SelectableItem::operator==; + static constexpr std::array s_items{Methods::getItem()...}; + static_assert(std::is_same_v); + }; + sofa::Data< RotationMethodsItems > d_rotationMethod; + + /** + * Selects the rotation method based on the current index + * + * @note This function is typically called by the GUI to update the active method. + */ + void selectRotationMethod() + { + const std::size_t selectedId = d_rotationMethod.getValue(); + + if (selectedId < NumberOfMethods) + { + doSelectRotationMethod(selectedId); + } + } + +private: + + /** + * Recursively selects the method at the given index + * @tparam Id The index of the method among the available methods + * @param selectedId The index of the method to select + */ + template + void doSelectRotationMethod(const std::size_t selectedId) + { + if (selectedId == Id) + { + using SelectedMethod = std::variant_alternative_t; + m_rotationComputer.template emplace(); + } + else + { + if constexpr (Id > 0) + { + doSelectRotationMethod(selectedId); + } + } + } +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/StablePolarDecomposition.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/StablePolarDecomposition.h new file mode 100644 index 00000000000..7c275648f27 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/StablePolarDecomposition.h @@ -0,0 +1,78 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct StablePolarDecomposition +{ + using RotationMatrix = sofa::type::Mat>; + + template + sofa::Coord_t computeCentroid(const std::array, NumberOfNodesInElement>& nodes) + { + sofa::Coord_t centroid; + for (const auto node : nodes) + { + centroid += node; + } + centroid /= static_cast>(NumberOfNodesInElement); + return centroid; + } + + template + void computeRotation(RotationMatrix& rotationMatrix, const RotationMatrix& initialRotationMatrix, + const std::array, NumberOfNodesInElement>& nodesPosition, + const std::array, NumberOfNodesInElement>& nodesRestPosition) + { + SOFA_UNUSED(initialRotationMatrix); + + const auto t = computeCentroid(nodesPosition); + const auto t0 = computeCentroid(nodesRestPosition); + + sofa::type::Mat> P(sofa::type::NOINIT); + sofa::type::Mat> 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); + + sofa::helper::Decompose>::polarDecomposition_stable(H, rotationMatrix); + } + + static constexpr sofa::helper::Item getItem() + { + return {"stable_polar", "Stable polar decomposition"}; + } +}; + +} 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 new file mode 100644 index 00000000000..dea8ce8055a --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/TriangleRotation.h @@ -0,0 +1,72 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct TriangleRotation +{ + using RotationMatrix = sofa::type::Mat>; + + template + void computeRotation(RotationMatrix& rotationMatrix, const RotationMatrix& initialRotationMatrix, + const std::array, NumberOfNodesInElement>& nodesPosition, + const std::array, NumberOfNodesInElement>& nodesRestPosition) + { + SOFA_UNUSED(nodesRestPosition); + + RotationMatrix currentRotation(sofa::type::NOINIT); + computeRotationFrom3Points(currentRotation, {nodesPosition[0], nodesPosition[1], nodesPosition[2]}); + + rotationMatrix = currentRotation.multTranspose(initialRotationMatrix); + } + + static constexpr sofa::helper::Item getItem() + { + return {"triangle", "Compute the rotation based on the Gram-Schmidt frame alignment"}; + } + +private: + + void computeRotationFrom3Points(RotationMatrix& rotationMatrix, + const std::array, 3>& nodesPosition) + { + using Coord = sofa::Coord_t; + + const Coord xAxis = (nodesPosition[1] - nodesPosition[0]).normalized(); + Coord yAxis = nodesPosition[2] - nodesPosition[0]; + const Coord zAxis = cross( xAxis, yAxis ).normalized(); + yAxis = cross( zAxis, xAxis ); //yAxis is a unit vector because zAxis and xAxis are orthogonal unit vectors + + rotationMatrix[0] = xAxis; + rotationMatrix[1] = yAxis; + rotationMatrix[2] = zAxis; + } +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h new file mode 100644 index 00000000000..2b49c360af6 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h @@ -0,0 +1,75 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +struct trait +{ + using DataVecCoord = sofa::DataVecDeriv_t; + using DataVecDeriv = sofa::DataVecDeriv_t; + using VecCoord = sofa::VecCoord_t; + using VecDeriv = sofa::VecDeriv_t; + using Coord = sofa::Coord_t; + using Deriv = sofa::Deriv_t; + using Real = sofa::Real_t; + + using FiniteElement = sofa::component::solidmechanics::fem::elastic::FiniteElement; + using TopologyElement = typename FiniteElement::TopologyElement; + using ReferenceCoord = typename FiniteElement::ReferenceCoord; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + static constexpr sofa::Size NumberOfDofsInElement = NumberOfNodesInElement * spatial_dimensions; + static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; + static constexpr sofa::Size NbQuadraturePoints = FiniteElement::quadraturePoints().size(); + + /// type of 2nd-order tensor for the elasticity tensor for isotropic materials + using ElasticityTensor = FullySymmetric4Tensor; + + /// the type of B in e = B d, if e is the strain, and d is the displacement + using StrainDisplacement = sofa::component::solidmechanics::fem::elastic::StrainDisplacement; + + /// the concatenation of the displacement of the element nodes in a single vector + using ElementDisplacement = sofa::type::Vec; + + /// tells how to compute the matrix-vector product of the stiffness matrix with a displacement + /// vector. It does not change the result, but it can have an impact on performances. + static constexpr MatrixVectorProductType matrixVectorProductType = + NbQuadraturePoints > 1 + ? MatrixVectorProductType::Dense + : MatrixVectorProductType::Factorization; + + /// the type of the element stiffness matrix + using ElementStiffness = sofa::component::solidmechanics::fem::elastic::FactorizedElementStiffness; + + using ElementForce = sofa::type::Vec>; +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp index db73328f65e..343b296dbba 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp @@ -27,6 +27,8 @@ namespace sofa::component::solidmechanics::fem::elastic { extern void registerBeamFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerElementCorotationalFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerElementLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* factory); extern void registerFastTetrahedralCorotationalForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedralFEMForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedralFEMForceFieldAndMass(sofa::core::ObjectFactory* factory); @@ -65,6 +67,8 @@ const char* getModuleVersion() void registerObjects(sofa::core::ObjectFactory* factory) { registerBeamFEMForceField(factory); + registerElementCorotationalFEMForceField(factory); + registerElementLinearSmallStrainFEMForceField(factory); registerFastTetrahedralCorotationalForceField(factory); registerHexahedralFEMForceField(factory); registerHexahedralFEMForceFieldAndMass(factory);