Skip to content

Commit

Permalink
[FEM.Elastic] Remove unused methods
Browse files Browse the repository at this point in the history
  • Loading branch information
alxbilger committed Jul 11, 2024
1 parent 0fd439a commit d0ea18c
Show file tree
Hide file tree
Showing 2 changed files with 0 additions and 89 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -359,8 +359,6 @@ class TetrahedronFEMForceField : public BaseLinearElasticityFEMForceField<DataTy

// Getting the stiffness matrix of index i
void getElementStiffnessMatrix(Real* stiffness, Index nodeIdx);
void getElementStiffnessMatrix(Real* stiffness, Tetrahedron& te);
virtual void computeMaterialStiffness(MaterialStiffness& materialMatrix, Index&a, Index&b, Index&c, Index&d);

protected:
void computeStrainDisplacement( StrainDisplacement &J, Coord a, Coord b, Coord c, Coord d );
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -280,57 +280,6 @@ inline void TetrahedronFEMForceField<DataTypes>::getElementStiffnessMatrix(Real*
}
}

template <class DataTypes>
inline void TetrahedronFEMForceField<DataTypes>::getElementStiffnessMatrix(Real* stiffness, Tetrahedron& te)
{
if (needUpdateTopology)
{
reinit();
needUpdateTopology = false;
}
const VecCoord X0=this->mstate->read(core::ConstVecCoordId::restPosition())->getValue();
Index a = te[0];
Index b = te[1];
Index c = te[2];
Index d = te[3];

Transformation R_0_1;
computeRotationLarge( R_0_1, (X0), a, b, c);

MaterialStiffness materialMatrix;
StrainDisplacement strainMatrix;
type::fixed_array<Coord,4> rotatedInitialElements;

rotatedInitialElements[0] = R_0_1*(X0)[a];
rotatedInitialElements[1] = R_0_1*(X0)[b];
rotatedInitialElements[2] = R_0_1*(X0)[c];
rotatedInitialElements[3] = R_0_1*(X0)[d];

rotatedInitialElements[1] -= rotatedInitialElements[0];
rotatedInitialElements[2] -= rotatedInitialElements[0];
rotatedInitialElements[3] -= rotatedInitialElements[0];
rotatedInitialElements[0] = Coord(0,0,0);

computeMaterialStiffness(materialMatrix,a,b,c,d);
computeStrainDisplacement(strainMatrix, rotatedInitialElements[0], rotatedInitialElements[1], rotatedInitialElements[2], rotatedInitialElements[3]);

Transformation Rot;
StiffnessMatrix JKJt,tmp;
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;

R_0_1.transpose();
computeStiffnessMatrix(JKJt, tmp, materialMatrix, strainMatrix, R_0_1);
for(int i=0; i<12; i++)
{
for(int j=0; j<12; j++)
stiffness[i*12+j]=tmp(i,j);
}
}


template<class DataTypes>
void TetrahedronFEMForceField<DataTypes>::computeMaterialStiffness(Index i, Index&a, Index&b, Index&c, Index&d)
{
Expand Down Expand Up @@ -376,42 +325,6 @@ void TetrahedronFEMForceField<DataTypes>::computeMaterialStiffness(Index i, Inde
materialsStiffnesses[i] /= (volumes6*6); // 36*Volume in the formula
}

template<class DataTypes>
void TetrahedronFEMForceField<DataTypes>::computeMaterialStiffness(MaterialStiffness& materialMatrix, Index&a, Index&b, Index&c, Index&d)
{
const Real youngModulus = this->d_youngModulus.getValue()[0];
const Real poissonRatio = this->d_poissonRatio.getValue();

materialMatrix[0][0] = materialMatrix[1][1] = materialMatrix[2][2] = 1;
materialMatrix[0][1] = materialMatrix[0][2] = materialMatrix[1][0] = materialMatrix[1][2] = materialMatrix[2][0] = materialMatrix[2][1] = poissonRatio/(1-poissonRatio);
materialMatrix[0][3] = materialMatrix[0][4] = materialMatrix[0][5] = 0;
materialMatrix[1][3] = materialMatrix[1][4] = materialMatrix[1][5] = 0;
materialMatrix[2][3] = materialMatrix[2][4] = materialMatrix[2][5] = 0;
materialMatrix[3][0] = materialMatrix[3][1] = materialMatrix[3][2] = materialMatrix[3][4] = materialMatrix[3][5] = 0;
materialMatrix[4][0] = materialMatrix[4][1] = materialMatrix[4][2] = materialMatrix[4][3] = materialMatrix[4][5] = 0;
materialMatrix[5][0] = materialMatrix[5][1] = materialMatrix[5][2] = materialMatrix[5][3] = materialMatrix[5][4] = 0;
materialMatrix[3][3] = materialMatrix[4][4] = materialMatrix[5][5] = (1-2*poissonRatio)/(2*(1-poissonRatio));
materialMatrix *= (youngModulus*(1-poissonRatio))/((1+poissonRatio)*(1-2*poissonRatio));

// divide by 36 times volumes of the element
const VecCoord X0=this->mstate->read(core::ConstVecCoordId::restPosition())->getValue();

Coord A = (X0)[b] - (X0)[a];
Coord B = (X0)[c] - (X0)[a];
Coord C = (X0)[d] - (X0)[a];
Coord AB = cross(A, B);
Real volumes6 = fabs( dot( AB, C ) );
m_restVolume += volumes6/6;
if (volumes6<0)
{
msg_error() << "Negative volume for tetra"<<a<<','<<b<<','<<c<<','<<d<<"> = "<<volumes6/6 ;
}
materialMatrix /= (volumes6*6);
}




template<class DataTypes>
inline void TetrahedronFEMForceField<DataTypes>::computeForce( Displacement &F, const Displacement &Depl, VoigtTensor &plasticStrain, const MaterialStiffness &K, const StrainDisplacement &J )
{
Expand Down

0 comments on commit d0ea18c

Please sign in to comment.