diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index 85447e27cdc..f086bbf1c90 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt @@ -11,11 +11,20 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseElementLinearFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseSourceTerm.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CauchyStressEvaluator.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.h + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/PressureSourceTerm.h + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/PressureSourceTerm.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/StressSourceTerm.h + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/StressSourceTerm.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VectorSourceTerm.h + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VectorSourceTerm.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.h @@ -70,9 +79,14 @@ set(SOURCE_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/init.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseElementLinearFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseSourceTerm.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BeamFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/CorotationalFEMForceField.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMSourceTermIntegrator.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/LinearSmallStrainFEMForceField.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/PressureSourceTerm.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/StressSourceTerm.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VectorSourceTerm.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FastTetrahedralCorotationalForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/FEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/HexahedralFEMForceField.cpp diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.cpp new file mode 100644 index 00000000000..b6cc5a93170 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.cpp @@ -0,0 +1,42 @@ +/****************************************************************************** +* 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 SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_BASE_SOURCE_TERM_CPP + +#include + +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h new file mode 100644 index 00000000000..6e3c6ed0cdc --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseSourceTerm.h @@ -0,0 +1,198 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_BASE_SOURCE_TERM_CPP) +#include +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @struct QuadratureContext + * @brief Everything the integrator knows at one quadrature point. + * + * Built once per quadrature point and handed to every integrated term. + * A source term reads from it and returns an integrand. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + * @tparam TElementType The type of finite element (e.g., sofa::geometry::Tetrahedron). + */ +template +struct QuadratureContext +{ + using DataTypes = TDataTypes; + using ElementType = TElementType; + using FiniteElement = sofa::fem::FiniteElement; + + using Real = sofa::Real_t; + using Coord = sofa::Coord_t; + using Deriv = sofa::Deriv_t; + + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; + + using Element = typename FiniteElement::TopologyElement; + using ShapeFunctions = sofa::type::Vec; + using GradientShapeFunctions = sofa::type::Mat; + using Jacobian = sofa::type::Mat; + + /// Node indices of the element being integrated, with which a term gathers its own nodal + /// degrees of freedom. + const Element& element; + + /// Shape function value at this quadrature point. + ShapeFunctions N; + + /// Reference-space gradients of the shape functions at this quadrature point. + GradientShapeFunctions gradientShapeFunctions; + + /// dx/dq of the reference-to-physical mapping, on the configuration the integrator chose. + Jacobian jacobian; + + /// \f$ |\det J| \f$, for information only: the integrator applies it, a term must not. + Real measure; + + /// Interpolated rest position at this quadrature point. + Coord restPosition; + + /// Interpolated displacement at this quadrature point. + Deriv displacement; +}; + +/** + * @brief Unit normal of an element, from the jacobian of its mapping. + * + * Defined only where the element spans one dimension less than the space it lives in: a surface + * element in 3D, an edge in 2D. Its orientation follows the node ordering of the element. + * + * @param jacobian dx/dq of the reference-to-physical mapping at the point of interest. + */ +template +sofa::type::Vec elementNormal( + const sofa::type::Mat& jacobian) +{ + static_assert(TopologicalDimension + 1 == spatial_dimensions, + "A normal is only defined for an element of codimension 1."); + + if constexpr (spatial_dimensions == 3) + { + return jacobian.col(0).cross(jacobian.col(1)).normalized(); + } + else + { + const sofa::type::Vec<2, Real> tangent = jacobian.col(0); + return sofa::type::Vec<2, Real>(tangent[1], -tangent[0]).normalized(); + } +} + +/** + * @class BaseSourceTerm + * @brief A source density whose value is determined by the geometry, not by the solution. + * + * The component calculates the integrand in evaluate() given a QuadratureContext. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + * @tparam TElementType The type of finite element (e.g., sofa::geometry::Tetrahedron). + */ +template +class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent +{ +public: + using DataTypes = TDataTypes; + using ElementType = TElementType; + + SOFA_CLASS(SOFA_TEMPLATE2(BaseSourceTerm, DataTypes, ElementType), + sofa::core::objectmodel::BaseComponent); + + using Deriv = sofa::Deriv_t; + using QuadratureContext_t = QuadratureContext; + + /** + * @brief Source density at one quadrature point, per unit physical measure. + * + * @param context Geometry of the quadrature point. + */ + virtual Deriv evaluate(const QuadratureContext_t& context) const = 0; + +protected: + + BaseSourceTerm() = default; + + /** + * @brief Value of a nodal property interpolated at the quadrature point. + * + * The property is gathered at the nodes of the element being integrated and combined with the + * shape functions evaluated at that point. + * + * @param property The component holding the nodal values. + * @param context Geometry of the quadrature point. + */ + template + static PropertyType interpolateProperty( + const sofa::core::BaseNodalProperty& property, + const QuadratureContext_t& context) + { + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + + sofa::helper::ReadAccessor>> propertyAccessor { + property.d_property}; + + std::array elementNodesProperty; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + elementNodesProperty[i] = property.getNodeProperty(context.element[i], propertyAccessor); + } + + return QuadratureContext_t::FiniteElement::Helper::evaluateValueInElement( + elementNodesProperty, context.N); + } +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_BASE_SOURCE_TERM_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API BaseSourceTerm; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.cpp new file mode 100644 index 00000000000..f7a908ed703 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.cpp @@ -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 * +******************************************************************************/ +#define SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP + +#include + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerFEMSourceTermIntegrator(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Consistent nodal load of a source term, integrated with the finite-element quadrature from a nodal source density field") + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + .add< FEMSourceTermIntegrator >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h new file mode 100644 index 00000000000..c401aa202c4 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.h @@ -0,0 +1,166 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP) +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class FEMSourceTermIntegrator + * @brief Integrates a source density into consistent nodal loads. + * + * A source term contributes \f$ \int_{\Omega} N_a \, r \, d\Omega \f$ to the right-hand side, where + * r is the density evaluated by a linked BaseSourceTerm (through l_constantSources) at + * each quadrature point. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + * @tparam TElementType The type of finite element (e.g., sofa::geometry::Tetrahedron). + */ +template +class FEMSourceTermIntegrator : + public sofa::core::behavior::ForceField, + public virtual sofa::core::behavior::TopologyAccessor +{ +public: + using DataTypes = TDataTypes; + using ElementType = TElementType; + SOFA_CLASS2(SOFA_TEMPLATE2(FEMSourceTermIntegrator, DataTypes, ElementType), + sofa::core::behavior::ForceField, + sofa::core::behavior::TopologyAccessor); + +protected: + using FiniteElement = sofa::fem::FiniteElement; + using Real = sofa::Real_t; + + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + +public: + + /** + * @brief Source terms integrated by this component. + * + * If left empty, the BaseSourceTerm components found in the current context are used. + */ + sofa::MultiLink, BaseSourceTerm, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_constantSources; + + /** + * @brief Initializes the component. + * + * This method performs several initialization steps: + * 1. Initializes the base force field. + * 2. Initializes the topology accessor. + * 3. Validates the linked source terms. + * 4. Integrates the source terms into the nodal force. + */ + void init() override; + + /** + * @brief Adds the nodal source term to the RHS vector. + * + * @param mparams Mechanical parameters for the computation. + * @param f The force vector to which the source term will be added. + * @param x The current positions. + * @param v The current velocities. + */ + void addForce( + const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& f, + const sofa::DataVecCoord_t& x, + const sofa::DataVecDeriv_t& v) override; + + /** + * @brief No-op. + */ + void addDForce(const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& df, + const sofa::DataVecDeriv_t& dx) override; + + /** + * @brief No-op. + */ + void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; + + using sofa::core::behavior::ForceField::getPotentialEnergy; + /** + * @brief Potential energy of the nodal load, \f$ V = -\sum_a F_a \cdot (x_a - x_{0,a}) \f$. + */ + SReal getPotentialEnergy(const sofa::core::MechanicalParams* mparams, + const sofa::DataVecCoord_t& x) const override; + + /** + * @brief Degree of the quadrature rule integrating the source terms. + */ + sofa::Data d_quadratureDegree; + +protected: + + /** + * @brief Default constructor. + */ + FEMSourceTermIntegrator(); + + /** + * @brief Ensures that valid source terms are linked, falling back to the current context. + */ + void validateSources(); + + /** + * @brief Runs the quadrature and accumulates every linked source term into m_constantForce. + * + * For each element and each quadrature point, a QuadratureContext is built and handed to every + * source term; the density it returns is weighted by \f$ w \, |\det J| \, N_a \f$ and scattered + * onto the element nodes. + */ + void assembleConstantForce(); + + /** + * @brief Nodal load of every term in l_constantSources. + */ + sofa::VecDeriv_t m_constantForce; +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_SOURCE_TERM_INTEGRATOR_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMSourceTermIntegrator; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl new file mode 100644 index 00000000000..1bc2ba9d2df --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMSourceTermIntegrator.inl @@ -0,0 +1,212 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +FEMSourceTermIntegrator::FEMSourceTermIntegrator() + : l_constantSources(initLink("constantSources", "Source terms of the weak form integrated by " + "this component. If empty, the ones found in the current context are used.")) + , d_quadratureDegree(initData(&d_quadratureDegree, static_cast(1), "quadratureDegree", + "Degree of the quadrature rule integrating the source terms.")) +{ + this->addUpdateCallback("reassembleConstantForce", {&d_quadratureDegree}, + [this](const sofa::core::DataTracker&) + { + if (!this->isComponentStateInvalid() && this->l_topology && this->mstate) + { + assembleConstantForce(); + } + + return this->getComponentState(); + }, {}); +} + +template +void FEMSourceTermIntegrator::init() +{ + sofa::core::behavior::ForceField::init(); + + if (!this->isComponentStateInvalid()) + { + sofa::core::behavior::TopologyAccessor::init(); + } + + if (!this->isComponentStateInvalid()) + { + this->validateSources(); + } + + if (!this->isComponentStateInvalid() && this->l_topology && this->mstate) + { + this->assembleConstantForce(); + } + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + } +} + +template +void FEMSourceTermIntegrator::validateSources() +{ + // Gather all BaseSourceTerm components in Context if empty + if (l_constantSources.empty()) + { + const auto sourcesInContext = this->getContext()->template getObjects >( + sofa::core::objectmodel::BaseContext::Local); + + for (const auto& source : sourcesInContext) + l_constantSources.add(source); + + msg_info_when(!sourcesInContext.empty(), this) << "No source term linked: the " + << sourcesInContext.size() << " one(s) found in the current context are used."; + } + + msg_warning_when(l_constantSources.empty(), this) + << "No source term linked, and none found in the current context '" + << this->getContext()->getName() << "'. This component has zero force contribution."; +} + +template +void FEMSourceTermIntegrator::assembleConstantForce() +{ + m_constantForce.assign(this->mstate->getSize(), sofa::Deriv_t{}); + + const auto restPositionsAccessor = this->mstate->readRestPositions(); + const auto positionsAccessor = this->mstate->readPositions(); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + const auto quadratureRule = FiniteElement::quadratureRule(d_quadratureDegree.getValue()); + + for (const auto& element : elements) + { + const std::array, NumberOfNodesInElement> elementNodesRestCoordinates = + extractNodesVectorFromGlobalVector(element, restPositionsAccessor.ref()); + const std::array, NumberOfNodesInElement> elementNodesCoordinates = + extractNodesVectorFromGlobalVector(element, positionsAccessor.ref()); + + std::array, NumberOfNodesInElement> elementNodesDisplacement; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + elementNodesDisplacement[i] = elementNodesCoordinates[i] - elementNodesRestCoordinates[i]; + } + + for (const auto& [quadraturePoint, weight] : quadratureRule) + { + const auto N = FiniteElement::shapeFunctions(quadraturePoint); + const auto dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + + const auto jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + elementNodesRestCoordinates, dN_dq_ref); + const auto measure = static_cast(sofa::type::absGeneralizedDeterminant(jacobian)); + + const auto restPosition = + FiniteElement::Helper::evaluateValueInElement(elementNodesRestCoordinates, N); + const auto displacement = + FiniteElement::Helper::evaluateValueInElement(elementNodesDisplacement, N); + + const QuadratureContext context{ + element, N, dN_dq_ref, jacobian, measure, restPosition, displacement}; + + const auto weightTimesMeasure = static_cast(weight) * measure; + + for (const auto& source : l_constantSources) + { + const auto density = source->evaluate(context); + + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + m_constantForce[element[i]] += density * (weightTimesMeasure * N[i]); + } + } + } + } +} + +template +void FEMSourceTermIntegrator::addForce(const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& f, + const sofa::DataVecCoord_t& x, + const sofa::DataVecDeriv_t& v) +{ + SOFA_UNUSED(mparams); + SOFA_UNUSED(x); + SOFA_UNUSED(v); + + if (this->isComponentStateInvalid()) + { + return; + } + + auto forceAccessor = sofa::helper::getWriteAccessor(f); + + for (sofa::Index i = 0; i < m_constantForce.size(); ++i) + { + forceAccessor[i] += m_constantForce[i]; + } +} + +template +void FEMSourceTermIntegrator::addDForce(const sofa::core::MechanicalParams* mparams, + sofa::DataVecDeriv_t& df, + const sofa::DataVecDeriv_t& dx) +{ + SOFA_UNUSED(mparams); + SOFA_UNUSED(df); + SOFA_UNUSED(dx); +} + +template +void FEMSourceTermIntegrator::buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) +{ + SOFA_UNUSED(matrix); +} + +template +SReal FEMSourceTermIntegrator::getPotentialEnergy(const sofa::core::MechanicalParams* mparams, + const sofa::DataVecCoord_t& x) const +{ + SOFA_UNUSED(mparams); + + if (this->isComponentStateInvalid()) + { + return 0.0; + } + + const sofa::helper::ReadAccessor positionAccessor = sofa::helper::getReadAccessor(x); + const auto restPositionAccessor = this->mstate->readRestPositions(); + + SReal energy = 0.0; + for (sofa::Index i = 0; i < m_constantForce.size(); ++i) + { + energy -= dot(m_constantForce[i], positionAccessor[i] - restPositionAccessor.ref()[i]); + } + return energy; +} + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.cpp new file mode 100644 index 00000000000..44fe9dc1c25 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.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 SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_PRESSURE_SOURCE_TERM_CPP + +#include + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerNodalPressure(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Definition of a nodal pressure (one scalar per dof).") + .add< NodalPressure >() + .add< NodalPressure >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalPressure; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalPressure; + +void registerPressureSourceTerm(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Traction obtained from a pressure prescribed at the nodes, acting along the normal of the element") + .add< PressureSourceTerm >() + .add< PressureSourceTerm >() + .add< PressureSourceTerm >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API PressureSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API PressureSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API PressureSourceTerm; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h new file mode 100644 index 00000000000..ddae6978a7b --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.h @@ -0,0 +1,115 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_PRESSURE_SOURCE_TERM_CPP) +#include +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class NodalPressure + * @brief A pressure prescribed at the nodes, one scalar per node. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + */ +template +class NodalPressure : public sofa::core::BaseNodalProperty> +{ +public: + using DataTypes = TDataTypes; + using Real = sofa::Real_t; + + SOFA_CLASS(SOFA_TEMPLATE(NodalPressure, DataTypes), + SOFA_TEMPLATE(sofa::core::BaseNodalProperty, sofa::Real_t)); + +protected: + + NodalPressure() : sofa::core::BaseNodalProperty(Real{}) {} +}; + +/** + * @class PressureSourceTerm + * @brief A traction \f$ p \, n \f$ built from a pressure prescribed at the nodes. + * + * The linked NodalPressure is interpolated at the quadrature point and multiplied by the unit + * normal of the element. A positive pressure acts along the normal, whose orientation follows the + * node ordering of the element. + * + * Only available on elements of codimension 1, the ones that have a normal. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + * @tparam TElementType The type of finite element (e.g., sofa::geometry::Triangle). + */ +template +class PressureSourceTerm : public BaseSourceTerm +{ +public: + using DataTypes = TDataTypes; + using ElementType = TElementType; + + SOFA_CLASS(SOFA_TEMPLATE2(PressureSourceTerm, DataTypes, ElementType), + SOFA_TEMPLATE2(BaseSourceTerm, DataTypes, ElementType)); + + using Deriv = sofa::Deriv_t; + using QuadratureContext_t = QuadratureContext; + using NodalPressure = ::sofa::component::solidmechanics::fem::elastic::NodalPressure; + + /** + * @brief Nodal values of the pressure this term integrates. + */ + sofa::SingleLink, NodalPressure, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_pressure; + + /** + * @brief Initializes the component and checks that a pressure is linked. + */ + void init() override; + + /** + * @brief The linked pressure interpolated at the quadrature point, times the unit normal. + */ + Deriv evaluate(const QuadratureContext_t& context) const override; + +protected: + + PressureSourceTerm(); +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_PRESSURE_SOURCE_TERM_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalPressure; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalPressure; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API PressureSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API PressureSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API PressureSourceTerm; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.inl new file mode 100644 index 00000000000..f04fb48a320 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/PressureSourceTerm.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 + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +PressureSourceTerm::PressureSourceTerm() + : l_pressure(initLink("pressure", "Nodal pressure integrated by this term.")) +{ +} + +template +void PressureSourceTerm::init() +{ + BaseSourceTerm::init(); + + if (this->isComponentStateInvalid()) + { + return; + } + + if (!l_pressure) + { + msg_error(this) << "The 'pressure' link must be set to a NodalPressure component. " + "Linked path: '" << l_pressure.getLinkedPath() << "'."; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; + } + + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); +} + +template +sofa::Deriv_t PressureSourceTerm::evaluate( + const QuadratureContext_t& context) const +{ + if (!l_pressure) + { + return Deriv{}; + } + + return elementNormal(context.jacobian) * this->interpolateProperty(*l_pressure, context); +} + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.cpp new file mode 100644 index 00000000000..f8d63ab23a9 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.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 SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_STRESS_SOURCE_TERM_CPP + +#include + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerNodalStress(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Definition of a nodal symmetric stress tensor (one tensor per dof).") + .add< NodalStress >() + .add< NodalStress >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalStress; + +void registerStressSourceTerm(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Traction obtained from a symmetric stress tensor prescribed at the nodes, contracted with the normal of the element") + .add< StressSourceTerm >() + .add< StressSourceTerm >() + .add< StressSourceTerm >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API StressSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API StressSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API StressSourceTerm; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.h new file mode 100644 index 00000000000..d967ea2c92f --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.h @@ -0,0 +1,125 @@ +/****************************************************************************** +* 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(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_STRESS_SOURCE_TERM_CPP) +#include +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/// The symmetric stress tensor prescribed at one node. +template +using StressTensor = + sofa::type::MatSym>; + +/** + * @class NodalStress + * @brief A symmetric stress tensor prescribed at the nodes, one tensor per node. + * + * A tensor is read from a scene as a full matrix, row by row: nine components in 3D, four in 2D, + * of which MatSym keeps the independent ones. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + */ +template +class NodalStress : public sofa::core::BaseNodalProperty> +{ +public: + using DataTypes = TDataTypes; + using Tensor = StressTensor; + + SOFA_CLASS(SOFA_TEMPLATE(NodalStress, DataTypes), + SOFA_TEMPLATE(sofa::core::BaseNodalProperty, StressTensor)); + +protected: + + NodalStress() : sofa::core::BaseNodalProperty(Tensor{}) {} +}; + +/** + * @class StressSourceTerm + * @brief A traction \f$ \sigma \, n \f$ built from a stress tensor prescribed at the nodes. + * + * The linked NodalStress is interpolated at the quadrature point and contracted with the unit + * normal of the element, whose orientation follows the node ordering. + * + * The tensor must be symmetric, which excludes the first Piola-Kirchhoff stress. + * + * Only available on elements of codimension 1, the ones that have a normal. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + * @tparam TElementType The type of finite element (e.g., sofa::geometry::Triangle). + */ +template +class StressSourceTerm : public BaseSourceTerm +{ +public: + using DataTypes = TDataTypes; + using ElementType = TElementType; + + SOFA_CLASS(SOFA_TEMPLATE2(StressSourceTerm, DataTypes, ElementType), + SOFA_TEMPLATE2(BaseSourceTerm, DataTypes, ElementType)); + + using Deriv = sofa::Deriv_t; + using QuadratureContext_t = QuadratureContext; + using NodalStress = ::sofa::component::solidmechanics::fem::elastic::NodalStress; + + /** + * @brief Nodal values of the stress tensor this term integrates. + */ + sofa::SingleLink, NodalStress, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_stress; + + /** + * @brief Initializes the component and checks that a stress is linked. + */ + void init() override; + + /** + * @brief The linked stress interpolated at the quadrature point, contracted with the normal. + */ + Deriv evaluate(const QuadratureContext_t& context) const override; + +protected: + + StressSourceTerm(); +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_STRESS_SOURCE_TERM_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalStress; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalStress; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API StressSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API StressSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API StressSourceTerm; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.inl new file mode 100644 index 00000000000..5fbb13d465c --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.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 + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +StressSourceTerm::StressSourceTerm() + : l_stress(initLink("stress", "Nodal stress tensor integrated by this term.")) +{ +} + +template +void StressSourceTerm::init() +{ + BaseSourceTerm::init(); + + if (this->isComponentStateInvalid()) + { + return; + } + + if (!l_stress) + { + msg_error(this) << "The 'stress' link must be set to a NodalStress component. " + "Linked path: '" << l_stress.getLinkedPath() << "'."; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; + } + + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); +} + +template +sofa::Deriv_t StressSourceTerm::evaluate( + const QuadratureContext_t& context) const +{ + if (!l_stress) + { + return Deriv{}; + } + + return this->interpolateProperty(*l_stress, context) * elementNormal(context.jacobian); +} + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.cpp new file mode 100644 index 00000000000..3c6aa474337 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.cpp @@ -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 * +******************************************************************************/ +#define SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_VECTOR_SOURCE_TERM_CPP + +#include + +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +void registerNodalSourceDensity(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Definition of a nodal source density (one vector per dof).") + .add< NodalSourceDensity >() + .add< NodalSourceDensity >() + .add< NodalSourceDensity >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalSourceDensity; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalSourceDensity; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalSourceDensity; + +void registerVectorSourceTerm(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Source density given as a vector at each node, per unit measure of the element") + .add< VectorSourceTerm >() + .add< VectorSourceTerm >() + .add< VectorSourceTerm >() + .add< VectorSourceTerm >() + .add< VectorSourceTerm >() + .add< VectorSourceTerm >() + .add< VectorSourceTerm >() + .add< VectorSourceTerm >() + .add< VectorSourceTerm >() + ); +} + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h new file mode 100644 index 00000000000..7fb83bcf25a --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.h @@ -0,0 +1,121 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_VECTOR_SOURCE_TERM_CPP) +#include +#include +#endif + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @class NodalSourceDensity + * @brief A source density prescribed at the nodes, one vector per node. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + */ +template +class NodalSourceDensity : public sofa::core::BaseNodalProperty> +{ +public: + using DataTypes = TDataTypes; + using Deriv = sofa::Deriv_t; + + SOFA_CLASS(SOFA_TEMPLATE(NodalSourceDensity, DataTypes), + SOFA_TEMPLATE(sofa::core::BaseNodalProperty, sofa::Deriv_t)); + +protected: + + NodalSourceDensity() : sofa::core::BaseNodalProperty(Deriv{}) {} +}; + +/** + * @class VectorSourceTerm + * @brief A source density given directly as a vector at each node. + * + * The linked NodalSourceDensity is interpolated at the quadrature point and returned unchanged. + * The measure the integrator divides it by is the one of the element it is attached to: a body + * force on a volume element, a traction on a boundary one. + * + * @tparam TDataTypes The data types used for positions, velocities, etc. (e.g., Vec3Types). + * @tparam TElementType The type of finite element (e.g., sofa::geometry::Tetrahedron). + */ +template +class VectorSourceTerm : public BaseSourceTerm +{ +public: + using DataTypes = TDataTypes; + using ElementType = TElementType; + + SOFA_CLASS(SOFA_TEMPLATE2(VectorSourceTerm, DataTypes, ElementType), + SOFA_TEMPLATE2(BaseSourceTerm, DataTypes, ElementType)); + + using Deriv = sofa::Deriv_t; + using QuadratureContext_t = QuadratureContext; + using NodalSourceDensity = + ::sofa::component::solidmechanics::fem::elastic::NodalSourceDensity; + + /** + * @brief Nodal values of the source density this term integrates. + */ + sofa::SingleLink, NodalSourceDensity, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_sourceDensity; + + /** + * @brief Initializes the component and checks that a source density is linked. + */ + void init() override; + + /** + * @brief The linked source density interpolated at the quadrature point. + */ + Deriv evaluate(const QuadratureContext_t& context) const override; + +protected: + + VectorSourceTerm(); +}; + +#if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_VECTOR_SOURCE_TERM_CPP) +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalSourceDensity; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalSourceDensity; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API NodalSourceDensity; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VectorSourceTerm; +#endif + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.inl new file mode 100644 index 00000000000..3acbcec062e --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VectorSourceTerm.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 + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +VectorSourceTerm::VectorSourceTerm() + : l_sourceDensity(initLink("sourceDensity", "Nodal source density integrated by this term.")) +{ +} + +template +void VectorSourceTerm::init() +{ + BaseSourceTerm::init(); + + if (this->isComponentStateInvalid()) + { + return; + } + + if (!l_sourceDensity) + { + msg_error(this) << "The 'sourceDensity' link must be set to a NodalSourceDensity " + "component. Linked path: '" << l_sourceDensity.getLinkedPath() << "'."; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; + } + + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); +} + +template +sofa::Deriv_t VectorSourceTerm::evaluate( + const QuadratureContext_t& context) const +{ + if (!l_sourceDensity) + { + return Deriv{}; + } + + return this->interpolateProperty(*l_sourceDensity, context); +} + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp index 87d14c740ff..23056027c0c 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/init.cpp @@ -28,19 +28,26 @@ namespace sofa::component::solidmechanics::fem::elastic extern void registerBeamFEMForceField(sofa::core::ObjectFactory* factory); extern void registerCorotationalFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerFEMSourceTermIntegrator(sofa::core::ObjectFactory* factory); extern void registerFastTetrahedralCorotationalForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedralFEMForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedralFEMForceFieldAndMass(sofa::core::ObjectFactory* factory); extern void registerHexahedronFEMForceField(sofa::core::ObjectFactory* factory); extern void registerHexahedronFEMForceFieldAndMass(sofa::core::ObjectFactory* factory); extern void registerLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerNodalPressure(sofa::core::ObjectFactory* factory); +extern void registerNodalSourceDensity(sofa::core::ObjectFactory* factory); +extern void registerNodalStress(sofa::core::ObjectFactory* factory); +extern void registerPressureSourceTerm(sofa::core::ObjectFactory* factory); extern void registerQuadBendingFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerStressSourceTerm(sofa::core::ObjectFactory* factory); extern void registerTetrahedralCorotationalFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTetrahedronFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangleFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangularAnisotropicFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangularFEMForceField(sofa::core::ObjectFactory* factory); extern void registerTriangularFEMForceFieldOptim(sofa::core::ObjectFactory* factory); +extern void registerVectorSourceTerm(sofa::core::ObjectFactory* factory); extern void registerVonMisesStress(sofa::core::ObjectFactory* factory); extern "C" { @@ -69,19 +76,26 @@ void registerObjects(sofa::core::ObjectFactory* factory) { registerBeamFEMForceField(factory); registerCorotationalFEMForceField(factory); + registerFEMSourceTermIntegrator(factory); registerFastTetrahedralCorotationalForceField(factory); registerHexahedralFEMForceField(factory); registerHexahedralFEMForceFieldAndMass(factory); registerHexahedronFEMForceField(factory); registerHexahedronFEMForceFieldAndMass(factory); registerLinearSmallStrainFEMForceField(factory); + registerNodalPressure(factory); + registerNodalSourceDensity(factory); + registerNodalStress(factory); + registerPressureSourceTerm(factory); registerQuadBendingFEMForceField(factory); + registerStressSourceTerm(factory); registerTetrahedralCorotationalFEMForceField(factory); registerTetrahedronFEMForceField(factory); registerTriangleFEMForceField(factory); registerTriangularAnisotropicFEMForceField(factory); registerTriangularFEMForceField(factory); registerTriangularFEMForceFieldOptim(factory); + registerVectorSourceTerm(factory); registerVonMisesStress(factory); } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt index d5e6a1bcb25..b8a38e68eb4 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt @@ -6,6 +6,7 @@ set(SOURCE_FILES BaseTetrahedronFEMForceField_test.h BeamFEMForceField_test.cpp FastTetrahedralCorotationalForceField_test.cpp + FEMSourceTermIntegrator_test.cpp HexahedronFEMForceField_test.cpp StrainDisplacement_test.cpp TetrahedralCorotationalFEMForceField_test.cpp diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp new file mode 100644 index 00000000000..2e76cacf63a --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/FEMSourceTermIntegrator_test.cpp @@ -0,0 +1,152 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace sofa +{ +using sofa::simpleapi::createObject; +using sofa::simpleapi::createRootNode; +using DataTypes = defaulttype::Vec2Types; +using Integrator = component::solidmechanics::fem::elastic::FEMSourceTermIntegrator; +using DOF = component::statecontainer::MechanicalObject; +using VecCoord = DataTypes::VecCoord; +using VecDeriv = DataTypes::VecDeriv; + +class FEMSourceTermIntegrator_test : public testing::BaseTest +{ +protected: + simulation::Simulation* m_simulation = nullptr; + simulation::Node::SPtr m_root; + + void doSetUp() override + { + m_simulation = sofa::simulation::getSimulation(); + } + + void doTearDown() override + { + if (m_root != nullptr) + sofa::simulation::node::unload(m_root); + } + + simulation::Node::SPtr makeMesh() + { + this->loadPlugins({"Sofa.Component.StateContainer", + "Sofa.Component.Topology.Container.Constant", "Sofa.Component.SolidMechanics.FEM.Elastic"}); + + auto root = createRootNode(m_simulation, "root"); + createObject(root, "MechanicalObject", {{"template", "Vec2"}, {"position", "0 0 1 0 1 1 0 1"}}); + createObject(root, "MeshTopology", {{"name", "mesh"}, {"triangles", "0 1 2 0 2 3"}}); + return root; + } + + static VecDeriv addForce(Integrator* integrator) + { + core::MechanicalParams mparams; + Data f; + f.setValue(VecDeriv(4)); + Data x; + Data v; + integrator->addForce(&mparams, f, x, v); + return f.getValue(); + } +}; + +// Splitting one VectorSourceTerm into several must not change the integrated force. +TEST_F(FEMSourceTermIntegrator_test, MultipleSourcesSumToOne) +{ + m_root = makeMesh(); + + createObject(m_root, "NodalSourceDensity", {{"name", "fullDensity"}, {"template", "Vec2"}, {"property", "300 -600"}}); + createObject(m_root, "VectorSourceTerm", {{"name", "full"}, {"template", "Vec2,Triangle"}, {"sourceDensity", "@fullDensity"}}); + auto* one = dynamic_cast(createObject(m_root, "FEMSourceTermIntegrator", + {{"name", "one"}, {"template", "Vec2,Triangle"}, {"topology", "@mesh"}, {"constantSources", "@full"}}).get()); + + createObject(m_root, "NodalSourceDensity", {{"name", "thirdDensity"}, {"template", "Vec2"}, {"property", "100 -200"}}); + createObject(m_root, "VectorSourceTerm", {{"name", "a"}, {"template", "Vec2,Triangle"}, {"sourceDensity", "@thirdDensity"}}); + createObject(m_root, "VectorSourceTerm", {{"name", "b"}, {"template", "Vec2,Triangle"}, {"sourceDensity", "@thirdDensity"}}); + createObject(m_root, "VectorSourceTerm", {{"name", "c"}, {"template", "Vec2,Triangle"}, {"sourceDensity", "@thirdDensity"}}); + auto* three = dynamic_cast(createObject(m_root, "FEMSourceTermIntegrator", + {{"name", "three"}, {"template", "Vec2,Triangle"}, {"topology", "@mesh"}, {"constantSources", "@a @b @c"}}).get()); + + simulation::node::initRoot(m_root.get()); + ASSERT_NE(one, nullptr); + ASSERT_NE(three, nullptr); + + const VecDeriv fOne = addForce(one); + const VecDeriv fThree = addForce(three); + ASSERT_EQ(fOne.size(), fThree.size()); + for (std::size_t i = 0; i < fOne.size(); ++i) + EXPECT_EQ(fOne[i], fThree[i]); +} + +// dE = -dX . F, the conservative-force identity getPotentialEnergy relies on. +TEST_F(FEMSourceTermIntegrator_test, PotentialEnergyMatchesWork) +{ + m_root = makeMesh(); + + createObject(m_root, "NodalSourceDensity", {{"name", "bodyForceDensity"}, {"template", "Vec2"}, {"property", "300 -600"}}); + createObject(m_root, "VectorSourceTerm", {{"name", "bodyForce"}, {"template", "Vec2,Triangle"}, {"sourceDensity", "@bodyForceDensity"}}); + auto* integrator = dynamic_cast(createObject(m_root, "FEMSourceTermIntegrator", + {{"name", "source"}, {"template", "Vec2,Triangle"}, {"topology", "@mesh"}, {"constantSources", "@bodyForce"}}).get()); + + simulation::node::initRoot(m_root.get()); + ASSERT_NE(integrator, nullptr); + + const VecDeriv f = addForce(integrator); + + VecCoord x0(4); + testing::copyFromData(x0, m_root->get()->readPositions()); + + VecCoord x1 = x0; + for (auto& xi : x1) + xi += DataTypes::Coord(0.01, -0.02); + + core::MechanicalParams mparams; + Data x0Data; + x0Data.setValue(x0); + Data x1Data; + x1Data.setValue(x1); + + const SReal e0 = integrator->getPotentialEnergy(&mparams, x0Data); + const SReal e1 = integrator->getPotentialEnergy(&mparams, x1Data); + + SReal work = 0; + for (std::size_t i = 0; i < f.size(); ++i) + work += dot(f[i], x1[i] - x0[i]); + + EXPECT_NEAR(e1 - e0, -work, 1e-9); +} + +} // namespace sofa diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h index a490fddd6a3..9c389914f9b 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h @@ -86,7 +86,7 @@ struct FiniteElementHelper const std::array& valuesAtNodes, const sofa::type::Vec& shapeFunctions) { - return std::inner_product(valuesAtNodes.begin(), valuesAtNodes.end(), shapeFunctions.begin(), T(0)); + return std::inner_product(valuesAtNodes.begin(), valuesAtNodes.end(), shapeFunctions.begin(), T{}); } }; diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h index fbffb5dbfed..53570c30d05 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_EDGE_CPP) #include @@ -55,12 +57,43 @@ struct FiniteElement return {{-static_cast(0.5)}, {static_cast(0.5)}}; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr sofa::type::Vec q0(static_cast(0)); - return { - std::make_pair(q0, static_cast(2)) - }; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point midpoint rule (default). + return std::array{ + std::make_pair(ReferenceCoord(static_cast(0)), static_cast(2)) + }; + } + else if constexpr (Degree <= 3) + { + // Degrees 2-3: 2-point Gauss-Legendre rule. + constexpr Real sqrt3 = 1.73205080757; + constexpr Real g = static_cast(1) / sqrt3; + return std::array{ + std::make_pair(ReferenceCoord(-g), static_cast(1)), + std::make_pair(ReferenceCoord( g), static_cast(1)) + }; + } + else + { + static_assert(Degree <= 3, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(sofa::Size degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: + case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h index 38339624e47..d855c9f8c72 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_HEXAHEDRON_CPP) #include @@ -96,24 +98,68 @@ struct FiniteElement return gradient; } + template static constexpr auto quadraturePoints() { - constexpr Real sqrt3 = 1.73205080757; //sqrt(3.) - constexpr Real sqrt3_1 = static_cast(1) / sqrt3; - constexpr Real one = static_cast(1); + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule. + return std::array{ + std::make_pair(ReferenceCoord(static_cast(0), static_cast(0), static_cast(0)), static_cast(8)) + }; + } + else if constexpr (Degree <= 3) + { + // Degrees 2-3: 2x2x2 Gauss-Legendre rule (default). + constexpr Real sqrt3 = 1.73205080757; //sqrt(3.) + constexpr Real sqrt3_1 = static_cast(1) / sqrt3; + constexpr Real one = static_cast(1); - constexpr std::array q { - std::pair{referenceElementNodes[0] * sqrt3_1, one}, - std::pair{referenceElementNodes[1] * sqrt3_1, one}, - std::pair{referenceElementNodes[2] * sqrt3_1, one}, - std::pair{referenceElementNodes[3] * sqrt3_1, one}, - std::pair{referenceElementNodes[4] * sqrt3_1, one}, - std::pair{referenceElementNodes[5] * sqrt3_1, one}, - std::pair{referenceElementNodes[6] * sqrt3_1, one}, - std::pair{referenceElementNodes[7] * sqrt3_1, one}, - }; + return std::array { + std::pair{referenceElementNodes[0] * sqrt3_1, one}, + std::pair{referenceElementNodes[1] * sqrt3_1, one}, + std::pair{referenceElementNodes[2] * sqrt3_1, one}, + std::pair{referenceElementNodes[3] * sqrt3_1, one}, + std::pair{referenceElementNodes[4] * sqrt3_1, one}, + std::pair{referenceElementNodes[5] * sqrt3_1, one}, + std::pair{referenceElementNodes[6] * sqrt3_1, one}, + std::pair{referenceElementNodes[7] * sqrt3_1, one}, + }; + } + else if constexpr (Degree <= 5) + { + // Degrees 4-5: 3x3x3 Gauss-Legendre rule. + constexpr Real g = 0.77459666924; //sqrt(3./5.) + constexpr std::array node{ -g, static_cast(0), g }; + constexpr std::array weight{ static_cast(5./9.), static_cast(8./9.), static_cast(5./9.) }; + + std::array q{}; + sofa::Size k = 0; + for (sofa::Size i = 0; i < 3; ++i) + for (sofa::Size j = 0; j < 3; ++j) + for (sofa::Size l = 0; l < 3; ++l) + q[k++] = std::make_pair(ReferenceCoord(node[i], node[j], node[l]), weight[i] * weight[j] * weight[l]); + return q; + } + else + { + static_assert(Degree <= 5, "FiniteElement: no quadrature rule for the requested degree"); + } + } - return q; + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(sofa::Size degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: + case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; } + case 4: + case 5: { static constexpr auto rule = quadraturePoints<5>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h index 4e37585eace..261eeae7036 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_QUAD_CPP) #include @@ -67,21 +69,62 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.) - constexpr Real sqrt6 = 2.44948974278; //sqrt(6.) - constexpr Real sqrt2 = 1.41421356237; //sqrt(2.) + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule. + return std::array{ + std::make_pair(ReferenceCoord(static_cast(0), static_cast(0)), static_cast(4)) + }; + } + else if constexpr (Degree <= 2) + { + // Degree 2: 3-point rule (default). + constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.) + constexpr Real sqrt6 = 2.44948974278; //sqrt(6.) + constexpr Real sqrt2 = 1.41421356237; //sqrt(2.) - constexpr sofa::type::Vec q0(sqrt2_3, 0.); - constexpr sofa::type::Vec q1(-1/sqrt6, -1./sqrt2); - constexpr sofa::type::Vec q2(-1/sqrt6, 1./sqrt2); + constexpr ReferenceCoord q0(sqrt2_3, 0.); + constexpr ReferenceCoord q1(-1/sqrt6, -1./sqrt2); + constexpr ReferenceCoord q2(-1/sqrt6, 1./sqrt2); - return { - std::make_pair(q0, 4./3.), - std::make_pair(q1, 4./3.), - std::make_pair(q2, 4./3.), - }; + return std::array{ + std::make_pair(q0, 4./3.), + std::make_pair(q1, 4./3.), + std::make_pair(q2, 4./3.) + }; + } + else if constexpr (Degree <= 3) + { + // Degree 3: 2x2 Gauss-Legendre rule. + constexpr Real sqrt3 = 1.73205080757; //sqrt(3.) + constexpr Real g = static_cast(1) / sqrt3; + return std::array{ + std::make_pair(ReferenceCoord(-g, -g), static_cast(1)), + std::make_pair(ReferenceCoord( g, -g), static_cast(1)), + std::make_pair(ReferenceCoord( g, g), static_cast(1)), + std::make_pair(ReferenceCoord(-g, g), static_cast(1)) + }; + } + else + { + static_assert(Degree <= 3, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(sofa::Size degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; } + case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h index ac14bb9c59b..bca61671e0a 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_TETAHEDRON_CPP) #include @@ -68,11 +70,45 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr sofa::type::Vec q0(1./4., 1./4., 1./4.); - constexpr std::array q { std::make_pair(q0, 1./6.) }; - return q; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule (default). + return std::array{ + std::make_pair(ReferenceCoord(1./4., 1./4., 1./4.), Real(1./6.)) + }; + } + else if constexpr (Degree <= 2) + { + // Degree 2: 4-point rule. + constexpr Real sqrt5 = 2.2360679774997896; + constexpr Real a = (5. - sqrt5) / 20.; + constexpr Real b = (5. + 3. * sqrt5) / 20.; + return std::array{ + std::make_pair(ReferenceCoord(a, a, a), Real(1./24.)), + std::make_pair(ReferenceCoord(b, a, a), Real(1./24.)), + std::make_pair(ReferenceCoord(a, b, a), Real(1./24.)), + std::make_pair(ReferenceCoord(a, a, b), Real(1./24.)) + }; + } + else + { + static_assert(Degree <= 2, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(sofa::Size degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h index 0d9659040de..84c9ec19a74 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include namespace sofa::fem { @@ -64,11 +66,41 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - return { - std::make_pair(sofa::type::Vec(1./3., 1./3.), 1./2.) - }; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule (default). + return std::array{ + std::make_pair(ReferenceCoord(1./3., 1./3.), Real(1./2.)) + }; + } + else if constexpr (Degree <= 2) + { + // Degree 2: 3-point interior rule. + return std::array{ + std::make_pair(ReferenceCoord(1./6., 1./6.), Real(1./6.)), + std::make_pair(ReferenceCoord(2./3., 1./6.), Real(1./6.)), + std::make_pair(ReferenceCoord(1./6., 2./3.), Real(1./6.)) + }; + } + else + { + static_assert(Degree <= 2, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(sofa::Size degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; diff --git a/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn b/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn new file mode 100644 index 00000000000..138e29f88f0 --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn @@ -0,0 +1,113 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Component/SolidMechanics/FEM/VectorSourceTerm.scn b/examples/Component/SolidMechanics/FEM/VectorSourceTerm.scn new file mode 100644 index 00000000000..daddebd6e4d --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/VectorSourceTerm.scn @@ -0,0 +1,292 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +