Partial rigidification of a deformable object

[AliceCatalano]

This discussion is created to have insights about the partial rigidification of a deformable ogject. With this I mean creating rigid ares in a mesh representing a deformable object. In my study case, for example, I have a deformable cube, on which I want to define 4 rigid planes, which orientation will be controlled by values coming from four IMUs.
The objective is to detect deformations with the IMU and seeing it represented in the simulation in the most realistic way possible. The deformations are imposed by simply pass as normal vector of the 4 rigid planes, the normal vector computed by the IMUs.

Untill now I succeded in the control part but not in the rigidification process and for this reason the simulation is not realistic at all. Deformations are happening even inside the constrained planes and the motion doesn't spread efficiently.

[alxbilger]

You should have a look on examples/Demos/rigidifiedSectionsInBeam.scn. The trick is to use a Rigid MechanicalObject to represent each of your planes, then use a RigidMapping.

[AliceCatalano]

Ah ok, sorry for the silly error in the code, I will go fix it immediately.
The problem with the original scn is still there. I just copied and pasted the code but still has the problem on the SparseLDLSolver. Should I try to change the solver?

[alxbilger]

The scene in v22.12 does not have a SparseLDLSolver, but a CGLinearSolver. If you want this scene to work with LDL, you need v23.06 (or use a MechanicalMatrixMapper but I don't recommend it)

[alxbilger]

To be more precise:
CGLinearSolver with template GraphScattered (this is the one in the scene v22.12) does not assemble the global matrix. It's a matrix-free solver. In opposite, SparseLDLSolver (this is the one in the scene v23.06) always assembles the global matrix. Therefore, there are two paths in the simulation depending on the solver you choose. Before v23.06, SOFA did not support all cases of mappings with assembling solvers. And it was the case of the rigidification. To make it work anyway, the trick was to use an additional component called MechanicalMatrixMapper. There are examples of how to use it in the SoftRobots plugin (before v23.06). But now, v23.06 supports the mappings and MechanicalMatrixMapper is no longer necessary.

[AliceCatalano]

Ok, thanks for the insight. I managed to stick the the v22.12 and implement the rigidification, but I still have some problems. I will attach here the new scene, and as you will be able to see, now I have the rigid frame in every node of my mesh, while the objective was to indicate 4 nodes of the mesh to make it delimit the planes on which I want to attach the rigid frame.

import Sofa
import os
import numpy as np
os.path.expanduser('~/Desktop/SOFA_v22.12.00_Linux/share/sofa/examples/Tutorials/Sofa_scene')
from imu_4planes import MyController

def createScene(root):
    root.gravity = [0, -10, 0]
    root.dt = 0.02
    root.animate=False
    root.addObject('VisualStyle', displayFlags="showBehaviorModels showWireframe")
    
    root.addObject('OglSceneFrame', style="Arrows", alignment="TopRight")
    root.addObject('MeshOBJLoader', name="Loaderk", filename="mesh/kidneyRM.obj")


    # Required Plugins
    root.addObject('RequiredPlugin', pluginName=["Sofa.Component.Constraint.Projective","Sofa.Component.Engine.Generate","Sofa.Component.Engine.Select","Sofa.Component.Engine.Transform",
                                                 "Sofa.Component.LinearSolver.Iterative","Sofa.Component.Mapping.Linear","Sofa.Component.Mapping.NonLinear","Sofa.Component.Mass",
                                                 "Sofa.Component.ODESolver.Backward","Sofa.Component.SolidMechanics.FEM.Elastic","Sofa.Component.StateContainer","Sofa.Component.Topology.Container.Dynamic",
                                                 "Sofa.Component.Topology.Container.Grid","Sofa.Component.Visual","Sofa.GL.Component.Rendering2D","Sofa.GL.Component.Rendering3D", "Sofa.Component.IO.Mesh",
                                                 "Sofa.Component.LinearSolver.Direct","Sofa.Component.SolidMechanics.Spring"])
    root.addObject('DefaultAnimationLoop')
    root.addObject('DefaultVisualManagerLoop')
    root.addObject('MeshGmshLoader', name="Loadercube", filename="mesh/cube.msh")

    rigid = root.addChild('Rigid')
    rigid.addObject('EulerImplicitSolver', name="odesolver", rayleighStiffness=0.1, rayleighMass=0.1)
    rigid.addObject('CGLinearSolver', name="CG Solver",template="GraphScattered", iterations=25, tolerance=1e-5, threshold=1e-10)
    rigid.addObject('TetrahedronSetTopologyContainer', name="topo", src="@../../Loadercube")
    rigid.addObject('MechanicalObject', src='@topo', template="Rigid3d", name="MechanicalModel", showObject=1, showObjectScale=0.01, showIndices=1, showIndicesScale=0.0001)
    
    square = root.addChild('Square')
    square.addObject('TetrahedronSetTopologyContainer', name="topo", src="@../../Loadercube")
    square.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="GeomAlgo")
    **square.addObject('MechanicalObject', src='@topo', template="Vec3d", name="MechanicalModel", position=[0, 0.035, 0.025], showObject=1, showObjectScale=0.01, showIndices=1, showIndicesScale=0.001)**
    square.addObject('MeshSpringForceField', name="Springs", stiffness=200, damping=10)
    square.addObject('DiagonalMass', name="Mass", massDensity=2.6)
    square.addObject('FixedConstraint', indices=[0, 1, 2, 3, 4, 5, 6, 7])

The line that is Bold is the line I thought would determinate the position of the rigid frame, but is not like this. I also tried to set up the RigidMapping, but honestly I didn't understand how to pass the indices in the right way, because in the example suggested it uses a matrix which shape i don't understand

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 
9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 
10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 
11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 
12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 
13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 
14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 
15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 
16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 
17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 
18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 
19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19

The scene is without visual model and mapping by intention because until now I just want to focus on defining just the planes I need. If you could help also with that would be amazing.

[alxbilger]

1. Define the 4 rigid DoFs

I remind the answer I gave earlier:

in the definition of your rigid mechanical object, you should define 4 positions. With rigids, a position has 7 dofs (3 for the translation, 4 for a quaternion). In total, the number of dofs is 4*7, but you defined only 3. You should get rid of src='@topo', since you don't want the rigids to be based on the topology.

So you should have something like

rigid.addObject('MechanicalObject', template="Rigid3d", position=[x0, y0, z0, qx0, qy0, qz0, qw0, x1, ...], name="MechanicalModel", showObject=1, showObjectScale=0.01, showIndices=1, showIndicesScale=0.0001)

And remove the TetrahedronSetTopologyContainer in the Node of the rigid DoFs. You don't need a topology because you have only 4 points and they are not linked together by any primitive.

2. Select the deformable DoFs mapped on the rigid DoFs

Use a BoxROI to select the deformable DoFs that you want to map on the rigid frames. You should have 4 BoxROIs (one per plane).

3. Map the deformable DoFs on the rigid DoFs

Add a RigidMapping. I agree that defining this component is not easy. The Data of interest is rigidIndexPerPoint. Its description is For each mapped point, the index of the Rigid it is mapped from. So you should work with the output of the BoxROIs, and define for each of the point of the plane, the index of the rigid frame. You may need to insert an additional Node with a SubsetMapping.

My explanation may be confusing but go step by step, and don't hesitate to ask questions again.

An alternative is to use STLIB, in particular, the rigidification module: https://github.com/SofaDefrost/STLIB/blob/master/python3/src/stlib3/physics/mixedmaterial/rigidification.py

[hugtalbot]

Hi @AliceCatalano

I finally had some time to build you a dedicated scene based on the rigidification example scene.

Here is the case you are looking at, a cube with rigidified faces:

# -*- coding: utf-8 -*-
"""
Templates to rigidify a deformable object.

The rigidification consist in mixing in a single object rigid and deformable parts.
The rigid and deformable parts are interacting together.

**Sofa Templates:**

.. autosummary::

    Rigidify

stlib3.physics.mixedmaterial.Rigidify
*******************************
.. autofunction:: Rigidify


Contributors:
        [email protected]
        [email protected]
"""

import Sofa
from splib3.numerics import Vec3, Quat, sdiv


def getBarycenter(selectedPoints):
    poscenter = [0., 0., 0.]
    if len(selectedPoints) != 0:
            poscenter = sdiv(sum(selectedPoints), float(len(selectedPoints)))
    return poscenter


