Skip to content

Commit

Permalink
[ForceField] add buildStiffnessMatrix to HyperReducedTetrahedralCorot…
Browse files Browse the repository at this point in the history 
…ationalFEMForceField (#92)

* [ForceField] add buildStiffnessMatrix to HyperReducedTetrahedralCorotationalFEMForceField

* Clean: remove remaining & normally useless addKToMatrin
  • Loading branch information
@VannesteFelix
VannesteFelix authored Dec 18, 2023
1 parent 7e9224e commit a798e13
Show file tree
Hide file tree
Showing 8 changed files with 37 additions and 281 deletions.
1 change: 0 additions & 1 deletion src/ModelOrderReduction/component/forcefield/HyperReducedHexahedronFEMForceField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,6 @@ class HyperReducedHexahedronFEMForceField : public virtual HexahedronFEMForceFie
using InheritForceField::getPotentialEnergy;
// getPotentialEnergy is implemented for polar method

void addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override;
void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /* matrix */) override;

void draw(const core::visual::VisualParams* vparams) override;
Expand Down
68 changes: 2 additions & 66 deletions src/ModelOrderReduction/component/forcefield/HyperReducedHexahedronFEMForceField.inl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -388,70 +388,6 @@ template<class DataTypes>
/////////////////////////////////////////////////


template<class DataTypes>
void HyperReducedHexahedronFEMForceField<DataTypes>::addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix)
{
// Build Matrix Block for this ForceField
int i,j,n1, n2, e;

typename VecElement::const_iterator it, it0;

Index node1, node2;

sofa::core::behavior::MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(this->mstate);

it0=this->getIndexedElements()->begin();
size_t nbElementsConsidered;
if (!d_performECSW.getValue()){
nbElementsConsidered = this->getIndexedElements()->size();
}
else
{
nbElementsConsidered = m_RIDsize;
}
for( unsigned int numElem = 0 ; numElem<nbElementsConsidered ;++numElem)
{
if (!d_performECSW.getValue()){
e = numElem;
}
else
{
e = reducedIntegrationDomain(numElem);
}
it = it0 + e;

const ElementStiffness &Ke = _elementStiffnesses.getValue()[e];
// const Transformation& Rt = _rotations[e];
// Transformation R; R.transpose(Rt);

Transformation Rot = this->getElementRotation(e);
Real kFactorTimesWeight;
if (!d_performECSW.getValue())
kFactorTimesWeight = (Real)mparams->kFactorIncludingRayleighDamping(this->rayleighStiffness.getValue());
else
kFactorTimesWeight = weights(e)*(Real)mparams->kFactorIncludingRayleighDamping(this->rayleighStiffness.getValue());


// find index of node 1
for (n1=0; n1<8; n1++)
{
node1 = (*it)[n1];
// find index of node 2
for (n2=0; n2<8; n2++)
{
node2 = (*it)[n2];

Mat33 tmp = Rot.multTranspose( Mat33(Coord(Ke[3*n1+0][3*n2+0],Ke[3*n1+0][3*n2+1],Ke[3*n1+0][3*n2+2]),
Coord(Ke[3*n1+1][3*n2+0],Ke[3*n1+1][3*n2+1],Ke[3*n1+1][3*n2+2]),
Coord(Ke[3*n1+2][3*n2+0],Ke[3*n1+2][3*n2+1],Ke[3*n1+2][3*n2+2])) ) * Rot;
for(i=0; i<3; i++)
for (j=0; j<3; j++)
r.matrix->add(r.offset+3*node1+i, r.offset+3*node2+j, - tmp[i][j]*kFactorTimesWeight);
}
}
}
}

