diff --git a/Sofa/Component/CMakeLists.txt b/Sofa/Component/CMakeLists.txt index 02e18a5d583..0a96063b030 100644 --- a/Sofa/Component/CMakeLists.txt +++ b/Sofa/Component/CMakeLists.txt @@ -15,10 +15,10 @@ sofa_add_subdirectory_modules(SOFACOMPONENT_TARGETS LinearSolver Mass ODESolver - Diffusion Mapping StateContainer SolidMechanics + Diffusion Constraint AnimationLoop MechanicalLoad diff --git a/Sofa/Component/Diffusion/CMakeLists.txt b/Sofa/Component/Diffusion/CMakeLists.txt index 0d73d1e6ab9..803d37b7b46 100644 --- a/Sofa/Component/Diffusion/CMakeLists.txt +++ b/Sofa/Component/Diffusion/CMakeLists.txt @@ -6,19 +6,22 @@ set(SOFACOMPONENTDIFFUSION_SOURCE_DIR "src/sofa/component/diffusion") set(HEADER_FILES ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/config.h.in ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/init.h + ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/DiffusionFEMForceField.h + ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/DiffusionFEMForceField.inl ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/TetrahedronDiffusionFEMForceField.h ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/TetrahedronDiffusionFEMForceField.inl ) set(SOURCE_FILES ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/init.cpp + ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/DiffusionFEMForceField.cpp ${SOFACOMPONENTDIFFUSION_SOURCE_DIR}/TetrahedronDiffusionFEMForceField.cpp ) -sofa_find_package(Sofa.Simulation.Core REQUIRED) +sofa_find_package(Sofa.Component.SolidMechanics.FEM.Elastic REQUIRED) add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES}) -target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core) +target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Component.SolidMechanics.FEM.Elastic) sofa_create_package_with_targets( PACKAGE_NAME ${PROJECT_NAME} diff --git a/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.cpp b/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.cpp new file mode 100644 index 00000000000..e52ef635a91 --- /dev/null +++ b/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.cpp @@ -0,0 +1,128 @@ +/****************************************************************************** +* 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_DIFFUSION_DIFFUSION_FEM_FORCE_FIELD_CPP + +#include + +#include +#include + +namespace sofa::component::diffusion +{ + +using sofa::defaulttype::Vec1Types; +using sofa::defaulttype::Vec2Types; +using sofa::defaulttype::Vec3Types; +using sofa::geometry::Edge; +using sofa::geometry::Triangle; +using sofa::geometry::Quad; +using sofa::geometry::Tetrahedron; +using sofa::geometry::Hexahedron; +using sofa::geometry::Prism; +using sofa::geometry::Pyramid; + +void registerDiffusionFEMForceField(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Generic isotropic FEM diffusion force field") + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + .add< DiffusionFEMForceField >() + ); +} + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +} // namespace sofa::component::diffusion diff --git a/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.h b/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.h new file mode 100644 index 00000000000..24307302609 --- /dev/null +++ b/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.h @@ -0,0 +1,163 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +#include +#include + +#if !defined(SOFA_COMPONENT_DIFFUSION_DIFFUSION_FEM_FORCE_FIELD_CPP) +#include +#endif + +namespace sofa::component::diffusion +{ +template +class DiffusionFEMForceField : + public sofa::component::solidmechanics::fem::elastic::FEMForceField +{ +public: + SOFA_CLASS( + SOFA_TEMPLATE3(DiffusionFEMForceField, DiffusionFieldTypes, SpatialTypes, ElementType), + SOFA_TEMPLATE3(sofa::component::solidmechanics::fem::elastic::FEMForceField, DiffusionFieldTypes, SpatialTypes, ElementType)); + +private: + using Inherited = sofa::component::solidmechanics::fem::elastic::FEMForceField; + using trait = sofa::component::solidmechanics::fem::elastic::trait; + + using Real = sofa::Real_t; // type of the diffused quantity + using GeometryReal = sofa::Real_t; // type of the geometry + + using ElementHessian = typename trait::ElementHessian; + using ElementGradient = typename trait::ElementGradient; + using ElementValue = typename trait::ElementDisplacement; + + using GeometryCoord = sofa::Coord_t; + + /// Number of components of the diffused quantity at each node (1 for a scalar field) + static constexpr sofa::Size FieldDimension = trait::field_dimensions; + /// Dimension of the space in which the mesh lives + static constexpr sofa::Size SpatialDimension = trait::spatial_dimensions; + static constexpr sofa::Size NumberOfNodes = trait::NumberOfNodesInElement; + + using FieldBlock = sofa::type::Mat; + +public: + DiffusionFEMForceField(); + + void init() override; + + void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override; + + using sofa::core::behavior::ForceField::getPotentialEnergy; + SReal getPotentialEnergy( + const sofa::core::MechanicalParams*, + const sofa::DataVecCoord_t& fieldValues) const override; + + using sofa::core::behavior::ForceField::addKToMatrix; + void addKToMatrix(sofa::linearalgebra::BaseMatrix* matrix, SReal kFact, unsigned& offset) override; + +protected: + void computeElementsForces( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForces, + const sofa::VecCoord_t& nodeValues) override; + + void computeElementsForcesDeriv( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams* mparams, + sofa::type::vector& elementForcesDeriv, + const sofa::VecDeriv_t& nodeDx) override; + + void validateSpatialState(); + +private: + void precomputeElementStiffness(); + + /// Gathers the nodal values of the diffused quantity (or of its increment) of one element + template + ElementValue computeElementValue( + const typename trait::TopologyElement& element, + const VecType& nodeValues) const; + + Real getElementDiffusionCoefficient(const sofa::helper::ReadAccessor > >& elementDiffusionCoefficient, sofa::Index elementId) const; + +public: + static constexpr Real defaultDiffusionCoefficient = 1; + + sofa::Data > d_elementDiffusionCoefficient; + sofa::Data > d_elementHessian; + + sofa::SingleLink, + sofa::core::objectmodel::BaseLink::FLAG_STRONGLINK> l_spatialState; +}; + +#if !defined(SOFA_COMPONENT_DIFFUSION_DIFFUSION_FEM_FORCE_FIELD_CPP) +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; + +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +extern template class SOFA_COMPONENT_DIFFUSION_API DiffusionFEMForceField; +#endif + +} // namespace sofa::component::diffusion diff --git a/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.inl b/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.inl new file mode 100644 index 00000000000..ef260205ae2 --- /dev/null +++ b/Sofa/Component/Diffusion/src/sofa/component/diffusion/DiffusionFEMForceField.inl @@ -0,0 +1,372 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +#include +#include + +#include +#include + +#include + +namespace sofa::component::diffusion +{ +template +DiffusionFEMForceField::DiffusionFEMForceField() + : d_elementDiffusionCoefficient(initData( + &d_elementDiffusionCoefficient, + sofa::type::vector(defaultDiffusionCoefficient), + "elementDiffusionCoefficient", + "List of isotropic diffusion coefficient per element. The last coefficient is used if the list does not match the number of elements.")) + , d_elementHessian(initData( + &d_elementHessian, + "elementStiffness", + "Precomputed diffusion stiffness matrix per element")) + , l_spatialState(initLink("spatialState", "Link to the state holding the spatial positions of the mesh")) +{ + this->addUpdateCallback("precomputeDiffusionStiffness", {&d_elementDiffusionCoefficient}, + [this](const sofa::core::DataTracker&) + { + precomputeElementStiffness(); + return this->getComponentState(); + }, + {}); +} + +template +void DiffusionFEMForceField::validateSpatialState() +{ + if (!l_spatialState.get()) + { + l_spatialState.set( + this->getContext()->template get >()); + + if (!l_spatialState) + { + msg_error() << "No compatible MechanicalState found in the current context. " + "Please set the link `spatialState` explicitly."; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } + else + { + msg_info() << "No link given for `spatialState`, but it has been found in the context"; + } + } +} + +template +void DiffusionFEMForceField::init() +{ + Inherited::init(); + + if (!this->isComponentStateInvalid()) + { + validateSpatialState(); + } + + if (!this->isComponentStateInvalid()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + precomputeElementStiffness(); + } +} + +template +void DiffusionFEMForceField::precomputeElementStiffness() +{ + if (!this->l_topology || !this->mstate || !this->l_spatialState) + return; + + if (this->isComponentStateInvalid()) + return; + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + + const auto elementDiffusion = sofa::helper::getReadAccessor(d_elementDiffusionCoefficient); + + const auto restPositionAccessor = this->l_spatialState->readRestPositions(); + + if (restPositionAccessor.size() != this->mstate->getSize()) + { + msg_error() << "The spatial state (" << restPositionAccessor.size() + << " nodes) and the diffused field state (" << this->mstate->getSize() + << " nodes) must have the same number of nodes"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; + } + + auto elementHessian = sofa::helper::getWriteOnlyAccessor(d_elementHessian); + elementHessian.resize(elements.size()); + + SCOPED_TIMER("precomputeDiffusionHessian"); + + sofa::helper::IotaView indices { static_cast(0ul), elements.size() }; + std::for_each(indices.begin(), indices.end(), + [&](const auto elementId) + { + const auto& element = elements[elementId]; + const Real diffusion = getElementDiffusionCoefficient(elementDiffusion, elementId); + + const std::array nodesCoordinates = + sofa::component::solidmechanics::fem::elastic::extractNodesVectorFromGlobalVector( + element, restPositionAccessor.ref()); + + ElementHessian K; + K.clear(); + + for (const auto& [quadraturePoint, weight] : trait::FiniteElement::quadraturePoints()) + { + // Gradients of the shape functions in the reference element + const sofa::type::Mat dN_dq_ref = + trait::FiniteElement::gradientShapeFunctions(quadraturePoint); + + // Jacobian of the geometric mapping (spatial dim x topological dim): geometry only + sofa::type::Mat jacobian; + jacobian.clear(); + for (sofa::Size i = 0; i < NumberOfNodes; ++i) + { + jacobian += sofa::type::dyad(nodesCoordinates[i], dN_dq_ref[i]); + } + + const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); + + const sofa::type::Mat J_inv = + sofa::type::inverse(jacobian); + + // Gradients of the shape functions in the spatial frame + sofa::type::Mat dN_dx(sofa::type::NOINIT); + for (sofa::Size i = 0; i < NumberOfNodes; ++i) + { + dN_dx[i] = J_inv.transposed() * dN_dq_ref[i]; + } + + const Real integrationFactor = diffusion * static_cast(weight * detJ); + + for (sofa::Index i = 0; i < NumberOfNodes; ++i) + { + for (sofa::Index j = 0; j < NumberOfNodes; ++j) + { + // scalar product of the spatial gradients: sum over SPATIAL dimensions + GeometryReal gradDot {}; + for (sofa::Index d = 0; d < SpatialDimension; ++d) + { + gradDot += dN_dx[i][d] * dN_dx[j][d]; + } + + const Real scalarStiffness = integrationFactor * static_cast(gradDot); + + // the same scalar stiffness applies to each component of the FIELD + for (sofa::Index c = 0; c < FieldDimension; ++c) + { + K(i * FieldDimension + c, j * FieldDimension + c) += scalarStiffness; + } + } + } + } + + elementHessian[elementId] = K; + }); +} + +template +template +auto DiffusionFEMForceField::computeElementValue( + const typename trait::TopologyElement& element, + const VecType& nodeValues) const -> ElementValue +{ + ElementValue values { sofa::type::NOINIT }; + + for (sofa::Size nodeId = 0; nodeId < NumberOfNodes; ++nodeId) + { + const auto globalNodeId = element[nodeId]; + + for (sofa::Size c = 0; c < FieldDimension; ++c) + { + values[nodeId * FieldDimension + c] = nodeValues[globalNodeId][c]; + } + } + + return values; +} + +template +auto DiffusionFEMForceField::getElementDiffusionCoefficient( + const sofa::helper::ReadAccessor>>& elementDiffusionCoefficient, + sofa::Index elementId) const -> Real +{ + if (elementId < elementDiffusionCoefficient.size()) + return elementDiffusionCoefficient[elementId]; + if (!elementDiffusionCoefficient.empty()) + return elementDiffusionCoefficient->back(); + return defaultDiffusionCoefficient; +} + +template +void DiffusionFEMForceField::computeElementsForces( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams*, + sofa::type::vector& elementForces, + const sofa::VecCoord_t& nodeValues) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + const auto elementStiffness = sofa::helper::getReadAccessor(d_elementHessian); + + for (std::size_t elementId = range.start; elementId < range.end; ++elementId) + { + const auto elementValue = computeElementValue(elements[elementId], nodeValues); + elementForces[elementId] = (elementStiffness[elementId] * elementValue); + } +} + +template +void DiffusionFEMForceField::computeElementsForcesDeriv( + const sofa::simulation::Range& range, + const sofa::core::MechanicalParams*, + sofa::type::vector& elementForcesDeriv, + const sofa::VecDeriv_t& nodeDx) +{ + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + const auto elementStiffness = sofa::helper::getReadAccessor(d_elementHessian); + + for (std::size_t elementId = range.start; elementId < range.end; ++elementId) + { + const auto elementDx = computeElementValue(elements[elementId], nodeDx); + + elementForcesDeriv[elementId] = (elementStiffness[elementId] * elementDx); + } +} + +template +void DiffusionFEMForceField::buildStiffnessMatrix( + sofa::core::behavior::StiffnessMatrix* matrix) +{ + if (this->isComponentStateInvalid() || !this->l_topology) + return; + + auto dfdx = matrix->getForceDerivativeIn(this->mstate) + .withRespectToPositionsIn(this->mstate); + + FieldBlock localMatrix { sofa::type::NOINIT }; + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + const auto elementStiffness = sofa::helper::getReadAccessor(d_elementHessian); + + if (elementStiffness.size() < elements.size()) + return; + + auto elementStiffnessIt = elementStiffness.begin(); + for (const auto& element : elements) + { + const auto& stiffnessMatrix = *elementStiffnessIt++; + + for (sofa::Index n1 = 0; n1 < NumberOfNodes; ++n1) + { + for (sofa::Index n2 = 0; n2 < NumberOfNodes; ++n2) + { + stiffnessMatrix.getsub( + FieldDimension * n1, + FieldDimension * n2, + localMatrix); + + dfdx( + element[n1] * FieldDimension, + element[n2] * FieldDimension) += -localMatrix; + } + } + } +} + +template +SReal DiffusionFEMForceField::getPotentialEnergy( + const sofa::core::MechanicalParams*, + const sofa::DataVecCoord_t& fieldValues) const +{ + if (this->isComponentStateInvalid() || !this->l_topology) + return 0; + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + const auto elementStiffness = sofa::helper::getReadAccessor(d_elementHessian); + const auto fieldAccessor = sofa::helper::getReadAccessor(fieldValues); + + if (elementStiffness.size() < elements.size()) + return 0; + + Real energy {}; + + for (std::size_t elementId = 0; elementId < elements.size(); ++elementId) + { + const auto elementValue = computeElementValue(elements[elementId], fieldAccessor.ref()); + + energy += elementValue * (elementStiffness[elementId] * elementValue); + } + + return static_cast(0.5 * energy); +} + +template +void DiffusionFEMForceField::addKToMatrix( + sofa::linearalgebra::BaseMatrix* matrix, + SReal kFact, + unsigned& offset) +{ + if (this->isComponentStateInvalid() || !this->l_topology) + return; + + FieldBlock localMatrix { sofa::type::NOINIT }; + + const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology); + const auto elementStiffness = sofa::helper::getReadAccessor(d_elementHessian); + + if (elementStiffness.size() < elements.size()) + return; + + auto elementStiffnessIt = elementStiffness.begin(); + for (const auto& element : elements) + { + const auto& stiffnessMatrix = *elementStiffnessIt++; + + for (sofa::Index n1 = 0; n1 < NumberOfNodes; ++n1) + { + for (sofa::Index n2 = 0; n2 < NumberOfNodes; ++n2) + { + stiffnessMatrix.getsub( + FieldDimension * n1, + FieldDimension * n2, + localMatrix); + + const auto value = + (-static_cast(kFact)) + * static_cast >(localMatrix); + + matrix->add( + offset + element[n1] * FieldDimension, + offset + element[n2] * FieldDimension, + value); + } + } + } +} + +} // namespace sofa::component::diffusion diff --git a/Sofa/Component/Diffusion/src/sofa/component/diffusion/init.cpp b/Sofa/Component/Diffusion/src/sofa/component/diffusion/init.cpp index a28a8777262..ae4d1b764bc 100644 --- a/Sofa/Component/Diffusion/src/sofa/component/diffusion/init.cpp +++ b/Sofa/Component/Diffusion/src/sofa/component/diffusion/init.cpp @@ -27,6 +27,7 @@ namespace sofa::component::diffusion { extern void registerTetrahedronDiffusionFEMForceField(sofa::core::ObjectFactory* factory); +extern void registerDiffusionFEMForceField(sofa::core::ObjectFactory* factory); extern "C" { SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule(); @@ -53,6 +54,7 @@ const char* getModuleVersion() void registerObjects(sofa::core::ObjectFactory* factory) { registerTetrahedronDiffusionFEMForceField(factory); + registerDiffusionFEMForceField(factory); } void init() diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h index 19a441f9b82..e5cd4f0b7fd 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h @@ -50,7 +50,7 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f void init() override; private: - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; using ElementHessian = typename trait::ElementHessian; using StrainDisplacement = typename trait::StrainDisplacement; using ElementDisplacement = typename trait::ElementDisplacement; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h index beaaf5995ef..0f0a61f55f5 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h @@ -129,17 +129,17 @@ struct RotationMethods : Rotati template class CorotationalFEMForceField : public BaseElementLinearFEMForceField, - public FEMForceField, + public FEMForceField, public CauchyStressEvaluator { public: SOFA_CLASS2( SOFA_TEMPLATE2(CorotationalFEMForceField, DataTypes, ElementType), SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType), - SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType)); + SOFA_TEMPLATE3(FEMForceField, DataTypes, DataTypes, ElementType)); private: - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; using ElementGradient = typename trait::ElementGradient; using ElementDisplacement = typename trait::ElementDisplacement; using RotationMatrix = sofa::type::Mat>; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl index fc5fc9fd809..b14850ffffc 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.inl @@ -50,7 +50,7 @@ template void CorotationalFEMForceField::init() { BaseElementLinearFEMForceField::init(); - FEMForceField::init(); + FEMForceField::init(); if (!this->isComponentStateInvalid()) { diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.cpp index 4ab52ad618e..fe1679fb524 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.cpp @@ -25,13 +25,60 @@ namespace sofa::component::solidmechanics::fem::elastic { -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +using sofa::defaulttype::Vec1Types; +using sofa::defaulttype::Vec2Types; +using sofa::defaulttype::Vec3Types; +using sofa::geometry::Edge; +using sofa::geometry::Triangle; +using sofa::geometry::Quad; +using sofa::geometry::Tetrahedron; +using sofa::geometry::Hexahedron; +using sofa::geometry::Prism; +using sofa::geometry::Pyramid; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h index c85babae37a..73f379f2120 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.h @@ -36,20 +36,27 @@ namespace sofa::component::solidmechanics::fem::elastic { -template +/** + * Base class for finite element force fields. + * + * @tparam DataTypes The type of the quantity to integrate in the finite elements + * @tparam SpatialTypes The type of the spatial coordinates + * @tparam ElementType The type of the finite element + */ +template class FEMForceField : public virtual sofa::core::behavior::ForceField, public virtual sofa::core::behavior::TopologyAccessor, public virtual sofa::simulation::TaskSchedulerUser { public: - SOFA_CLASS3(SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType), + SOFA_CLASS3(SOFA_TEMPLATE3(FEMForceField, DataTypes, SpatialTypes, ElementType), sofa::core::behavior::ForceField, sofa::core::behavior::TopologyAccessor, sofa::simulation::TaskSchedulerUser); protected: - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; using ElementGradient = typename trait::ElementGradient; public: @@ -72,7 +79,7 @@ class FEMForceField : sofa::Data d_computeForceStrategy; sofa::Data d_computeForceDerivStrategy; - sofa::Data> d_elementSpace; + sofa::Data> d_elementSpace; protected: @@ -137,15 +144,50 @@ class FEMForceField : }; #if !defined(SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_FEM_FORCEFIELD_CPP) -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; -extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API FEMForceField; #endif } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl index 44d6dc85e07..36859cafca7 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl @@ -27,11 +27,11 @@ namespace sofa::component::solidmechanics::fem::elastic { -template -FEMForceField::FEMForceField() +template +FEMForceField::FEMForceField() : d_computeForceStrategy(initData(&d_computeForceStrategy, "computeForceStrategy", std::string("The compute strategy used to compute the forces.\n" + ComputeStrategy::dataDescription()).c_str())) , d_computeForceDerivStrategy(initData(&d_computeForceDerivStrategy, "computeForceDerivStrategy", std::string("The compute strategy used to compute the forces derivatives.\n" + ComputeStrategy::dataDescription()).c_str())) - , d_elementSpace(initData(&d_elementSpace, static_cast>(0.125), "elementSpace", "When rendering, the space between elements")) + , d_elementSpace(initData(&d_elementSpace, static_cast>(0.125), "elementSpace", "When rendering, the space between elements")) { d_elementSpace.setGroup("Visualization"); @@ -39,8 +39,8 @@ FEMForceField::FEMForceField() d_computeForceDerivStrategy.setGroup("Multithreading"); } -template -void FEMForceField::init() +template +void FEMForceField::init() { sofa::core::behavior::ForceField::init(); @@ -55,16 +55,18 @@ void FEMForceField::init() } } -template -void FEMForceField::addForce( +template +void FEMForceField::addForce( const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& f, const sofa::DataVecCoord_t& x, const sofa::DataVecDeriv_t& v) { - SOFA_UNUSED(mparams); SOFA_UNUSED(v); + if (this->isComponentStateInvalid()) + return; + auto positionAccessor = sofa::helper::getReadAccessor(x); if (this->l_topology == nullptr) return; @@ -85,8 +87,8 @@ void FEMForceField::addForce( dispatchElementForcesToNodes(elements, forceAccessor.wref()); } -template -void FEMForceField::computeElementsForces( +template +void FEMForceField::computeElementsForces( const sofa::core::MechanicalParams* mparams, sofa::type::vector& f, const sofa::VecCoord_t& x) @@ -104,8 +106,8 @@ void FEMForceField::computeElementsForces( }); } -template -void FEMForceField::dispatchElementForcesToNodes( +template +void FEMForceField::dispatchElementForcesToNodes( const sofa::type::vector& elements, sofa::VecDeriv_t& nodeForces) { @@ -119,16 +121,16 @@ void FEMForceField::dispatchElementForcesToNodes( for (sofa::Size j = 0; j < trait::NumberOfNodesInElement; ++j) { auto& nodeForce = nodeForces[element[j]]; - for (sofa::Size k = 0; k < trait::spatial_dimensions; ++k) + for (sofa::Size k = 0; k < trait::field_dimensions; ++k) { - nodeForce[k] -= elementForce[j * trait::spatial_dimensions + k]; + nodeForce[k] -= elementForce[j * trait::field_dimensions + k]; } } } } -template -void FEMForceField::addDForce( +template +void FEMForceField::addDForce( const sofa::core::MechanicalParams* mparams, sofa::DataVecDeriv_t& df, const sofa::DataVecDeriv_t& dx) @@ -157,13 +159,12 @@ void FEMForceField::addDForce( dispatchElementForcesDerivToNodes(mparams, elements, dfAccessor.wref()); } -template -void FEMForceField::computeElementsForcesDeriv( +template +void FEMForceField::computeElementsForcesDeriv( const sofa::core::MechanicalParams* mparams, sofa::type::vector& df, const sofa::VecDeriv_t& dx) { SCOPED_TIMER("ElementForcesDeriv"); - beforeElementForceDeriv(mparams); this->beforeElementForceDeriv(mparams); @@ -178,8 +179,8 @@ void FEMForceField::computeElementsForcesDeriv( }); } -template -void FEMForceField::dispatchElementForcesDerivToNodes(const sofa::core::MechanicalParams* mparams, +template +void FEMForceField::dispatchElementForcesDerivToNodes(const sofa::core::MechanicalParams* mparams, const sofa::type::vector& elements, sofa::VecDeriv_t& nodeForcesDeriv) { @@ -198,17 +199,17 @@ void FEMForceField::dispatchElementForcesDerivToNodes(co const auto nodeId = element[i]; auto& df = nodeForcesDeriv[nodeId]; - for (sofa::Size dim = 0; dim < trait::spatial_dimensions; ++dim) + for (sofa::Size dim = 0; dim < trait::field_dimensions; ++dim) { - df[dim] -= kFactor * elementDForce[i * trait::spatial_dimensions + dim]; + df[dim] -= kFactor * elementDForce[i * trait::field_dimensions + dim]; } } } } -template -sofa::simulation::ForEachExecutionPolicy FEMForceField::getExecutionPolicy( +template +sofa::simulation::ForEachExecutionPolicy FEMForceField::getExecutionPolicy( const sofa::Data& strategy) const { auto computeForceStrategyAccessor = sofa::helper::getReadAccessor(strategy); @@ -219,8 +220,8 @@ sofa::simulation::ForEachExecutionPolicy FEMForceField:: : sofa::simulation::ForEachExecutionPolicy::SEQUENTIAL; } -template -void FEMForceField::draw(const sofa::core::visual::VisualParams* vparams) +template +void FEMForceField::draw(const sofa::core::visual::VisualParams* vparams) { if (!vparams->displayFlags().getShowForceFields()) return; @@ -239,8 +240,8 @@ void FEMForceField::draw(const sofa::core::visual::Visua m_drawMesh.drawAllElements(vparams->drawTool(), x, this->l_topology.get()); } -template -void FEMForceField::computeBBox(const core::ExecParams* params, bool onlyVisible) +template +void FEMForceField::computeBBox(const core::ExecParams* params, bool onlyVisible) { SOFA_UNUSED(params); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h index 2beaf66e1b3..0ad3366190e 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h @@ -37,18 +37,18 @@ namespace sofa::component::solidmechanics::fem::elastic template class LinearSmallStrainFEMForceField : public BaseElementLinearFEMForceField, - public FEMForceField, + public FEMForceField, public CauchyStressEvaluator { public: SOFA_CLASS3( SOFA_TEMPLATE2(LinearSmallStrainFEMForceField, DataTypes, ElementType), SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType), - SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType), + SOFA_TEMPLATE3(FEMForceField, DataTypes, DataTypes, ElementType), CauchyStressEvaluator); private: - using trait = typename FEMForceField::trait; + using trait = typename FEMForceField::trait; using ElementHessian = typename trait::ElementHessian; using ElementDisplacement = typename trait::ElementDisplacement; using StrainDisplacement = typename trait::StrainDisplacement; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl index 1a4d22e427f..84c5faa3e47 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.inl @@ -31,7 +31,7 @@ template void LinearSmallStrainFEMForceField::init() { BaseElementLinearFEMForceField::init(); - FEMForceField::init(); + FEMForceField::init(); if (!this->isComponentStateInvalid()) { diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h index 97fa1d979e4..e4563978f45 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/ElementStiffnessMatrix.h @@ -45,7 +45,7 @@ auto integrate( static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; - typename trait::ElementHessian K; + typename trait::ElementHessian K; auto C = elasticityTensor.toVoigtMatSym().toMat(); diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h index 91b94764ba5..5d9b67612c4 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/trait.h @@ -30,10 +30,11 @@ namespace sofa::component::solidmechanics::fem::elastic { -template +template struct trait { - using DataVecCoord = sofa::DataVecDeriv_t; + // ---- Types of the field (displacement, temperature, concentration, ...) ---- + using DataVecCoord = sofa::DataVecCoord_t; using DataVecDeriv = sofa::DataVecDeriv_t; using VecCoord = sofa::VecCoord_t; using VecDeriv = sofa::VecDeriv_t; @@ -41,23 +42,35 @@ struct trait using Deriv = sofa::Deriv_t; using Real = sofa::Real_t; - using FiniteElement = sofa::fem::FiniteElement; + // ---- Types of the geometry ---- + using GeometryCoord = sofa::Coord_t; + using GeometryReal = sofa::Real_t; + + using FiniteElement = sofa::fem::FiniteElement; using TopologyElement = typename FiniteElement::TopologyElement; using ReferenceCoord = typename FiniteElement::ReferenceCoord; - static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + /// Dimension of the space in which the mesh lives + static constexpr sofa::Size spatial_dimensions = SpatialTypes::coord_total_size; + /// Number of components of the field at each node (1 for a scalar field) + static constexpr sofa::Size field_dimensions = DataTypes::coord_total_size; + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; - static constexpr sofa::Size NumberOfDofsInElement = NumberOfNodesInElement * spatial_dimensions; + static constexpr sofa::Size NumberOfDofsInElement = NumberOfNodesInElement * field_dimensions; static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; static constexpr sofa::Size NbQuadraturePoints = FiniteElement::quadraturePoints().size(); - /// type of 2nd-order tensor for the elasticity tensor for isotropic materials - using ElasticityTensor = sofa::type::FullySymmetric4Tensor; + static_assert(field_dimensions >= 1, "The field must have at least one component"); + static_assert(TopologicalDimension <= spatial_dimensions, + "The element cannot have a higher dimension than the space it is embedded in"); + + /// 4th-order tensor for isotropic elasticity: lives in physical space + using ElasticityTensor = sofa::type::FullySymmetric4Tensor; /// the type of B in e = B d, if e is the strain, and d is the displacement - using StrainDisplacement = sofa::component::solidmechanics::fem::elastic::StrainDisplacement; + using StrainDisplacement = sofa::component::solidmechanics::fem::elastic::StrainDisplacement; - /// the concatenation of the displacement of the element nodes in a single vector + /// the concatenation of the nodal values (displacement, temperature, ...) of the element nodes in a single vector using ElementDisplacement = sofa::type::Vec; /// the type of the element hessian matrix diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.h index e1e536b6370..6d4ef440de2 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.h @@ -40,12 +40,12 @@ struct StiffnessFilter; template class HyperelasticityFEMForceField : - public sofa::component::solidmechanics::fem::elastic::FEMForceField + public sofa::component::solidmechanics::fem::elastic::FEMForceField { public: SOFA_CLASS( SOFA_TEMPLATE2(HyperelasticityFEMForceField, TDataTypes, TElementType), - SOFA_TEMPLATE2(sofa::component::solidmechanics::fem::elastic::FEMForceField, TDataTypes, TElementType)); + SOFA_TEMPLATE3(sofa::component::solidmechanics::fem::elastic::FEMForceField, TDataTypes, TDataTypes, TElementType)); using DataTypes = TDataTypes; @@ -102,7 +102,7 @@ class HyperelasticityFEMForceField : HyperelasticityFEMForceField(); - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; using ElementGradient = typename trait::ElementGradient; void validateMaterial(); diff --git a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.inl index f9969e26f31..96b3030bdd6 100644 --- a/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.inl +++ b/Sofa/Component/SolidMechanics/FEM/HyperElastic/src/sofa/component/solidmechanics/fem/hyperelastic/HyperelasticityFEMForceField.inl @@ -43,7 +43,7 @@ HyperelasticityFEMForceField::HyperelasticityFEMForceF template void HyperelasticityFEMForceField::init() { - sofa::component::solidmechanics::fem::elastic::FEMForceField::init(); + sofa::component::solidmechanics::fem::elastic::FEMForceField::init(); if (!this->isComponentStateInvalid()) { diff --git a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaCorotationalFEMForceField.inl b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaCorotationalFEMForceField.inl index 982ad644a51..f4509b793fb 100644 --- a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaCorotationalFEMForceField.inl +++ b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaCorotationalFEMForceField.inl @@ -43,7 +43,7 @@ void CudaCorotationalFEMForceField::init() template void CudaCorotationalFEMForceField::uploadStiffnessAndConnectivity() { - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; if (!this->l_topology) return; @@ -152,7 +152,7 @@ void CudaCorotationalFEMForceField::uploadStiffnessAndCo template void CudaCorotationalFEMForceField::uploadRotations() { - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; constexpr auto dim = trait::spatial_dimensions; const auto& rotations = this->m_rotations; @@ -176,7 +176,7 @@ void CudaCorotationalFEMForceField::uploadRotations() template void CudaCorotationalFEMForceField::uploadInitialRotationsTransposed() { - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; constexpr auto dim = trait::spatial_dimensions; constexpr auto nNodes = trait::NumberOfNodesInElement; @@ -206,7 +206,7 @@ void CudaCorotationalFEMForceField::uploadInitialRotatio template void CudaCorotationalFEMForceField::downloadRotations() { - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; constexpr auto dim = trait::spatial_dimensions; if (!m_gpuRotationsUploaded) return; @@ -251,7 +251,7 @@ void CudaCorotationalFEMForceField::addForce( return; } - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; constexpr auto nNodes = trait::NumberOfNodesInElement; constexpr auto dim = trait::spatial_dimensions; @@ -314,7 +314,7 @@ void CudaCorotationalFEMForceField::addDForce( return; } - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; constexpr auto nNodes = trait::NumberOfNodesInElement; constexpr auto dim = trait::spatial_dimensions; diff --git a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaLinearSmallStrainFEMForceField.inl b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaLinearSmallStrainFEMForceField.inl index 967bb65ccf5..e9e0ab6eace 100644 --- a/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaLinearSmallStrainFEMForceField.inl +++ b/applications/plugins/SofaCUDA/Component/src/SofaCUDA/component/solidmechanics/fem/elastic/CudaLinearSmallStrainFEMForceField.inl @@ -42,7 +42,7 @@ void CudaLinearSmallStrainFEMForceField::init() template void CudaLinearSmallStrainFEMForceField::uploadStiffnessAndConnectivity() { - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; if (!this->l_topology) return; @@ -160,7 +160,7 @@ void CudaLinearSmallStrainFEMForceField::addForce( return; } - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; constexpr auto nNodes = trait::NumberOfNodesInElement; constexpr auto dim = trait::spatial_dimensions; @@ -203,7 +203,7 @@ void CudaLinearSmallStrainFEMForceField::addDForce( return; } - using trait = sofa::component::solidmechanics::fem::elastic::trait; + using trait = sofa::component::solidmechanics::fem::elastic::trait; constexpr auto nNodes = trait::NumberOfNodesInElement; constexpr auto dim = trait::spatial_dimensions; diff --git a/examples/Component/Diffusion/DiffusionFEMForceField.scn b/examples/Component/Diffusion/DiffusionFEMForceField.scn new file mode 100644 index 00000000000..40f7c2d0630 --- /dev/null +++ b/examples/Component/Diffusion/DiffusionFEMForceField.scn @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +