Skip to content

[FEM.Elastic] Cache stiffness matrices before computations #6069

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
fredroy wants to merge 2 commits into
base: master
Choose a base branch
from
+28 −9
Open

[FEM.Elastic] Cache stiffness matrices before computations #6069

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4d81580
fix: getExecutionPolicy did not use the provided data
fredroy Apr 3, 2026
941aec2
Improve parallel performance of element force derivative computation
fredroy Apr 3, 2026
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
12 changes: 12 additions & 0 deletions ...EM/Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,10 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f
using StrainDisplacement = typename trait::StrainDisplacement;
using Real = typename trait::Real;

public:
using AssembledStiffnessMatrix = sofa::type::Mat<
trait::NumberOfDofsInElement, trait::NumberOfDofsInElement, Real>;

protected:

BaseElementLinearFEMForceField();
Expand All @@ -70,6 +74,14 @@ class BaseElementLinearFEMForceField : public sofa::component::solidmechanics::f
* List of precomputed element stiffness matrices
*/
sofa::Data<sofa::type::vector<ElementStiffness> > d_elementStiffness;

/**
* Contiguous buffer of assembled stiffness matrices (one per element).
* Extracted from d_elementStiffness for cache-friendly access in the hot path.
* This avoids loading the full FactorizedElementStiffness struct (~15 KB per element)
* when only the assembled matrix (~4.6 KB) is needed.
*/
sofa::type::vector<AssembledStiffnessMatrix> m_assembledStiffnessMatrices;
};

#if !defined(ELASTICITY_COMPONENT_BASE_ELEMENT_LINEAR_FEM_FORCEFIELD_CPP)
Expand Down
7 changes: 7 additions & 0 deletions .../Elastic/src/sofa/component/solidmechanics/fem/elastic/BaseElementLinearFEMForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,13 @@ void BaseElementLinearFEMForceField<DataTypes, ElementType>::precomputeElementSt
const std::array<sofa::Coord_t<DataTypes>, trait::NumberOfNodesInElement> nodesCoordinates = extractNodesVectorFromGlobalVector(element, restPositionAccessor.ref());
elementStiffness[elementId] = integrate<DataTypes, ElementType, trait::matrixVectorProductType>(nodesCoordinates, elasticityTensor);
});

// Extract assembled matrices into a contiguous buffer for cache-friendly access
m_assembledStiffnessMatrices.resize(elements.size());
for (std::size_t i = 0; i < elements.size(); ++i)
{
m_assembledStiffnessMatrices[i] = elementStiffness[i].getAssembledMatrix();
}
}

}
8 changes: 4 additions & 4 deletions ...lastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsFo
{
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
auto restPositionAccessor = this->mstate->readRestPositions();
auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
const auto& assembledMatrices = this->m_assembledStiffnessMatrices;

for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
{
Expand All @@ -112,7 +112,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsFo
transformedDisplacement = elementRotation.multTranspose(elementNodesCoordinates[j] - t) - (restElementNodesCoordinates[j] - t0);
}

const auto& stiffnessMatrix = elementStiffness[elementId];
const auto& stiffnessMatrix = assembledMatrices[elementId];

auto& elementForce = elementForces[elementId];
elementForce = stiffnessMatrix * displacement;
Expand All @@ -134,7 +134,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsFo
const sofa::VecDeriv_t<DataTypes>& nodeDx)
{
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
const auto& assembledMatrices = this->m_assembledStiffnessMatrices;

for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
{
Expand All @@ -150,7 +150,7 @@ void ElementCorotationalFEMForceField<DataTypes, ElementType>::computeElementsFo
rotated_dx = elementRotation.multTranspose(nodeDx[element[n]]);
}

const auto& stiffnessMatrix = elementStiffness[elementId];
const auto& stiffnessMatrix = assembledMatrices[elementId];

auto& df = elementForcesDeriv[elementId];
df = stiffnessMatrix * element_dx;
Expand Down
8 changes: 4 additions & 4 deletions ...c/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,12 +49,12 @@ void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::computeEleme
{
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
auto restPositionAccessor = this->mstate->readRestPositions();
auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
const auto& assembledMatrices = this->m_assembledStiffnessMatrices;

for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
{
const auto& element = elements[elementId];
const auto& stiffnessMatrix = elementStiffness[elementId];
const auto& stiffnessMatrix = assembledMatrices[elementId];

typename trait::ElementDisplacement displacement{ sofa::type::NOINIT };

Expand All @@ -79,12 +79,12 @@ void ElementLinearSmallStrainFEMForceField<DataTypes, ElementType>::computeEleme
const sofa::VecDeriv_t<DataTypes>& nodeDx)
{
const auto& elements = trait::FiniteElement::getElementSequence(*this->l_topology);
auto elementStiffness = sofa::helper::getReadAccessor(this->d_elementStiffness);
const auto& assembledMatrices = this->m_assembledStiffnessMatrices;

for (std::size_t elementId = range.start; elementId < range.end; ++elementId)
{
const auto& element = elements[elementId];
const auto& stiffnessMatrix = elementStiffness[elementId];
const auto& stiffnessMatrix = assembledMatrices[elementId];

const std::array<sofa::Coord_t<DataTypes>, trait::NumberOfNodesInElement> elementNodesDx =
extractNodesVectorFromGlobalVector(element, nodeDx);
Expand Down
2 changes: 1 addition & 1 deletion ...olidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FEMForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,7 @@ template <class DataTypes, class ElementType>
sofa::simulation::ForEachExecutionPolicy FEMForceField<DataTypes, ElementType>::getExecutionPolicy(
const sofa::Data<ComputeStrategy>& strategy) const
{
auto computeForceStrategyAccessor = sofa::helper::getReadAccessor(d_computeForceStrategy);
auto computeForceStrategyAccessor = sofa::helper::getReadAccessor(strategy);
const auto& computeForceStrategy = computeForceStrategyAccessor->key();

return (computeForceStrategy == parallelComputeStrategy)
Expand Down
Loading