SOFA scene with imposed displacements and Von Mises stress computation in 2D

[pvillard]

I am trying to write a simple SOFA scene including imposed displacements to a given geometry and export the results with the Von Mises stress for each node.

My scene is currently the following:

<?xml version='1.0'?>
<Node name="root" dt="0.001" gravity="0 0 0" > 
    <VisualStyle displayFlags="hideVisualModels hideBehaviorModels hideCollisionModels hideMappings showForceFields showInteractionForceFields showWireframe" />	
    <BackgroundSetting color="1 1 1 1" />    
        <Node name="object">
            <EulerImplicit name="cg_odesolver" printLog="false"  rayleighStiffness="0.01" rayleighMass="0.01" />
            <CGLinearSolver />
            <MeshObjLoader name="meshLoader" filename="object.obj" />
            <TriangleSetTopologyContainer name="Container" src="@meshLoader" />
            <TriangleSetGeometryAlgorithms name="GeomAlgo"/>
            <MechanicalObject name="MO" showIndicesScale="0.01" showIndices="false" />
            <DiagonalMass massDensity="1050" /> <!-- 1.05 g/cm^3 according to https://bionumbers.hms.harvard.edu/files/Density%20and%20mass%20of%20each%20organ-tissue.pdf -->
            <TriangularFEMForceField name="FEM" youngModulus="100000" poissonRatio="0.3" method="large" /> <!-- data from the literature -->
            <MeshObjLoader name="displ" filename="displ.obj"/>
            <MeshObjLoader name="indices" filename="indices.obj"/>
            <AffineMovementConstraint indices="@indices.position" meshIndices = "@indices.position" translation="@displ.position"   drawConstrainedPoints="1"/>
        </Node> 
</Node>

1) To impose the displacements, I have a list of indices given by the file “indices.obj”, I have the 2D displacement vectors given by the file “displ.obj”. I thought the command to use would be AffineMovementConstraint but it doesn’t seem to work as the translation should only be a 3D vector and not a list of 2D vectors as I have. Anyway, by using that, I have continuous displacements that are never-ending. What is the proper way to impose displacement on a list of nodes within a SOFA scene?

2) I didn’t find how to export the stresses from the WriteState export command. Is there a way to do it without modifying SOFA code and recompiling it?

Thank you in advance for any answer.

Best,

Pierre

[Hugo]

Hey @pvillard

1) The field indices of the AffineMovementConstraint is a list of integers corresponding to the indices of the constraint DOFs. I therefore don’t think that providing the field position of the loader as the input of indices is appropriate. If you want to constraint all points, you can create a box bounding the object and use the output indices.
As you noticed, the field translation is a Vector3 (you can make this translation in plane, i.e. 2D) and the rotation is a RotationMatrix. This one transformation will be applied on all selected indices.

Do you want to constraint a different displacement for each point?

2) WriteState indeed exports the MechanicalState. To export any data, I would recommand to use the VTKExporter.

Best,

Hugo

[pvillard]

Hi Hugo,

Thank you very much for your (fast) answer.

1) Instead of trying to assign my 974 imposed displacements with 1 line (passing the indices with an obj file threw the attribute position), I wrote 974 lines as follow:
<AffineMovementConstraint indices="i" meshIndices = "i" translation="-x -y 0" drawConstrainedPoints="1"/> where the vertex of index i has the displacement (x,y). I still have no warning and no error and the displayed vertices with the AffineMovementConstraint are the good ones. I have the same problem in the simulation as before.

2) I managed to export the nodal forces with:
<VTKExporter name="vTKExporter" listening="true" pointsDataFields="MO.force" .... How can I get the stress instead?
Otherwise, is it possible to directly visualize Von-Mises stresses within runSofa?

Thank you,

Pierre

[Noura]

Hello Pierre,

I think that you should be able to visualize the stresses using the following parameters from “TriangularFEMForceField.h”

Data<bool> showStressValue;
Data<bool> showStressVector;

Try this:

<TriangularFEMForceField name="FEM" youngModulus="100000" poissonRatio="0.3" method="large"  showStressValue="true" showStressVector="false" /> 