def Rigidify(targetObject, sourceObject, groupIndices, frames=None, name=None):
        """ Transform a deformable object into a mixed one containing both rigid and deformable parts.

            :param targetObject: parent node where to attach the final object.
            :param sourceObject: node containing the deformable object. The object should be following
                                 the ElasticMaterialObject template.
            :param list groupIndices: array of array indices to rigidify. The length of the array should be equal to the number
                                      of rigid component.
            :param list frames: array of frames. The length of the array should be equal to the number
                                of rigid component. The orientation are given in eulerAngles (in degree) by passing
                                three values or using a quaternion by passing four values.
                                [[rx,ry,rz], [qx,qy,qz,w]]
                                User can also specify the position of the frame by passing six values (position and orientation in degree)
                                or seven values (position and quaternion).
                                [[x,y,z,rx,ry,rz], [x,y,z,qx,qy,qz,w]]
                                If the position is not specified, the position of the rigids will be the barycenter of the region to rigidify.
            :param str name: specify the name of the Rigidified object, is none provided use the name of the SOurceObject.
        """
        if frames is None:
            frames = [[0., 0., 0.]]*len(groupIndices)

        assert len(groupIndices) == len(frames), "size mismatch."

        if name is None:
            name = sourceObject.name

        # sourceObject.reinit()
        ero = targetObject.addChild(name)

        allPositions = sourceObject.container.position.value
        
        allIndices =list(range(len(allPositions)))

        rigids = []
        indicesMap = []

        def mfilter(si, ai, pts):
                tmp = []
                for i in ai:
                        if i in si:
                                tmp.append(pts[i])
                return tmp

        # get all the points from the source.
        selectedIndices = []
        for i in range(len(groupIndices)):
                selectedPoints = mfilter(groupIndices[i], allIndices, allPositions)
                if len(frames[i]) == 3:
                        orientation = Quat.createFromEuler(frames[i], inDegree=True)
                        poscenter = getBarycenter(selectedPoints)
                elif len(frames[i]) == 4:
                        orientation = frames[i]
                        poscenter = getBarycenter(selectedPoints)
                elif len(frames[i]) == 6:
                        orientation = Quat.createFromEuler([frames[i][3], frames[i][4], frames[i][5]], inDegree=True)
                        poscenter = [frames[i][0], frames[i][1], frames[i][2]]
                elif len(frames[i]) == 7:
                        orientation = [frames[i][3], frames[i][4], frames[i][5], frames[i][6]]
                        poscenter = [frames[i][0], frames[i][1], frames[i][2]]
                else:
                        Sofa.msg_error("Do not understand the size of a frame.")

                rigids.append(poscenter + list(orientation))

                selectedIndices += map(lambda x: x, groupIndices[i])
                indicesMap += [i] * len(groupIndices[i])

        otherIndices = list(filter(lambda x: x not in selectedIndices, allIndices))
        Kd = {v: None for k, v in enumerate(allIndices)}
        Kd.update({v: [0, k] for k, v in enumerate(otherIndices)})
        Kd.update({v: [1, k] for k, v in enumerate(selectedIndices)})
        indexPairs = [v for kv in Kd.values() for v in kv]

        freeParticules = ero.addChild("DeformableParts")
        freeParticules.addObject("MechanicalObject", template="Vec3", name="dofs",
                                    position=[allPositions[i] for i in otherIndices])

        rigidParts = ero.addChild("RigidParts")
        rigidParts.addObject("MechanicalObject", template="Rigid3", name="dofs", reserve=len(rigids), position=rigids)

        rigidifiedParticules = rigidParts.addChild("RigidifiedParticules")
        rigidifiedParticules.addObject("MechanicalObject", template="Vec3", name="dofs",
                                          position=[allPositions[i] for i in selectedIndices])
        rigidifiedParticules.addObject("RigidMapping", name="mapping", globalToLocalCoords=True, rigidIndexPerPoint=indicesMap)

        if "solver" in sourceObject.objects:
            sourceObject.removeObject(sourceObject.solver)
        if "integration" in sourceObject.objects:
            sourceObject.removeObject(sourceObject.integration)
        if "correction" in sourceObject.objects:
            sourceObject.removeObject(sourceObject.correction)

        sourceObject.addObject("SubsetMultiMapping", name="mapping", template="Vec3,Vec3",
                              input=[freeParticules.dofs.getLinkPath(),rigidifiedParticules.dofs.getLinkPath()],
                              output=sourceObject.dofs.getLinkPath(),
                              indexPairs=indexPairs)

        rigidifiedParticules.addChild(sourceObject)
        freeParticules.addChild(sourceObject)
        return ero


