diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index fddd6033a0d..85447e27cdc 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt @@ -13,6 +13,7 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/BaseLinearElasticityFEMForceField.inl ${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}/LinearSmallStrainFEMForceField.h @@ -45,6 +46,8 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceField.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceFieldOptim.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceFieldOptim.inl + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VonMisesStress.h + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VonMisesStress.inl ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/ComputeStrategy.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/ElementStiffnessMatrix.h @@ -84,6 +87,7 @@ set(SOURCE_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularAnisotropicFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceField.cpp ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/TriangularFEMForceFieldOptim.cpp + ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/VonMisesStress.cpp ) sofa_find_package(Sofa.Simulation.Core REQUIRED) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h new file mode 100644 index 00000000000..b890687e46b --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CauchyStressEvaluator.h @@ -0,0 +1,45 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include + +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +class CauchyStressEvaluator : public virtual core::objectmodel::BaseComponent +{ +public: + SOFA_ABSTRACT_CLASS(CauchyStressEvaluator, BaseComponent) + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + using StressVoigtVector = sofa::type::Vec, sofa::Real_t>; + using DeformationGradient = sofa::type::Mat>; + + virtual StressVoigtVector computeStress(const DeformationGradient& F, sofa::Size elementId) = 0; +}; + +} 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 899edf4c084..45bc796d644 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 @@ -24,7 +24,7 @@ #include #include #include - +#include #include #if !defined(ELASTICITY_COMPONENT_ELEMENT_LINEAR_SMALL_STRAIN_FEM_FORCE_FIELD_CPP) @@ -37,13 +37,15 @@ namespace sofa::component::solidmechanics::fem::elastic template class LinearSmallStrainFEMForceField : public BaseElementLinearFEMForceField, - public FEMForceField + public FEMForceField, + public CauchyStressEvaluator { public: - SOFA_CLASS2( + SOFA_CLASS3( SOFA_TEMPLATE2(LinearSmallStrainFEMForceField, DataTypes, ElementType), SOFA_TEMPLATE2(BaseElementLinearFEMForceField, DataTypes, ElementType), - SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType)); + SOFA_TEMPLATE2(FEMForceField, DataTypes, ElementType), + CauchyStressEvaluator); private: using trait = typename FEMForceField::trait; @@ -51,6 +53,8 @@ class LinearSmallStrainFEMForceField : using ElementDisplacement = typename trait::ElementDisplacement; using StrainDisplacement = typename trait::StrainDisplacement; using ElementGradient = typename trait::ElementGradient; + using DeformationGradient = typename trait::DeformationGradient; + using StressVoigtVector = typename trait::StressVoigtVector; public: void init() override; @@ -64,6 +68,8 @@ class LinearSmallStrainFEMForceField : // almost deprecated, but here for compatibility with unit tests void addKToMatrix(sofa::linearalgebra::BaseMatrix* matrix, SReal kFact, unsigned& offset) override; + StressVoigtVector computeStress(const DeformationGradient& F, sofa::Size elementId) override; + protected: void computeElementsForces( 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 6fea5667a99..caeaaef83f8 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 @@ -180,4 +180,32 @@ void LinearSmallStrainFEMForceField::addKToMatrix( } } +template +auto LinearSmallStrainFEMForceField::computeStress( + const DeformationGradient& F, sofa::Size elementId) -> StressVoigtVector +{ + const auto strainTensor = static_cast>(1)/2 * (F + F.transposed()) - DeformationGradient::Identity(); + + sofa::type::Vec, Real_t > strainVoigt; + for (sofa::Size i = 0; i < type::NumberOfIndependentElements; ++i) + { + const auto [p, q] = type::toTensorIndices(i); + strainVoigt[i] = strainTensor(p, q); + } + + const auto youngModulus = this->getYoungModulusInElement(elementId); + const auto poissonRatio = this->getPoissonRatioInElement(elementId); + + LameLambda> lambda { 0 }; + LameMu> mu { 0 }; + + sofa::component::solidmechanics::fem::elastic::toLameParameters>( + YoungModulus>(youngModulus), PoissonRatio>(poissonRatio), + lambda, mu); + + const auto elasticityTensor = makeIsotropicElasticityTensor>(mu, lambda); + + return elasticityTensor.toVoigtMatSym().toMat() * strainVoigt; +} + } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.cpp new file mode 100644 index 00000000000..6eb5fe5af53 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.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 * +******************************************************************************/ +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API VonMisesStress; + +void registerVonMisesStress(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Compute and draw von Mises stress based on a local least-square projection in each element") + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + .add< VonMisesStress >() + ); +} + +} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h new file mode 100644 index 00000000000..6ddcb89e844 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.h @@ -0,0 +1,137 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +/** + * @brief Component that computes and visualizes the von Mises stress field on a FEM mesh. + * + * This component evaluates the Cauchy stress at each quadrature point of the elements + * using a linked @see CauchyStressEvaluator. Since stress values are naturally calculated + * at quadrature points (e.g., Gauss points), but visualization typically requires nodal + * values, this component performs a least-squares projection. + * + * The projection maps the stress tensor components from the quadrature points to the + * element nodes by solving a local system involving the Gram matrix (the mass matrix + * of the shape functions). + * + * Once the stress tensor is projected to the nodes, the deviatoric part and the + * resulting von Mises scalar are computed for visualization via a color map. + * + * @tparam DataTypes The SOFA DataTypes (e.g. Vec3d) + * @tparam ElementType The geometric element type (e.g. Tetrahedron, Hexahedron) + */ +template +class VonMisesStress : public core::behavior::TopologyAccessor, public core::behavior::SingleStateAccessor +{ +public: + SOFA_CLASS2( + SOFA_TEMPLATE2(VonMisesStress, DataTypes, ElementType), + core::behavior::TopologyAccessor, + core::behavior::SingleStateAccessor); + +protected: + using FiniteElement = sofa::fem::FiniteElement; + + static constexpr sofa::Size spatial_dimensions = DataTypes::spatial_dimensions; + static constexpr sofa::Size NumberOfNodesInElement = ElementType::NumberOfNodes; + static constexpr sofa::Size NumberOfDofsInElement = NumberOfNodesInElement * spatial_dimensions; + static constexpr sofa::Size TopologicalDimension = FiniteElement::TopologicalDimension; + static constexpr sofa::Size NumberOfQuadraturePoints = FiniteElement::quadraturePoints().size(); + + // a stress tensor represented as a vector using the Voigt mapping + using StressVoigtVector = sofa::type::Vec, sofa::Real_t>; + + using DeformationGradient = sofa::type::Mat>; + +public: + + void init() override; + void draw(const core::visual::VisualParams* vparams) override; + void computeBBox(const core::ExecParams* params, bool onlyVisible) override; + + // A stress value for each node in an element + using LocalStressValues = std::array, NumberOfNodesInElement>; + Data> d_nodalStress; + + Data d_colorMap; + Data d_lighting; + + sofa::SingleLink, + sofa::BaseLink::FLAG_STOREPATH | sofa::BaseLink::FLAG_STRONGLINK> l_stressEvaluator; + + VonMisesStress(); + + void handleEvent(core::objectmodel::Event*) override; + +protected: + + void validateStressEvaluatorLink(); + + using ElementGramMatrix = sofa::type::Mat>; + sofa::type::vector m_elementInverseGramMatrices; + + struct PrecomputedData + { + sofa::type::Mat> jacobian { sofa::type::NOINIT }; + sofa::type::Mat> jacobianInv { sofa::type::NOINIT }; + Real_t detJacobian {}; + sofa::type::Mat> dN_dQ { sofa::type::NOINIT }; + }; + + sofa::type::vector> m_precomputedData; + + void precomputeData(); + + /** + * @brief Computes the inverse of the Gram matrix (integrated N^T * N) for each element. + * This matrix is used for the least-square projection of values from quadrature points to nodes. + */ + void calculateElementInverseGramMatrices(const auto& elements, sofa::type::vector& inverseGramMatrices); + + /** + * @brief Projects values evaluated at quadrature points to nodal values using least-squares. + * @param elementId Index of the element + * @param valuesAtQuadraturePoints Array of values evaluated at each quadrature point of the element + * @return Array of projected values at each node of the element + */ + std::array projectQuadraturePointValuesToNodes( + sofa::Size elementId, + const std::array& valuesAtQuadraturePoints) const; + + static StressVoigtVector deviatoricStress(const StressVoigtVector& sigma); + static Real_t vonMisesStress(const StressVoigtVector& deviatoricStress); + + core::visual::DrawElementColoredMesh m_renderer; +}; + +} // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl new file mode 100644 index 00000000000..1070d6758d4 --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/VonMisesStress.inl @@ -0,0 +1,369 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include +#include +#include +#include + +namespace sofa::component::solidmechanics::fem::elastic +{ + +template +VonMisesStress::VonMisesStress() + : d_nodalStress(initData(&d_nodalStress, sofa::type::vector{}, "nodalStress", + "Local nodal von Mises stress values")) + , d_colorMap(initData(&d_colorMap, sofa::helper::ColorMap(), "colorMap", "Color map")) + , d_lighting(initData(&d_lighting, true, "lighting", "If true, light is simulated on the mesh. Otherwise, no lighting effect.")) + , l_stressEvaluator(initLink("stressEvaluator", "The component in charge of evaluating the Cauchy stress.")) +{ + // This component must receive events + f_listening.setValue(true); +} + +template +void VonMisesStress::init() +{ + core::behavior::SingleStateAccessor::init(); + + if (!this->isComponentStateInvalid()) + { + this->validateTopology(); + } + + if (!this->isComponentStateInvalid()) + { + validateStressEvaluatorLink(); + } + + if (!this->isComponentStateInvalid()) + { + auto nodalStress = sofa::helper::getWriteOnlyAccessor(d_nodalStress); + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + nodalStress->resize(elements.size()); + + this->precomputeData(); + this->calculateElementInverseGramMatrices(elements, m_elementInverseGramMatrices); + } +} + +template +void VonMisesStress::handleEvent(core::objectmodel::Event* event) +{ + if (simulation::AnimateEndEvent::checkEventType(event)) + { + SCOPED_TIMER("vonMisesStress"); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + const auto nbElements = elements.size(); + + auto nodalStress = sofa::helper::getWriteOnlyAccessor(d_nodalStress); + nodalStress->clear(); + nodalStress->resize(nbElements); + + auto positionAccessor = this->mstate->readPositions(); + + helper::IotaView indices{static_cast(0ul), nbElements}; + + std::for_each( + indices.begin(), indices.end(), + [&](const auto elementId) + { + const auto& element = elements[elementId]; + + std::array, NumberOfNodesInElement> nodeCoordinatesInElement; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + nodeCoordinatesInElement[i] = positionAccessor[element[i]]; + } + + static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); + + std::array stressAtQuadraturePoints; + for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) + { + const auto& dN_dq_ref = gradients[q]; + const auto J_q = FiniteElement::Helper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); + + // Deformation Gradient F = J_curr * J_rest_inv + const DeformationGradient F = J_q * m_precomputedData[elementId][q].jacobianInv; + stressAtQuadraturePoints[q] = l_stressEvaluator->computeStress(F, elementId); + } + + // Least-square projection from quadrature points to nodes + const auto nodalStressInElement = projectQuadraturePointValuesToNodes(elementId, stressAtQuadraturePoints); + + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + nodalStress[elementId][i] = vonMisesStress(deviatoricStress(nodalStressInElement[i])); + } + + }); + } +} + +template +void VonMisesStress::validateStressEvaluatorLink() +{ + if (l_stressEvaluator.empty()) + { + msg_info() << "Link to a valid stress evaluator should be set to ensure right behavior. The first " + "stress evaluator found in current context will be used."; + l_stressEvaluator.set(this->getContext()->template get>()); + } + + if (l_stressEvaluator == nullptr) + { + msg_error() << "No stress evaluator component found at path: '" << this->l_stressEvaluator.getLinkedPath() + << "', nor in current context: " << this->getContext()->name + << ". Object must have a stress evaluator. " + << "The list of available stress evaluator components is: " + << sofa::core::ObjectFactory::getInstance() + ->listClassesDerivedFrom>(); + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + } +} + +template +void VonMisesStress::precomputeData() +{ + if (this->l_topology == nullptr) return; + if (this->mstate == nullptr) return; + + auto restPositionAccessor = this->mstate->readRestPositions(); + const auto& restPosition = restPositionAccessor.ref(); + + const auto& elements = FiniteElement::getElementSequence(*this->l_topology); + m_precomputedData.resize(elements.size()); + + static constexpr auto gradients = sofa::fem::FiniteElementHelper::gradientShapeFunctionAtQuadraturePoints(); + + for (std::size_t i = 0; i < elements.size(); ++i) + { + const auto& element = elements[i]; + std::array, NumberOfNodesInElement> nodeCoordinatesInElement; + for (sofa::Size n = 0; n < NumberOfNodesInElement; ++n) + nodeCoordinatesInElement[n] = restPosition[element[n]]; + + for (std::size_t j = 0; j < NumberOfQuadraturePoints; ++j) + { + const auto& dN_dq_ref = gradients[j]; + PrecomputedData& data = m_precomputedData[i][j]; + data.jacobian = sofa::fem::FiniteElementHelper::jacobianFromReferenceToPhysical(nodeCoordinatesInElement, dN_dq_ref); + data.jacobianInv = sofa::type::inverse(data.jacobian); + data.detJacobian = sofa::type::absGeneralizedDeterminant(data.jacobian); + + for (sofa::Size n = 0; n < NumberOfNodesInElement; ++n) + { + data.dN_dQ[n] = data.jacobianInv.multTranspose(dN_dq_ref[n]); + } + } + } +} + +template +void VonMisesStress::calculateElementInverseGramMatrices( + const auto& elements, sofa::type::vector& elementMassMatrices) +{ + const auto nbElements = elements.size(); + elementMassMatrices.resize(nbElements); + + auto positionAccessor = this->mstate->readPositions(); + + SCOPED_TIMER("elementMassMatrix"); + helper::IotaView indices{static_cast(0ul), nbElements}; + std::for_each( + indices.begin(), indices.end(), + [&](const auto elementId) + { + const auto& element = elements[elementId]; + auto& elementMassMatrix = elementMassMatrices[elementId]; + + std::array, NumberOfNodesInElement> nodeCoordinatesInElement; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + nodeCoordinatesInElement[i] = positionAccessor[element[i]]; + } + + for (const auto& [quadraturePoint, weight] : FiniteElement::quadraturePoints()) + { + // gradient of shape functions in the reference element evaluated at the quadrature + // point + const sofa::type::Mat> + dN_dq_ref = FiniteElement::gradientShapeFunctions(quadraturePoint); + + // jacobian of the mapping from the reference space to the physical space, evaluated + // at the quadrature point + sofa::type::Mat> + jacobian = FiniteElement::Helper::jacobianFromReferenceToPhysical( + nodeCoordinatesInElement, dN_dq_ref); + + const auto detJ = sofa::type::absGeneralizedDeterminant(jacobian); + + // shape functions in the reference element evaluated at the quadrature point + const auto N = FiniteElement::shapeFunctions(quadraturePoint); + + const auto NT_N = sofa::type::dyad(N, N); + + elementMassMatrix += (weight * detJ) * NT_N; + } + + sofa::type::invertMatrix(elementMassMatrix, elementMassMatrix); + }); +} + +template +auto VonMisesStress::projectQuadraturePointValuesToNodes( + sofa::Size elementId, + const std::array& valuesAtQuadraturePoints) const +-> std::array +{ + using Real = sofa::Real_t; + static constexpr auto quadraturePoints = FiniteElement::quadraturePoints(); + static constexpr auto voigtSize = sofa::type::NumberOfIndependentElements; + + std::array projectedNodalValues; + for (auto& val : projectedNodalValues) val.clear(); + + // Project each Voigt component independently: M * x_i = b_i + for (sofa::Size i = 0; i < voigtSize; ++i) + { + sofa::type::Vec b; + b.clear(); + + for (sofa::Size q = 0; q < NumberOfQuadraturePoints; ++q) + { + const auto& weight = quadraturePoints[q].second; + const auto detJ = m_precomputedData[elementId][q].detJacobian; + const auto N = FiniteElement::shapeFunctions(quadraturePoints[q].first); + + const Real val = valuesAtQuadraturePoints[q][i] * weight * detJ; + for (sofa::Size node = 0; node < NumberOfNodesInElement; ++node) + { + b[node] += N[node] * val; + } + } + + const auto x = m_elementInverseGramMatrices[elementId] * b; + for (sofa::Size node = 0; node < NumberOfNodesInElement; ++node) + { + projectedNodalValues[node][i] = x[node]; + } + } + + return projectedNodalValues; +} + +template +auto VonMisesStress::deviatoricStress(const StressVoigtVector& sigma) -> StressVoigtVector +{ + StressVoigtVector s = sigma; + + Real_t trace {}; + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + trace += sigma[type::tensorToVoigtIndex(i, i)]; + } + + static constexpr Real_t dim_inv = static_cast>(1.0 / spatial_dimensions); + + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + s[type::tensorToVoigtIndex(i, i)] -= dim_inv * trace; + } + + return s; +} + +template +Real_t +VonMisesStress::vonMisesStress(const StressVoigtVector& deviatoricStress) +{ + sofa::type::Mat> s; + for (sofa::Size i = 0; i < spatial_dimensions; ++i) + { + for (sofa::Size j = 0; j < spatial_dimensions; ++j) + { + s[i][j] = deviatoricStress[type::tensorToVoigtIndex(i, j)]; + } + } + + return std::sqrt(static_cast>(3) / 2 * sofa::type::trace(s.multTranspose(s))); +} + +template +void VonMisesStress::draw(const core::visual::VisualParams* vparams) +{ + const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); + + vparams->drawTool()->setLightingEnabled(d_lighting.getValue()); + + const auto nodalStress = d_nodalStress.getValue(); + if (nodalStress.empty()) + return; + + const auto positions = this->mstate->readPositions(); + const auto& colorMap = d_colorMap.getValue(); + + Real_t minStress = std::numeric_limits>::max(); + Real_t maxStress = std::numeric_limits>::lowest(); + + for (const auto& elementStress : nodalStress) + { + for (const auto& stress : elementStress) + { + minStress = std::min(minStress, stress); + maxStress = std::max(maxStress, stress); + } + } + + const auto evaluator = colorMap.getEvaluator(minStress, maxStress); + + sofa::type::vector> nodesColors; + for (const auto& elementStress : nodalStress) + { + std::array nodesColorsInElement; + for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i) + { + nodesColorsInElement[i] = (evaluator(elementStress[i])); + } + nodesColors.push_back(nodesColorsInElement); + } + + m_renderer.drawAllElements(vparams->drawTool(), positions.ref(), this->l_topology.get(), nodesColors); +} + +template +void VonMisesStress::computeBBox(const core::ExecParams* params, bool onlyVisible) +{ + SOFA_UNUSED(params); + SOFA_UNUSED(onlyVisible); + + if (!this->mstate) return; + + const auto bbox = this->mstate->computeBBox(); //this may compute twice the mstate bbox, but there is no way to determine if the bbox has already been computed + this->f_bbox.setValue(std::move(bbox)); +} + +} // namespace sofa::component::solidmechanics::fem::elastic 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 ccc0b9f98a6..91b94764ba5 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 @@ -67,6 +67,10 @@ struct trait sofa::Real_t>; using ElementGradient = sofa::type::Vec>; + + using DeformationGradient = sofa::type::Mat>; + + using StressVoigtVector = sofa::type::Vec, sofa::Real_t>; }; } 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 2c37e45f155..87d14c740ff 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 @@ -41,6 +41,7 @@ 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 registerVonMisesStress(sofa::core::ObjectFactory* factory); extern "C" { SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule(); @@ -81,6 +82,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) registerTriangularAnisotropicFEMForceField(factory); registerTriangularFEMForceField(factory); registerTriangularFEMForceFieldOptim(factory); + registerVonMisesStress(factory); } void init() diff --git a/Sofa/GL/src/sofa/gl/DrawToolGL.cpp b/Sofa/GL/src/sofa/gl/DrawToolGL.cpp index fbbe4087f13..7f2fbd705d4 100644 --- a/Sofa/GL/src/sofa/gl/DrawToolGL.cpp +++ b/Sofa/GL/src/sofa/gl/DrawToolGL.cpp @@ -979,6 +979,10 @@ void DrawToolGL::drawQuads(const std::vector &points, const type::RGBAColo void DrawToolGL::drawQuads(const std::vector &points, const std::vector& colors) { + glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE); + glEnable(GL_COLOR_MATERIAL); + setMaterial(colors[0]); + glBegin(GL_QUADS); { for (std::size_t i=0; i &points, const std::vector + static sofa::type::RGBAColor getColor(const ColorContainer& container, std::size_t elementId, std::size_t localId, std::size_t globalId) + { + if constexpr (std::is_same_v) + { + // Global approach + return container[globalId]; + } + else + { + // Local approach + return container[elementId][localId]; + } + } + template PositionType applyElementSpace(const PositionType& position, const PositionType& elementCenter) const { @@ -92,6 +107,60 @@ struct BaseDrawColoredMesh template struct DrawElementColoredMesh{}; +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Edge; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Triangle; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Quad; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + template<> struct SOFA_CORE_API DrawElementColoredMesh : public BaseDrawColoredMesh> @@ -139,7 +208,12 @@ struct SOFA_CORE_API DrawElementColoredMesh const auto p = this->applyElementSpace(position[vertexId], center); renderedPoints.push_back(p); - colors.push_back(nodesColors[vertexId]); + + // Find local ID of the vertex in the Tetrahedron to support local color mapping + std::size_t localId = 0; + for (; localId < 4; ++localId) if (element[localId] == vertexId) break; + + colors.push_back(this->getColor(nodesColors, i, localId, vertexId)); } } } @@ -149,6 +223,99 @@ struct SOFA_CORE_API DrawElementColoredMesh }; +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Hexahedron; + friend BaseDrawColoredMesh; + static constexpr std::size_t NumberQuadsInHexahedron = 6; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + const auto& elements = topology->getHexahedra(); + const auto& facets = topology->getQuads(); + + if(facets.empty()) + { + msg_error_once("DrawElementColoredMesh") << "Drawing hexahedra needs the associated quads in the topology."; + return; + } + + renderedPoints.clear(); + colors.clear(); + + for (auto i : elementIndices) + { + const auto& element = elements[i]; + const auto& facetsInElement = topology->getQuadsInHexahedron(i); + const auto center = this->elementCenter(position, element); + + for (std::size_t j = 0; j < NumberQuadsInHexahedron; ++j) + { + const auto faceId = facetsInElement[j]; + for (const auto vertexId : facets[faceId]) + { + const auto p = this->applyElementSpace(position[vertexId], center); + + renderedPoints.push_back(p); + + // Find local ID of the vertex in the Hexahedron to support local color mapping + std::size_t localId = 0; + for (; localId < 8; ++localId) if (element[localId] == vertexId) break; + + colors.push_back(this->getColor(nodesColors, i, localId, vertexId)); + } + } + } + + drawTool->drawQuads(renderedPoints, colors); + } + +}; + +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Prism; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + +template<> +struct SOFA_CORE_API DrawElementColoredMesh + : public BaseDrawColoredMesh> +{ + using ElementType = sofa::geometry::Pyramid; + friend BaseDrawColoredMesh; + + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& nodesColors) + { + } +}; + class SOFA_CORE_API DrawColoredMesh { public: @@ -159,6 +326,7 @@ class SOFA_CORE_API DrawColoredMesh void draw(sofa::helper::visual::DrawTool* drawTool, const PositionContainer& position, const ColorContainer& nodesColors, sofa::core::topology::BaseMeshTopology* topology) { std::get>(m_meshes).drawAllElements(drawTool, position, topology, nodesColors); + std::get>(m_meshes).drawAllElements(drawTool, position, topology, nodesColors); } private: @@ -166,8 +334,8 @@ class SOFA_CORE_API DrawColoredMesh // DrawElementColoredMesh, // DrawElementColoredMesh, // DrawElementColoredMesh, - DrawElementColoredMesh - // DrawElementColoredMesh, + DrawElementColoredMesh, + DrawElementColoredMesh // DrawElementColoredMesh, // DrawElementColoredMesh > m_meshes; diff --git a/examples/Component/SolidMechanics/FEM/VonMisesStress.scn b/examples/Component/SolidMechanics/FEM/VonMisesStress.scn new file mode 100644 index 00000000000..afb3186a691 --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/VonMisesStress.scn @@ -0,0 +1,37 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +