Skip to content

Commit a798e13

Browse files
VannesteFelixVannesteFelix
authored
[ForceField] add buildStiffnessMatrix to HyperReducedTetrahedralCorotationalFEMForceField (#92)
* [ForceField] add buildStiffnessMatrix to HyperReducedTetrahedralCorotationalFEMForceField * Clean: remove remaining & normally useless addKToMatrin
1 parent 7e9224e commit a798e13

8 files changed

+37
-281
lines changed

‎src/ModelOrderReduction/component/forcefield/HyperReducedHexahedronFEMForceField.h‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -143,7 +143,6 @@ class HyperReducedHexahedronFEMForceField : public virtual HexahedronFEMForceFie
143143
using InheritForceField::getPotentialEnergy;
144144
// getPotentialEnergy is implemented for polar method
145145

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

149148
void draw(const core::visual::VisualParams* vparams) override;

‎src/ModelOrderReduction/component/forcefield/HyperReducedHexahedronFEMForceField.inl‎

Lines changed: 2 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -388,70 +388,6 @@ template<class DataTypes>
388388
/////////////////////////////////////////////////
389389

390390

391-
template<class DataTypes>
392-
void HyperReducedHexahedronFEMForceField<DataTypes>::addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix)
393-
{
394-
// Build Matrix Block for this ForceField
395-
int i,j,n1, n2, e;
396-
397-
typename VecElement::const_iterator it, it0;
398-
399-
Index node1, node2;
400-
401-
sofa::core::behavior::MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(this->mstate);
402-
403-
it0=this->getIndexedElements()->begin();
404-
size_t nbElementsConsidered;
405-
if (!d_performECSW.getValue()){
406-
nbElementsConsidered = this->getIndexedElements()->size();
407-
}
408-
else
409-
{
410-
nbElementsConsidered = m_RIDsize;
411-
}
412-
for( unsigned int numElem = 0 ; numElem<nbElementsConsidered ;++numElem)
413-
{
414-
if (!d_performECSW.getValue()){
415-
e = numElem;
416-
}
417-
else
418-
{
419-
e = reducedIntegrationDomain(numElem);
420-
}
421-
it = it0 + e;
422-
423-
const ElementStiffness &Ke = _elementStiffnesses.getValue()[e];
424-
// const Transformation& Rt = _rotations[e];
425-
// Transformation R; R.transpose(Rt);
426-
427-
Transformation Rot = this->getElementRotation(e);
428-
Real kFactorTimesWeight;
429-
if (!d_performECSW.getValue())
430-
kFactorTimesWeight = (Real)mparams->kFactorIncludingRayleighDamping(this->rayleighStiffness.getValue());
431-
else
432-
kFactorTimesWeight = weights(e)*(Real)mparams->kFactorIncludingRayleighDamping(this->rayleighStiffness.getValue());
433-
434-
435-
// find index of node 1
436-
for (n1=0; n1<8; n1++)
437-
{
438-
node1 = (*it)[n1];
439-
// find index of node 2
440-
for (n2=0; n2<8; n2++)
441-
{
442-
node2 = (*it)[n2];
443-
444-
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]),
445-
Coord(Ke[3*n1+1][3*n2+0],Ke[3*n1+1][3*n2+1],Ke[3*n1+1][3*n2+2]),
446-
Coord(Ke[3*n1+2][3*n2+0],Ke[3*n1+2][3*n2+1],Ke[3*n1+2][3*n2+2])) ) * Rot;
447-
for(i=0; i<3; i++)
448-
for (j=0; j<3; j++)
449-
r.matrix->add(r.offset+3*node1+i, r.offset+3*node2+j, - tmp[i][j]*kFactorTimesWeight);
450-
}
451-
}
452-
}
453-
}
454-
455391
template <class DataTypes>
456392
void HyperReducedHexahedronFEMForceField<DataTypes>::buildStiffnessMatrix(
457393
core::behavior::StiffnessMatrix* matrix)
@@ -472,14 +408,14 @@ void HyperReducedHexahedronFEMForceField<DataTypes>::buildStiffnessMatrix(
472408
typename VecElement::const_iterator it, it0;
473409

474410
it0 = indexedElements->begin();
475-
int nbElementsConsidered;
411+
std::size_t nbElementsConsidered;
476412

477413
if (!d_performECSW.getValue())
478414
nbElementsConsidered = indexedElements->size();
479415
else
480416
nbElementsConsidered = m_RIDsize;
481417

482-
for( unsigned int numElem = 0 ; numElem<nbElementsConsidered ;++numElem)
418+
for(std::size_t numElem = 0 ; numElem<nbElementsConsidered ;++numElem)
483419
{
484420
if (!d_performECSW.getValue()){
485421
e = numElem;

‎src/ModelOrderReduction/component/forcefield/HyperReducedRestShapeSpringsForceField.h‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -103,9 +103,6 @@ class HyperReducedRestShapeSpringsForceField : public virtual RestShapeSpringsFo
103103

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

106-
/// Brings ForceField contribution to the global system stiffness matrix.
107-
virtual void addKToMatrix(const core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix ) override;
108-
109106
void buildStiffnessMatrix(core::behavior::StiffnessMatrix* /* matrix */) override;
110107

111108
virtual void draw(const core::visual::VisualParams* vparams) override;

‎src/ModelOrderReduction/component/forcefield/HyperReducedRestShapeSpringsForceField.inl‎

Lines changed: 0 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -437,83 +437,6 @@ void HyperReducedRestShapeSpringsForceField<DataTypes>::draw(const VisualParams
437437
vparams->drawTool()->restoreLastState();
438438
}
439439

440-
template<class DataTypes>
441-
void HyperReducedRestShapeSpringsForceField<DataTypes>::addKToMatrix(const MechanicalParams* mparams, const MultiMatrixAccessor* matrix )
442-
{
443-
// remove to be able to build in parallel
444-
// const VecIndex& indices = points.getValue();
445-
// const VecReal& k = stiffness.getValue();
446-
MultiMatrixAccessor::MatrixRef mref = matrix->getMatrix(this->mstate);
447-
BaseMatrix* mat = mref.matrix;
448-
unsigned int offset = mref.offset;
449-
Real kFact = (Real)mparams->kFactorIncludingRayleighDamping(this->rayleighStiffness.getValue());
450-
451-
const int N = Coord::total_size;
452-
453-
unsigned int curIndex = 0;
454-
455-
const VecReal &k = d_stiffness.getValue();
456-
if (k.size()!= m_indices.size() )
457-
{
458-
const Real k0 = k[0];
459-
if (d_performECSW.getValue()){
460-
for(unsigned int index = 0 ; index <m_RIDsize ; index++)
461-
{
462-
curIndex = m_indices[reducedIntegrationDomain(index)];
463-
464-
for(int i = 0; i < N; i++)
465-
{
466-
mat->add(offset + N * curIndex + i, offset + N * curIndex + i, -kFact * k0 * weights(reducedIntegrationDomain(index)));
467-
}
468-
469-
}
470-
}
471-
else
472-
{
473-
for (unsigned int index = 0; index < m_indices.size(); index++)
474-
{
475-
curIndex = m_indices[index];
476-
477-
for(int i = 0; i < N; i++)
478-
{
479-
mat->add(offset + N * curIndex + i, offset + N * curIndex + i, -kFact * k0);
480-
}
481-
msg_info() << "1st: kfact is :" << kFact << ". Contrib is: " << -kFact * k0;
482-
}
483-
}
484-
485-
}
486-
else
487-
{
488-
if (d_performECSW.getValue()){
489-
for(unsigned int index = 0 ; index <m_RIDsize ; index++)
490-
{
491-
curIndex = m_indices[reducedIntegrationDomain(index)];
492-
493-
for(int i = 0; i < N; i++)
494-
{
495-
mat->add(offset + N * curIndex + i, offset + N * curIndex + i, -kFact * k[reducedIntegrationDomain(index)] * weights(reducedIntegrationDomain(index)));
496-
}
497-
498-
}
499-
}
500-
else
501-
{
502-
for (unsigned int index = 0; index < m_indices.size(); index++)
503-
{
504-
curIndex = m_indices[index];
505-
506-
for(int i = 0; i < N; i++)
507-
{
508-
mat->add(offset + N * curIndex + i, offset + N * curIndex + i, -kFact * k[index]);
509-
}
510-
msg_info() << "2nd: kfact is :" << kFact << ". Contrib is " << -kFact * k[index];
511-
}
512-
}
513-
514-
}
515-
}
516-
517440
template <class DataTypes>
518441
void HyperReducedRestShapeSpringsForceField<DataTypes>::buildStiffnessMatrix(
519442
core::behavior::StiffnessMatrix* matrix)
@@ -536,8 +459,6 @@ void HyperReducedRestShapeSpringsForceField<DataTypes>::buildStiffnessMatrix(
536459
nbIndicesConsidered = m_indices.size();
537460

538461
for (sofa::Index index = 0; index < nbIndicesConsidered ; index++)
539-
540-
// for (const auto index : m_indices)
541462
{
542463
if (!d_performECSW.getValue())
543464
curIndex = index;

‎src/ModelOrderReduction/component/forcefield/HyperReducedTetrahedralCorotationalFEMForceField.h‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -130,7 +130,6 @@ class HyperReducedTetrahedralCorotationalFEMForceField : public virtual Tetrahed
130130
virtual void addForce(const core::MechanicalParams* mparams, DataVecDeriv& d_f, const DataVecCoord& d_x, const DataVecDeriv& d_v) override;
131131
virtual void addDForce(const core::MechanicalParams* mparams, DataVecDeriv& d_df, const DataVecDeriv& d_dx) override;
132132

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

136135
void draw(const core::visual::VisualParams* vparams) override;

‎src/ModelOrderReduction/component/forcefield/HyperReducedTetrahedralCorotationalFEMForceField.inl‎

Lines changed: 35 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#pragma once
1818

1919
#include <ModelOrderReduction/component/forcefield/HyperReducedTetrahedralCorotationalFEMForceField.h>
20+
#include <sofa/core/behavior/BaseLocalForceFieldMatrix.h>
2021
#include <sofa/core/visual/VisualParams.h>
2122
#include <sofa/core/MechanicalParams.h>
2223
#include <sofa/component/topology/container/grid/GridTopology.h>
@@ -686,82 +687,62 @@ void HyperReducedTetrahedralCorotationalFEMForceField<DataTypes>::draw(const cor
686687
vparams->drawTool()->restoreLastState();
687688
}
688689

689-
690-
template<class DataTypes>
691-
void HyperReducedTetrahedralCorotationalFEMForceField<DataTypes>::addKToMatrix(sofa::linearalgebra::BaseMatrix *mat, SReal k, unsigned int &offset)
690+
template <class DataTypes>
691+
void HyperReducedTetrahedralCorotationalFEMForceField<DataTypes>::
692+
buildStiffnessMatrix(core::behavior::StiffnessMatrix* matrix)
692693
{
693-
// Build Matrix Block for this ForceField
694-
unsigned int i,j,n1, n2, row, column, ROW, COLUMN;
695-
696-
Transformation Rot;
697-
StiffnessMatrix JKJt,tmp;
694+
StiffnessMatrix JKJt, RJKJtRt;
695+
sofa::type::Mat<3, 3, Real> localMatrix(type::NOINIT);
698696

699697
const type::vector<typename TetrahedralCorotationalFEMForceField<DataTypes>::TetrahedronInformation>& tetrahedronInf = tetrahedronInfo.getValue();
700-
701-
Index noeud1, noeud2;
702-
703-
Rot[0][0]=Rot[1][1]=Rot[2][2]=1;
704-
Rot[0][1]=Rot[0][2]=0;
705-
Rot[1][0]=Rot[1][2]=0;
706-
Rot[2][0]=Rot[2][1]=0;
707-
const sofa::core::topology::BaseMeshTopology::SeqTetrahedra& tetras = _topology->getTetrahedra();
698+
const sofa::core::topology::BaseMeshTopology::SeqTetrahedra& tetrahedra = _topology->getTetrahedra();
708699

709700
unsigned int nbElementsConsidered;
701+
710702
if (!d_performECSW.getValue())
711703
nbElementsConsidered = _topology->getNbTetrahedra();
712704
else
713705
nbElementsConsidered = m_RIDsize;
714-
SReal kTimesWeight;
715706

716-
unsigned int IT;
717-
for(unsigned int tetNum=0 ; tetNum < nbElementsConsidered ; ++tetNum)
707+
708+
static constexpr Transformation identity = []
709+
{
710+
Transformation i;
711+
i.identity();
712+
return i;
713+
}();
714+
715+
716+
auto dfdx = matrix->getForceDerivativeIn(this->mstate)
717+
.withRespectToPositionsIn(this->mstate);
718+
719+
std::size_t IT;
720+
for(std::size_t tetNum=0 ; tetNum < nbElementsConsidered ; ++tetNum)
718721
{
719722
if (!d_performECSW.getValue())
720723
IT = tetNum;
721724
else
722725
IT = reducedIntegrationDomain(tetNum);
723726

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

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

739-
for(i=0; i<3; i++)
733+
static constexpr auto S = DataTypes::deriv_total_size; // size of node blocks
734+
for (sofa::Size n1 = 0; n1 < core::topology::BaseMeshTopology::Tetrahedron::size(); ++n1)
735+
{
736+
for (sofa::Size n2 = 0; n2 < core::topology::BaseMeshTopology::Tetrahedron::size(); ++n2)
740737
{
741-
ROW = offset+3*noeud1+i;
742-
row = 3*n1+i;
743-
// find index of node 2
744-
for (n2=0; n2<4; n2++)
745-
{
746-
noeud2 = t[n2];
747-
748-
for (j=0; j<3; j++)
749-
{
750-
COLUMN = offset+3*noeud2+j;
751-
column = 3*n2+j;
752-
mat->add(ROW, COLUMN, - tmp[row][column]*kTimesWeight);
753-
}
754-
}
738+
RJKJtRt.getsub(S * n1, S * n2, localMatrix); //extract the submatrix corresponding to the coupling of nodes n1 and n2
739+
if (!d_performECSW.getValue())
740+
dfdx(S * tetra[n1], S * tetra[n2]) += -localMatrix;
741+
else
742+
dfdx(S * tetra[n1], S * tetra[n2]) += -localMatrix*weights(IT);
755743
}
756744
}
757745
}
758746
}
759747

760-
template <class DataTypes>
761-
void HyperReducedTetrahedralCorotationalFEMForceField<DataTypes>::
762-
buildStiffnessMatrix(core::behavior::StiffnessMatrix* matrix)
763-
{
764-
core::behavior::ForceField<DataTypes>::buildStiffnessMatrix(matrix);
765-
}
766-
767748
} // namespace sofa::component::solidmechanics::fem::elastic

‎src/ModelOrderReduction/component/forcefield/HyperReducedTetrahedronHyperelasticityFEMForceField.h‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -118,7 +118,6 @@ public :
118118

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

124123
void draw(const core::visual::VisualParams* vparams) override;

0 commit comments

Comments
 (0)