Issue with TopologyDataHandler initialization and visualization in co-rotational FEM model

[sodasoda520]

Hello,

I'm working on a soft body simulation using co-rotational FEM in SOFA, and I've encountered an interesting issue regarding initialization and visualization.

Current Situation:

1. I have a soft body model using TetrahedralCorotationalFEMForceField
2. The model includes rigid parts created using Rigidify
3. The scene setup includes both mechanical and visual models

Issue Description:
When I run my scene, I face a dilemma:

Case 1 (Without initialization):

def rigidify(self):
    deformableObject = self.ElasticBody.MechanicalModel
    # No self.ElasticBody.init() here
    # ... rest of the rigidify code

• No error messages
• But the model is not visible in the scene, despite setting displayFlags="showVisual showBehavior showForceFields"

Case 2 (With initialization):

def rigidify(self):
    deformableObject = self.ElasticBody.MechanicalModel
    self.ElasticBody.init() # Adding this line
    # ... rest of the rigidify code

• Model becomes visible
• But generates these errors:

[ERROR] [UniformMass(mass)] TopologyData: indices already has a TopologyDataHandler. createTopologyHandler should only be called once at init of the TopologyData.
[WARNING] [UniformMass(mass)] totalMass value overriding the value of the attribute vertexMass.
[ERROR] [TetrahedralCorotationalFEMForceField(corotationalFEM)] TopologyData: tetrahedronInfo already has a TopologyDataHandler. createTopologyHandler should only be called once at init of the TopologyData.

Questions:

1. What is the correct way to handle initialization in this case?
2. Why is the model not visible without explicit initialization?
3. How can I properly initialize the topology without causing the TopologyDataHandler errors?

Relevant Code:

# -*- coding: utf-8 -*-
# Implementation of a soft body model using co-rotational FEM method
import Sofa.Core
import math
import numpy as np
from mixedmaterial import Rigidify
from components import addOrientedBoxRoi
from splib3.numerics import vec3
from splib3.numerics.quat import Quat
from tutorial import *
from splib3.animation import animate
from stlib3.scene import Scene
from fixingbox import FixingBox

def ElasticBody(name="ElasticBody", rotation=[0, 0, 0], translation=[0, 0, 0], color=[0.82, 0.41, 0.12]):
    # Create an elastic body using co-rotational FEM
    self = Sofa.Core.Node(name)
    mechanicalmodel = self.addChild("MechanicalModel")
    # Load mesh and setup topology
    mechanicalmodel.addObject('MeshGmshLoader',
                              name='loader',
                              rotation=rotation,
                              translation=translation,
                              filename="mesh/dense.msh")
    mechanicalmodel.addObject('TetrahedronSetTopologyContainer',
                              src='@loader',
                              name='container')
    mechanicalmodel.addObject('TetrahedronSetGeometryAlgorithms',
                              template="Vec3d")
    mechanicalmodel.addObject('MechanicalObject',
                              name='dofs',
                              template="Vec3d",
                              position=mechanicalmodel.loader.position.getLinkPath(),
                              showObject=False,
                              showObjectScale=5.0)
    mechanicalmodel.addObject('UniformMass',
                              name="mass",
                              totalMass=0.00355)
    mechanicalmodel.addObject('TetrahedralCorotationalFEMForceField',
                              name="corotationalFEM",
                              template="Vec3d",
                              method="large",
                              poissonRatio=0.34,
                              youngModulus=2500)
    visualmodel = self.addChild("VisualModel")
    visualmodel.addObject('MeshGmshLoader',
                           name='loader',
                           filename="mesh/dense.msh",
                           rotation=rotation,
                           translation=translation)
    visualmodel.addObject('OglModel',
                           name='renderer',
                           src='@loader',
                           color=color)
    visualmodel.addObject('IdentityMapping',
                           input='@../MechanicalModel/dofs',
                           output='@renderer')
    return self

