Skip to content

Lazy VonMisesStress computation #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
damienmarchal wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Lazy VonMisesStress computation #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
23e44cc
Move the template specialization #includes out of the namespaces.
damienmarchal Aug 9, 2025
513299b
FIXUP
damienmarchal Aug 9, 2025
8a95d01
Attempt to move the template specialization out of the common header …
damienmarchal Aug 9, 2025
aee7b43
Fixes needed to compile with clang15
damienmarchal Aug 9, 2025
19405ad
Use the DDG update mechanism to compute the vonMises stree only when …
damienmarchal Aug 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/Elasticity/component/ElementLinearSmallStrainFEMForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ namespace elasticity
template <class DataTypes, class ElementType>
void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::selectFEMTypes()
{
this->addLinearFEMType<ElementType>();
this->template addLinearFEMType<ElementType>();
}

} // namespace elasticity
2 changes: 1 addition & 1 deletion src/Elasticity/component/LinearSmallStrainFEMForceField.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ class LinearSmallStrainFEMForceField : public sofa::core::behavior::ForceField<D
SOFA_TEMPLATE(sofa::core::behavior::ForceField, DataTypes));

private:
using Inherit1::mstate;
using DataVecCoord = sofa::DataVecDeriv_t<DataTypes>;
using DataVecDeriv = sofa::DataVecDeriv_t<DataTypes>;
using Real = sofa::Real_t<DataTypes>;
Expand All @@ -30,7 +31,6 @@ class LinearSmallStrainFEMForceField : public sofa::core::behavior::ForceField<D

sofa::Data<Real> d_poissonRatio;
sofa::Data<Real> d_youngModulus;
sofa::Data<bool> d_computeVonMisesStress;
sofa::Data<sofa::type::vector<Real> > d_vonMisesStressValues;

