Skip to content
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

Jump to bottom

[FEM] Add method computeVonMisesStress in TetrahedralCorotationalForceField and option to draw them #4945

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

[FEM] Add method computeVonMisesStress in TetrahedralCorotationalForceField and option to draw them #4945

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
2022a0d
[FEM] Add method to computeVonMisesStress in TetrahedralCorotational…
epernod Jun 14, 2024
744ead0
[FEM] Update TetrahedralCoro von mises drawing logic
epernod Jul 15, 2024
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
13 changes: 13 additions & 0 deletions ...stic/src/sofa/component/solidmechanics/fem/elastic/TetrahedralCorotationalFEMForceField.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
#include <sofa/type/Vec.h>
#include <sofa/type/Mat.h>
#include <sofa/helper/map.h>
#include <sofa/helper/ColorMap.h>

// corotational tetrahedron from
// @InProceedings{NPF05,
Expand Down Expand Up @@ -102,6 +103,7 @@ class TetrahedralCorotationalFEMForceField : public BaseLinearElasticityFEMForce
StrainDisplacementTransposed strainDisplacementTransposedMatrix;
/// large displacement method
type::fixed_array<Coord,4> rotatedInitialElements;
type::Mat<4, 4, Real> elemShapeFun;
Transformation rotation;
/// polar method
Transformation initialTransformation;
Expand Down Expand Up @@ -190,12 +192,22 @@ class TetrahedralCorotationalFEMForceField : public BaseLinearElasticityFEMForce
Data<sofa::type::RGBAColor> d_drawColor3; ///< draw color for faces 3
Data<sofa::type::RGBAColor> d_drawColor4; ///< draw color for faces 4
Data<std::map < std::string, sofa::type::vector<double> > > _volumeGraph;

Data<bool> d_computeVonMisesStress;
Data<type::vector<Real> > d_vonMisesPerElement; ///< von Mises Stress per element
Data<type::vector<Real> > d_vonMisesPerNode; ///< von Mises Stress per node

sofa::helper::ColorMap* m_VonMisesColorMap;
Comment on lines +198 to +200
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Data<type::vector<Real> > d_vonMisesPerNode; ///< von Mises Stress per node
sofa::helper::ColorMap* m_VonMisesColorMap;
Data<type::vector<Real> > d_vonMisesPerNode; ///< von Mises Stress per node
sofa::helper::ColorMap* m_vonMisesColorMap;

I guess this is a copy-paste from the other file 🫢but it is nice to have consistency

Real prevMaxStress = -1.0;

using Inherit1::l_topology;


protected:
TetrahedralCorotationalFEMForceField();

~TetrahedralCorotationalFEMForceField() override;

/// Pointer to the topology container. Will be set by link @sa l_topology
sofa::core::topology::BaseMeshTopology* m_topology;
public:
Expand Down Expand Up @@ -235,6 +247,7 @@ class TetrahedralCorotationalFEMForceField : public BaseLinearElasticityFEMForce

void computeBBox(const core::ExecParams* params, bool onlyVisible) override;

void computeVonMisesStress();