template <class DataTypes>
void HyperReducedHexahedronFEMForceField<DataTypes>::buildStiffnessMatrix(
core::behavior::StiffnessMatrix* matrix)
Expand All @@ -472,14 +408,14 @@ void HyperReducedHexahedronFEMForceField<DataTypes>::buildStiffnessMatrix(
typename VecElement::const_iterator it, it0;

it0 = indexedElements->begin();
int nbElementsConsidered;
std::size_t nbElementsConsidered;

if (!d_performECSW.getValue())
nbElementsConsidered = indexedElements->size();
else
nbElementsConsidered = m_RIDsize;

for( unsigned int numElem = 0 ; numElem<nbElementsConsidered ;++numElem)
for(std::size_t numElem = 0 ; numElem<nbElementsConsidered ;++numElem)
{
if (!d_performECSW.getValue()){
e = numElem;
Expand Down
3 changes: 0 additions & 3 deletions src/ModelOrderReduction/component/forcefield/HyperReducedRestShapeSpringsForceField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,9 +103,6 @@ class HyperReducedRestShapeSpringsForceField : public virtual RestShapeSpringsFo

virtual void addDForce(const core::MechanicalParams* mparams, DataVecDeriv& df, const DataVecDeriv& dx) override;

/// Brings ForceField contribution to the global system stiffness matrix.
virtual void addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix ) override;

void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /* matrix */) override;

virtual void draw(const core::visual::VisualParams* vparams) override;
Expand Down
79 changes: 0 additions & 79 deletions src/ModelOrderReduction/component/forcefield/HyperReducedRestShapeSpringsForceField.inl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -437,83 +437,6 @@ void HyperReducedRestShapeSpringsForceField<DataTypes>::draw(const VisualParams
vparams->drawTool()->restoreLastState();
}

template<class DataTypes>
void HyperReducedRestShapeSpringsForceField<DataTypes>::addKToMatrix(const MechanicalParams* mparams, const MultiMatrixAccessor* matrix )
{
// remove to be able to build in parallel
// const VecIndex& indices = points.getValue();
// const VecReal& k = stiffness.getValue();
MultiMatrixAccessor::MatrixRef mref = matrix->getMatrix(this->mstate);
BaseMatrix* mat = mref.matrix;
unsigned int offset = mref.offset;
Real kFact = (Real)mparams->kFactorIncludingRayleighDamping(this->rayleighStiffness.getValue());

const int N = Coord::total_size;

unsigned int curIndex = 0;

const VecReal &k = d_stiffness.getValue();
if (k.size()!= m_indices.size() )
{
const Real k0 = k[0];
if (d_performECSW.getValue()){
for(unsigned int index = 0 ; index <m_RIDsize ; index++)
{
curIndex = m_indices[reducedIntegrationDomain(index)];

for(int i = 0; i < N; i++)
{
mat->add(offset + N * curIndex + i, offset + N * curIndex + i, -kFact * k0 * weights(reducedIntegrationDomain(index)));
}

}
}
else
{
for (unsigned int index = 0; index < m_indices.size(); index++)
{
curIndex = m_indices[index];

for(int i = 0; i < N; i++)
{
mat->add(offset + N * curIndex + i, offset + N * curIndex + i, -kFact * k0);
}
msg_info() << "1st: kfact is :" << kFact << ". Contrib is: " << -kFact * k0;
}
}

}
else
{
if (d_performECSW.getValue()){
for(unsigned int index = 0 ; index <m_RIDsize ; index++)
{
curIndex = m_indices[reducedIntegrationDomain(index)];

for(int i = 0; i < N; i++)
{
mat->add(offset + N * curIndex + i, offset + N * curIndex + i, -kFact * k[reducedIntegrationDomain(index)] * weights(reducedIntegrationDomain(index)));
}

}
}
else
{
for (unsigned int index = 0; index < m_indices.size(); index++)
{
curIndex = m_indices[index];

for(int i = 0; i < N; i++)
{
mat->add(offset + N * curIndex + i, offset + N * curIndex + i, -kFact * k[index]);
}
msg_info() << "2nd: kfact is :" << kFact << ". Contrib is " << -kFact * k[index];
}
}

}
}

template <class DataTypes>
void HyperReducedRestShapeSpringsForceField<DataTypes>::buildStiffnessMatrix(
core::behavior::StiffnessMatrix* matrix)
Expand All @@ -536,8 +459,6 @@ void HyperReducedRestShapeSpringsForceField<DataTypes>::buildStiffnessMatrix(
nbIndicesConsidered = m_indices.size();

for (sofa::Index index = 0; index < nbIndicesConsidered ; index++)

// for (const auto index : m_indices)
{
if (!d_performECSW.getValue())
curIndex = index;
Expand Down
1 change: 0 additions & 1 deletion ...delOrderReduction/component/forcefield/HyperReducedTetrahedralCorotationalFEMForceField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,6 @@ class HyperReducedTetrahedralCorotationalFEMForceField : public virtual Tetrahed
virtual void addForce(const core::MechanicalParams* mparams, DataVecDeriv& d_f, const DataVecCoord& d_x, const DataVecDeriv& d_v) override;
virtual void addDForce(const core::MechanicalParams* mparams, DataVecDeriv& d_df, const DataVecDeriv& d_dx) override;

virtual void addKToMatrix(sofa::linearalgebra::BaseMatrix *m, SReal kFactor, unsigned int &offset) override;
void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /* matrix */) override;

void draw(const core::visual::VisualParams* vparams) override;
Expand Down
89 changes: 35 additions & 54 deletions ...lOrderReduction/component/forcefield/HyperReducedTetrahedralCorotationalFEMForceField.inl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#pragma once

#include <ModelOrderReduction/component/forcefield/HyperReducedTetrahedralCorotationalFEMForceField.h>
#include <sofa/core/behavior/BaseLocalForceFieldMatrix.h>
#include <sofa/core/visual/VisualParams.h>
#include <sofa/core/MechanicalParams.h>
#include <sofa/component/topology/container/grid/GridTopology.h>
Expand Down Expand Up @@ -686,82 +687,62 @@ void HyperReducedTetrahedralCorotationalFEMForceField<DataTypes>::draw(const cor
vparams->drawTool()->restoreLastState();
}


template<class DataTypes>
void HyperReducedTetrahedralCorotationalFEMForceField<DataTypes>::addKToMatrix(sofa::linearalgebra::BaseMatrix *mat, SReal k, unsigned int &offset)
template <class DataTypes>
void HyperReducedTetrahedralCorotationalFEMForceField<DataTypes>::
buildStiffnessMatrix(core::behavior::StiffnessMatrix* matrix)
{
// Build Matrix Block for this ForceField
unsigned int i,j,n1, n2, row, column, ROW, COLUMN;

Transformation Rot;
StiffnessMatrix JKJt,tmp;
StiffnessMatrix JKJt, RJKJtRt;
sofa::type::Mat<3, 3, Real> localMatrix(type::NOINIT);

const type::vector<typename TetrahedralCorotationalFEMForceField<DataTypes>::TetrahedronInformation>& tetrahedronInf = tetrahedronInfo.getValue();

Index noeud1, noeud2;

Rot[0][0]=Rot[1][1]=Rot[2][2]=1;
Rot[0][1]=Rot[0][2]=0;
Rot[1][0]=Rot[1][2]=0;
Rot[2][0]=Rot[2][1]=0;
const sofa::core::topology::BaseMeshTopology::SeqTetrahedra& tetras = _topology->getTetrahedra();
const sofa::core::topology::BaseMeshTopology::SeqTetrahedra& tetrahedra = _topology->getTetrahedra();

unsigned int nbElementsConsidered;

if (!d_performECSW.getValue())
nbElementsConsidered = _topology->getNbTetrahedra();
else
nbElementsConsidered = m_RIDsize;
SReal kTimesWeight;

unsigned int IT;
for(unsigned int tetNum=0 ; tetNum < nbElementsConsidered ; ++tetNum)

static constexpr Transformation identity = []
{
Transformation i;
i.identity();
return i;
}();


auto dfdx = matrix->getForceDerivativeIn(this->mstate)
.withRespectToPositionsIn(this->mstate);

std::size_t IT;
for(std::size_t tetNum=0 ; tetNum < nbElementsConsidered ; ++tetNum)
{
if (!d_performECSW.getValue())
IT = tetNum;
else
IT = reducedIntegrationDomain(tetNum);

if (method == SMALL)
this->computeStiffnessMatrix(JKJt,tmp,tetrahedronInf[IT].materialMatrix,tetrahedronInf[IT].strainDisplacementTransposedMatrix,Rot);
else
this->computeStiffnessMatrix(JKJt,tmp,tetrahedronInf[IT].materialMatrix,tetrahedronInf[IT].strainDisplacementTransposedMatrix,tetrahedronInf[IT].rotation);
const core::topology::BaseMeshTopology::Tetrahedron t=tetras[IT];
if (!d_performECSW.getValue())
kTimesWeight = k;
else
kTimesWeight = k*weights(IT);
const auto& rotation = method == SMALL ? identity : tetrahedronInf[IT].rotation;
this->computeStiffnessMatrix(JKJt, RJKJtRt, tetrahedronInf[IT].materialMatrix,
tetrahedronInf[IT].strainDisplacementTransposedMatrix, rotation);

// find index of node 1
for (n1=0; n1<4; n1++)
{
noeud1 = t[n1];
const core::topology::BaseMeshTopology::Tetrahedron tetra = tetrahedra[IT];

for(i=0; i<3; i++)
static constexpr auto S = DataTypes::deriv_total_size; // size of node blocks
for (sofa::Size n1 = 0; n1 < core::topology::BaseMeshTopology::Tetrahedron::size(); ++n1)
{
for (sofa::Size n2 = 0; n2 < core::topology::BaseMeshTopology::Tetrahedron::size(); ++n2)
{
ROW = offset+3*noeud1+i;
row = 3*n1+i;
// find index of node 2
for (n2=0; n2<4; n2++)
{
noeud2 = t[n2];

for (j=0; j<3; j++)
{
COLUMN = offset+3*noeud2+j;
column = 3*n2+j;
mat->add(ROW, COLUMN, - tmp[row][column]*kTimesWeight);
}
}
RJKJtRt.getsub(S * n1, S * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2
if (!d_performECSW.getValue())
dfdx(S * tetra[n1], S * tetra[n2]) += -localMatrix;
else
dfdx(S * tetra[n1], S * tetra[n2]) += -localMatrix*weights(IT);
}
}
}
}

template <class DataTypes>
void HyperReducedTetrahedralCorotationalFEMForceField<DataTypes>::
buildStiffnessMatrix(core::behavior::StiffnessMatrix* matrix)
{
core::behavior::ForceField<DataTypes>::buildStiffnessMatrix(matrix);
}

} // namespace sofa::component::solidmechanics::fem::elastic
1 change: 0 additions & 1 deletion ...OrderReduction/component/forcefield/HyperReducedTetrahedronHyperelasticityFEMForceField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,6 @@ public :

virtual void addForce(const core::MechanicalParams* mparams /* PARAMS FIRST */, DataVecDeriv& d_f, const DataVecCoord& d_x, const DataVecDeriv& d_v) override;
virtual void addDForce(const core::MechanicalParams* mparams /* PARAMS FIRST */, DataVecDeriv& d_df, const DataVecDeriv& d_dx) override;
virtual void addKToMatrix(sofa::linearalgebra::BaseMatrix *mat, SReal k, unsigned int &offset) override;
void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /* matrix */) override;

void draw(const core::visual::VisualParams* vparams) override;
Expand Down
Loading

0 comments on commit a798e13

Please sign in to comment.