You can refer to the function “computePrincipalStress” from “TriangularFEMForceField.h” to verify if this is what you need.

Otherwise, for Von Mises stresses, I’m not sure if they are implemented for the triangular FF. In the case of a tetrahedron FF you can visualize them for a FF assigned on a given volume mesh by:

<TetrahedronFEMForceField template="Vec3d" name="FEM"  printLog="1"  method="svd"  poissonRatio="0.48"  youngModulus="16.34" computeVonMisesStress="1"  isToPrint="false" showVonMisesStressPerNode="false" showStressColorMap="blue 1 0 0 1  1 0.5 0.5 1" updateStiffness="true" /> 

Hope this helps!

Noura

[pvillard]

Hi Noura,

Thank you very much for your help. I have missed these boolean data in TriangularFEMForceField. Now, I have a stress visualization within runSOFA.

Do you know by any chance how I can export these data? I know I have to go through a VTKExporter and I guess from a field of FEM (the name of my TriangularFEMForceField) but I didn’t figure out one that works

Thank you in advance,

Pierre

[Noura]

Hi Pierre,

I have no clue how to export the stress values directly via the “VTKExporter” component. Maybe it is possible !
Also, I’m not sure if it is possible to get the stresses directly in the scene since that there is no data field called “stress” in the “TriangularFEMForceField” (FF).

Nevertheless, the stress values can be accessed and exported to an output stream by slightly modifying Sofa code.

I’ll explain in detail how to do it, but it is really simple.

1- There is a parameter called f_computePrincipalStress in “TriangularFEMForceField.h”

2- The previous parameter is not activated by default in the constructor of the FF in “TriangularFEMForceField.inl”

,f_computePrincipalStress(initData(&f_computePrincipalStress,false,"computePrincipalStress","Compute principal stress for each triangle"))

3- In your scene, assign “true” to the attribute “computePrincipalStress” in the FF
<TriangularFEMForceField name="FEM" youngModulus="10000" poissonRatio="0.3" method="large" showStressValue="true" computePrincipalStress = "true" />

4- Note that a handy class called TriangleInformation is defined in “TriangularFEMForceField.h”. The latter contains many information related to the FF (including the stresses).
Further, this class has an output stream operator operator<< which enables writing/exporting the values to an output file for instance.

5- Then, you have to seek where to export them in Sofa code. An option would be to do that in the “AddForce” function where the boolean is examined, there you have

if (f_computePrincipalStress.getValue())
...
for(unsigned int i=0; i<nbTriangles; ++i)
            computePrincipalStress(i, triangleInf[i].stress);
...

You may replace it by

...
std::ofstream file_principal_stresses ("stress.csv",std::ios::app );
if (f_computePrincipalStress.getValue())
for(unsigned int i=0; i<nbTriangles; ++i){
            computePrincipalStress(i, triangleInf[i].stress);
            file_principal_stresses<<triangleInf[i].stress<<",";
}
...

Finally, that was a quick way to get them, BUT I would not do it that way to avoid unnecessary writing to a stream each time “AddForce” is called!
Instead, you can define an additional parameter in “TriangularFEMForceField.h” to be activated only when writing to a file is desired

Data<bool> f_exportStressValuesToFile;

Then, you initialize it in “TriangularFEMForceField.inl” as all other attributes:

TriangularFEMForceField<DataTypes>::TriangularFEMForceField()
...
,f_exportStressValuesToFile(initData(&f_exportStressValuesToFile,false,"exportPrincipalStress","Export principal stress for each triangle"))
...

You include exportPrincipalStress = “true” in the FF component in the scene, and you check if it is true somewhere in the code to write the stresses or not.
An additional parameter to hold the name of the file can be added in the same manner.

Hope that this solves your issue in case you couldn’t figure out how to do it in a simpler way, otherwise I would be glad to learn how 🙂

Noura

[Hugo]

Hey @pvillard

@noura is right, there is no data allowing to access the stress (for now). It can therefore not be exported (VTK or other) as is. But the change proposed by Noura would be nice.

@noura have you already implemented it? It would be worth sharing in a PR.