def createScene(rootNode):
        """
        """
        from stlib3.scene import MainHeader
        from stlib3.physics.deformable import ElasticMaterialObject
        from splib3.objectmodel import setData

        rootNode.addObject("DefaultAnimationLoop")
        rootNode.addObject("VisualStyle")
        rootNode.VisualStyle.displayFlags = "showBehavior"
        rootNode.addObject("RegularGridTopology", name="grid", nx="8", ny="8", nz="8", xmin="-1", xmax="1", ymin="-1", ymax="1", zmin="-1", zmax="1", drawQuads=False)

        plugins = rootNode.addChild("LoadedPlugins")
        plugins.addObject('RequiredPlugin', name='Sofa.Component.Mass') # Needed to use components [UniformMass]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.SolidMechanics.FEM.Elastic') # Needed to use components [TetrahedronFEMForceField]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.StateContainer') # Needed to use components [MechanicalObject]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.Topology.Container.Dynamic') # Needed to use components [TetrahedronSetGeometryAlgorithms,TetrahedronSetTopologyContainer,TetrahedronSetTopologyModifier]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.Topology.Container.Grid') # Needed to use components [RegularGridTopology]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.Topology.Mapping') # Needed to use components [Hexa2TetraTopologicalMapping]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.Visual') # Needed to use components [VisualStyle]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.Constraint.Projective') # Needed to use components [FixedConstraint]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.LinearSolver.Iterative') # Needed to use components [CGLinearSolver]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.Mapping.Linear') # Needed to use components [SubsetMultiMapping]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.Mapping.NonLinear') # Needed to use components [RigidMapping]  
        plugins.addObject('RequiredPlugin', name='Sofa.Component.ODESolver.Backward') # Needed to use components [EulerImplicitSolver]  

        ###########################################################
        # Define the deformable object (cube)
        tetraTopo = rootNode.addChild("TetrahedralTopology")
        tetraTopo.addObject('TetrahedronSetTopologyContainer', position="@../grid.position",name='container')
        tetraTopo.addObject('TetrahedronSetTopologyModifier', name='modifier')
        tetraTopo.addObject('Hexa2TetraTopologicalMapping', name='mapping')
        
        deformableCube = tetraTopo.addChild("DeformableCube")
        deformableCube.addObject('EulerImplicitSolver', name='integration')
        deformableCube.addObject('SparseLDLSolver', name="solver", template='CompressedRowSparseMatrixd')
        deformableCube.addObject('TetrahedronSetTopologyContainer', src="@../container",name='container')
        deformableCube.addObject('TetrahedronSetGeometryAlgorithms',  template="Vec3d", name="GeomAlgo", drawTetrahedra=False, drawScaleTetrahedra=0.5)
        deformableCube.addObject('MechanicalObject', template='Vec3', position="@../../grid.position",name='dofs', showObject=False, showObjectScale=10)
        deformableCube.addObject('UniformMass', totalMass=1., name='mass')
        deformableCube.addObject('TetrahedronFEMForceField', template='Vec3',
                                         method='large', name='forcefield',
                                         poissonRatio=0.4, youngModulus=5e3)

        rootNode.init()

        ###########################################################
        # Rigidification of the elasticobject for given indices with given frameOrientations.
        o = Rigidify(tetraTopo,
                     deformableCube,
                     name="RigidifiedStructure",
                     frames=[[0, 0, 1], [1, 0, 0]],                            # HERE you define the rigid frames located in the middle of the faces
                     groupIndices=[[475,473,483,484], [223, 231,287,295]])     # HERE you define which are the DOFs of the deformable associated to each rigid frame
        ###########################################################

        # Activate some rendering on the rigidified object.
        setData(o.RigidParts.dofs, showObject=True, showObjectScale=1, drawMode=2)
        setData(o.RigidParts.RigidifiedParticules.dofs, showObject=True, showObjectScale=0.1,
                drawMode=1, showColor=[1., 1., 0., 1.])
        setData(o.DeformableParts.dofs, showObject=True, showObjectScale=0.1, drawMode=2)
        o.RigidParts.addObject("FixedConstraint", indices=0)

        simulationNode = rootNode.addChild("Simulation")
        simulationNode.addObject("EulerImplicitSolver")
        simulationNode.addObject("CGLinearSolver", iterations=25, tolerance=1e-5, threshold=1e-5)
        simulationNode.addChild(o)
        return rootNode

All you need to run this scene (with runSofa) is to have the PYTHONPATH pointing to the SofaPython3 sites-packages and the STLIB sites-packages