void init() override;
Expand Down
41 changes: 25 additions & 16 deletions src/Elasticity/component/LinearSmallStrainFEMForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ LinearSmallStrainFEMForceField<DataTypes>::LinearSmallStrainFEMForceField()
: l_topology(initLink("topology", "Link to a topology containing elements"))
, d_poissonRatio(initData(&d_poissonRatio, static_cast<Real>(0.45), "poissonRatio", "Poisson's ratio"))
, d_youngModulus(initData(&d_youngModulus, static_cast<Real>(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";
Expand All @@ -31,6 +30,31 @@ void LinearSmallStrainFEMForceField<DataTypes>::init()
{
sofa::core::behavior::ForceField<DataTypes>::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();
Expand Down Expand Up @@ -74,21 +98,6 @@ void LinearSmallStrainFEMForceField<DataTypes>::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 <class DataTypes>
Expand Down
7 changes: 7 additions & 0 deletions src/Elasticity/finiteelement/FiniteElement[all].h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
#pragma once

#include <Elasticity/finiteelement/FiniteElement[Edge].h>
#include <Elasticity/finiteelement/FiniteElement[Hexahedron].h>
#include <Elasticity/finiteelement/FiniteElement[Quad].h>
#include <Elasticity/finiteelement/FiniteElement[Tetrahedron].h>
#include <Elasticity/finiteelement/FiniteElement[Triangle].h>
2 changes: 1 addition & 1 deletion src/Elasticity/impl/CorotationalFEM.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#pragma once

#include <Elasticity/impl/LinearFEM.h>
#include <Elasticity/impl/LinearFEM[all].h>

namespace elasticity
{
Expand Down
9 changes: 2 additions & 7 deletions src/Elasticity/impl/LinearFEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,15 +1,10 @@
#define ELASTICITY_LINEARFEM_CPP

#include <Elasticity/impl/LinearFEM.inl>
#include <Elasticity/config.h>
#include <Elasticity/finiteelement/FiniteElement[all].h>
#include <Elasticity/impl/LinearFEM.inl>
#include <sofa/defaulttype/VecTypes.h>

#include <Elasticity/finiteelement/FiniteElement[Edge].h>
#include <Elasticity/finiteelement/FiniteElement[Hexahedron].h>
#include <Elasticity/finiteelement/FiniteElement[Quad].h>
#include <Elasticity/finiteelement/FiniteElement[Tetrahedron].h>
#include <Elasticity/finiteelement/FiniteElement[Triangle].h>

namespace elasticity
{

Expand Down
17 changes: 0 additions & 17 deletions src/Elasticity/impl/LinearFEM.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,21 +99,4 @@ class LinearFEM : public BaseLinearFEM<DataTypes>
const sofa::type::vector<ElementStiffness>& stiffnessMatrices() const;
};

#if !defined(ELASTICITY_LINEARFEM_CPP)
#include <Elasticity/finiteelement/FiniteElement[Edge].h>
#include <Elasticity/finiteelement/FiniteElement[Hexahedron].h>
#include <Elasticity/finiteelement/FiniteElement[Quad].h>
#include <Elasticity/finiteelement/FiniteElement[Tetrahedron].h>
#include <Elasticity/finiteelement/FiniteElement[Triangle].h>

extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec2Types, sofa::geometry::Edge>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Edge>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec2Types, sofa::geometry::Triangle>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Triangle>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec2Types, sofa::geometry::Quad>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
#endif
}
20 changes: 20 additions & 0 deletions src/Elasticity/impl/LinearFEM[all].h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
#pragma once

#include <Elasticity/config.h>
#include <Elasticity/impl/LinearFEM.h>
#include <Elasticity/finiteelement/FiniteElement[all].h>

namespace elasticity
{
#if !defined(ELASTICITY_LINEARFEM_CPP)
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec2Types, sofa::geometry::Edge>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Edge>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec2Types, sofa::geometry::Triangle>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Triangle>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec2Types, sofa::geometry::Quad>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
extern template class ELASTICITY_API LinearFEM<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
#endif
}
2 changes: 1 addition & 1 deletion src/Elasticity/impl/MatrixTools.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
namespace elasticity
{

template <size_t N, typename real>
template <sofa::Size N, typename real>
real determinantSquareMatrix(const sofa::type::Mat<N, N, real>& mat)
{
if constexpr (N == 1)
Expand Down
12 changes: 6 additions & 6 deletions tests/CorotationalFEM_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,12 @@ TEST(CorotationalFEM, deformationGradientTetra)
{
using FE = FiniteElement<sofa::geometry::Tetrahedron, sofa::defaulttype::Vec3Types>;

constexpr std::array<sofa::type::Vec3, 4> tetraNodesCoordinates({
{0, 0, 0},
{1, 0, 0},
{0, 1, 0},
{0, 0, 1}
});
constexpr std::array<sofa::type::Vec3, 4> 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<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>::computeCentroid(
Expand Down
20 changes: 10 additions & 10 deletions tests/HexahedronLinearSmallStrainFEMForceField_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,16 +19,16 @@ TEST(HexahedronLinearSmallStrainFEMForceField, jacobian)
using Force = ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
using FE = FiniteElement<sofa::geometry::Hexahedron, sofa::defaulttype::Vec3Types>;

constexpr std::array<sofa::type::Vec3, 8> 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<sofa::type::Vec3, 8> 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())
{
Expand Down
12 changes: 6 additions & 6 deletions tests/TetrahedronLinearSmallStrainFEMForceField_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,12 +135,12 @@ TEST(TET4LinearSmallStrainFEMForceField, jacobian)
{
using FE = FiniteElement<sofa::geometry::Tetrahedron, sofa::defaulttype::Vec3Types>;

constexpr std::array<sofa::type::Vec3, 4> tetraNodesCoordinates({
{0, 0, 0},
{1, 0, 0},
{0, 1, 0},
{0, 0, 1}
});
constexpr std::array<sofa::type::Vec3, 4> 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);
Expand Down