def SOFT(name="SOFT"):  
    def rigidify(self):
        #self.ElasticBody.init()
        """
        Without this line of code, there will be no error after importing, but it cannot be rendered, even though showVisual showBehavior showForceFields has been set up.
Adding this line will result in an error:
            [ERROR]   [UniformMass(mass)] TopologyData: indices already has a TopologyDataHandler. createTopologyHandler should only be called once at init of the TopologyData.
            [WARNING] [UniformMass(mass)] totalMass value overriding the value of the attribute vertexMass.
            vertexMass = totalMass / nb_dofs.
            To remove this warning you need to set either totalMass or vertexMass data field, but not both.
            [ERROR]   [TetrahedralCorotationalFEMForceField(corotationalFEM)] TopologyData: tetrahedronInfo already has a TopologyDataHandler. createTopologyHandler should only be called once at init of the TopologyData.
        """
        deformableObject = self.ElasticBody.MechanicalModel
        box = addOrientedBoxRoi(
            self,
            position=deformableObject.dofs.position.value,
            name="RigidBoxROI",
            translation=[20, 2.5, 0.125],
            scale=[1.2, 6, 0.4],
            drawBoxes=True
        )
        box.init()
        rigidIndices = [ind for ind in box.indices.value]
        rigidFrame = [20, 2.5, 0.125, 0.0, 0.0, 0.0, 1.0]
        Rigidify(
            self,
            sourceObject=deformableObject,
            groupIndices=[rigidIndices],
            frames=[rigidFrame],
            name="RigidifiedStructure"
        )

    def addFixedConstraint(self):
        self.fixedConstraint = FixingBox(
            self,
            self.ElasticBody.MechanicalModel,
            scale=[1.2, 6, 0.4],
            translation=[0, 2.5, 0.125]
        )
        self.fixedConstraint.BoxROI.drawBoxes = True

    self = Sofa.Core.Node(name)
    elasticBody = ElasticBody(color=[0.82, 0.41, 0.12])
    self.addChild(elasticBody)
    self.ElasticBody = elasticBody
    rigidify(self)
    addFixedConstraint(self)
    return self

def createScene(rootNode):
    rootNode.addObject("RequiredPlugin", pluginName=[
        "ArticulatedSystemPlugin",
        "Sofa.Component.AnimationLoop",
        "Sofa.Component.Constraint.Lagrangian.Correction",
        "Sofa.Component.Constraint.Lagrangian.Solver",
        "Sofa.Component.Constraint.Projective",
        "Sofa.Component.Engine.Select",
        "Sofa.Component.IO.Mesh",
        "Sofa.Component.LinearSolver.Direct",
        "Sofa.Component.Mass",
        "Sofa.Component.ODESolver.Backward",
        "Sofa.Component.SolidMechanics.FEM.Elastic",
        "Sofa.Component.SolidMechanics.Spring",
        "Sofa.Component.StateContainer",
        "Sofa.Component.Topology.Container.Constant",
        "Sofa.Component.Topology.Container.Dynamic",
        "Sofa.Component.Visual",
        "Sofa.GL.Component.Rendering3D",
        "Sofa.GUI.Component",
        "Sofa.Component.Mapping.Linear",
        "Sofa.Component.Mapping.NonLinear"
    ])
    rootNode.gravity = [0.0, -9810, 0.0]
    rootNode.dt = 0.05
    rootNode.addObject('DefaultAnimationLoop')
    rootNode.addObject('DefaultVisualManagerLoop')
    rootNode.addObject('VisualStyle', displayFlags="showVisual showBehavior showForceFields")
    rootNode.addObject('FreeMotionAnimationLoop')
    rootNode.addObject('GenericConstraintSolver',
                       maxIterations=50,
                       tolerance=1e-5)
    simulation = rootNode.addChild('Simulation')
    simulation.addObject('EulerImplicitSolver',
                         name='odesolver',
                         rayleighStiffness=0.1,
                         rayleighMass=0.1)
    simulation.addObject('SparseLDLSolver',
                         name='precond',
                         template="CompressedRowSparseMatrix")
    simulation.addObject('GenericConstraintCorrection')
    modelling = rootNode.addChild('Modelling')
    soft = modelling.addChild(SOFT())
    simulation.addChild(soft)
    return rootNode
[I can provide more code if needed]

Environment:

• SOFA version: v23.06
• Python version: 3.8
• OS: Windows 10

Any insights or suggestions would be greatly appreciated.

Thank you!

[bakpaul]

Hello !

If you are launching your scene with runSofa you should never call init explicitly because you object will be initialized after the call to createScene. And during the createScene execution, the rest of the tree will not have been initialized.

Your version of SOFA is quite old, this issue might have already been fixed in newer version: try upgrading !