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 index 457b3cec534..63a516b3b31 100644 --- 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 @@ -44,10 +44,7 @@ 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. + * @brief Geometry of one 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). @@ -72,36 +69,32 @@ struct QuadratureContext 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. + /// Node indices of the element. const Element& element; - /// Shape function value at this quadrature point. + /// Shape function values. ShapeFunctions N; - /// Reference-space gradients of the shape functions at this quadrature point. + /// Reference-space gradients of the shape functions. GradientShapeFunctions gradientShapeFunctions; - /// dx/dq of the reference-to-physical mapping, on the configuration the integrator chose. + /// dx/dq of the reference-to-physical mapping. Jacobian jacobian; - /// \f$ |\det J| \f$, for information only: the integrator applies it, a term must not. + /// \f$ |\det J| \f$ Real measure; - /// Interpolated rest position at this quadrature point. + /// Interpolated rest position. Coord restPosition; - /// Interpolated displacement at this quadrature point. + /// Interpolated displacement. 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. + * @brief Unit normal of an element of codimension 1, oriented by its node ordering. * - * @param jacobian dx/dq of the reference-to-physical mapping at the point of interest. + * @param jacobian dx/dq of the reference-to-physical mapping. */ template sofa::type::Vec elementNormal( @@ -126,10 +119,55 @@ sofa::type::Vec elementNormal( } /** - * @class BaseSourceTerm - * @brief A source density whose value is determined by the geometry, not by the solution. + * @brief Derivative of the area normal of an element with respect to a node position: + * + * \f[ \frac{\partial (|\det J| \, n)}{\partial x} + * = \sum_i \frac{\partial (|\det J| \, n)}{\partial t_i} \frac{\partial N}{\partial q_i}, \f] + * + * with \f$ t_i \f$ the columns of J and N the shape function of the node. * - * The component calculates the integrand in evaluate() given a QuadratureContext. + * @param jacobian dx/dq of the reference-to-physical mapping. + * @param gradientOfShapeFunction dN/dq of the node. + */ +template +sofa::type::Mat elementAreaNormalDerivative( + const sofa::type::Mat& jacobian, + const sofa::type::Vec& gradientOfShapeFunction) +{ + static_assert(TopologicalDimension + 1 == spatial_dimensions, + "An area normal is only defined for an element of codimension 1."); + + using Derivative = sofa::type::Mat; + + if constexpr (spatial_dimensions == 3) + { + const auto tangent0 = jacobian.col(0); + const auto tangent1 = jacobian.col(1); + + const Derivative skewTangent0 { + {Real{0}, -tangent0[2], tangent0[1]}, + {tangent0[2], Real{0}, -tangent0[0]}, + {-tangent0[1], tangent0[0], Real{0}}}; + + const Derivative skewTangent1 { + {Real{0}, -tangent1[2], tangent1[1]}, + {tangent1[2], Real{0}, -tangent1[0]}, + {-tangent1[1], tangent1[0], Real{0}}}; + + return skewTangent0 * gradientOfShapeFunction[1] + - skewTangent1 * gradientOfShapeFunction[0]; + } + else + { + const Derivative rotation {{Real{0}, Real{1}}, {Real{-1}, Real{0}}}; + + return rotation * gradientOfShapeFunction[0]; + } +} + +/** + * @class BaseSourceTerm + * @brief A source term of the weak form. * * @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). @@ -144,9 +182,15 @@ class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent SOFA_CLASS(SOFA_TEMPLATE2(BaseSourceTerm, DataTypes, ElementType), sofa::core::objectmodel::BaseComponent); + using Real = sofa::Real_t; using Deriv = sofa::Deriv_t; using QuadratureContext_t = QuadratureContext; + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + + /// Derivative of a nodal force with respect to a node position. + using SourceDerivative = sofa::type::Mat; + /** * @brief Source density at one quadrature point, per unit physical measure. * @@ -155,27 +199,31 @@ class BaseSourceTerm : public sofa::core::objectmodel::BaseComponent virtual Deriv evaluate(const QuadratureContext_t& context) const = 0; /** - * @brief The Data this term's integrand reads, beyond the geometry it is handed. - * - * Access to this Data is enabled so that a callback on FEMSourceTermIntegrator can be activated. - * They are tracked and the cached integral is rebuilt whenever one of them is written. + * @brief Data the integrand depends on; used for callbacks. */ virtual sofa::type::vector integrandInputs() const { return {}; } + /** + * @brief Derivative of the measure-weighted density with respect to a node position, + * \f$ \partial (|\det J| \, r) / \partial x_b \f$, with r the source density. + * + * @param context Geometry of the quadrature point. + * @param node Index of the node in the element. + */ + virtual SourceDerivative evaluateStiffness(const QuadratureContext_t& context, + sofa::Size node) 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. + * @brief Nodal property interpolated at the quadrature point. * - * @param property The component holding the nodal values. + * @param property The nodal values. * @param context Geometry of the quadrature point. */ template 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 index 6f230e46aa5..dcee7fc3096 100644 --- 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 @@ -37,11 +37,9 @@ 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. + * @brief Integrates source terms into consistent nodal loads, + * \f$ \int_{\Omega} N_a \, r \, d\Omega \f$, with \f$ N_a \f$ the shape function of node a and r + * the source density. * * @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). @@ -68,30 +66,27 @@ class FEMSourceTermIntegrator : public: /** - * @brief Source terms integrated by this component. - * - * If left empty, the BaseSourceTerm components found in the current context are used. + * @brief Source terms integrated once on the rest configuration. */ sofa::MultiLink, BaseSourceTerm, sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_constantSources; + /** + * @brief Source terms integrated on the current configuration. + */ + sofa::MultiLink, BaseSourceTerm, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_nonConstantSources; + /** * @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. Tracks the Data the constant source terms read. - * 5. Integrates the source terms into the nodal force. */ void init() override; /** - * @brief Adds the nodal source term to the RHS vector. + * @brief Adds the nodal loads of the source terms. * * @param mparams Mechanical parameters for the computation. - * @param f The force vector to which the source term will be added. + * @param f The force vector, to which the source term will be added. * @param x The current positions. * @param v The current velocities. */ @@ -102,14 +97,14 @@ class FEMSourceTermIntegrator : const sofa::DataVecDeriv_t& v) override; /** - * @brief No-op. + * @brief Applies the tangent of the non-constant source terms to a displacement increment. */ void addDForce(const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& df, const sofa::DataVecDeriv_t& dx) override; /** - * @brief No-op. + * @brief Assembles the tangent of the non-constant source terms. */ void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; @@ -125,11 +120,13 @@ class FEMSourceTermIntegrator : */ sofa::Data d_quadratureDegree; -protected: - /** - * @brief Default constructor. + * @brief Whether the non-symmetric tangent of the non-constant source terms is assembled. */ + sofa::Data d_useTangentStiffness; + +protected: + FEMSourceTermIntegrator(); /** @@ -143,16 +140,21 @@ class FEMSourceTermIntegrator : void doUpdateInternal() override; /** - * @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. + * @brief Integrates the constant source terms on the rest configuration. */ void assembleConstantForce(); /** - * @brief Nodal load of every term in l_constantSources. + * @brief Integrates the non-constant source terms. + * + * @param x The configuration the terms are integrated on. + * @param f The force vector the nodal forces are accumulated into. + */ + void assembleNonConstantForce(const sofa::VecCoord_t& x, + sofa::VecDeriv_t& f) const; + + /** + * @brief Nodal load of the constant source terms. */ sofa::VecDeriv_t m_constantForce; }; 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 index 264bea419e6..a25cd6876b5 100644 --- 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 @@ -22,16 +22,22 @@ #pragma once #include #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.")) + : l_constantSources(initLink("constantSources", "Source terms of the weak form integrated once " + "on the rest configuration.")) + , l_nonConstantSources(initLink("nonConstantSources", "Source terms of the weak form " + "integrated at every step on the current configuration. If both source links are " + "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.")) + , d_useTangentStiffness(initData(&d_useTangentStiffness, true, "useTangentStiffness", + "Whether the (non-symmetric) tangent of the non-constant source terms is assembled.")) { this->addUpdateCallback("reassembleConstantForce", {&d_quadratureDegree}, [this](const sofa::core::DataTracker&) @@ -94,20 +100,21 @@ void FEMSourceTermIntegrator::doUpdateInternal() template void FEMSourceTermIntegrator::validateSources() { - // Gather all BaseSourceTerm components in Context if empty - if (l_constantSources.empty()) + // Gather all BaseSourceTerm components in Context if both links are empty + if (l_constantSources.empty() && l_nonConstantSources.empty()) { const auto sourcesInContext = this->getContext()->template getObjects >( sofa::core::objectmodel::BaseContext::Local); for (const auto& source : sourcesInContext) - l_constantSources.add(source); + l_nonConstantSources.add(source); msg_info_when(!sourcesInContext.empty(), this) << "No source term linked: the " - << sourcesInContext.size() << " one(s) found in the current context are used."; + << sourcesInContext.size() << " one(s) found in the current context are integrated on " + "the current configuration."; } - msg_warning_when(l_constantSources.empty(), this) + msg_warning_when(l_constantSources.empty() && l_nonConstantSources.empty(), this) << "No source term linked, and none found in the current context '" << this->getContext()->getName() << "'. This component has zero force contribution."; } @@ -168,6 +175,65 @@ void FEMSourceTermIntegrator::assembleConstantForce() } } +template +void FEMSourceTermIntegrator::assembleNonConstantForce( + const sofa::VecCoord_t& x, sofa::VecDeriv_t& f) const +{ + if (l_nonConstantSources.empty()) + { + return; + } + + const auto restPositionsAccessor = this->mstate->readRestPositions(); + + 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, x); + + 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( + elementNodesCoordinates, 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_nonConstantSources) + { + const auto density = source->evaluate(context); + + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + f[element[i]] += density * (weightTimesMeasure * N[i]); + } + } + } + } +} + template void FEMSourceTermIntegrator::addForce(const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& f, @@ -175,7 +241,6 @@ void FEMSourceTermIntegrator::addForce(const sofa::core: const sofa::DataVecDeriv_t& v) { SOFA_UNUSED(mparams); - SOFA_UNUSED(x); SOFA_UNUSED(v); if (this->isComponentStateInvalid()) @@ -189,6 +254,9 @@ void FEMSourceTermIntegrator::addForce(const sofa::core: { forceAccessor[i] += m_constantForce[i]; } + + const sofa::helper::ReadAccessor positionAccessor = sofa::helper::getReadAccessor(x); + assembleNonConstantForce(positionAccessor.ref(), forceAccessor.wref()); } template @@ -196,15 +264,143 @@ void FEMSourceTermIntegrator::addDForce(const sofa::core sofa::DataVecDeriv_t& df, const sofa::DataVecDeriv_t& dx) { - SOFA_UNUSED(mparams); - SOFA_UNUSED(df); - SOFA_UNUSED(dx); + if (this->isComponentStateInvalid() || l_nonConstantSources.empty() + || !d_useTangentStiffness.getValue()) + { + return; + } + + const auto kFactor = static_cast(sofa::core::mechanicalparams::kFactorIncludingRayleighDamping( + mparams, this->rayleighStiffness.getValue())); + + auto forceDerivativeAccessor = sofa::helper::getWriteAccessor(df); + const sofa::helper::ReadAccessor displacementDerivativeAccessor = sofa::helper::getReadAccessor(dx); + + 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( + elementNodesCoordinates, 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}; + + sofa::Deriv_t contraction{}; + for (const auto& source : l_nonConstantSources) + { + for (sofa::Size b = 0; b < NumberOfNodesInElement; ++b) + { + contraction += source->evaluateStiffness(context, b) + * displacementDerivativeAccessor[element[b]]; + } + } + + for (sofa::Size a = 0; a < NumberOfNodesInElement; ++a) + { + forceDerivativeAccessor[element[a]] += + contraction * (kFactor * static_cast(weight) * N[a]); + } + } + } } template void FEMSourceTermIntegrator::buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) { - SOFA_UNUSED(matrix); + if (this->isComponentStateInvalid() || l_nonConstantSources.empty() + || !d_useTangentStiffness.getValue()) + { + return; + } + + auto dfdx = matrix->getForceDerivativeIn(this->mstate).withRespectToPositionsIn(this->mstate); + + 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( + elementNodesCoordinates, 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}; + + std::array::SourceDerivative, + NumberOfNodesInElement> sourceDerivative{}; + + for (const auto& source : l_nonConstantSources) + { + for (sofa::Size b = 0; b < NumberOfNodesInElement; ++b) + { + sourceDerivative[b] += source->evaluateStiffness(context, b); + } + } + + for (sofa::Size a = 0; a < NumberOfNodesInElement; ++a) + { + const auto weightTimesShapeFunction = static_cast(weight) * N[a]; + + for (sofa::Size b = 0; b < NumberOfNodesInElement; ++b) + { + dfdx(element[a] * spatial_dimensions, element[b] * spatial_dimensions) + += sourceDerivative[b] * weightTimesShapeFunction; + } + } + } + } } template 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 index 6ba2126bff6..f92a0808543 100644 --- 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 @@ -59,11 +59,7 @@ class NodalPressure : public sofa::core::BaseNodalProperty SOFA_CLASS(SOFA_TEMPLATE2(PressureSourceTerm, DataTypes, ElementType), SOFA_TEMPLATE2(BaseSourceTerm, DataTypes, ElementType)); + using Real = sofa::Real_t; using Deriv = sofa::Deriv_t; using QuadratureContext_t = QuadratureContext; + using SourceDerivative = typename BaseSourceTerm::SourceDerivative; using NodalPressure = ::sofa::component::solidmechanics::fem::elastic::NodalPressure; /** @@ -89,20 +87,26 @@ class PressureSourceTerm : public BaseSourceTerm sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_pressure; /** - * @brief Initializes the component and checks that a pressure is linked. + * @brief Initializes the component. */ void init() override; /** - * @brief The linked pressure interpolated at the quadrature point, times the unit normal. + * @brief The pressure interpolated at the quadrature point, times the unit normal. */ Deriv evaluate(const QuadratureContext_t& context) const override; /** - * @brief The nodal values this term interpolates, for FEMSourceTermIntegrator to track. + * @brief The nodal values of the pressure. */ sofa::type::vector integrandInputs() const override; + /** + * @brief The pressure times the derivative of the area normal of the element. + */ + SourceDerivative evaluateStiffness(const QuadratureContext_t& context, + sofa::Size node) const override; + protected: PressureSourceTerm(); 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 index 4947aa792a8..9caffcbe7de 100644 --- 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 @@ -76,4 +76,17 @@ PressureSourceTerm::integrandInputs() const return {&l_pressure->d_property}; } +template +auto PressureSourceTerm::evaluateStiffness( + const QuadratureContext_t& context, sofa::Size node) const -> SourceDerivative +{ + if (!l_pressure) + { + return SourceDerivative{}; + } + + return elementAreaNormalDerivative(context.jacobian, context.gradientShapeFunctions[node]) + * 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.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/StressSourceTerm.h index 5ea8c00fb91..28c363f6f41 100644 --- 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 @@ -120,6 +120,7 @@ class BaseStressSourceTerm : public BaseSourceTerm using Deriv = sofa::Deriv_t; using QuadratureContext_t = QuadratureContext; + using SourceDerivative = typename BaseSourceTerm::SourceDerivative; /** * @brief Nodal values of the stress tensor this term integrates. @@ -142,6 +143,12 @@ class BaseStressSourceTerm : public BaseSourceTerm */ sofa::type::vector integrandInputs() const override; + /** + * @brief The stress tensor applied to the derivative of the area normal of the element. + */ + SourceDerivative evaluateStiffness(const QuadratureContext_t& context, + sofa::Size node) const override; + protected: BaseStressSourceTerm(); 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 index d6a322d2c9c..aedb23fb3dc 100644 --- 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 @@ -77,4 +77,17 @@ BaseStressSourceTerm::integrandInputs() co return {&l_stress->d_property}; } +template +auto BaseStressSourceTerm::evaluateStiffness( + const QuadratureContext_t& context, sofa::Size node) const -> SourceDerivative +{ + if (!l_stress) + { + return SourceDerivative{}; + } + + return this->interpolateProperty(*l_stress, context) + * elementAreaNormalDerivative(context.jacobian, context.gradientShapeFunctions[node]); +} + } // 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 index 424b9f7e527..c0eb99e2675 100644 --- 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 @@ -59,10 +59,6 @@ class NodalSourceDensity : public sofa::core::BaseNodalProperty using Deriv = sofa::Deriv_t; using QuadratureContext_t = QuadratureContext; + using SourceDerivative = typename BaseSourceTerm::SourceDerivative; using NodalSourceDensity = ::sofa::component::solidmechanics::fem::elastic::NodalSourceDensity; @@ -88,20 +85,29 @@ class VectorSourceTerm : public BaseSourceTerm sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_sourceDensity; /** - * @brief Initializes the component and checks that a source density is linked. + * @brief Initializes the component. */ void init() override; /** - * @brief The linked source density interpolated at the quadrature point. + * @brief The source density interpolated at the quadrature point. */ Deriv evaluate(const QuadratureContext_t& context) const override; /** - * @brief The nodal values this term interpolates, for FEMSourceTermIntegrator to track. + * @brief The nodal values of the source density. */ sofa::type::vector integrandInputs() const override; + /** + * @brief The source density times the derivative of the measure: + * + * \f[ \frac{\partial (|\det J| \, r)}{\partial x_b} + * = |\det J| \; r \otimes \left( J^{-T} \frac{\partial N_b}{\partial q} \right) \f] + */ + SourceDerivative evaluateStiffness(const QuadratureContext_t& context, + sofa::Size node) const override; + protected: VectorSourceTerm(); 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 index 1cead57afb4..a8b3c2b3593 100644 --- 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 @@ -76,4 +76,19 @@ VectorSourceTerm::integrandInputs() const return {&l_sourceDensity->d_property}; } +template +auto VectorSourceTerm::evaluateStiffness( + const QuadratureContext_t& context, sofa::Size node) const -> SourceDerivative +{ + if (!l_sourceDensity) + { + return SourceDerivative{}; + } + + const auto measureDerivative = sofa::type::inverse(context.jacobian).transposed() + * context.gradientShapeFunctions[node] * context.measure; + + return sofa::type::dyad(this->interpolateProperty(*l_sourceDensity, context), measureDerivative); +} + } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn new file mode 100644 index 00000000000..81d0718da1b --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/PressureSourceTerm.scn @@ -0,0 +1,165 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn b/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn index fb3723ee104..61fa3584ca3 100644 --- a/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn +++ b/examples/Component/SolidMechanics/FEM/StressSourceTerm.scn @@ -1,13 +1,10 @@ - @@ -30,20 +27,97 @@ - - + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + @@ -62,7 +136,7 @@ - + - + @@ -95,7 +169,7 @@ - +