[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

[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

[Author]

Posts