From 23e44cca40f9740d3d6c9ca1af99600fe806cf6e Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Sat, 9 Aug 2025 07:54:36 +0200 Subject: [PATCH 1/5] Move the template specialization #includes out of the namespaces. (otherwise it produced nested elasticity namespaces) --- src/Elasticity/impl/LinearFEM.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Elasticity/impl/LinearFEM.h b/src/Elasticity/impl/LinearFEM.h index 0bc6472..e5c32ff 100644 --- a/src/Elasticity/impl/LinearFEM.h +++ b/src/Elasticity/impl/LinearFEM.h @@ -9,6 +9,12 @@ #include +#include +#include +#include +#include +#include + namespace elasticity { @@ -100,12 +106,6 @@ class LinearFEM : public BaseLinearFEM }; #if !defined(ELASTICITY_LINEARFEM_CPP) -#include -#include -#include -#include -#include - extern template class ELASTICITY_API LinearFEM; extern template class ELASTICITY_API LinearFEM; extern template class ELASTICITY_API LinearFEM; From 513299b37690b5e77c396d6bf5bd1bee73ab558d Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Sat, 9 Aug 2025 07:57:30 +0200 Subject: [PATCH 2/5] FIXUP --- src/Elasticity/impl/LinearFEM.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Elasticity/impl/LinearFEM.h b/src/Elasticity/impl/LinearFEM.h index e5c32ff..eebe7a0 100644 --- a/src/Elasticity/impl/LinearFEM.h +++ b/src/Elasticity/impl/LinearFEM.h @@ -9,11 +9,13 @@ #include +#if !defined(ELASTICITY_LINEARFEM_CPP) #include #include #include #include #include +#endif namespace elasticity { From 8a95d015ad58fbac143cb1c5344bff4b68b41b2a Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Sat, 9 Aug 2025 08:11:24 +0200 Subject: [PATCH 3/5] Attempt to move the template specialization out of the common header so it is no more hardcoded. --- CMakeLists.txt | 2 ++ .../finiteelement/FiniteElement[all].h | 7 +++++++ src/Elasticity/impl/CorotationalFEM.h | 2 +- src/Elasticity/impl/LinearFEM.cpp | 9 ++------- src/Elasticity/impl/LinearFEM.h | 19 ------------------ src/Elasticity/impl/LinearFEM[all].h | 20 +++++++++++++++++++ 6 files changed, 32 insertions(+), 27 deletions(-) create mode 100644 src/Elasticity/finiteelement/FiniteElement[all].h create mode 100644 src/Elasticity/impl/LinearFEM[all].h diff --git a/CMakeLists.txt b/CMakeLists.txt index a648be7..a1cb502 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,6 +9,7 @@ set(HEADER_FILES src/Elasticity/init.h src/Elasticity/finiteelement/FiniteElement.h + src/Elasticity/finiteelement/FiniteElement[all].h src/Elasticity/finiteelement/FiniteElement[Edge].h src/Elasticity/finiteelement/FiniteElement[Hexahedron].h src/Elasticity/finiteelement/FiniteElement[Quad].h @@ -26,6 +27,7 @@ set(HEADER_FILES src/Elasticity/impl/CorotationalFEM.inl src/Elasticity/impl/ElasticityTensor.h src/Elasticity/impl/ElementStiffnessMatrix.h + src/Elasticity/impl/LinearFEM[all].h src/Elasticity/impl/LinearFEM.h src/Elasticity/impl/LinearFEM.inl src/Elasticity/impl/MatrixTools.h diff --git a/src/Elasticity/finiteelement/FiniteElement[all].h b/src/Elasticity/finiteelement/FiniteElement[all].h new file mode 100644 index 0000000..8209736 --- /dev/null +++ b/src/Elasticity/finiteelement/FiniteElement[all].h @@ -0,0 +1,7 @@ +#pragma once + +#include +#include +#include +#include +#include diff --git a/src/Elasticity/impl/CorotationalFEM.h b/src/Elasticity/impl/CorotationalFEM.h index 6c95574..a7f6f77 100644 --- a/src/Elasticity/impl/CorotationalFEM.h +++ b/src/Elasticity/impl/CorotationalFEM.h @@ -1,6 +1,6 @@ #pragma once -#include +#include namespace elasticity { diff --git a/src/Elasticity/impl/LinearFEM.cpp b/src/Elasticity/impl/LinearFEM.cpp index dd025f6..45bef2d 100644 --- a/src/Elasticity/impl/LinearFEM.cpp +++ b/src/Elasticity/impl/LinearFEM.cpp @@ -1,15 +1,10 @@ #define ELASTICITY_LINEARFEM_CPP -#include #include +#include +#include #include -#include -#include -#include -#include -#include - namespace elasticity { diff --git a/src/Elasticity/impl/LinearFEM.h b/src/Elasticity/impl/LinearFEM.h index eebe7a0..eaa2d91 100644 --- a/src/Elasticity/impl/LinearFEM.h +++ b/src/Elasticity/impl/LinearFEM.h @@ -9,14 +9,6 @@ #include -#if !defined(ELASTICITY_LINEARFEM_CPP) -#include -#include -#include -#include -#include -#endif - namespace elasticity { @@ -107,15 +99,4 @@ class LinearFEM : public BaseLinearFEM const sofa::type::vector& stiffnessMatrices() const; }; -#if !defined(ELASTICITY_LINEARFEM_CPP) -extern template class ELASTICITY_API LinearFEM; -extern template class ELASTICITY_API LinearFEM; -extern template class ELASTICITY_API LinearFEM; -extern template class ELASTICITY_API LinearFEM; -extern template class ELASTICITY_API LinearFEM; -extern template class ELASTICITY_API LinearFEM; -extern template class ELASTICITY_API LinearFEM; -extern template class ELASTICITY_API LinearFEM; -extern template class ELASTICITY_API LinearFEM; -#endif } diff --git a/src/Elasticity/impl/LinearFEM[all].h b/src/Elasticity/impl/LinearFEM[all].h new file mode 100644 index 0000000..351e680 --- /dev/null +++ b/src/Elasticity/impl/LinearFEM[all].h @@ -0,0 +1,20 @@ +#pragma once + +#include +#include +#include + +namespace elasticity +{ +#if !defined(ELASTICITY_LINEARFEM_CPP) +extern template class ELASTICITY_API LinearFEM; +extern template class ELASTICITY_API LinearFEM; +extern template class ELASTICITY_API LinearFEM; +extern template class ELASTICITY_API LinearFEM; +extern template class ELASTICITY_API LinearFEM; +extern template class ELASTICITY_API LinearFEM; +extern template class ELASTICITY_API LinearFEM; +extern template class ELASTICITY_API LinearFEM; +extern template class ELASTICITY_API LinearFEM; +#endif +} From aee7b435f99352d9d4da6aacb74697c803377256 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Sat, 9 Aug 2025 08:20:46 +0200 Subject: [PATCH 4/5] Fixes needed to compile with clang15 --- .../ElementLinearSmallStrainFEMForceField.inl | 2 +- src/Elasticity/impl/MatrixTools.h | 2 +- tests/CorotationalFEM_test.cpp | 12 +++++------ ...ronLinearSmallStrainFEMForceField_test.cpp | 20 +++++++++---------- ...ronLinearSmallStrainFEMForceField_test.cpp | 12 +++++------ 5 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/Elasticity/component/ElementLinearSmallStrainFEMForceField.inl b/src/Elasticity/component/ElementLinearSmallStrainFEMForceField.inl index d77d690..f4d664c 100644 --- a/src/Elasticity/component/ElementLinearSmallStrainFEMForceField.inl +++ b/src/Elasticity/component/ElementLinearSmallStrainFEMForceField.inl @@ -10,7 +10,7 @@ namespace elasticity template void ElementLinearSmallStrainFEMForceField::selectFEMTypes() { - this->addLinearFEMType(); + this->template addLinearFEMType(); } } // namespace elasticity diff --git a/src/Elasticity/impl/MatrixTools.h b/src/Elasticity/impl/MatrixTools.h index 0d3540a..7e28a6f 100644 --- a/src/Elasticity/impl/MatrixTools.h +++ b/src/Elasticity/impl/MatrixTools.h @@ -7,7 +7,7 @@ namespace elasticity { -template +template real determinantSquareMatrix(const sofa::type::Mat& mat) { if constexpr (N == 1) diff --git a/tests/CorotationalFEM_test.cpp b/tests/CorotationalFEM_test.cpp index 9d02b96..7fa4820 100644 --- a/tests/CorotationalFEM_test.cpp +++ b/tests/CorotationalFEM_test.cpp @@ -24,12 +24,12 @@ TEST(CorotationalFEM, deformationGradientTetra) { using FE = FiniteElement; - constexpr std::array tetraNodesCoordinates({ - {0, 0, 0}, - {1, 0, 0}, - {0, 1, 0}, - {0, 0, 1} - }); + constexpr std::array tetraNodesCoordinates{{ + {0_sreal, 0_sreal, 0_sreal}, + {1_sreal, 0_sreal, 0_sreal}, + {0_sreal, 1_sreal, 0_sreal}, + {0_sreal, 0_sreal, 1_sreal} + }}; const auto q = FE::quadraturePoints(); static const auto centroid = CorotationalFEM::computeCentroid( diff --git a/tests/HexahedronLinearSmallStrainFEMForceField_test.cpp b/tests/HexahedronLinearSmallStrainFEMForceField_test.cpp index ec196fc..9e05578 100644 --- a/tests/HexahedronLinearSmallStrainFEMForceField_test.cpp +++ b/tests/HexahedronLinearSmallStrainFEMForceField_test.cpp @@ -19,16 +19,16 @@ TEST(HexahedronLinearSmallStrainFEMForceField, jacobian) using Force = ElementLinearSmallStrainFEMForceField; using FE = FiniteElement; - constexpr std::array hexaNodesCoordinates({ - {-1, -1, -1}, - {1, -1, -1}, - {1, 1, -1}, - {-1, 1, -1}, - {-1, -1, 1}, - {1, -1, 1}, - {1, 1, 1}, - {-1, 1, 1} - }); + constexpr std::array hexaNodesCoordinates{{ + {-1_sreal, -1_sreal, -1_sreal}, + {1_sreal, -1_sreal, -1_sreal}, + {1_sreal, 1_sreal, -1_sreal}, + {-1_sreal, 1_sreal, -1_sreal}, + {-1_sreal, -1_sreal, 1_sreal}, + {1_sreal, -1_sreal, 1_sreal}, + {1_sreal, 1_sreal, 1_sreal}, + {-1_sreal, 1_sreal, 1_sreal} + }}; for (const auto& [q, w] : FE::quadraturePoints()) { diff --git a/tests/TetrahedronLinearSmallStrainFEMForceField_test.cpp b/tests/TetrahedronLinearSmallStrainFEMForceField_test.cpp index 8dd34b4..489c5d1 100644 --- a/tests/TetrahedronLinearSmallStrainFEMForceField_test.cpp +++ b/tests/TetrahedronLinearSmallStrainFEMForceField_test.cpp @@ -135,12 +135,12 @@ TEST(TET4LinearSmallStrainFEMForceField, jacobian) { using FE = FiniteElement; - constexpr std::array tetraNodesCoordinates({ - {0, 0, 0}, - {1, 0, 0}, - {0, 1, 0}, - {0, 0, 1} - }); + constexpr std::array tetraNodesCoordinates{{ + {0_sreal, 0_sreal, 0_sreal}, + {1_sreal, 0_sreal, 0_sreal}, + {0_sreal, 1_sreal, 0_sreal}, + {0_sreal, 0_sreal, 1_sreal} + }}; const auto q = FE::quadraturePoints(); const auto dN_dq_ref = FE::gradientShapeFunctions(q[0].first); From 19405ad8e55a8445a8af80f386f2bc2113ede728 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Sat, 9 Aug 2025 09:04:48 +0200 Subject: [PATCH 5/5] Use the DDG update mechanism to compute the vonMises stree only when needed. (And consequently removes the d_computeVonMises data) NB: Having a mstate pointer in the callback is not really satisfactory, any better suggestion welcome. I have a vague remembering that a long time ago Bruno Marques brings to link the callback mechanisme which could solve the problem. Signed-off-by: Damien Marchal --- .../LinearSmallStrainFEMForceField.h | 2 +- .../LinearSmallStrainFEMForceField.inl | 41 +++++++++++-------- 2 files changed, 26 insertions(+), 17 deletions(-) diff --git a/src/Elasticity/component/LinearSmallStrainFEMForceField.h b/src/Elasticity/component/LinearSmallStrainFEMForceField.h index b78bd36..fa32da2 100644 --- a/src/Elasticity/component/LinearSmallStrainFEMForceField.h +++ b/src/Elasticity/component/LinearSmallStrainFEMForceField.h @@ -17,6 +17,7 @@ class LinearSmallStrainFEMForceField : public sofa::core::behavior::ForceField; using DataVecDeriv = sofa::DataVecDeriv_t; using Real = sofa::Real_t; @@ -30,7 +31,6 @@ class LinearSmallStrainFEMForceField : public sofa::core::behavior::ForceField d_poissonRatio; sofa::Data d_youngModulus; - sofa::Data d_computeVonMisesStress; sofa::Data > d_vonMisesStressValues; void init() override; diff --git a/src/Elasticity/component/LinearSmallStrainFEMForceField.inl b/src/Elasticity/component/LinearSmallStrainFEMForceField.inl index 7e6c769..26be219 100644 --- a/src/Elasticity/component/LinearSmallStrainFEMForceField.inl +++ b/src/Elasticity/component/LinearSmallStrainFEMForceField.inl @@ -16,7 +16,6 @@ LinearSmallStrainFEMForceField::LinearSmallStrainFEMForceField() : l_topology(initLink("topology", "Link to a topology containing elements")) , d_poissonRatio(initData(&d_poissonRatio, static_cast(0.45), "poissonRatio", "Poisson's ratio")) , d_youngModulus(initData(&d_youngModulus, static_cast(1e6), "youngModulus", "Young's modulus")) - , d_computeVonMisesStress(initData(&d_computeVonMisesStress, true, "computeVonMisesStress", "Compute Von Mises stress")) , d_vonMisesStressValues(initData(&d_vonMisesStressValues, "vonMisesStressValues", "Von Mises stress values")) { static std::string groupName = "Mechanical parameter"; @@ -31,6 +30,31 @@ void LinearSmallStrainFEMForceField::init() { sofa::core::behavior::ForceField::init(); + // We have a mstate here + auto position = mstate->write (sofa::core::vec_id::write_access::position); + auto restPosition = mstate->write (sofa::core::vec_id::write_access::restPosition); + sofa::core::objectmodel::Base::addUpdateCallback("computeVonMisesStress", {position,restPosition}, [this,position](const sofa::core::DataTracker& ) + { + auto positionAccessor = mstate->readPositions(); + auto restPositionAccessor = mstate->readRestPositions(); + + std::cout << "COMPUTING THE VON MISE STRESS FOR [" << this->getPathName() << "] AS POSITION OR REST POSITION CHANGED " << position->getCounter() << std::endl; + + m_vonMisesStressContainer.clear(); + m_vonMisesStressContainer.resize(positionAccessor.size()); + + for (const auto& finiteElement : m_finiteElements) + { + finiteElement->computeVonMisesStress(m_vonMisesStressContainer, positionAccessor.ref(), restPositionAccessor.ref()); + } + + m_vonMisesStressContainer.accept(); + + d_vonMisesStressValues.setValue(m_vonMisesStressContainer.getStressValues()); + return sofa::core::objectmodel::ComponentState::Valid; + }, {&d_vonMisesStressValues}); + + if (!this->isComponentStateInvalid()) { this->validateTopology(); @@ -74,21 +98,6 @@ void LinearSmallStrainFEMForceField::addForce( { finiteElement->addForce(forceAccessor.wref(), positionAccessor.ref(), restPositionAccessor.ref()); } - - if (d_computeVonMisesStress.getValue()) - { - m_vonMisesStressContainer.clear(); - m_vonMisesStressContainer.resize(positionAccessor.size()); - - for (const auto& finiteElement : m_finiteElements) - { - finiteElement->computeVonMisesStress(m_vonMisesStressContainer, positionAccessor.ref(), restPositionAccessor.ref()); - } - - m_vonMisesStressContainer.accept(); - - d_vonMisesStressValues.setValue(m_vonMisesStressContainer.getStressValues()); - } } template