Hugo

[Noura]

Hi @hugo,

I haven’t implemented it. I only had a look at the code to propose a quick solution. But I’ll try to dedicate some time to implement it (because I have to upgrade to the latest SOFA version first).

Noura

[Rayan]

Hi Pierre,

I don’t know if it could help you, but I’ve wrote a component, linked to a mechanical object and which returns the current strain applied on this mechanical state. I am using deformation gradient. If you want to get the stress, maybe you can easily link it to the strain, depending on the model used.

Hope it will help you,

Rayan

[pvillard]

Hi Rayan,

Thank you for your answer. I am not sure if I will use it but can you share it just in case?

Best,

Pierre

[Hugo]

Hi @rayan

Do not hesitate to keep Pierre updated about your work 😉
Cheers

Hugo

[Rayan]

Hi,

Sorry for the wait 🙂

Here is the component:

#ifndef SOFA_COMPONENT_FORCEFIELD_CURRENTSTRAIN_INL
#define SOFA_COMPONENT_FORCEFIELD_CURRENTSTRAIN_INL

#include "CurrentStrain.h"
#include <sofa/core/behavior/ForceField.inl>
#include <SofaBaseTopology/TopologyData.inl>
#include <sofa/defaulttype/VecTypes.h>
#include <SofaBaseMechanics/MechanicalObject.h>
#include <sofa/core/ObjectFactory.h>
#include <fstream> // for reading the file
#include <iostream> //for debugging
#include <Eigen/Eigenvalues>
#include <unsupported/Eigen/MatrixFunctions>