protected:
/** Method to create @sa TetrahedronInformation when a new tetrahedron is created.
Expand Down
251 changes: 214 additions & 37 deletions ...ic/src/sofa/component/solidmechanics/fem/elastic/TetrahedralCorotationalFEMForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,10 @@ TetrahedralCorotationalFEMForceField<DataTypes>::TetrahedralCorotationalFEMForce
, d_drawColor2(initData(&d_drawColor2, sofa::type::RGBAColor(0.0f, 0.5f, 1.0f, 1.0f), "drawColor2", " draw color for faces 2"))
, d_drawColor3(initData(&d_drawColor3, sofa::type::RGBAColor(0.0f, 1.0f, 1.0f, 1.0f), "drawColor3", " draw color for faces 3"))
, d_drawColor4(initData(&d_drawColor4, sofa::type::RGBAColor(0.5f, 1.0f, 1.0f, 1.0f), "drawColor4", " draw color for faces 4"))
, d_computeVonMisesStress(initData(&d_computeVonMisesStress, false, "computeVonMisesStress", "compute and display von Mises stress: 0: no computations, 1: using corotational strain, 2: using full Green strain. Set listening=1"))
, d_vonMisesPerElement(initData(&d_vonMisesPerElement, "vonMisesPerElement", "von Mises Stress per element"))
, d_vonMisesPerNode(initData(&d_vonMisesPerNode, "vonMisesPerNode", "von Mises Stress per node"))
, m_VonMisesColorMap(nullptr)
{
this->addAlias(&d_assembling, "assembling");

Expand All @@ -96,6 +100,14 @@ TetrahedralCorotationalFEMForceField<DataTypes>::TetrahedralCorotationalFEMForce

}


template <class DataTypes>
TetrahedralCorotationalFEMForceField<DataTypes>::~TetrahedralCorotationalFEMForceField()
{
if (m_VonMisesColorMap != nullptr)
delete m_VonMisesColorMap;
}

template <class DataTypes>
void TetrahedralCorotationalFEMForceField<DataTypes>::init()
{
Expand All @@ -117,6 +129,11 @@ void TetrahedralCorotationalFEMForceField<DataTypes>::init()
}

reinit(); // compute per-element stiffness matrices and other precomputed values

if (m_VonMisesColorMap == nullptr)
{
m_VonMisesColorMap = new helper::ColorMap(256, std::string("Blue to Red"));
}
}


Expand Down Expand Up @@ -195,6 +212,9 @@ void TetrahedralCorotationalFEMForceField<DataTypes>::addForce(const core::Mecha
}
}
d_f.endEdit();

if (d_computeVonMisesStress.getValue())
computeVonMisesStress();
}

template<class DataTypes>
Expand Down Expand Up @@ -874,7 +894,19 @@ void TetrahedralCorotationalFEMForceField<DataTypes>::initLarge(int i, Index&a,

computeStrainDisplacement( tetrahedronInf[i].strainDisplacementTransposedMatrix,tetrahedronInf[i].rotatedInitialElements[0], tetrahedronInf[i].rotatedInitialElements[1],tetrahedronInf[i].rotatedInitialElements[2],tetrahedronInf[i].rotatedInitialElements[3] );

d_tetrahedronInfo.endEdit();

type::Mat<4, 4, Real> matVert;
const core::topology::BaseMeshTopology::Tetrahedron tetra = { a, b, c, d };
for (Index k = 0; k < 4; k++) {
Index ix = tetra[k];
matVert[k][0] = 1.0;
for (Index l = 1; l < 4; l++)
matVert[k][l] = X0[ix][l - 1];
}

const bool canInvert = type::invertMatrix(tetrahedronInf[i].elemShapeFun, matVert);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const bool canInvert = type::invertMatrix(tetrahedronInf[i].elemShapeFun, matVert);
[[maybe_unused]] const bool canInvert = type::invertMatrix(tetrahedronInf[i].elemShapeFun, matVert);


tetrahedronInfo.endEdit();
}

template<class DataTypes>
Expand Down Expand Up @@ -1197,12 +1229,119 @@ void TetrahedralCorotationalFEMForceField<DataTypes>::computeBBox(const core::Ex
}


template<class DataTypes>
void TetrahedralCorotationalFEMForceField<DataTypes>::computeVonMisesStress()
{
typename core::behavior::MechanicalState<DataTypes>* mechanicalObject;
this->getContext()->get(mechanicalObject);
const VecCoord& X = mechanicalObject->read(core::ConstVecCoordId::position())->getValue();

const sofa::core::topology::BaseMeshTopology::SeqTetrahedra& tetras = m_topology->getTetrahedra();
const type::vector<typename TetrahedralCorotationalFEMForceField<DataTypes>::TetrahedronInformation>& tetrahedronInf = tetrahedronInfo.getValue();

helper::WriteAccessor<Data<type::vector<Real> > > vME = d_vonMisesPerElement;
vME.resize(tetras.size());

for (Size i = 0; i < m_topology->getNbTetrahedra(); ++i)
{
type::Vec<6, Real> vStrain;
type::Mat<3, 3, Real> gradU;
Transformation R_0_2;
Displacement D;

const auto& tetra = tetras[i];
const TetrahedronInformation& tetraInfo = tetrahedronInf[i];

computeRotationLarge(R_0_2, X, tetra[0], tetra[1], tetra[2]);
//tetraInfo.rotation.transpose(R_0_2);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
//tetraInfo.rotation.transpose(R_0_2);


// positions of the deformed and displaced Tetrahedron in its frame
type::fixed_array<Coord, 4> deforme;
for (int j = 0; j < 4; ++j)
deforme[j] = R_0_2 * X[tetra[j]];

deforme[1][0] -= deforme[0][0];
deforme[2][0] -= deforme[0][0];
deforme[2][1] -= deforme[0][1];
deforme[3] -= deforme[0];

// displacement
D[0] = 0;
D[1] = 0;
D[2] = 0;
D[3] = tetraInfo.rotatedInitialElements[1][0] - deforme[1][0];
D[4] = 0;
D[5] = 0;
D[6] = tetraInfo.rotatedInitialElements[2][0] - deforme[2][0];
D[7] = tetraInfo.rotatedInitialElements[2][1] - deforme[2][1];
D[8] = 0;
D[9] = tetraInfo.rotatedInitialElements[3][0] - deforme[3][0];
D[10] = tetraInfo.rotatedInitialElements[3][1] - deforme[3][1];
D[11] = tetraInfo.rotatedInitialElements[3][2] - deforme[3][2];

const type::Mat<4, 4, Real>& shf = tetraInfo.elemShapeFun;

/// compute gradU
for (Index k = 0; k < 3; k++) {
for (Index l = 0; l < 3; l++) {
gradU[k][l] = 0.0;
for (Index m = 0; m < 4; m++)
gradU[k][l] += shf[l + 1][m] * D[3 * m + k];
}
}

type::Mat<3, 3, Real> strain = Real(0.5) * (gradU + gradU.transposed());

for (Index j = 0; j < 3; j++)
vStrain[j] = strain[j][j];
vStrain[3] = strain[1][2];
vStrain[4] = strain[0][2];
vStrain[5] = strain[0][1];

Real lambda = tetraInfo.materialMatrix[0][1];
Real mu = tetraInfo.materialMatrix[3][3];

/// stress
type::VecNoInit<6, Real> s; // VoigtTensor

Real traceStrain = 0.0;
for (Index k = 0; k < 3; k++) {
traceStrain += vStrain[k];
s[k] = vStrain[k] * 2 * mu;
}

for (Index k = 3; k < 6; k++)
s[k] = vStrain[k] * 2 * mu;

for (Index k = 0; k < 3; k++)
s[k] += lambda * traceStrain;

vME[i] = helper::rsqrt(s[0] * s[0] + s[1] * s[1] + s[2] * s[2] - s[0] * s[1] - s[1] * s[2] - s[2] * s[0] + 3 * s[3] * s[3] + 3 * s[4] * s[4] + 3 * s[5] * s[5]);
if (vME[i] < 1e-10)
vME[i] = 0.0;
}

helper::WriteAccessor<Data<type::vector<Real> > > vMN = d_vonMisesPerNode;
vMN.resize(X.size());

/// compute the values of vonMises stress in nodes
for (Index dof = 0; dof < X.size(); dof++) {
core::topology::BaseMeshTopology::TetrahedraAroundVertex tetrasAroundDOF = m_topology->getTetrahedraAroundVertex(dof);

vMN[dof] = 0.0;
for (size_t at = 0; at < tetrasAroundDOF.size(); at++)
vMN[dof] += vME[tetrasAroundDOF[at]];
if (!tetrasAroundDOF.empty())
vMN[dof] /= Real(tetrasAroundDOF.size());
}
}


template<class DataTypes>
void TetrahedralCorotationalFEMForceField<DataTypes>::draw(const core::visual::VisualParams* vparams)
{
if (!vparams->displayFlags().getShowForceFields()) return;
if (!this->mstate) return;
if (!d_drawing.getValue()) return;

const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle();

Expand All @@ -1211,46 +1350,84 @@ void TetrahedralCorotationalFEMForceField<DataTypes>::draw(const core::visual::V
if (vparams->displayFlags().getShowWireFrame())
vparams->drawTool()->setPolygonMode(0,true);


std::vector< type::Vec3 > points[4];
for(Size i = 0 ; i<m_topology->getNbTetrahedra(); ++i)
if (d_computeVonMisesStress.getValue())
{
const core::topology::BaseMeshTopology::Tetrahedron t=m_topology->getTetrahedron(i);

const auto& [a, b, c, d] = t.array();
Coord center = (x[a]+x[b]+x[c]+x[d])*0.125;
Coord pa = (x[a]+center)*(Real)0.666667;
Coord pb = (x[b]+center)*(Real)0.666667;
Coord pc = (x[c]+center)*(Real)0.666667;
Coord pd = (x[d]+center)*(Real)0.666667;

points[0].push_back(pa);
points[0].push_back(pb);
points[0].push_back(pc);

points[1].push_back(pb);
points[1].push_back(pc);
points[1].push_back(pd);

points[2].push_back(pc);
points[2].push_back(pd);
points[2].push_back(pa);

points[3].push_back(pd);
points[3].push_back(pa);
points[3].push_back(pb);
Real minVM = (Real)1e20, maxVM = (Real)-1e20;
Real minVMN = (Real)1e20, maxVMN = (Real)-1e20;
helper::ReadAccessor<Data<type::vector<Real> > > vM = d_vonMisesPerElement;
helper::ReadAccessor<Data<type::vector<Real> > > vMN = d_vonMisesPerNode;

if (m_VonMisesColorMap == nullptr)
return;

for (size_t i = 0; i < vM.size(); i++)
{
minVM = (vM[i] < minVM) ? vM[i] : minVM;
maxVM = (vM[i] > maxVM) ? vM[i] : maxVM;
}
if (maxVM < prevMaxStress)
{
maxVM = prevMaxStress;
}
for (size_t i = 0; i < vMN.size(); i++)
{
minVMN = (vMN[i] < minVMN) ? vMN[i] : minVMN;
maxVMN = (vMN[i] > maxVMN) ? vMN[i] : maxVMN;
}

// Draw nodes (if node option enabled)
if (vMN.size() == x.size())
{
std::vector<sofa::type::RGBAColor> nodeColors(x.size());
std::vector<type::Vec3> pts(x.size());
helper::ColorMap::evaluator<Real> evalColor = m_VonMisesColorMap->getEvaluator(minVMN, maxVMN);
for (size_t nd = 0; nd < x.size(); nd++) {
pts[nd] = x[nd];
nodeColors[nd] = evalColor(vMN[nd]);
}
vparams->drawTool()->drawPoints(pts, 10, nodeColors);
}
}

if (d_drawing.getValue())
{
std::vector< type::Vec3 > points[4];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could reserve the memory beforehand as the size seems computable

for (Size i = 0; i < m_topology->getNbTetrahedra(); ++i)
{
const core::topology::BaseMeshTopology::Tetrahedron t = m_topology->getTetrahedron(i);

const auto& [a, b, c, d] = t.array();
Coord center = (x[a] + x[b] + x[c] + x[d]) * 0.125;
Coord pa = (x[a] + center) * (Real)0.666667;
Coord pb = (x[b] + center) * (Real)0.666667;
Coord pc = (x[c] + center) * (Real)0.666667;
Coord pd = (x[d] + center) * (Real)0.666667;

points[0].push_back(pa);
points[0].push_back(pb);
points[0].push_back(pc);

points[1].push_back(pb);
points[1].push_back(pc);
points[1].push_back(pd);

points[2].push_back(pc);
points[2].push_back(pd);
points[2].push_back(pa);

points[3].push_back(pd);
points[3].push_back(pa);
points[3].push_back(pb);
}

vparams->drawTool()->drawTriangles(points[0], d_drawColor1.getValue());
vparams->drawTool()->drawTriangles(points[1], d_drawColor2.getValue());
vparams->drawTool()->drawTriangles(points[2], d_drawColor3.getValue());
vparams->drawTool()->drawTriangles(points[3], d_drawColor4.getValue());

vparams->drawTool()->drawTriangles(points[0], d_drawColor1.getValue());
vparams->drawTool()->drawTriangles(points[1], d_drawColor2.getValue());
vparams->drawTool()->drawTriangles(points[2], d_drawColor3.getValue());
vparams->drawTool()->drawTriangles(points[3], d_drawColor4.getValue());
}

if (vparams->displayFlags().getShowWireFrame())
vparams->drawTool()->setPolygonMode(0,false);



}


Expand Down
4 changes: 2 additions & 2 deletions examples/Component/SolidMechanics/FEM/TetrahedralCorotationalFEMForceField.scn
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@
<Hexa2TetraTopologicalMapping input="@grid" output="@Tetra_topo" />

<DiagonalMass massDensity="0.2" />
<TetrahedralCorotationalFEMForceField name="CFEM" youngModulus="1000" poissonRatio="0.3" method="large" />
<TetrahedralCorotationalFEMForceField name="CFEM" youngModulus="1000" poissonRatio="0.3" method="large" computeVonMisesStress="1"/>

<BoxROI template="Vec3" name="box_roi" box="-6 -6 -1 30 6 0.1" drawBoxes="1" />
<FixedProjectiveConstraint template="Vec3" indices="@box_roi.indices" />
Expand All @@ -67,7 +67,7 @@
<Hexa2TetraTopologicalMapping input="@grid" output="@Tetra_topo" />

<DiagonalMass massDensity="0.2" />
<TetrahedralCorotationalFEMForceField name="CFEM" youngModulus="1000" poissonRatio="0.3" method="polar" />
<TetrahedralCorotationalFEMForceField name="CFEM" youngModulus="1000" poissonRatio="0.3" method="polar" computeVonMisesStress="1"/>

<BoxROI template="Vec3" name="box_roi" box="-6 -6 -1 30 6 0.1" drawBoxes="1" />
<FixedProjectiveConstraint template="Vec3" indices="@box_roi.indices" />
Expand Down
Loading