[jnbrunet]

It takes care of propagating the mapped force vector into the system (top level mechanical object) force vector.

[jnbrunet]

Hey @jjcasmar,

Have a look at the StaticSolver which does a series of Newton-Raphson iterations. It should give you hints on how to achieve such state updates within the same time step.

Let us know if you hit specific difficulties and don’t hesitate to show your work to the community, this is a component that could be very useful to others.

Jean-Nicolas

[jnbrunet]

Hi @jjcasmar,

This is usually done using mechanical operators. In your case, you could create a component that does this:

sofa::core::MechanicalParams mparams (sofa::core::ExecParams::defaultInstance());
sofa::simulation::common::VectorOperations vop( &mparams, this->getContext());
sofa::simulation::common::MechanicalOperations mop( &mparams, this->getContext());
sofa::core::behavior::MultiVecDeriv F( &vop, sofa::core::VecDerivId::force());

// Accumulate the force vectors
// 1. Clear the force vector (F := 0)
// 2. Go down in the current context tree calling addForce on every forcefields
// 3. Go up from the current context tree leaves calling applyJT on every mechanical mappings
mop.computeForce(F);

Jean-Nicolas

[jnbrunet]

Hey @nguyenvanpho,

My bad, the fix was in SOFA, not SP3. Hence got to your SOFA repository and do the git pull there.

Also, make sure your SP3 origin is the one from sofa-framework, and not SofaDefrost, ie:

$ cd $SP3_SRC
$ git remote -vv
(..)
origin git@github.com:sofa-framework/plugin.SofaPython3.git
(..)

J-N

[jnbrunet]

Hey @nguyenvanpho,

The fix has been merge in the master branch. Hence, you can do

$ cd $SP3_SRC
$ git checkout master
$ git pull
$ cd $SP3_BLD
$ cmake -DCMAKE_PREFIX_PATH=$SOFA_BLD/install/lib/cmake $SP3_SRC

J-N

[jnbrunet]

Hey @jacqueline,

Awesome. We will now investigate your error with the fetch option, this isn’t normal and looks like a bug, thank you for reporting it.

J-N

[jnbrunet]

Hey @nguyenvanpho,

Simply to let you know that the merge has been made. You can now update your SP3 repository.

J-N

[jnbrunet]

Hey @jacqueline,

Yeah you got it at a bad time, there was a change made in SOFA a couple of days ago and not propagated to SofaPython3. We merge the fix in SofaPython3 yesterday morning.

It should be fine now, simply git pull the SofaPython3 repository.

JN

[jnbrunet]

Hey @antonio-8196,

Sorry for the delay. Have you solved your problem yet?

I cannot really help you for your problem on Windows.

For Linux, it seems you are missing the  libpython2.7.so.1.0 file. I’m not sure why as you seems to have python 2.7 installed. Can you give me the output of:

$ sudo updatedb
$ locate libpython2.7.so

[jnbrunet]

Hey @nguyenvanpho,

That’s the inconvenient of following so closely a development branch, there are often bug introduced! This problem was introduced recently (last week) in SOFA and makes external plugins fail to find the modules of SOFA.

There’s a fix being made right now in SOFA, see this pull request. Hopefully it should be merged during our next SOFA developer meeting, next Wednesday. If you cannot wait until then, you can try to checkout the fix branch (in SOFA, not SP3):

git remote add fredroy https://github.com/fredroy/sofa.git
git fetch fredroy
git checkout -b fix_outoftree_compil

Sorry for this. We are working hard right now to make SOFA and SP3 stable and release a first version very soon.

[jnbrunet]

Hi @nguyenvanpho,

Looks like I cannot access your image. Could you paste the output somewhere else?

JN

[jnbrunet]

Hi @nguyenvanpho,

