Strain Visualization in a Color Map.

[Pasquale94]

Hello everyone,

I am working on two simulations of soft robots, one using TetrahedronFEMForceField and one using TetrahedronHyperelasticityFEMForceField , in particular the Ogden model.
I want to visualize the strains along my structures, in a colour map,while the robot is inflating.
I already seen how to create a colour map, but for now I am able to create only for displacements with the command EvalPointsDistance, and I would like to do this for the strains, I checked also in the forum and the more relate topics talk about stresses not strains.
Did anyone risolved this issue?
Thanks in advance

[Pasquale94]

Hello, I saw a little bit the codes TetrahedronFEMForceField.h and TetrahedronHyperelasticityFEMForceField.h.

In the first one there is this function:

void computeStrainDisplacement( StrainDisplacement &J, Coord a, Coord b, Coord c, Coord d );

that is defined in the TetraHedronFEMForceField.inl. Anf I don’t understand how (or if it is possible) to visualize the strain from this function.

instead for the TetrahedronHyperelasticityFEMForceField.h, there is the class :

class TetrahedronRestInformation : public fem::StrainInformation<DataTypes>
    {
    public:
        /// shape vector at the rest configuration
        Coord m_shapeVector[4];
        /// fiber direction in rest configuration
        Coord m_fiberDirection;
        /// rest volume
        Real m_restVolume;
        /// current tetrahedron volume
        Real m_volScale;
        Real m_volume;
        /// volume/ restVolume
        MatrixSym m_SPKTensorGeneral;
        /// deformation gradient = gradPhi
        Matrix3 m_deformationGradient;
        /// right Cauchy-Green deformation tensor C (gradPhi^T gradPhi)
        Real m_strainEnergy;

        /// Output stream
        inline friend ostream& operator<< ( ostream& os, const TetrahedronRestInformation& /*eri*/ ) {  return os;  }
        /// Input stream
        inline friend istream& operator>> ( istream& in, TetrahedronRestInformation& /*eri*/ ) { return in; }

        TetrahedronRestInformation() {}
    };

but I don’t now how to use it. I wrote simply like this :

        finger.createObject('TetrahedronHyperelasticityFEMForceField',name = 'FEM', materialName="Ogden", ParameterSet=str(k0_) + " "+str(mu_) +' '+ str(alpha1_), AnisotropyDirections=" ", tetrahedronInfo = 'true')        

But the scene in sofa crashes.

[jnbrunet]

Hey @Pasquale94,

Unfortunately, I am not aware of any way to do this at the current state of SOFA.

The “tetrahedronInfo” data field you used is supposed to be a vector<TetrahedronRestInformation> having one entry per tetrahedral element. Normally, to set or to get this data vector, the method ostream& operator<< and istream& operator>> in TetrahedronRestInformation should be implemented. But as you can see, they are doing absolutely nothing here. This is probably why your scene crash.

This question comes up quite often on the forum, so if you are comfortable enough in C++, you could create a plugin containing a single component that does the visualization. This would be very useful for many here. If you (or someone that read this) are up to the challenge, here is a small draft I just did that could help you start:

#include <sofa/core/topology/BaseMeshTopology.h>
#include <sofa/core/objectmodel/BaseObject.h>
#include <SofaMiscFem/TetrahedronHyperelasticityFEMForceField.h>
#include <sofa/core/visual/VisualParams.h>
#include <sofa/helper/accessor.h>

class StrainVisualization : public sofa::core::objectmodel::BaseObject {
public:
    using Vec3Types = sofa::defaulttype::Vec3Types;
    using HyperelasticForcefield = sofa::component:forcefield::TetrahedronHyperelasticityFEMForceField<Vec3Types>; 
    template <typename ObjectType>
    using Link = SingleLink<StrainVisualization, ObjectType, BaseLink::FLAG_STRONGLINK>;
    
    StrainVisualization () 
    :  d_forcefield(initLink("forcefield", "The TetrahedronHyperelasticityFEMForceField component that will be used for the strain map."))
    ,  d_topology_container(initLink("topology_container", "The topology that contains the tetrahedrons."))
    ,  d_mechanical_state(initLink("mechanical_state", "The mechanical object that contains the mesh position."))
    {}

    void init() {
        if (! d_forcefield.get()) {
            auto ff = this->getContext()->template get<HyperelasticForcefield>();
            if (ff)
                d_forcefield.set(ff);
            else
                msg_error() << "No TetrahedronHyperelasticityFEMForceField provided or found in the current context.";
        }
        
        if (d_forcefield.get()) {
            HyperelasticForcefield * forcefield = d_forcefield.get();
            d_topology_container.set(forcefield->getContext()->template get<sofa::core::topology::BaseMeshTopology>());
            d_mechanical_state.set(forcefield->getContext()->template get<sofa::core::behavior::MechanicalState<Vec3Types>);
        }
    }
    