namespace sofa
{
namespace component
{
namespace forcefield
{
using namespace sofa::defaulttype;
using namespace	sofa::component::topology;
using namespace core::topology;

template<class DataTypes> CurrentStrain<DataTypes>::CurrentStrain(sofa::core::behavior::MechanicalState<DataTypes> *ff)
    : core::objectmodel::BaseObject()
    , d_currentStrain(initData(&d_currentStrain,"CurrentStrain","Current strain per Element"))
    , l_mstate(initLink("mstate", "input MechanicalState"), ff)
    , m_topology(nullptr)
    , l_topology(initLink("topology", "link to the topology container"))
{

}

template<class DataTypes>
void CurrentStrain<DataTypes>::init()
{
    std::cout<<"I'm in init()" <<std::endl;

    if (!l_mstate)
    {
        msg_error(this) << "currentStrain not linked to any mechanical state" ;
        return;
    }

    if (l_topology.empty())
    {
        msg_info() << "link to Topology container should be set to ensure right behavior. First Topology found in current context will be used.";
        l_topology.set(this->getContext()->getMeshTopologyLink());
    }

    m_topology = l_topology.get();
    msg_info() << "Topology path used: '" << l_topology.getLinkedPath() << "'";

    if (m_topology == nullptr)
    {
        msg_error() << "No topology component found at path: " << l_topology.getLinkedPath() << ", nor in current context: " << this->getContext()->name;
        sofa::core::objectmodel::BaseObject::d_componentstate.setValue(sofa::core::objectmodel::ComponentState::Invalid);
        return;
    }

    unsigned int i,j;
    unsigned int nbTetrahedra=m_topology->getNbTetrahedra();

    typename DataTypes::Real volume;
    typename DataTypes::Coord point[4];
    const typename DataTypes::VecCoord restPosition = l_mstate->read(sofa::core::ConstVecCoordId::restPosition())->getValue();

    for(i=0; i<nbTetrahedra; i++ )
    {

        m_shapeVector.resize(nbTetrahedra);

        const vector< Tetrahedron > &tetrahedronArray= m_topology->getTetrahedra();

        ///describe the indices of the 4 tetrahedron vertices
        const Tetrahedron &t= tetrahedronArray[i];

        for(j=0;j<4;++j)
        {
            point[j]=(restPosition)[t[j]];
        }
        /// compute 6 times the rest volume
        volume=dot(cross(point[2]-point[0],point[3]-point[0]),point[1]-point[0]);

            for(j=0;j<4;++j) {
                if (!(j%2))
                    m_shapeVector[i].push_back(-cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/ volume);
                else
                    m_shapeVector[i].push_back(cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/ volume);
            }
    }

    BaseData* pos_init = l_mstate->read(sofa::core::ConstVecCoordId::position())->getData();
    c_callback.addInput(pos_init);
    c_callback.addCallback(std::bind(&CurrentStrain::strainplaceholder,this));

    }

template<class DataTypes>
void CurrentStrain<DataTypes>::reinit()
{
    this->init();
}

template <class DataTypes>
void CurrentStrain<DataTypes>::update()
{

}

template <class DataTypes>
void CurrentStrain<DataTypes>::strain(const core::MechanicalParams*, const DataVecCoord& d_x)
{

    //std::cout<<"I'm in strain" <<std::endl;

    const VecCoord& x = d_x.getValue();

    unsigned int i=0,j=0,k=0,l=0;
    unsigned int nbTetrahedra=m_topology->getNbTetrahedra();

    assert(this->mstate); // ask why

    Coord dp[3],x0,sv;

    helper::vector<double>* currentStrain = d_currentStrain.beginEdit();
    currentStrain->resize(nbTetrahedra);

    for(i=0; i<nbTetrahedra; i++ )
    {

        Matrix3 deformationGradient;

        const Tetrahedron &ta= m_topology->getTetrahedron(i);

        x0=x[ta[0]];
        dp[0]=x[ta[1]]-x0;

        sv=m_shapeVector[i][1];

        //std::cout<<"m_shapeVector vaut : "<< m_shapeVector <<std::endl;

        for (k=0;k<3;++k) {
                for (l=0;l<3;++l) {
                        deformationGradient[k][l]=dp[0][k]*sv[l];
                }
        }

        for (j=1;j<3;++j) {
                dp[j]=x[ta[j+1]]-x0;
                sv=m_shapeVector[i][j+1];
                for (k=0;k<3;++k) {
                        for (l=0;l<3;++l) {
                                deformationGradient[k][l]+=dp[j][k]*sv[l];
                        }
                }
        }

        Matrix3 Id,C;

        for (k=0;k<3;++k) {
            for (l=k;l<3;++l) {
                C(k,l) = (deformationGradient(0,k)*deformationGradient(0,l)+
                         deformationGradient(1,k)*deformationGradient(1,l)+
                         deformationGradient(2,k)*deformationGradient(2,l));
                C(l,k) = C(k,l);
            }
        }

        //Create identity
        Id = Id*0;
        Id(0,0) = 1;
        Id(1,1) = 1;
        Id(2,2) = 1;

        Matrix3 E = 1/2.0*(C-Id);

        currentStrain->operator[](i) = sqrt(pow(E(0,0),2) + pow(E(1,1),2) + pow(E(2,2),2) + pow(E(1,2),2) + pow(E(0,1),2) + pow(E(0,2),2));

//        std::cout<<"la force courante est:  " << d_currentStrain<<std::endl;

    }

    d_currentStrain.endEdit();

}

template<class DataTypes>
void CurrentStrain<DataTypes>::strainplaceholder()
{

    DataVecCoord d_pos;
    VecCoord &pos = *d_pos.beginEdit();
    pos =  this->l_mstate->read(core::ConstVecCoordId::position())->getValue();

    d_pos.setValue(pos);

    strain(core::MechanicalParams::defaultInstance(), d_pos);

    d_pos.endEdit();
}

}   //namespace forcefield
}   //namespace component
}   //namespace sofa

#endif

It’s quiet raw but don’t hesitate if you have questions 🙂

Regards,

Rayan

[bani]

Hi Rayan,
Could you please also tell how to implement it using sofapython ?

Thanks
Vani

[Hugo]

hey @vanisri

Rayan is not the area anymore, he was an intern.
You can get inspired from the C++ and translate it into python.

All you need is to have an access to the topology and the mechanical object (mstate).
Then compute at init the shapeVectors.

When you want it, compute the strain as it is done by the function strainplaceholder (itself calling the strain()) function in Rayan’s code.

Best wishes,

Hugo

[Author]

Posts