Could you try to update both your SOFA and SofaPython3 source code from their respective master branches (i.e. git checkout master && git pull). Try cleaning any previous builds from both of them (i.e. rm -rf build/*) and restart the compilations and installations (ninja install).

[jnbrunet]

Hey @antonio-8196,

Compiling SOFA and its plugins can be difficult, especially when some of the plugins depend on external libraries such as libpython, and I know it can be sometime quite frustrating.

Note that there are two different plugins for SOFA:
1. SofaPython which comes with SOFA by default and is only compatible with Python. 2.X
2. SofaPython3 which comes from another repository and is only compatible with Python 3.x

Make sure you do not mix them as they cannot be compiled or loaded at the same time in SOFA.

From your last output, it looks like you are loading the SofaRobots plugin. Did you compiled it or downloaded it? If you successfully compiled it, it means you have everything needed to also run it. Can you go in your build directory and give us the output of grep -i libpython CMakeCache.txt

J-N

[jnbrunet]

Hey @antonio-8196,

Indeed it can be confusing. When you are building SofaPython3 within SOFA (using the SOFA_EXTERNAL_DIRECTORIES CMake option), they share the same build folder and the install command has to be executed there, so you did good here.

The alternative would be to do an “out-of-tree” build where SofaPython3 has its own build folder and compilation process. In this case you would have something like

.
├── SOFA
│   ├── src
|   └── build
└── SofaPython3
    ├── src
    └── build

Jean-Nicolas

[jnbrunet]

Hey @antonio-8196,

CMake seems to indicate that you are giving him a path which is not the one of SOFA’s source code, it can’t find the file CMakeLists.txt (which should be at the root of the SOFA’s source).

If your directories looks like :

My-SOFA
├── src
│   ├── (...)
|   └── CMakeLists.txt
└── build

Then you should go into My-SOFA/build and launch cmake with

My-SOFA/build $ cmake [OPTIONS] ../src

JN

[jnbrunet]

Hey @psomers,

Sorry your issue on github passed under the radar.

I just replied to it.

For others, here is the link to the github issue.

J-N

[jnbrunet]

Hey @rob,

The Hexa2TetraTopologicalMapping is there simply to create a tetrahedral topology from an hexahedral one (the regular grid). Hence you should not need it in your scene if you already have tetrahedrons.

You should however have a tetrahedral topology container (TetrahedronSetTopologyContainer). Make sure it is filled with tetrahedrons by double clicking on it and looking at the data field tetrahedra in the Property tab.

Also, open up the properties of the TetrahedronHyperelasticityFEMForceField (again by double clicking on it) and look at the data field *topology* in the Links tab. Does it point to your tetrahedron container? If not, you can force it with, for example:

root.meca.createObject('TetrahedronSetTopologyContainer', name='Container')
(...)
root.meca.createObject('TetrahedronHyperelasticityFEMForceField', 
    (...),
    topology='@./Container'
)

J-N

[jnbrunet]

Hey @rob,

Does runSofa displays any error or warning before the crash? Can you paste its log output here?

The following python scene works well on my side:

import Sofa

poisson_ratio = 0.45
young_modulus = 1e4
mu = young_modulus / (2. * (1. + poisson_ratio))
lm = young_modulus * poisson_ratio / ((1. + poisson_ratio) * (1. - 2.*poisson_ratio))
    
def createScene(root):
    root.createObject('VisualStyle', displayFlags=['showForceFields', 'showBehaviorModels'])
    root.createChild('meca')
    root.meca.createObject('EulerImplicitSolver')
    root.meca.createObject('CGLinearSolver')
    
    root.meca.createObject('RegularGridTopology', name='grid', min=[0,0,0], max=[1,1,2.7], n=[3,3,8])
    root.meca.createObject('MechanicalObject', name='mo')
    root.meca.createObject('TetrahedronSetTopologyContainer', name='Container')
    root.meca.createObject('TetrahedronSetTopologyModifier')
    root.meca.createObject('TetrahedronSetTopologyAlgorithms')
    root.meca.createObject('Hexa2TetraTopologicalMapping', input="@./grid", output="@./Container" )
    root.meca.createObject('TetrahedronHyperelasticityFEMForceField', 
        ParameterSet=[mu, lm],  
        materialName="StVenantKirchhoff", 
        AnisotropyDirections=[0, 0, 1]
    )
    root.meca.createObject('BoxROI', drawBoxes=True, box=[-0.05, -0.05, -0.05, 1.05, 1.05, 0.05], name="fixed_box")
    root.meca.createObject('FixedConstraint', indices="@fixed_box.indices")
    root.meca.createObject('BoxROI', drawBoxes=True, box=[-0.05, -0.05, +2.60, 1.05, 1.05, 2.7], name="traction_box")
    root.meca.createObject('TriangleSetGeometryAlgorithms')
    root.meca.createObject('TrianglePressureForceField', 
        pressure=[0, -2e2, 0], 
        triangleList='@traction_box.triangleIndices', 
        showForces=True
    )

[jnbrunet]

Hey @rob,

The python code looks fine. You will probably have to set up correctly the material parameters (mu and lambda) so that it makes senses with the magnitude of your applied external forces . If you had found a Young modulus and Poisson ratio that worked well with the corotational model, start with these values using the equation I gave you to get the equivalent Lamé parameters.

Remove the anisotropy to make sure the crash doesn’t come from that.

Hyperelastic materials might also be harder to reach convergence, you can play your time step to hopefully find what is going on.

J-N

[jnbrunet]

Hey @rob,

The TetrahedronFEMForceField implements a linear isotropic material, which means that it offers the same resistance to strain in every directions. It also means that it can only represent small and linear displacements. Since rotations are non-linear rigid motions, it will generate wrong elastic forces if your object is rotated. This is why you have an over-expansion.

When you set the method to “large”, it uses a corotational approach, which is still isotropic linear elasticity, but tries to remove the rotation before computing the elastic forces).

If you want anisotropy (more resistance to stretch in some given direction), or if you want a true non-linear material that is not sensible to rigid motions, you might have a look at the TetrahedronHyperelasticityFEMForceField component. For example,

<TetrahedronHyperelasticityFEMForceField ParameterSet="3448.2759 31034.483" AnisotropyDirections="0 1 0" materialName="StVenantKirchhoff"/>

Here the ParameterSet here uses, in order, the mu and lambda parameters instead of the Young modulus and Poisson ratio used in the TetrahedronFEMForceField. You can use this formulae to pass from one parameter set to the other:

mu     = young_modulus / (2. * (1. + poisson_ratio))
lambda = young_modulus * poisson_ratio / ((1. + poisson_ratio) * (1. - 2.*poisson_ratio))

Other examples can be found in the file “examples/Components/forcefield/TetrahedronHyperelasticityFEMForceField.scn”

J-N

[jnbrunet]

Yeah I completely understand your issue. Like I said, I have no clue on where the migration towards SofaPython3 of STLIB and SofaRobots stands right now. I am sure the devs are eager to do the swap to SP3, but they might be running low on available time.

We will have to wait for @damien-marchaluniv-lille1-fr or @bmarques or someone in the Defrost team to give us some feedback.

[jnbrunet]

Hey Paul,

I never tried, by I don’t see anything that would prevent changing the color of an OglModel during runtime.

However, as far as I know, OglModel is using a texture map, so you would have to provide it with an image containing the different colors of your nodes, and assign the texture coordinates to each (or some) nodes during runtime. The texture image must be provided by a filename and is loaded during the init phase, so you probably won’t be able to change it during runtime.

In your example, you can have a look at the texture coordinates of your OglModel nodes with:

print(visual.VisualModel.texcoords.value)

and you can change the texture coordinates of a subset of nodes with, for example,

blue_nodes = [0, 5, 6, 7] # indices of nodes that must be blue
blue_coordinates = [0.7, 0.9] # coordinates of a blue pixel in the texture image
with visual.VisualModel.texcoords.writeableArray() as wa:
    wa[blue_nodes] = blue_coordinates    

J-N

Edit:
Don’t forget that the color of a triangle is interpolated from the color of its nodes.
If you want to change the color of triangles, I guess you could do something like

triangles = visual.VisualModel.triangles.value
blue_triangles = [0, 4, 9] # indices of triangles that must be blue
blue_nodes = triangles[blue_triangles]
blue_coordinates = [0.7, 0.9] # coordinates of a blue pixel in the texture image
with visual.VisualModel.texcoords.writeableArray() as wa:
    wa[blue_nodes] = blue_coordinates    

[jnbrunet]

Hey Paul,

Looks like the OglModel needs a special visual initialization: Sofa.Simulation.initVisual(marker). The following scene works for me, it adds a sphere at each time step.

import Sofa
import Sofa.Simulation

def Marker(rootNode, name, pose):

    marker = rootNode.addChild(name)
    marker.addObject(
        "MechanicalObject", name="MechanicalModel", template="Rigid3", position=pose
    )

    visual = marker.addChild("VisualModel")

    visual.loader = visual.addObject(
        "MeshObjLoader",
        handleSeams="1",
        name="VisualMeshLoader",
        scale=1,
        filename="mesh/ball.obj",
    )

    visual.addObject(
        "OglModel", src="@VisualMeshLoader", name="VisualModel", color=[1, 1, 1]
    )

    visual.addObject(
        "RigidMapping",
        input="@../MechanicalModel",
        name="MM->VM mapping",
        output="@VisualModel",
    )
    print(f"Created {visual} with name {name}!")
    
    marker.init()
    Sofa.Simulation.initVisual(marker)

    return marker
    
class Controller(Sofa.Core.Controller):

    def __init__(self):
        Sofa.Core.Controller.__init__(self)
        self.count = 0
        
    def onAnimateBeginEvent(self, e):
        Marker(self.getContext(), f'marker_{self.count}', [self.count*2, 0, 0, 0, 0, 0, 1])
        self.count+=1
    
def createScene(root):
    root.bbox = [[-1, -1, -1], [1, 1, 1]]
    root.addObject(Controller())

J-N

[jnbrunet]

Hey @jlorenze,

Unfortunately, SofaPython(2) is indeed incompatible with SofaPython3 as it would be a lot of work to make sure there isn’t mix-ups between the two when they are both loaded in memory, and I’m not even sure python’s libraries allow it.

For the migration SoftRobot and STLIB to SofaPython3, maybe @damien-marchaluniv-lille1-fr or @bmarques would have an idea of the current progression state.

Converting a scene from SofaPython(2) to SofaPython3 doesn’t look too difficult. You can have a look at this PR for an example.

Small note: The SofaPython3 repository from SofaDefrost Github has been deprecated as its development is now part of SOFA consortium, make sure you are using the new repository here.

J-N

[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

[Author]

Posts