    void draw(const sofa::core::visual::VisualParams* vparams) {
        using Vector3 = sofa::core::visual::DrawTool::Vector3;
        using Color = sofa::core::visual::DrawTool::RGBAColor;
        
        if (! d_forcefield.get() || !d_topology_container.get() || !d_mechanical_state.get())
            return;
        
        HyperelasticForcefield * forcefield = d_forcefield.get();
        auto * tetrahedrons_info_data = dynamic_cast< HyperelasticForcefield::TetrahedronData<sofa::helper::vector<TetrahedronRestInformation>> * >(
            forcefield->findData("tetrahedronInfo")
        );
        
        const sofa::helper::vector<TetrahedronRestInformation> & tetrahedrons_info = sofa::helper::getReadAccessor(tetrahedrons_info_data);
        if (tetrahedrons_info.empty())
            return;
            
        sofa::helper::vector<Vector3> triangles_points;
        sofa::helper::vector<Color> triangles_points_color;
        
        sofa::helper::vector<HyperelasticForcefield::Matrix3> nodal_strain(d_mechanical_state->size(), HyperelasticForcefield::Matrix3(0.));
        sofa::helper::vector<Color> nodal_color(d_mechanical_state->size());
        
        const auto number_of_tetrahedrons = tetrahedrons_info.size();
        const auto & x = d_mechanical_state->read(sofa::core::ConstVecCoordId::position())->getValue();
        const auto Id = HyperelasticForcefield::Matrix3::Identity();
        
        // Step 1: accumulate the strain on each node
        for (std::size_t tetra_id = 0; tetra_id < number_of_tetrahedrons; ++tetra_id) {
            const auto & node_indices = d_topology_container->getTetrahedron(tetra_id);
            const auto & tetra_info = tetrahedrons_info[tetra_id];
            const auto & F = tetra_info.m_deformationGradient; // Deformation tensor F at the center of the tetra
            const auto E = 0.5(F*F.transposed() - Id); // Strain tensor E at the center of the tetra
        
            nodal_strain[node_indices[0]] += E/4;
        }
        
        // Step 2: fill-in the color of the nodes from their strain
        for (std::size_t node_id = 0; node_id < nodal_color.size(); ++node_id) {
            const auto & E = nodal_strain[node_id];
            nodal_color[node_id] = color_ramp(E); // You will have to implement the color_ramp here
        }
        
        // Step 3: fill-in the triangle nodes vector
        for (std::size_t tetra_id = 0; tetra_id < number_of_tetrahedrons; ++tetra_id) {
            const auto & node_indices = d_topology_container->getTetrahedron(tetra_id);
            
            const auto & node_0 = node_indices[0],
                         node_1 = node_indices[1],
                         node_2 = node_indices[2],
                         node_3 = node_indices[3];
        
            // Triangle 1
            triangles_points.push_back(x[node_0]);
            triangles_points.push_back(x[node_1]);
            triangles_points.push_back(x[node_2]);
            triangles_points_color.push_back(nodal_color[node_0]);
            triangles_points_color.push_back(nodal_color[node_1]);
            triangles_points_color.push_back(nodal_color[node_2]);
            
            // Triangle 2
            triangles_points.push_back(x[node_1]);
            triangles_points.push_back(x[node_2]);
            triangles_points.push_back(x[node_3]);
            triangles_points_color.push_back(nodal_color[node_1]);
            triangles_points_color.push_back(nodal_color[node_2]);
            triangles_points_color.push_back(nodal_color[node_3]);
            
            // Triangle 3
            triangles_points.push_back(x[node_2]);
            triangles_points.push_back(x[node_3]);
            triangles_points.push_back(x[node_0]);
            triangles_points_color.push_back(nodal_color[node_2]);
            triangles_points_color.push_back(nodal_color[node_3]);
            triangles_points_color.push_back(nodal_color[node_0]);
            
            // Triangle 4
            triangles_points.push_back(x[node_3]);
            triangles_points.push_back(x[node_0]);
            triangles_points.push_back(x[node_1]);
            triangles_points_color.push_back(nodal_color[node_3]);
            triangles_points_color.push_back(nodal_color[node_0]);
            triangles_points_color.push_back(nodal_color[node_1]);
        }
        
        vparams->drawTool()->drawTriangles(triangles_points, triangles_points_color);
    }

    // Data fields
    Link<HyperelasticForcefield> d_forcefield;
    Link<sofa::core::topology::BaseMeshTopology> d_topology_container;
    Link<sofa::core::behavior::MechanicalState<Vec3Types>> d_mechanical_state;
}

J-N

[Pasquale94]

Hello @jnbrunet,

Yeah I am able to code in C++, but I don’t know pretty much of SOFA API, because is few months that i use it, and now I am too busy to dedicate time to this.
Anyway, in Sofa there is already a command that create the colour map, so I think that the function have to be similar to ‘computevonMisesStresses’ and ‘StressperNode’ attributes that there are present in the TetraHedronForceField class.

So I think that would be nice that for both forcefields, Linear Elastic and Hyperelastic one, there would be accessible datafields for stress and strains, not only for the visualization but also to export them, in a file and moreover in a ROS topic, that would great for example for the control (in dynamics) of soft robots, for grasping and in-hand manipulation research-field.

Thank you anyway for the information.

Pasquale.

[Author]

Posts