diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/HyperElastic/CMakeLists.txt index 12a1833d4ad..8cec2c82e59 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/CMakeLists.txt @@ -6,6 +6,12 @@ set(SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR "src/sofa/component/so set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/config.h.in ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/init.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/HyperelasticityFEMForceField.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/HyperelasticityFEMForceField.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/HyperelasticMaterial.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/HyperelasticMaterial.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/PK2HyperelasticMaterial.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/PK2HyperelasticMaterial.inl ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/StandardTetrahedralFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/StandardTetrahedralFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/TetrahedronHyperelasticityFEMDrawing.h @@ -22,19 +28,45 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/material/StableNeoHookean.h ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/material/STVenantKirchhoff.h ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/material/VerondaWestman.h + + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/impl/MajorSymmetric4Tensor.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/impl/Strain.h + + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/LinearMechanicalParametersComponent.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/LinearMechanicalParametersComponent.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/MooneyRivlinMaterial.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/MooneyRivlinMaterial.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/NeoHookeanMaterial.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/NeoHookeanMaterial.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/OgdenMaterial.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/OgdenMaterial.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/StVenantKirchhoffMaterial.h + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/StVenantKirchhoffMaterial.inl ) set(SOURCE_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/init.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/PK2HyperelasticMaterial.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/HyperelasticityFEMForceField.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/HyperelasticMaterial.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/StandardTetrahedralFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/TetrahedronHyperelasticityFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/material/PlasticMaterial.cpp + + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/LinearMechanicalParametersComponent.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/MooneyRivlinMaterial.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/NeoHookeanMaterial.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/OgdenMaterial.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMHYPERELASTIC_SOURCE_DIR}/hyperelasticmaterials/StVenantKirchhoffMaterial.cpp ) -sofa_find_package(Sofa.Simulation.Core REQUIRED) +sofa_find_package(Sofa.Simulation.Core Sofa.Component.SolidMechanics.FEM.Elastic REQUIRED) add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES}) -target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core) +target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core Sofa.Component.SolidMechanics.FEM.Elastic) sofa_create_package_with_targets( PACKAGE_NAME ${PROJECT_NAME} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/Sofa.Component.SolidMechanics.FEM.HyperElasticConfig.cmake.in b/Sofa/Component/SolidMechanics/FEM/HyperElastic/Sofa.Component.SolidMechanics.FEM.HyperElasticConfig.cmake.in index 5e94811a2f4..e174907f90b 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/Sofa.Component.SolidMechanics.FEM.HyperElasticConfig.cmake.in +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/Sofa.Component.SolidMechanics.FEM.HyperElasticConfig.cmake.in @@ -5,6 +5,7 @@ find_package(Sofa.Config QUIET REQUIRED) sofa_find_package(Sofa.Simulation.Core QUIET REQUIRED) +sofa_find_package(Sofa.Component.SolidMechanics.FEM.Elastic QUITE REQUIRED) sofa_find_package(Sofa.Component.Topology.Container.Grid QUIET REQUIRED) if(NOT TARGET @PROJECT_NAME@) diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.cpp new file mode 100644 index 00000000000..a9cfa931fe4 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.cpp @@ -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 * +******************************************************************************/ +#define ELASTICITY_COMPONENT_HYPERELASTIC_MATERIAL_CPP + +#include +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticMaterial; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.h new file mode 100644 index 00000000000..ba253a4bc77 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.h @@ -0,0 +1,82 @@ +/****************************************************************************** +* 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 + +#if !defined(ELASTICITY_COMPONENT_HYPERELASTIC_MATERIAL_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template +class HyperelasticMaterial : public virtual sofa::core::objectmodel::BaseObject +{ +public: + SOFA_ABSTRACT_CLASS(HyperelasticMaterial, sofa::core::objectmodel::BaseObject); + using DataTypes = TDataTypes; + +protected: + using Real = sofa::Real_t; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + using DeformationGradient = sofa::type::Mat; + using RightCauchyGreenTensor = sofa::type::Mat; + using StressTensor = sofa::type::Mat; + using ElasticityTensor = sofa::type::FullySymmetric4Tensor; + using TangentModulus = MajorSymmetric4Tensor; + +public: + void init() override; + + /** + * Computes the First Piola-Kirchhoff stress tensor for a given deformation gradient. + * + * It corresponds to the derivative of the strain energy density function w.r.t. deformation + * gradient. + */ + virtual StressTensor firstPiolaKirchhoffStress(Strain& strain) = 0; + + /** + * Compute the jacobian of the first Piola-Kirchhoff stress tensor with respect to the + * deformation gradient. + * + * It is called the material tangent modulus. + */ + virtual TangentModulus materialTangentModulus(Strain& strain) = 0; + +}; + +#if !defined(ELASTICITY_COMPONENT_HYPERELASTIC_MATERIAL_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticMaterial; +#endif + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.inl new file mode 100644 index 00000000000..225513ef0c0 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticMaterial.inl @@ -0,0 +1,39 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +void HyperelasticMaterial::init() +{ + sofa::core::objectmodel::BaseObject::init(); + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + } +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.cpp new file mode 100644 index 00000000000..046a8988450 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.cpp @@ -0,0 +1,57 @@ +/****************************************************************************** +* 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_HYPERLASTICITY_FEM_FORCE_FIELD_CPP + +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +void registerHyperelasticityFEMForceField(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Hyperelasticity") + .add< HyperelasticityFEMForceField >() + .add< HyperelasticityFEMForceField >() + .add< HyperelasticityFEMForceField >() + .add< HyperelasticityFEMForceField >() + .add< HyperelasticityFEMForceField >() + .add< HyperelasticityFEMForceField >() + .add< HyperelasticityFEMForceField >() + .add< HyperelasticityFEMForceField >() + .add< HyperelasticityFEMForceField >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.h new file mode 100644 index 00000000000..6056675d182 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.h @@ -0,0 +1,164 @@ +/****************************************************************************** +* 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_ELEMENT_HYPERLASTICITY_FEM_FORCE_FIELD_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template +class HyperelasticityFEMForceField : + public sofa::component::solidmechanics::fem::elastic::FEMForceField +{ +public: + SOFA_CLASS( + SOFA_TEMPLATE2(HyperelasticityFEMForceField, TDataTypes, TElementType), + SOFA_TEMPLATE2(sofa::component::solidmechanics::fem::elastic::FEMForceField, TDataTypes, TElementType)); + + using DataTypes = TDataTypes; + +private: + 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::fem::FiniteElement; + + static constexpr sofa::Size spatial_dimensions = TDataTypes::spatial_dimensions; + static constexpr sofa::Size NumberOfNodesInElement = TElementType::NumberOfNodes; + static constexpr sofa::Size NumberOfDofsInElement = NumberOfNodesInElement * spatial_dimensions; + static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; + static constexpr sofa::Size NumberOfQuadraturePoints = FiniteElement::quadraturePoints().size(); + + using DeformationGradient = sofa::type::Mat; + + /// the type of the element stiffness matrix + using ElementStiffness = sofa::type::Mat< + TElementType::NumberOfNodes * DataTypes::spatial_dimensions, + TElementType::NumberOfNodes * DataTypes::spatial_dimensions, + sofa::Real_t + >; + +public: + void init() override; + + void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; + + using Inherit1::getPotentialEnergy; + SReal getPotentialEnergy(const sofa::core::MechanicalParams*, const DataVecCoord& x) const override; + + using Inherit1::addKToMatrix; + // almost deprecated, but here for compatibility with unit tests + void addKToMatrix(sofa::linearalgebra::BaseMatrix* matrix, SReal kFact, unsigned& offset) override; + + sofa::SingleLink, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_material; + +protected: + using trait = sofa::component::solidmechanics::fem::elastic::trait; + using ElementGradient = typename trait::ElementGradient; + + void validateMaterial(); + + bool m_isHessianValid; + + void computeHessian(const VecCoord& coordinates); + + /** + * List of precomputed element stiffness matrices + */ + sofa::type::vector m_elementStiffness; + + const VecCoord* m_coordinates{ nullptr }; + + DeformationGradient computeDeformationGradient( + const sofa::type::Mat& J_q, + const sofa::type::Mat& J_Q_inv); + + struct PrecomputedData + { + // jacobian of the mapping from the reference space to the rest physical space, evaluated at the + // quadrature point + sofa::type::Mat jacobian { sofa::type::NOINIT }; + + // inverse of the jacobian of the mapping from the reference space to the rest physical space, + // evaluated at the quadrature point + sofa::type::Mat jacobianInv { sofa::type::NOINIT }; + + Real detJacobian {}; + + // gradient of the shape functions in the physical element evaluated at the quadrature point + sofa::type::Mat dN_dQ { sofa::type::NOINIT }; + }; + + /** + * Data can be precomputed from the rest configuration and used later in computations for the + * current configuration. A piece of precomputed data is stored for each quadrature point in + * each element. + */ + sofa::type::vector> m_precomputedData; + + void precomputeData(); + + 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 beforeElementForceDeriv(const sofa::core::MechanicalParams* mparams) override; + + void computeElementsForcesDeriv( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& df, + const sofa::VecDeriv_t& dx) override; +}; + +#if !defined(ELASTICITY_COMPONENT_HYPERLASTICITY_FEM_FORCE_FIELD_CPP) +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API HyperelasticityFEMForceField; +#endif + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.inl new file mode 100644 index 00000000000..58a4fc1e627 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.inl @@ -0,0 +1,406 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +void HyperelasticityFEMForceField::init() +{ + sofa::component::solidmechanics::fem::elastic::FEMForceField::init(); + + if (!this->isComponentStateInvalid()) + { + this->validateMaterial(); + } + + if (!this->isComponentStateInvalid()) + { + precomputeData(); + } + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + } +} + +template +void HyperelasticityFEMForceField::buildStiffnessMatrix( + sofa::core::behavior::StiffnessMatrix* matrix) +{ + if (this->isComponentStateInvalid()) + return; + + auto dfdx = matrix->getForceDerivativeIn(this->mstate) + .withRespectToPositionsIn(this->mstate); + + if (!m_isHessianValid) + { + computeHessian(*m_coordinates); + } + + sofa::type::Mat localMatrix(sofa::type::NOINIT); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + auto elementStiffnessIt = m_elementStiffness.begin(); + for (const auto& element : elements) + { + const auto& stiffnessMatrix = *elementStiffnessIt++; + + for (sofa::Index n1 = 0; n1 < NumberOfNodesInElement; ++n1) + { + for (sofa::Index n2 = 0; n2 < NumberOfNodesInElement; ++n2) + { + stiffnessMatrix.getsub(spatial_dimensions * n1, spatial_dimensions * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2 + dfdx(element[n1] * spatial_dimensions, element[n2] * spatial_dimensions) += -localMatrix; + } + } + } +} + +template +SReal HyperelasticityFEMForceField::getPotentialEnergy( + const sofa::core::MechanicalParams*, const DataVecCoord& x) const +{ + return 0; +} + +template +void HyperelasticityFEMForceField::addKToMatrix( + sofa::linearalgebra::BaseMatrix* matrix, SReal kFact, unsigned& offset) +{ + if (this->isComponentStateInvalid()) + return; + + if (!m_isHessianValid) + { + computeHessian(*m_coordinates); + } + + using LocalMatType = sofa::type::Mat; + LocalMatType localMatrix{sofa::type::NOINIT}; + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + auto elementStiffnessIt = m_elementStiffness.begin(); + for (const auto& element : elements) + { + const auto& stiffnessMatrix = *elementStiffnessIt++; + + for (sofa::Index n1 = 0; n1 < NumberOfNodesInElement; ++n1) + { + for (sofa::Index n2 = 0; n2 < NumberOfNodesInElement; ++n2) + { + stiffnessMatrix.getsub(spatial_dimensions * n1, 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] * spatial_dimensions, + offset + element[n2] * spatial_dimensions, value); + } + } + } +} + +template +void HyperelasticityFEMForceField::validateMaterial() +{ + if (l_material.empty()) + { + msg_info() << "Link to a valid material should be set to ensure right behavior. First " + "material found in current context will be used."; + l_material.set(this->getContext()->template get>()); + } + + if (l_material == nullptr) + { + msg_error() << "No material component found at path: '" << this->l_material.getLinkedPath() + << "', nor in current context: " << this->getContext()->name + << ". Object must have a material. " + << "The list of available material components is: " + << sofa::core::ObjectFactory::getInstance() + ->listClassesDerivedFrom>(); + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } +} + +template +void HyperelasticityFEMForceField::computeHessian(const VecCoord& coordinates) +{ + if (this->l_topology == nullptr) return; + if (l_material == nullptr) return; + + SCOPED_TIMER("ComputeHessian"); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + + if (m_precomputedData.size() != elements.size()) + { + precomputeData(); + } + + m_elementStiffness.clear(); + m_elementStiffness.resize(elements.size()); + + auto elementStiffnessIt = m_elementStiffness.begin(); + + std::size_t elementIndex = 0; + for (const auto& element : elements) + { + SCOPED_TIMER_TR("Element"); + const std::array elementNodesCoordinates = sofa::component::solidmechanics::fem::elastic::extractNodesVectorFromGlobalVector(element, this->mstate->readPositions().ref()); + + ElementStiffness& K = *elementStiffnessIt++; + K.clear(); + + static constexpr auto quadraturePoints = FiniteElement::quadraturePoints(); + static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); + + for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) + { + const auto& weight = quadraturePoints[q].second; + const PrecomputedData& precomputedData = m_precomputedData[elementIndex][q]; + + // gradient of shape functions in the reference element evaluated at the quadrature point + const sofa::type::Mat& dN_dq_ref = gradients[q]; + + // jacobian of the mapping from the reference space to the physical space, evaluated at the + // quadrature point + const auto J_q = sofa::fem::FiniteElementHelper::jacobianFromReferenceToPhysical(elementNodesCoordinates, dN_dq_ref); + + const auto detJ_Q = precomputedData.detJacobian; + + const sofa::type::Mat& J_Q_inv = precomputedData.jacobianInv; + + // gradient of the shape functions in the physical element evaluated at the quadrature point + const sofa::type::Mat& dN_dQ = precomputedData.dN_dQ; + + const DeformationGradient F = computeDeformationGradient(J_q, J_Q_inv); + + Strain strain(deformationGradient, F); + + // derivative of first Piola-Kirchhoff stress tensor with respect to deformation gradient + const auto dPdF = [&]() + { + SCOPED_TIMER_VARNAME_TR(dPdFTimer, "dPdF"); + return l_material->materialTangentModulus(strain); + }(); + + const auto factor = detJ_Q * weight; + + SCOPED_TIMER_VARNAME_TR(tensorTimer, "Tensor"); + for (sofa::Size element_i = 0; element_i < NumberOfNodesInElement; ++element_i) + { + for (sofa::Size element_j = 0; element_j < NumberOfNodesInElement; ++element_j) + { + for (sofa::Size dimension_i = 0; dimension_i < spatial_dimensions; ++dimension_i) + { + for (sofa::Size dimension_j = 0; dimension_j < spatial_dimensions; ++dimension_j) + { + auto& k = K(element_i * spatial_dimensions + dimension_i, element_j * spatial_dimensions + dimension_j); + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + for (sofa::Size j = 0; j < spatial_dimensions; ++j) + { + k += factor * dN_dQ[element_i][i] * dPdF(dimension_i, i, dimension_j, j) * dN_dQ[element_j][j]; + } + } + } + } + } + } + } + ++elementIndex; + } + + m_isHessianValid = true; +} + +template +auto HyperelasticityFEMForceField::computeDeformationGradient( + const sofa::type::Mat& J_q, + const sofa::type::Mat& J_Q_inv) -> DeformationGradient +{ + return J_q * J_Q_inv; +} + +template +void HyperelasticityFEMForceField::precomputeData() +{ + if (this->l_topology == nullptr) return; + + auto restPositionAccessor = this->mstate->readRestPositions(); + const auto& restPosition = restPositionAccessor.ref(); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + m_precomputedData.resize(elements.size()); + + static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); + + for (std::size_t i = 0; i < elements.size(); ++i) + { + const auto& element = elements[i]; + const std::array elementNodesRestCoordinates = sofa::component::solidmechanics::fem::elastic::extractNodesVectorFromGlobalVector(element, restPosition); + + for (std::size_t j = 0; j < NumberOfQuadraturePoints; ++j) + { + // gradient of shape functions in the reference element evaluated at the quadrature point + const sofa::type::Mat& dN_dq_ref = gradients[j]; + + PrecomputedData& data = m_precomputedData[i][j]; + data.jacobian = sofa::fem::FiniteElementHelper::jacobianFromReferenceToPhysical(elementNodesRestCoordinates, dN_dq_ref); + data.jacobianInv = sofa::type::inverse(data.jacobian); + data.detJacobian = sofa::type::absGeneralizedDeterminant(data.jacobian); + + for (sofa::Size n = 0; n < NumberOfNodesInElement; ++n) + { + data.dN_dQ[n] = data.jacobianInv.multTranspose(dN_dq_ref[n]); + } + } + } +} +template +void HyperelasticityFEMForceField::beforeElementForce( + const sofa::core::MechanicalParams* mparams, sofa::type::vector& f, + const sofa::VecCoord_t& x) +{ + //store coordinates to use it later when computing the Hessian + m_coordinates = &x; + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + if (m_precomputedData.size() != elements.size()) + { + precomputeData(); + } + + m_isHessianValid = false; + + //reset force vector + for (auto& elF : f) + { + elF.clear(); + } +} + +template +void HyperelasticityFEMForceField::computeElementsForces( + const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, + sofa::type::vector& f, const sofa::VecCoord_t& x) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + + static constexpr auto quadraturePoints = FiniteElement::quadraturePoints(); + static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); + + for (std::size_t elementId = range.start; elementId < range.end; ++elementId) + { + const auto element = elements[elementId]; + const std::array elementNodesCoordinates = sofa::component::solidmechanics::fem::elastic::extractNodesVectorFromGlobalVector(element, x); + + //access to the force vector in the element. This is the value to compute in this iteration + //force assembly at the DoF level is done in a later function. + auto& elementForce = f[elementId]; + + for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) + { + const auto& weight = quadraturePoints[q].second; + const PrecomputedData& precomputedData = m_precomputedData[elementId][q]; + + // gradient of shape functions in the reference element evaluated at the quadrature point + const sofa::type::Mat& dN_dq_ref = gradients[q]; + + // jacobian of the mapping from the reference space to the physical space, evaluated at the + // quadrature point + const auto J_q = sofa::fem::FiniteElementHelper::jacobianFromReferenceToPhysical(elementNodesCoordinates, dN_dq_ref); + + const auto detJ_Q = precomputedData.detJacobian; + + const sofa::type::Mat& J_Q_inv = precomputedData.jacobianInv; + + // gradient of the shape functions in the physical element evaluated at the quadrature point + const sofa::type::Mat& dN_dQ = precomputedData.dN_dQ; + + const DeformationGradient F = computeDeformationGradient(J_q, J_Q_inv); + + Strain strain(deformationGradient, F); + const auto P = l_material->firstPiolaKirchhoffStress(strain); + + for (sofa::Size i = 0; i < trait::NumberOfNodesInElement; ++i) + { + const auto f_q = (detJ_Q * weight) * P * dN_dQ[i]; + for (sofa::Size j = 0; j < trait::spatial_dimensions; ++j) + { + elementForce[i * trait::spatial_dimensions + j] += f_q[j]; + } + } + } + } +} + +template +void HyperelasticityFEMForceField::beforeElementForceDeriv( + const sofa::core::MechanicalParams* mparams) +{ + if (!m_isHessianValid) + { + computeHessian(*m_coordinates); + } +} + +template +void HyperelasticityFEMForceField::computeElementsForcesDeriv( + const sofa::simulation::Range& range, const sofa::core::MechanicalParams* mparams, + sofa::type::vector& df, const sofa::VecDeriv_t& dx) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + + for (std::size_t elementId = range.start; elementId < range.end; ++elementId) + { + const auto& element = elements[elementId]; + const auto& stiffnessMatrix = m_elementStiffness[elementId]; + + const std::array, trait::NumberOfNodesInElement> elementNodesDx = + sofa::component::solidmechanics::fem::elastic::extractNodesVectorFromGlobalVector(element, dx); + + sofa::type::Vec> element_dx(sofa::type::NOINIT); + for (sofa::Size nodeId = 0; nodeId < trait::NumberOfNodesInElement; ++nodeId) + { + const auto& nodeDx = elementNodesDx[nodeId]; + for (sofa::Size dim = 0; dim < trait::spatial_dimensions; ++dim) + { + element_dx[nodeId * trait::spatial_dimensions + dim] = nodeDx[dim]; + } + } + + df[elementId] = stiffnessMatrix * element_dx; + } +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.cpp new file mode 100644 index 00000000000..3b6b715d009 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.cpp @@ -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 * +******************************************************************************/ +#define ELASTICITY_COMPONENT_PK2HYPERELASTIC_MATERIAL_CPP +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API PK2HyperelasticMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API PK2HyperelasticMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API PK2HyperelasticMaterial; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.h new file mode 100644 index 00000000000..f404918752e --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.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::hyperelastic +{ + +/** + * A hyperelastic material defined by its second Piola-Kirchhoff stress tensor and its Lagrangian + * elasticity tensor. + */ +template +class PK2HyperelasticMaterial : public HyperelasticMaterial +{ +public: + SOFA_CLASS(PK2HyperelasticMaterial, HyperelasticMaterial); + using DataTypes = TDataTypes; + +protected: + using DeformationGradient = typename HyperelasticMaterial::DeformationGradient; + using RightCauchyGreenTensor = typename HyperelasticMaterial::RightCauchyGreenTensor; + using StressTensor = typename HyperelasticMaterial::StressTensor; + using ElasticityTensor = typename HyperelasticMaterial::ElasticityTensor; + using TangentModulus = typename HyperelasticMaterial::TangentModulus; + using HyperelasticMaterial::spatial_dimensions; + +public: + StressTensor firstPiolaKirchhoffStress(Strain& strain) final; + TangentModulus materialTangentModulus(Strain& strain) final; + +protected: + virtual StressTensor secondPiolaKirchhoffStress(Strain& strain) = 0; + virtual ElasticityTensor elasticityTensor(Strain& strain) = 0; +}; + +#if !defined(ELASTICITY_COMPONENT_HYPERELASTIC_MATERIAL_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API PK2HyperelasticMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API PK2HyperelasticMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API PK2HyperelasticMaterial; +#endif + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.inl new file mode 100644 index 00000000000..eaa7790cd3a --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/PK2HyperelasticMaterial.inl @@ -0,0 +1,67 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +auto PK2HyperelasticMaterial::firstPiolaKirchhoffStress(Strain& strain) -> StressTensor +{ + const auto& F = strain.deformationGradient(); + const auto S = secondPiolaKirchhoffStress(strain); + return F * S; +} + +template +auto PK2HyperelasticMaterial::materialTangentModulus(Strain& strain) -> TangentModulus +{ + using Real = sofa::Real_t; + SCOPED_TIMER_TR("tangentModulus"); + + const auto& F = strain.deformationGradient(); + const auto C = elasticityTensor(strain); + const auto S = secondPiolaKirchhoffStress(strain); + + const auto A = TangentModulus([&F, &C, &S](sofa::Index i, sofa::Index j, sofa::Index k, sofa::Index l) + { + auto A_ijkl = sofa::component::solidmechanics::fem::elastic::kroneckerDelta(i,k) * S(l, j); + for (std::size_t q = 0; q < spatial_dimensions; ++q) + { + for (std::size_t r = 0; r < spatial_dimensions; ++r) + { + A_ijkl += F(i, q) * C(q, j, l, r) * F(k, r); + } + } + return A_ijkl; + }); + + return A; +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.cpp new file mode 100644 index 00000000000..be2ec4c7c6e --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.cpp @@ -0,0 +1,44 @@ +/****************************************************************************** +* 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_MATERIAL_INCOMPRESSIBLEMOONEYRIVLIN_CPP + +#include +#include + +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +void registerIncompressibleMooneyRivlinMaterial(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Incompressible Mooney-Rivlin material") + .add< IncompressibleMooneyRivlinMaterial >() + .add< IncompressibleMooneyRivlinMaterial >() + .add< IncompressibleMooneyRivlinMaterial >(true)); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API IncompressibleMooneyRivlinMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API IncompressibleMooneyRivlinMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API IncompressibleMooneyRivlinMaterial; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.h new file mode 100644 index 00000000000..097051d85ba --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.h @@ -0,0 +1,69 @@ +/****************************************************************************** +* 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 + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_INCOMPRESSIBLEMOONEYRIVLIN_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template +class IncompressibleMooneyRivlinMaterial : + public PK2HyperelasticMaterial +{ +public: + SOFA_CLASS(IncompressibleMooneyRivlinMaterial, PK2HyperelasticMaterial); + +private: + using Real = sofa::Real_t; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + using DeformationGradient = typename PK2HyperelasticMaterial::DeformationGradient; + using RightCauchyGreenTensor = typename PK2HyperelasticMaterial::RightCauchyGreenTensor; + using StressTensor = typename PK2HyperelasticMaterial::StressTensor; + using ElasticityTensor = typename PK2HyperelasticMaterial::ElasticityTensor; + +public: + sofa::Data m_mu10; + sofa::Data m_mu01; + + StressTensor secondPiolaKirchhoffStress(Strain& strain) override; + + ElasticityTensor elasticityTensor(Strain& strain) override; + +protected: + IncompressibleMooneyRivlinMaterial(); +}; + + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_INCOMPRESSIBLEMOONEYRIVLIN_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API IncompressibleMooneyRivlinMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API IncompressibleMooneyRivlinMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API IncompressibleMooneyRivlinMaterial; +#endif +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.inl new file mode 100644 index 00000000000..59a884b3cdb --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/IncompressibleMooneyRivlinMaterial.inl @@ -0,0 +1,68 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +IncompressibleMooneyRivlinMaterial::IncompressibleMooneyRivlinMaterial() + : m_mu10(initData(&m_mu10, static_cast(1e3), "mu10", + "Material constant associated to the first invariant")) + , m_mu01(initData(&m_mu01, static_cast(1e3), "mu01", + "Material constant associated to the second invariant")) +{} + +template +auto IncompressibleMooneyRivlinMaterial::secondPiolaKirchhoffStress(Strain& strain) -> StressTensor +{ + static constexpr auto& I = Strain::identity; + const auto& C = strain.getRightCauchyGreenTensor(); + + const auto invariant1 = strain.getInvariant1(); + + const auto mu10 = m_mu10.getValue(); + const auto mu01 = m_mu01.getValue(); + + return static_cast(2) * mu10 * I + static_cast(2) * mu01 * (invariant1 * I - C); +} + +template +auto IncompressibleMooneyRivlinMaterial::elasticityTensor(Strain& strain) -> ElasticityTensor +{ + auto delta = [](auto i, auto j){ return sofa::component::solidmechanics::fem::elastic::kroneckerDelta(i, j); }; + + return ElasticityTensor( + [mu01 = m_mu01.getValue(), &delta](sofa::Index i, sofa::Index j, sofa::Index k, sofa::Index l) + { + return 2 * mu01 * ( + 2 * delta(i, j) * delta(k, l) + - delta(i, k) * delta(j, l) + - delta(i, l) * delta(j, k) + ); + }); +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.cpp new file mode 100644 index 00000000000..83d947706c1 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.cpp @@ -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 * +******************************************************************************/ +#define ELASTICITY_COMPONENT_LINEAR_MECHANICAL_PARAMETERS_CPP +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API LinearMechanicalParametersComponent; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API LinearMechanicalParametersComponent; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API LinearMechanicalParametersComponent; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.h new file mode 100644 index 00000000000..68ee4eb8ee0 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.h @@ -0,0 +1,59 @@ +/****************************************************************************** +* 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 + +#if !defined(ELASTICITY_COMPONENT_LINEAR_MECHANICAL_PARAMETERS_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template +class LinearMechanicalParametersComponent : public virtual sofa::core::objectmodel::BaseObject +{ + using Real = sofa::Real_t; + +public: + SOFA_CLASS(LinearMechanicalParametersComponent, sofa::core::objectmodel::BaseObject); + + sofa::Data d_poissonRatio; + sofa::Data d_youngModulus; + +protected: + LinearMechanicalParametersComponent(); + + void setLameCoefficients(); + + // Lamé's coefficients + Real m_lambda, m_mu; +}; + +#if !defined(ELASTICITY_COMPONENT_LINEAR_MECHANICAL_PARAMETERS_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API LinearMechanicalParametersComponent; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API LinearMechanicalParametersComponent; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API LinearMechanicalParametersComponent; +#endif +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.inl new file mode 100644 index 00000000000..e28697f9ef1 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/LinearMechanicalParametersComponent.inl @@ -0,0 +1,60 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +LinearMechanicalParametersComponent::LinearMechanicalParametersComponent() +: d_poissonRatio(initData(&d_poissonRatio, static_cast(0.45), "poissonRatio", + "Poisson's ratio: represents the material's ability to undergo deformation in directions orthogonal to the applied stress")) +, d_youngModulus(initData(&d_youngModulus, static_cast(1e6), "youngModulus", + "Young's modulus: describes the material's stiffness")) +{ + this->addUpdateCallback("toLameCoefficients", {&this->d_youngModulus, &this->d_poissonRatio}, + [this](const sofa::core::DataTracker& ) + { + setLameCoefficients(); + return this->getComponentState(); + }, {}); + setLameCoefficients(); +} + +template +void LinearMechanicalParametersComponent::setLameCoefficients() +{ + sofa::component::solidmechanics::fem::elastic::LameLambda lambdaStrong { 0 }; + sofa::component::solidmechanics::fem::elastic::LameMu muStrong { 0 }; + + sofa::component::solidmechanics::fem::elastic::toLameParameters( + sofa::component::solidmechanics::fem::elastic::YoungModulus(this->d_youngModulus.getValue()), + sofa::component::solidmechanics::fem::elastic::PoissonRatio(this->d_poissonRatio.getValue()), + lambdaStrong, muStrong); + + m_lambda = lambdaStrong.get(); + m_mu = muStrong.get(); +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.cpp new file mode 100644 index 00000000000..bfe6206634b --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.cpp @@ -0,0 +1,44 @@ +/****************************************************************************** +* 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_MATERIAL_MOONEYRIVLIN_CPP + +#include +#include + +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +void registerMooneyRivlinMaterial(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Mooney-Rivlin material") + .add< MooneyRivlinMaterial >() + .add< MooneyRivlinMaterial >() + .add< MooneyRivlinMaterial >(true)); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API MooneyRivlinMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API MooneyRivlinMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API MooneyRivlinMaterial; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.h new file mode 100644 index 00000000000..655d8dd286e --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.h @@ -0,0 +1,70 @@ +/****************************************************************************** +* 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 + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_MOONEYRIVLIN_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template +class MooneyRivlinMaterial : + public PK2HyperelasticMaterial +{ +public: + SOFA_CLASS(MooneyRivlinMaterial, PK2HyperelasticMaterial); + +private: + using Real = sofa::Real_t; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + using DeformationGradient = typename PK2HyperelasticMaterial::DeformationGradient; + using RightCauchyGreenTensor = typename PK2HyperelasticMaterial::RightCauchyGreenTensor; + using StressTensor = typename PK2HyperelasticMaterial::StressTensor; + using ElasticityTensor = typename PK2HyperelasticMaterial::ElasticityTensor; + +public: + sofa::Data m_mu10; + sofa::Data m_mu01; + sofa::Data m_bulkModulus; + + StressTensor secondPiolaKirchhoffStress(Strain& strain) override; + + ElasticityTensor elasticityTensor(Strain& strain) override; + +protected: + MooneyRivlinMaterial(); +}; + + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_MOONEYRIVLIN_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API MooneyRivlinMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API MooneyRivlinMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API MooneyRivlinMaterial; +#endif +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.inl new file mode 100644 index 00000000000..707eae01137 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/MooneyRivlinMaterial.inl @@ -0,0 +1,132 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +MooneyRivlinMaterial::MooneyRivlinMaterial() + : m_mu10(initData(&m_mu10, static_cast(1e3), "mu10", + "Material constant associated to the first invariant")) + , m_mu01(initData(&m_mu01, static_cast(1e3), "mu01", + "Material constant associated to the second invariant")) + , m_bulkModulus(initData(&m_bulkModulus, static_cast(1e2), "bulkModulus", "Bulk modulus")) +{} + +template +auto MooneyRivlinMaterial::secondPiolaKirchhoffStress(Strain& strain) -> StressTensor +{ + const auto& C = strain.getRightCauchyGreenTensor(); + const auto C_1 = sofa::type::inverse(C); + + const auto J = strain.getDeterminantDeformationGradient(); + assert(J > 0); + + const auto S_isochoric = [this, J, &strain, &C_1, &C]() + { + static constexpr auto& I = Strain::identity; + static constexpr Real dim_1 = static_cast(1) / static_cast(spatial_dimensions); + + const auto mu10 = m_mu10.getValue(); + const auto mu01 = m_mu01.getValue(); + + const auto invariant1 = strain.getInvariant1(); + const auto invariant2 = strain.getInvariant2(); + + const auto S_mu_10 = pow(J, -static_cast(2) * dim_1) * (I - dim_1 * invariant1 * C_1); + const auto S_mu_01 = pow(J, -static_cast(4) * dim_1) * (invariant1 * I - C - static_cast(2) * dim_1 * invariant2 * C_1); + + return static_cast(2) * (mu10 * S_mu_10 + mu01 * S_mu_01); + }(); + + const auto S_volumetric = [this, J, &C_1]() + { + const auto bulk = m_bulkModulus.getValue(); + return bulk * log(J) * C_1; + }(); + + return S_isochoric + S_volumetric; +} + +template +auto MooneyRivlinMaterial::elasticityTensor(Strain& strain) -> ElasticityTensor +{ + static constexpr Real dim_1 = static_cast(1) / static_cast(spatial_dimensions); + const auto& C = strain.getRightCauchyGreenTensor(); + const auto J = strain.getDeterminantDeformationGradient(); + const auto logJ = log(J); + const auto J_2dim = pow(J, -2 * dim_1); + const auto J_4dim = pow(J, -4 * dim_1); + const auto C_1 = sofa::type::inverse(C); + const auto I1 = strain.getInvariant1(); + const auto I2 = strain.getInvariant2(); + + const auto mu01 = m_mu01.getValue(); + const auto mu10 = m_mu10.getValue(); + const auto bulk = m_bulkModulus.getValue(); + + auto delta = [](auto i, auto j){ return sofa::component::solidmechanics::fem::elastic::kroneckerDelta(i, j); }; + + return ElasticityTensor( + [&](sofa::Index i, sofa::Index j, sofa::Index k, sofa::Index l) + { + //derivative of C^{-1} with respect to C + const auto dC_1dC = -static_cast(0.5) * (C_1(i, k) * C_1(l, j) + C_1(i, l) * C_1(k, j)); + + // the derivative of S_mu_10 with respect to C + // each term has both minor and major symmetries + const Real dS_mu_10dC = -dim_1 * J_2dim * ( + delta(i, j) * C_1(k, l) + delta(k, l) * C_1(i, j) + - dim_1 * C_1(k, l) * I1 * C_1(i, j) + + I1 * dC_1dC + ); + + // The derivative of S_mu_01 with respect to C + // the terms have been grouped to highlight the minor and major symmetries + const Real dS_mu_01dC = J_4dim * ( + - 2 * dim_1 * ( + C_1(k, l) * (I1 * delta(i, j) - C(i, j)) + + C_1(i, j) * (I1 * delta(k, l) - C(k, l)) + + dC_1dC * I2 + ) + + delta(i, j) * delta(k, l) + - static_cast(0.5) * ( + delta(i, k) * delta(j, l) + + delta(i, l) * delta(j, k)) + + 4 * dim_1 * dim_1 * C_1(k, l) * C_1(i, j) * I2 + ); + + const Real dS_isochoric_dC = 2 * (mu10 * dS_mu_10dC + mu01 * dS_mu_01dC); + + // the derivative of S_volumetric with respect to C + // this term has both minor and major symmetries + const Real dS_volumetric_dC = bulk * (C_1(l, k) * C_1(i, j) / 2 + logJ * dC_1dC); + + return 2 * (dS_isochoric_dC + dS_volumetric_dC); + }); +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.cpp new file mode 100644 index 00000000000..c9c13f247d3 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.cpp @@ -0,0 +1,44 @@ +/****************************************************************************** +* 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_MATERIAL_NEOHOOKEANMATERIAL_CPP + +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +void registerNeoHookeanMaterial(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Neo-Hookean material") + .add< NeoHookeanMaterial >() + .add< NeoHookeanMaterial >() + .add< NeoHookeanMaterial >(true)); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API NeoHookeanMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API NeoHookeanMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API NeoHookeanMaterial; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.h new file mode 100644 index 00000000000..d981c2b7f52 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.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 + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_NEOHOOKEANMATERIAL_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +/** + * Formulation from Bonet, J. and R. D. Wood (2008). Nonlinear continuum mechanics for finite + * element analysis. Cambridge university press. + */ +template +class NeoHookeanMaterial : + public PK2HyperelasticMaterial, + public LinearMechanicalParametersComponent +{ +public: + SOFA_CLASS2(NeoHookeanMaterial, PK2HyperelasticMaterial, LinearMechanicalParametersComponent); + +private: + using Real = sofa::Real_t; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + using DeformationGradient = typename PK2HyperelasticMaterial::DeformationGradient; + using RightCauchyGreenTensor = typename PK2HyperelasticMaterial::RightCauchyGreenTensor; + using StressTensor = typename PK2HyperelasticMaterial::StressTensor; + using ElasticityTensor = typename PK2HyperelasticMaterial::ElasticityTensor; + + using LinearMechanicalParametersComponent::m_lambda; + using LinearMechanicalParametersComponent::m_mu; + +public: + StressTensor secondPiolaKirchhoffStress(Strain& strain) override; + + ElasticityTensor elasticityTensor(Strain& strain) override; +}; + + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_NEOHOOKEANMATERIAL_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API NeoHookeanMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API NeoHookeanMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API NeoHookeanMaterial; +#endif +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.inl new file mode 100644 index 00000000000..113158f3e93 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/NeoHookeanMaterial.inl @@ -0,0 +1,58 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +auto NeoHookeanMaterial::secondPiolaKirchhoffStress(Strain& strain) -> StressTensor +{ + static constexpr auto& I = Strain::identity; + const auto& C = strain.getRightCauchyGreenTensor(); + + const DeformationGradient C_1 = sofa::type::inverse(C); + const Real J = strain.getDeterminantDeformationGradient(); + + return m_mu * (I - C_1) + m_lambda * std::log(J) * C_1; +} + +template +auto NeoHookeanMaterial::elasticityTensor(Strain& strain) -> ElasticityTensor +{ + const auto& C = strain.getRightCauchyGreenTensor(); + const RightCauchyGreenTensor C_1 = sofa::type::inverse(C); + const Real J = strain.getDeterminantDeformationGradient(); + const Real logJ = std::log(J); + + return ElasticityTensor( + [mu = m_mu, lambda = m_lambda, &C_1, logJ](sofa::Index i, sofa::Index j, sofa::Index k, sofa::Index l) + { + return (mu - lambda * logJ) * (C_1(i, k) * C_1(l, j) + C_1(i, l) * C_1(k, j)) + + lambda * C_1(l, k) * C_1(i, j); + }); +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.cpp new file mode 100644 index 00000000000..6440da0ce51 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.cpp @@ -0,0 +1,44 @@ +/****************************************************************************** +* 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_MATERIAL_OGDEN_CPP + +#include +#include + +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +void registerOgdenMaterial(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Ogden material") + .add< OgdenMaterial >() + .add< OgdenMaterial >() + .add< OgdenMaterial >(true)); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API OgdenMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API OgdenMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API OgdenMaterial; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.h new file mode 100644 index 00000000000..cefeb525d4c --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.h @@ -0,0 +1,71 @@ +/****************************************************************************** +* 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 + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_OGDEN_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +template +class OgdenMaterial : + public PK2HyperelasticMaterial +{ +public: + SOFA_CLASS(OgdenMaterial, PK2HyperelasticMaterial); + +private: + using Real = sofa::Real_t; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + using DeformationGradient = typename PK2HyperelasticMaterial::DeformationGradient; + using RightCauchyGreenTensor = typename PK2HyperelasticMaterial::RightCauchyGreenTensor; + using StressTensor = typename PK2HyperelasticMaterial::StressTensor; + using ElasticityTensor = typename PK2HyperelasticMaterial::ElasticityTensor; + +public: + // Ogden material constants + sofa::Data m_mu; + sofa::Data m_alpha; + sofa::Data m_kappa; + + StressTensor secondPiolaKirchhoffStress(Strain& strain) override; + + ElasticityTensor elasticityTensor(Strain& strain) override; + +protected: + OgdenMaterial(); +}; + + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_OGDEN_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API OgdenMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API OgdenMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API OgdenMaterial; +#endif +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl new file mode 100644 index 00000000000..348453461ac --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/OgdenMaterial.inl @@ -0,0 +1,223 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +OgdenMaterial::OgdenMaterial() + : m_mu(initData(&m_mu, static_cast(1e2), "mu", + "Material constant relevant to the shear modulus")) + , m_alpha(initData(&m_alpha, static_cast(1.5), "alpha", + "Material constant exponent related to strain-stiffening")) + , m_kappa(initData(&m_kappa, static_cast(1e3), "kappa", + "Material constant related to the bulk modulus")) +{} + +template +auto OgdenMaterial::secondPiolaKirchhoffStress(Strain& strain) -> StressTensor +{ + using EigenMatrix = Eigen::Matrix; + + const auto& C = strain.getRightCauchyGreenTensor(); + const auto C_1 = sofa::type::inverse(C); + + const auto J = strain.getDeterminantDeformationGradient(); + assert(J > 0); + + const auto S_isochoric = [this, J, &strain, &C_1, &C] () + { + const Real mu = m_mu.getValue(); + const Real alpha = m_alpha.getValue(); + + const Real FJ = pow(J, -alpha / static_cast(3)); + + EigenMatrix CEigen; + for (sofa::Index m = 0; m < spatial_dimensions; ++m) + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + CEigen(m, n) = C(m, n); + + // The SelfAdjointEigenSolver is normally the suitable choice for symmetric matrices + // However, it has produced wrong results in some previous tests + // Disable temporarilly until fixed /*Eigen::SelfAdjointEigenSolver*/ + Eigen::EigenSolver EigenProblemSolver(CEigen, true); + if (EigenProblemSolver.info() != Eigen::Success) + dmsg_warning("OgdenMaterial") << "EigenSolver iterations failed to converge"; + + const auto eigenvectors = EigenProblemSolver.eigenvectors().real().eval(); + const auto eigenvalues = EigenProblemSolver.eigenvalues().real().eval(); + + // Precompute trace of C^(alpha/2) and eigenvalue powers: lambda^(alpha/2 - 1) + const Real aBy2 = static_cast(0.5) * alpha; + const Real aBy2Minus1 = aBy2 - static_cast(1); + Real trCaBy2{static_cast(0)}; + typename DataTypes::Coord eigenvaluePowers{}; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + { + trCaBy2 += pow(eigenvalues(n), aBy2); + eigenvaluePowers(n) = pow(eigenvalues(n), aBy2Minus1); + } + + const Real coefS1 = FJ * mu / (alpha * static_cast(2)); + const Real coefS2 = -FJ * mu / (alpha * static_cast(6)) * trCaBy2; + + StressTensor S{}; + + for (sofa::Index i = 0; i < spatial_dimensions; ++i) + { + for (sofa::Index j = 0; j < spatial_dimensions; ++j) + { + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + S(i,j) += coefS1 * eigenvaluePowers(n) * eigenvectors(i,n) * eigenvectors(j,n); + S(i,j) += coefS2 * C_1(j, i); + } + } + + return static_cast(2) * S; + }(); + + const auto S_volumetric = [this, J, &C_1]() + { + const Real kappa = m_kappa.getValue(); + return kappa * log(J) * C_1; + }(); + + return S_isochoric + S_volumetric; +} + +template +auto OgdenMaterial::elasticityTensor(Strain& strain) -> ElasticityTensor +{ + using EigenMatrix = Eigen::Matrix; + + const Real mu = m_mu.getValue(); + const Real alpha = m_alpha.getValue(); + const Real kappa = m_kappa.getValue(); + + const auto& C = strain.getRightCauchyGreenTensor(); + const auto C_1 = sofa::type::inverse(C); + + const Real J = strain.getDeterminantDeformationGradient(); + assert(J > 0); + const Real FJ = pow(J, -alpha / static_cast(3)); + + EigenMatrix CEigen; + for (sofa::Index m = 0; m < spatial_dimensions; ++m) + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + CEigen(m, n) = C(m, n); + + // Disable temporarilly until fixed /*Eigen::SelfAdjointEigenSolver*/ + Eigen::EigenSolver EigenProblemSolver(CEigen, true); + if (EigenProblemSolver.info() != Eigen::Success) + dmsg_warning("OgdenMaterial") << "EigenSolver iterations failed to converge"; + const auto eigenvectors = EigenProblemSolver.eigenvectors().real().eval(); + const auto eigenvalues = EigenProblemSolver.eigenvalues().real().eval(); + + // Precompute trace of C^(alpha/2) and eigenvalue powers: lambda^(a/2 - 1) and lambda^(a/2 - 2) + const Real aBy2 = static_cast(0.5) * alpha; + const Real aBy2Minus1 = aBy2 - static_cast(1); + Real trCaBy2{static_cast(0)}; + typename DataTypes::Coord eigenvaluePowers1{}, eigenvaluePowers2{}; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + { + trCaBy2 += pow(eigenvalues(n), aBy2); + eigenvaluePowers1(n) = pow(eigenvalues(n), aBy2Minus1); + eigenvaluePowers2(n) = pow(eigenvalues(n), aBy2Minus1 - static_cast(1)); + } + + const Real coefSpectral = FJ * mu / alpha * aBy2Minus1; + const Real coefRotational = FJ * mu / alpha; + + // Precompute C^(alpha1/2 - 1) from eigenbasis: V * D * V^T; D_i = lambda_i^(alpha1/2 - 1) + RightCauchyGreenTensor D, EigenBasis; + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + { + for (sofa::Index m = 0; m < spatial_dimensions; ++m) + EigenBasis(n,m) = eigenvectors(n,m); + D(n,n) = eigenvaluePowers1(n); + } + const RightCauchyGreenTensor CaBy2Minus1(EigenBasis*D*(EigenBasis.transposed())); + + return ElasticityTensor( + [&](sofa::Index i, sofa::Index j, sofa::Index k, sofa::Index l) + { + // Derivative of S_isochoric with respect to C + const Real T_isochoric = [&]() + { + Real sum{static_cast(0)}; + + // Spectral decomposition contributions; derivatives of eigenvalues and eigenvectors + for (sofa::Index n = 0; n < spatial_dimensions; ++n) + { + sum += coefSpectral * eigenvaluePowers2(n) * eigenvectors(i,n) * eigenvectors(j,n) + * eigenvectors(k,n) * eigenvectors(l,n); + + for (sofa::Index m = 0; m < spatial_dimensions; ++m) + { + // This sum is taken over only for m != n + if (m == n) continue; + + // Eigenvalue multiplicity (val1 == val2) causes an indeterminate form; + // In that case, the limit of the expression is used + const Real eigenvalueMultiplier = + std::fabs(eigenvalues[n] - eigenvalues[m]) < std::numeric_limits::epsilon() + ? aBy2Minus1 * eigenvaluePowers2(n) + : (eigenvaluePowers1(n) - eigenvaluePowers1(m)) / (eigenvalues(n) - eigenvalues(m)); + + sum += static_cast(0.5) * coefRotational * eigenvalueMultiplier * + (eigenvectors(i,n) * eigenvectors(j,m) * eigenvectors(k,m) * eigenvectors(l,n) + + eigenvectors(i,n) * eigenvectors(j,m) * eigenvectors(k,n) * eigenvectors(l,m)); + } + } + + // Remaining terms grouped to highlight respective contribution from + // differentiating coefS1, coefS2 and C^(-1) from S_isochoric + sum += FJ * mu * ( + // Contribution from derivative of F(J) + trCaBy2 / static_cast(18) * C_1(j,i) * C_1(l,k) - + // Contribution from derivatives of C^(alpha/2 - 1) and tr(C^(alpha/2)) + (CaBy2Minus1(j,i) * C_1(l,k) + C_1(j,i) * CaBy2Minus1(k,l)) + / static_cast(6) + + // Contribution from derivative of C^(-1) + trCaBy2 / (static_cast(6) * alpha) + * (C_1(j,k) * C_1(l,i) + C_1(j,l) * C_1(k,i)) + ); + + return sum; + }(); + + // Derivative of S_volumetric with respect to C + const Real T_volumetric = static_cast(0.5) * kappa * C_1(l,k) * C_1(j,i) + - static_cast(0.5) * kappa * log(J) + * (C_1(j,k) * C_1(l,i) + C_1(j,l) * C_1(k,i)); + + return static_cast(2) * (T_isochoric + T_volumetric); + }); +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.cpp new file mode 100644 index 00000000000..32abfcb114b --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.cpp @@ -0,0 +1,44 @@ +/****************************************************************************** +* 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_MATERIAL_ST_VENANT_KIRCHHOFF_MATERIAL_CPP + +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +void registerStVenantKirchhoffMaterial(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Saint Venant-Kirchhoff material model for hyperelastic materials") + .add< StVenantKirchhoffMaterial >() + .add< StVenantKirchhoffMaterial >() + .add< StVenantKirchhoffMaterial >(true)); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API StVenantKirchhoffMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API StVenantKirchhoffMaterial; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API StVenantKirchhoffMaterial; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.h new file mode 100644 index 00000000000..20d9cc37af9 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.h @@ -0,0 +1,86 @@ +/****************************************************************************** +* 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_MATERIAL_STVENANTKIRCHHOFFMATERIAL_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::hyperelastic +{ + +/** + * @class StVenantKirchhoffMaterial + * @brief Represents the St. Venant-Kirchhoff material model for hyperelastic materials. + * + * The St. Venant-Kirchhoff material is a simple model used to describe the stress-strain behavior + * of isotropic hyperelastic materials. It is based on the linear elastic theory but is extended + * to handle large deformations using the Green-Lagrange strain tensor and the Second + * Piola-Kirchhoff stress tensor. + * + * This material model is only valid for cases where the deformation is relatively small, despite + * being nonlinear in nature. Its application is mainly limited due to its inability to correctly + * predict behavior under significant strain, as it does not accurately represent material + * nonlinearity under large deformation. + * + * The material model is defined using two parameters: + * - The Young's modulus: Describes the material's stiffness. + * - The Poisson's ratio: Represents the material's ability to undergo deformation in directions + * orthogonal to the applied stress. + */ +template +class StVenantKirchhoffMaterial: + public PK2HyperelasticMaterial, + public LinearMechanicalParametersComponent +{ +public: + SOFA_CLASS2(StVenantKirchhoffMaterial, PK2HyperelasticMaterial, LinearMechanicalParametersComponent); + +private: + using Real = sofa::Real_t; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + using DeformationGradient = typename PK2HyperelasticMaterial::DeformationGradient; + using RightCauchyGreenTensor = typename PK2HyperelasticMaterial::RightCauchyGreenTensor; + using StressTensor = typename PK2HyperelasticMaterial::StressTensor; + using ElasticityTensor = typename PK2HyperelasticMaterial::ElasticityTensor; + + using LinearMechanicalParametersComponent::m_lambda; + using LinearMechanicalParametersComponent::m_mu; + +protected: + StressTensor secondPiolaKirchhoffStress(Strain& strain) override; + + ElasticityTensor elasticityTensor(Strain& strain) override; +}; + +#if !defined(ELASTICITY_COMPONENT_MATERIAL_STVENANTKIRCHHOFFMATERIAL_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API StVenantKirchhoffMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API StVenantKirchhoffMaterial; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_HYPERELASTIC_API StVenantKirchhoffMaterial; +#endif +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.inl new file mode 100644 index 00000000000..983999373a2 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/hyperelasticmaterials/StVenantKirchhoffMaterial.inl @@ -0,0 +1,58 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +template +auto StVenantKirchhoffMaterial::secondPiolaKirchhoffStress(Strain& strain) +-> StressTensor +{ + static const auto& I = sofa::type::Mat::Identity(); + + // Green-Lagrangian strain tensor + const auto& E = strain.getGreenLagrangeTensor(); + + // Second Piola-Kirchhoff stress tensor + return m_lambda * sofa::type::trace(E) * I + static_cast(2) * m_mu * E; +} + +template +auto StVenantKirchhoffMaterial::elasticityTensor(Strain& strain) -> ElasticityTensor +{ + SCOPED_TIMER_TR("elasticityTensor"); + SOFA_UNUSED(strain); + + sofa::component::solidmechanics::fem::elastic::LameLambda lambda { m_lambda }; + sofa::component::solidmechanics::fem::elastic::LameMu mu { m_mu }; + + return sofa::component::solidmechanics::fem::elastic::makeIsotropicElasticityTensor(mu, lambda); +} + +} // namespace elasticity diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/impl/MajorSymmetric4Tensor.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/impl/MajorSymmetric4Tensor.h new file mode 100644 index 00000000000..d06dd92ebfe --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/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 sofa::component::solidmechanics::fem::hyperelastic +{ + +/** + * 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/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/impl/Strain.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/impl/Strain.h new file mode 100644 index 00000000000..f7a6c52dde0 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/impl/Strain.h @@ -0,0 +1,273 @@ +/****************************************************************************** +* 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::hyperelastic +{ + +/** + * @struct DeformationGradientTag + * @brief Tag for deformation gradient initialization + */ +inline struct DeformationGradientTag {} deformationGradient; + +/** + * @struct RightCauchyGreenTensorTag + * @brief Tag for right Cauchy-Green tensor initialization + */ +inline struct RightCauchyGreenTensorTag {} rightCauchyGreenTensor; + +/** + * @class Strain + * @brief Class to compute and manage strain measures in elasticity simulations + * + * This class handles the computation of strain measures (deformation gradient, right Cauchy-Green tensor, Green-Lagrange tensor, and strain invariants) + * using an optional pattern for lazy computation. It supports initialization from either the deformation gradient (F) or the right Cauchy-Green tensor (C). + * + * @tparam DataTypes The data type of the simulation. Must provide: + * - `spatial_dimensions`: Spatial dimension of the problem (2 or 3) + * - `Real_t`: Scalar type for numerical operations + * + * @note The class follows the following invariant rules: + * - When initialized with F: Right Cauchy-Green tensor (C) is computed as C = F^T * F + * - When initialized with C: Deformation gradient (F) is not computed + * - All computed values are stored in `std::optional` to avoid unnecessary computations + * - Invariants are defined as: + * I1 = tr(C) + * I2 = (I1^2 - tr(C^2)) / 2 + * I3 = det(C) = (det(F))^2 + */ +template +struct Strain +{ + /** + * @brief Spatial dimension of the problem (2 or 3) + */ + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + using Real = sofa::Real_t; + + /** + * @brief Type for deformation gradient tensor (spatial x spatial) + */ + using DeformationGradient = sofa::type::Mat; + + /** + * @brief Type for right Cauchy-Green tensor (spatial x spatial) + */ + using RightCauchyGreenTensor = sofa::type::Mat; + + /** + * @brief Type for Green-Lagrange tensor (spatial x spatial) + */ + using GreenLagrangeTensor = sofa::type::Mat; + + static constexpr sofa::type::IdentityMatrix identity {}; + + /** + * @brief Constructor initializing from deformation gradient + * + * @param tag Tag for deformation gradient initialization + * @param F Deformation gradient tensor (F) + */ + constexpr Strain(DeformationGradientTag, const DeformationGradient& F) : m_deformationGradient(F) {} + + /** + * @brief Constructor initializing from right Cauchy-Green tensor + * + * @param tag Tag for right Cauchy-Green tensor initialization + * @param C Right Cauchy-Green tensor (C) + */ + constexpr Strain(RightCauchyGreenTensorTag, const RightCauchyGreenTensor& C) : m_rightCauchyGreenTensor(C) {} + + const DeformationGradient& deformationGradient() const + { + if (!m_deformationGradient.has_value()) + { + throw std::logic_error("Cannot access deformation gradient. It has not been provided"); + } + return *m_deformationGradient; + } + + void setDeformationGradient(const DeformationGradient& F) + { + reset(); + m_deformationGradient = F; + } + + Real getDeterminantDeformationGradient() + { + if (!m_determinantF) + { + computeDeterminant(); + } + return *m_determinantF; + } + + const RightCauchyGreenTensor& getRightCauchyGreenTensor() + { + if (!m_rightCauchyGreenTensor) + { + computeRightCauchyGreenTensor(); + } + return *m_rightCauchyGreenTensor; + } + + void setRightCauchyGreenTensor(const RightCauchyGreenTensor& C) + { + reset(); + m_rightCauchyGreenTensor = C; + } + + const GreenLagrangeTensor& getGreenLagrangeTensor() + { + if (!m_greenLagrangeTensor) + { + computeGreenLagrangeTensor(); + } + return *m_greenLagrangeTensor; + } + + Real getInvariant1() + { + if (!m_invariant1) + { + computeInvariant1(); + } + return *m_invariant1; + } + + Real getInvariant2() + { + if (!m_invariant2) + { + computeInvariant2(); + } + return *m_invariant2; + } + + Real getInvariant3() + { + if (!m_invariant3) + { + computeInvariant3(); + } + return *m_invariant3; + } + + /** + * @brief Reset all computed values + * + * Clears all optional values (deformation gradient, C, invariants, tensors) + */ + void reset() + { + m_deformationGradient.reset(); + m_rightCauchyGreenTensor.reset(); + m_determinantF.reset(); + m_greenLagrangeTensor.reset(); + m_invariant1.reset(); + m_invariant2.reset(); + m_invariant3.reset(); + } + +protected: + + std::optional m_deformationGradient; + std::optional m_rightCauchyGreenTensor; + + std::optional m_determinantF; + std::optional m_greenLagrangeTensor; + std::optional m_invariant1; + std::optional m_invariant2; + std::optional m_invariant3; + + void computeDeterminant() + { + if (m_deformationGradient.has_value()) + { + m_determinantF = sofa::type::determinant(*m_deformationGradient); + } + else if (m_rightCauchyGreenTensor.has_value()) + { + m_determinantF = sqrt(sofa::type::determinant(*m_rightCauchyGreenTensor)); + } + else + { + throw std::logic_error("Cannot compute determinant. Neither deformation gradient nor right Cauchy-Green tensor is provided"); + } + } + + void computeRightCauchyGreenTensor() + { + if (m_deformationGradient.has_value()) + { + const auto& F = *m_deformationGradient; + auto& C = m_rightCauchyGreenTensor.emplace(); + 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) + { + C(i, j) += F(k, i) * F(k, j); + } + } + } + } + else + { + throw std::logic_error("Cannot compute right Cauchy-Green tensor. Deformation gradient is not provided"); + } + } + + void computeGreenLagrangeTensor() + { + const auto& C = getRightCauchyGreenTensor(); + m_greenLagrangeTensor = static_cast(0.5) * (C - identity); + } + + void computeInvariant1() + { + const auto& C = getRightCauchyGreenTensor(); + m_invariant1 = sofa::type::trace(C); + } + + void computeInvariant2() + { + const auto I1 = getInvariant1(); + const auto& C = getRightCauchyGreenTensor(); + const auto trC2 = sofa::type::trace(C * C); + m_invariant2 = static_cast(0.5) * (I1 * I1 - trC2); + } + + void computeInvariant3() + { + const auto J = getDeterminantDeformationGradient(); + m_invariant3 = J * J; + } +}; + +} diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/init.cpp b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/init.cpp index 10d9577fb07..f0e2aa665e9 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/init.cpp +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/init.cpp @@ -25,7 +25,12 @@ namespace sofa::component::solidmechanics::fem::hyperelastic { + extern void registerHyperelasticityFEMForceField(sofa::core::ObjectFactory* factory); + extern void registerMooneyRivlinMaterial(sofa::core::ObjectFactory* factory); + extern void registerNeoHookeanMaterial(sofa::core::ObjectFactory* factory); + extern void registerOgdenMaterial(sofa::core::ObjectFactory* factory); extern void registerPlasticMaterial(sofa::core::ObjectFactory* factory); + extern void registerStVenantKirchhoffMaterial(sofa::core::ObjectFactory* factory); extern void registerStandardTetrahedralFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTetrahedronHyperelasticityFEMForceField(sofa::core::ObjectFactory* factory); @@ -53,7 +58,12 @@ const char* getModuleVersion() void registerObjects(sofa::core::ObjectFactory* factory) { + registerHyperelasticityFEMForceField(factory); + registerMooneyRivlinMaterial(factory); + registerNeoHookeanMaterial(factory); + registerOgdenMaterial(factory); registerPlasticMaterial(factory); + registerStVenantKirchhoffMaterial(factory); registerStandardTetrahedralFEMForceField(factory); registerTetrahedronHyperelasticityFEMForceField(factory); } diff --git a/Sofa/framework/Type/CMakeLists.txt b/Sofa/framework/Type/CMakeLists.txt index 9d8596f2809..822634784d8 100644 --- a/Sofa/framework/Type/CMakeLists.txt +++ b/Sofa/framework/Type/CMakeLists.txt @@ -12,6 +12,7 @@ set(HEADER_FILES ${SOFATYPESRC_ROOT}/Frame.h ${SOFATYPESRC_ROOT}/FullySymmetric4Tensor.h ${SOFATYPESRC_ROOT}/hardening.h + ${SOFATYPESRC_ROOT}/IdentityMatrix.h ${SOFATYPESRC_ROOT}/Mat.h ${SOFATYPESRC_ROOT}/MatSym.h ${SOFATYPESRC_ROOT}/Mat_solve_Cholesky.h diff --git a/Sofa/framework/Type/src/sofa/type/IdentityMatrix.h b/Sofa/framework/Type/src/sofa/type/IdentityMatrix.h new file mode 100644 index 00000000000..fd7c16227ab --- /dev/null +++ b/Sofa/framework/Type/src/sofa/type/IdentityMatrix.h @@ -0,0 +1,109 @@ +#pragma once + +#include + +namespace sofa::type +{ + +struct IdentityMatrix +{}; + +template +struct ScaledIdentityMatrix +{ + real scale; +}; + +template +constexpr Mat operator+(const IdentityMatrix& I, const Mat& M) +{ + Mat res(M); + for (sofa::Size i = 0; i < N; ++i) + { + res[i][i] += static_cast(1); + } + return res; +} + +template +constexpr Mat operator+(const Mat& M, const IdentityMatrix& I) +{ + return I + M; +} + +template +constexpr Mat operator-(const IdentityMatrix& I, const Mat& M) +{ + Mat res(-M); + for (sofa::Size i = 0; i < N; ++i) + { + res[i][i] += static_cast(1); + } + return res; +} + +template +constexpr Mat operator-(const Mat& M, const IdentityMatrix& I) +{ + Mat res(M); + for (sofa::Size i = 0; i < N; ++i) + { + res[i][i] -= static_cast(1); + } + return res; +} + +template +constexpr ScaledIdentityMatrix operator*(const IdentityMatrix& I, real s) +{ + return { s }; +} + +template +constexpr ScaledIdentityMatrix operator*(real s, const IdentityMatrix& I) +{ + return { s }; +} + +template +constexpr Mat operator-(const ScaledIdentityMatrix& I, const Mat& M) +{ + Mat res(-M); + for (sofa::Size i = 0; i < N; ++i) + { + res[i][i] += I.scale; + } + return res; +} + +template +constexpr Mat operator-(const Mat& M, const ScaledIdentityMatrix& I) +{ + Mat res(M); + for (sofa::Size i = 0; i < N; ++i) + { + res[i][i] -= I.scale; + } + return res; +} + +template +constexpr Mat operator+(const ScaledIdentityMatrix& I, const Mat& M) +{ + Mat res(M); + for (sofa::Size i = 0; i < N; ++i) + { + res[i][i] += I.scale; + } + return res; +} + +template +constexpr Mat operator+(const Mat& M, const ScaledIdentityMatrix& I) +{ + return I + M; +} + + + +} diff --git a/examples/Validation/cantilever_beam/CantileverBeam_ElementFEMForceField.xml b/examples/Validation/cantilever_beam/CantileverBeam_ElementFEMForceField.xml index 157de731557..7b950015a83 100644 --- a/examples/Validation/cantilever_beam/CantileverBeam_ElementFEMForceField.xml +++ b/examples/Validation/cantilever_beam/CantileverBeam_ElementFEMForceField.xml @@ -10,6 +10,7 @@ + diff --git a/examples/Validation/cantilever_beam/hexahedron/svk/parallel/assembled/HyperelasticityFEMForceField.scn b/examples/Validation/cantilever_beam/hexahedron/svk/parallel/assembled/HyperelasticityFEMForceField.scn new file mode 100644 index 00000000000..48cb648f56b --- /dev/null +++ b/examples/Validation/cantilever_beam/hexahedron/svk/parallel/assembled/HyperelasticityFEMForceField.scn @@ -0,0 +1,12 @@ + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/hexahedron/svk/parallel/matrixfree/HyperelasticityFEMForceField.scn b/examples/Validation/cantilever_beam/hexahedron/svk/parallel/matrixfree/HyperelasticityFEMForceField.scn new file mode 100644 index 00000000000..75bbfd8107b --- /dev/null +++ b/examples/Validation/cantilever_beam/hexahedron/svk/parallel/matrixfree/HyperelasticityFEMForceField.scn @@ -0,0 +1,11 @@ + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/hexahedron/svk/sequential/assembled/HyperelasticityFEMForceField.scn b/examples/Validation/cantilever_beam/hexahedron/svk/sequential/assembled/HyperelasticityFEMForceField.scn new file mode 100644 index 00000000000..4f8e213cf4f --- /dev/null +++ b/examples/Validation/cantilever_beam/hexahedron/svk/sequential/assembled/HyperelasticityFEMForceField.scn @@ -0,0 +1,13 @@ + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/hexahedron/svk/sequential/matrixfree/HyperelasticityFEMForceField.scn b/examples/Validation/cantilever_beam/hexahedron/svk/sequential/matrixfree/HyperelasticityFEMForceField.scn new file mode 100644 index 00000000000..26d221c61b8 --- /dev/null +++ b/examples/Validation/cantilever_beam/hexahedron/svk/sequential/matrixfree/HyperelasticityFEMForceField.scn @@ -0,0 +1,11 @@ + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/tetrahedron/svk/parallel/assembled/HyperelasticityFEMForceField.scn b/examples/Validation/cantilever_beam/tetrahedron/svk/parallel/assembled/HyperelasticityFEMForceField.scn new file mode 100644 index 00000000000..2d7cbda3f67 --- /dev/null +++ b/examples/Validation/cantilever_beam/tetrahedron/svk/parallel/assembled/HyperelasticityFEMForceField.scn @@ -0,0 +1,20 @@ + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/tetrahedron/svk/parallel/matrixfree/HyperelasticityFEMForceField.scn b/examples/Validation/cantilever_beam/tetrahedron/svk/parallel/matrixfree/HyperelasticityFEMForceField.scn new file mode 100644 index 00000000000..8234a43ee78 --- /dev/null +++ b/examples/Validation/cantilever_beam/tetrahedron/svk/parallel/matrixfree/HyperelasticityFEMForceField.scn @@ -0,0 +1,18 @@ + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/tetrahedron/svk/sequential/assembled/HyperelasticityFEMForceField.scn b/examples/Validation/cantilever_beam/tetrahedron/svk/sequential/assembled/HyperelasticityFEMForceField.scn new file mode 100644 index 00000000000..1228b5a5f02 --- /dev/null +++ b/examples/Validation/cantilever_beam/tetrahedron/svk/sequential/assembled/HyperelasticityFEMForceField.scn @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/tetrahedron/svk/sequential/matrixfree/HyperelasticityFEMForceField.scn b/examples/Validation/cantilever_beam/tetrahedron/svk/sequential/matrixfree/HyperelasticityFEMForceField.scn new file mode 100644 index 00000000000..02886d334dd --- /dev/null +++ b/examples/Validation/cantilever_beam/tetrahedron/svk/sequential/matrixfree/HyperelasticityFEMForceField.scn @@ -0,0 +1,18 @@ + + + + + + + + + + + + + + + + +