SOFA CUDA , collision implementation #5275

anyCli0 · 2025-02-17T13:50:20Z

anyCli0
Feb 17, 2025

Hi, I m trying to use SofaCUDA plugin to spead up my simulation. There s no so much documentation about the usage of this plugin, so probably I poorly understood it. I have understood that some components are not implemented in CUDA: for example Lagrange multiplayer collision is not implemented, as Linear Solver LDL, while for example Penalty method and CG solver are used. I would like to base my simulation on lagrange multiplayer for the collision and on LDLSolver because they are more precised. For these reasons I try to use some components on CPU while others on GPU, writing the following scripts. However the effect obtain is strange, the mechanicalObject mesh of the cylinder deforms really badly. Can you help me? is it possible and does it have sense to implement something like that?
thank you in advance

import Sofa
import Sofa.Simulation
import Sofa.Simulation, Sofa.Gui, Sofa.Core, SofaRuntime
import numpy as np
import os, time
from Sofa.constants import *

class ClosureController(Sofa.Core.Controller):
    def __init__(self, *args, **kwargs):
        super().__init__(*args,**kwargs)
        self.engines_clip1 = kwargs["engines_clip1"]
        self.engines_clip2 = kwargs["engines_clip2"]
        self.init_clip1 = kwargs["init_clip1"]
        self.init_clip2 = kwargs["init_clip2"]
        self.target_clip1 = kwargs["target_clip1"]
        self.target_clip2 = kwargs["target_clip2"]
        self.velocity = kwargs["velocity"]

        self.listening = True
        self.name = "ClosureController"
        print("Python::__init__::" + str(self.name.value))
        self.axe=2
        self.accumulated_time  = 0
        
    def rotate(self,engines,velocity, direction):
        for engine in engines:
            engine.rotation[self.axe] += direction* velocity
        
    def  onAnimateEndEvent(self, event): # begin event, sito
        dt = 0.02  # Il passo temporale che hai impostato
        
        # Aggiungi dt al tempo accumulato
        self.accumulated_time += dt
        print(f" Time: {self.accumulated_time:.2f} seconds")

        self.rotate(self.engines_clip1, self.velocity, -1) 
        self.rotate(self.engines_clip2, self.velocity, 1)   
        print("Rotating clips")

        # Salva la rotazione corrente dopo la rotazione
        self.current_rotation_clip1 = self.init_clip1[self.axe] + self.engines_clip1[0].rotation[self.axe]  
        self.current_rotation_clip2 = self.init_clip2[self.axe] + self.engines_clip2[0].rotation[self.axe] 

        # Ora possiamo fare il confronto con il target
        print(f"Current rotation clip1: {self.current_rotation_clip1}, Target: {self.target_clip1[self.axe]}")
        print(f"Current rotation clip2: {self.current_rotation_clip2}, Target: {self.target_clip2[self.axe]}")

        # Confronta la rotazione corrente con il target
        if abs(self.current_rotation_clip1) > abs(self.target_clip1[self.axe]):
            print("Rotation in progress for clip1")
        else:
            print("Target reached or no movement for clip1")
           
        if abs(self.current_rotation_clip2) > abs(self.target_clip2[self.axe]):
            print("Rotation in progress for clip2")
        else:
            self.velocity=0
            self.root.animate= False 
        print(f"Accumulated time: {self.accumulated_time:.2f} seconds")

def createScene(root):
    root.dt = 0.02
    root.name = 'root'
    root.gravity = [0.0, -9810, 0.0]
    root.addObject('VisualStyle', displayFlags='showVisualModels showForceFields showBehaviorModels')
    root.addObject('FreeMotionAnimationLoop', parallelCollisionDetectionAndFreeMotion=1)
    root.addObject('GenericConstraintSolver', maxIt= 1000, tolerance=1e-4, computeConstraintForces=1)
    root.addObject('CollisionPipeline', name="CollisionPipeline", verbose='0')
    root.addObject('ParallelBruteForceBroadPhase')
    root.addObject('ParallelBVHNarrowPhase')
   
    root.addObject('NewProximityIntersection', alarmDistance = 0.3, contactDistance = 0.1)
    root.addObject('DefaultContactManager', name = "ContactResponse", response = "FrictionContactConstraint" )
    
    # cylinder
    cylinder=root.addChild("Cylinder")
    cylinder.addObject("EulerImplicitSolver", name="cg_odesolver", rayleighStiffness= "0.1", rayleighMass="0.1") 
    cylinder.addObject("CGLinearSolver", name="linear solver", iterations="25",tolerance="1e-09", threshold="1e-09" )
    cylinder.addObject("MeshGmshLoader", name="meshLoader", filename="./cilindri/TETRA/tetra0.5.msh", translation=[0, -8, 0])
    
    cylinder.addObject("MechanicalObject", name="dofs", src="@meshLoader", template="CudaVec3f")
    cylinder.addObject("TetrahedronSetTopologyContainer", name="topo", src="@meshLoader")
    cylinder.addObject("TetrahedronSetGeometryAlgorithms", template="CudaVec3f", name="GeomAlgo")
    cylinder.addObject("DiagonalMass", template="CudaVec3f,CudaVec3f",name="DiagonalMass", totalMass="7e-5" )
    cylinder.addObject("TetrahedronFEMForceField", template="CudaVec3f",name="FEM", method="large", poissonRatio="0.45", youngModulus="216", computeGlobalMatrix="0" )
   
    cylinder.addObject('BoxROI', name='box_ind1', box="-2.5 -10.5 -10 2.5 -5.5 -9.8",  drawBoxes="1")
    cylinder.addObject('BoxROI', name='box_ind2', box="-2.5 -10.5 10 2.5 -5.5 9.8",  drawBoxes="1") 
    cylinder.addObject("FixedProjectiveConstraint", name="FixedProjectiveConstraint", template="CudaVec3f", indices= '@box_ind1.indices')
    cylinder.addObject("FixedProjectiveConstraint", name="FixedProjectiveConstraint", template="CudaVec3f", indices= '@box_ind2.indices')
    cylinder.addObject("UncoupledConstraintCorrection")
    # Collision cylinder
    surf=cylinder.addChild("Surf")
    Loader=surf.addObject("MeshSTLLoader", filename="./cilindri/SHELL/shell0.5.stl", translation=[0, -8, 0])
    surf.addObject("MechanicalObject", name="collision", position="@[-1].position", template="Vec3d" )
    surf.addObject("TriangleSetTopologyContainer", name="topology",src=Loader.getLinkPath())
    surf.addObject("TriangleCollisionModel", name="CollisionModel")
    surf.addObject("LineCollisionModel", name="CollisionModel")
    surf.addObject("PointCollisionModel", name="CollisionModel")
    surf.addObject("BarycentricMapping", name="collision mapping", input="@../dofs", output="@collision")

    # clip
    engines_clip1=[]
    engines_clip2=[]

    clip1=root.addChild("clip1")
    clip1.addObject('MeshObjLoader', name='loader1', filename='./clip/Laz45.750p1newcoll.obj',rotation= [0, 0, +15])
    clip1.addObject('MeshTopology', src='@loader1')
    te1=clip1.addObject('TransformEngine', name ='te1', input_position='@loader1.position')
    engines_clip1.append(te1)
    clip1.addObject('MechanicalObject', name="mo1", position=te1.output_position.getLinkPath(), template="Vec3d")
    clip1.addObject('TriangleCollisionModel', group='1', selfCollision="0")
    clip1.addObject('LineCollisionModel', group='1', selfCollision="0")
    clip1.addObject('PointCollisionModel', group='1', selfCollision="0")
    clip2=root.addChild("clip2")
    clip2.addObject('MeshObjLoader', name='loader2', filename='./clip/Laz45.750p2newcoll.obj', rotation= [0, 0, -15])
    clip2.addObject('MeshTopology', src='@loader2')
    te2=clip2.addObject('TransformEngine', name ='te2', input_position='@loader2.position')
    engines_clip2.append(te2)
    clip2.addObject('MechanicalObject', name="mo2", position=te2.output_position.getLinkPath(), template="Vec3d")
    clip2.addObject('TriangleCollisionModel', group='1', selfCollision="0")
    clip2.addObject('LineCollisionModel', group='1', selfCollision="0")
    clip2.addObject('PointCollisionModel', group='1', selfCollision="0")

    root.addObject(ClosureController(engines_clip1=engines_clip1, engines_clip2=engines_clip2, init_clip1= [0, 0, +15], init_clip2= [0, 0, -15], target_clip1=[0, 0, +3], target_clip2=[0, 0, -3], velocity=0.06))
USE_GUI = True
''
def main():
    import SofaRuntime
    import Sofa.Gui


    # aggiunta plugins
    plugins = ["MultiThreading", "SofaCUDA","SofaGeneralEngine", "Sofa.GL.Component.Rendering3D","Sofa.Component.Collision.Geometry", "Sofa.Component.Topology.Container.Dynamic", "Sofa.Component.ODESolver.Backward", "Sofa.Component.LinearSolver.Iterative", 'Sofa.Component.LinearSolver.Direct',
              "Sofa.Component.IO.Mesh", "Sofa.Component.Collision.Response.Contact", 'Sofa.Component.Collision.Detection.Intersection', "Sofa.Component.SolidMechanics.FEM.HyperElastic", "Sofa.Component.Collision.Detection.Algorithm", "Sofa.Component.AnimationLoop", "Sofa.Component.Constraint.Lagrangian.Solver", "Sofa.Component.Constraint.Lagrangian.Correction"]
    for p in plugins:
        SofaRuntime.importPlugin(p)

    root = Sofa.Core.Node("root") 
    createScene(root)
    Sofa.Simulation.init(root) 

    if not USE_GUI:
        for iteration in range(10):
            Sofa.Simulation.animate(root, root.dt.value)
    else:
        Sofa.Gui.GUIManager.Init("myscene", "qglviewer")
        Sofa.Gui.GUIManager.createGUI(root, __file__)
        Sofa.Gui.GUIManager.SetDimension(1080, 1080)
        Sofa.Gui.GUIManager.MainLoop(root)
        Sofa.Gui.GUIManager.closeGUI()

if __name__ == '__main__':
    main()

zhanglepy · 2025-02-25T08:19:01Z

zhanglepy
Feb 25, 2025

I am also very interested in the Sofa CUDA plugin！

0 replies

alxbilger · 2025-03-04T10:49:02Z

alxbilger
Mar 4, 2025
Maintainer

I recommend that you make a minimal example showing the bad behavior. Is the bad behavior due to collision? Does it work without collision?
It should be possible to have some components on CPU and others on GPU.

Some comments:

I see that you have a CGLinearSolver, but you mentioned you wanted a direct solver. Does it work better?
UncoupledConstraintCorrection is a rough approximation. Try with LinearSolverConstraintCorrection.

1 reply

anyCli0 Mar 6, 2025
Author

Hi, sorry for the late reply. I used CG Linear Solver because I saw that it is implemented in CUDA, the usage of the iterative solver needs UncoupledConstraintCorrection because LinearSolverConstraintCorrection and GenericConstraintCorrection don't accept the default template of the solver (i.e. GraphScattered). the other templates bring to instability. I moved to the Direct solver usage. Here below i report 2 simulations, the first one with LDL Sparse Solver, Lagrange Multiplayer based collision, totally in CPU. The second is the equivalent script but run in GPU. The first simulation goes really slowly but the collision works, the second one brings to instability. I attached screenshot of both scripts and results

--> first simulation (CPU)
import Sofa
import Sofa.Simulation
import Sofa.Simulation, Sofa.Gui, Sofa.Core, SofaRuntime
import numpy as np
import os, time
from Sofa.constants import *

class ClosureController(Sofa.Core.Controller):
def init(self, *args, **kwargs):
super().init(*args,**kwargs)
self.engines_clip1 = kwargs["engines_clip1"]
self.engines_clip2 = kwargs["engines_clip2"]
self.init_clip1 = kwargs["init_clip1"]
self.init_clip2 = kwargs["init_clip2"]
self.target_clip1 = kwargs["target_clip1"]
self.target_clip2 = kwargs["target_clip2"]
self.velocity = kwargs["velocity"]

    self.listening = True
    self.name = "ClosureController"
    print("Python::__init__::" + str(self.name.value))
    self.axe=2
    self.accumulated_time  = 0
    
def rotate(self,engines,velocity, direction):
    for engine in engines:
        engine.rotation[self.axe] += direction* velocity
    
def  onAnimateEndEvent(self, event): # begin event, sito
    dt = 0.02  # Il passo temporale che hai impostato
    
    # Aggiungi dt al tempo accumulato
    self.accumulated_time += dt
    print(f" Time: {self.accumulated_time:.2f} seconds")

    self.rotate(self.engines_clip1, self.velocity, -1) 
    self.rotate(self.engines_clip2, self.velocity, 1)   
    print("Rotating clips")

    # Salva la rotazione corrente dopo la rotazione
    self.current_rotation_clip1 = self.init_clip1[self.axe] + self.engines_clip1[0].rotation[self.axe]  
    self.current_rotation_clip2 = self.init_clip2[self.axe] + self.engines_clip2[0].rotation[self.axe] 

    # Ora possiamo fare il confronto con il target
    print(f"Current rotation clip1: {self.current_rotation_clip1}, Target: {self.target_clip1[self.axe]}")
    print(f"Current rotation clip2: {self.current_rotation_clip2}, Target: {self.target_clip2[self.axe]}")

    # Confronta la rotazione corrente con il target
    if abs(self.current_rotation_clip1) > abs(self.target_clip1[self.axe]):
        print("Rotation in progress for clip1")
    else:
        print("Target reached or no movement for clip1")
       
    if abs(self.current_rotation_clip2) > abs(self.target_clip2[self.axe]):
        print("Rotation in progress for clip2")
    else:
        self.velocity=0
        self.root.animate= False 
    print(f"Accumulated time: {self.accumulated_time:.2f} seconds")

def createScene(root):
root.name = 'root'
root.gravity = [0.0, -9810, 0.0]
root.addObject('VisualStyle', displayFlags='showVisualModels showForceFields showBehaviorModels')
root.addObject('FreeMotionAnimationLoop')
root.addObject('GenericConstraintSolver', maxIt= 1000, tolerance=1e-4)
root.addObject('CollisionPipeline', name="CollisionPipeline", verbose='0')
root.addObject('BruteForceBroadPhase')
root.addObject('BVHNarrowPhase')

root.addObject('NewProximityIntersection', alarmDistance = 0.3, contactDistance = 0.1)
root.addObject('CollisionResponse', name = "ContactResponse", response = "FrictionContactConstraint" )

# cylinder
cylinder=root.addChild("Cylinder")
cylinder.addObject("EulerImplicitSolver", name="cg_odesolver", rayleighStiffness= "0.1", rayleighMass="0.1") 
# cylinder.addObject("CGLinearSolver", name="linear solver",  iterations="25",tolerance="1e-09", threshold="1e-09")
cylinder.addObject('SparseLDLSolver', template="CompressedRowSparseMatrixd")
cylinder.addObject("MeshGmshLoader", name="meshLoader", filename="./cilindri/TETRA/tetra1.0.msh", translation=[0, -8, 0])

cylinder.addObject("MechanicalObject", name="dofs", src="@meshLoader", template="Vec3d")
cylinder.addObject("TetrahedronSetTopologyContainer", name="topo", src="@meshLoader")
cylinder.addObject("TetrahedronSetGeometryAlgorithms", template="Vec3d", name="GeomAlgo")
cylinder.addObject("DiagonalMass", template="Vec3d,Vec3d",name="DiagonalMass", totalMass="7e-5" )
cylinder.addObject("TetrahedronFEMForceField", template="Vec3d",name="FEM", method="large", poissonRatio="0.45", youngModulus="216", computeGlobalMatrix="0" )

cylinder.addObject('BoxROI', name='box_ind1', box="-2.5 -10.5 -10 2.5 -5.5 -9.8",  drawBoxes="1")
cylinder.addObject('BoxROI', name='box_ind2', box="-2.5 -10.5 10 2.5 -5.5 9.8",  drawBoxes="1") 
cylinder.addObject("FixedConstraint", name="FixedProjectiveConstraint", template="Vec3d", indices= '@box_ind1.indices')
cylinder.addObject("FixedConstraint", name="FixedProjectiveConstraint", template="Vec3d", indices= '@box_ind2.indices')
cylinder.addObject("GenericConstraintCorrection")
# Collision cylinder
surf=cylinder.addChild("Surf")
Loader=surf.addObject("MeshSTLLoader", filename="./cilindri/SHELL/shell1.0.stl", translation=[0, -8, 0])
surf.addObject("MechanicalObject", name="collision", position="@[-1].position", template="Vec3d" )
surf.addObject("TriangleSetTopologyContainer", name="topology",src=Loader.getLinkPath())
surf.addObject("TriangleCollisionModel", name="CollisionModel")
surf.addObject("LineCollisionModel", name="CollisionModel")
# surf.addObject("PointCollisionModel", name="CollisionModel")
surf.addObject("BarycentricMapping", name="collision mapping", input="@../dofs", output="@collision")

# # clip
engines_clip1=[]
engines_clip2=[]

clip1=root.addChild("clip1")
clip1.addObject('MeshObjLoader', name='loader1', filename='./clip/Laz45.750p1newcoll.obj',rotation= [0, 0, +15])
clip1.addObject('MeshTopology', src='@loader1')
te1=clip1.addObject('TransformEngine', name ='te1', input_position='@loader1.position')
engines_clip1.append(te1)
clip1.addObject('MechanicalObject', name="mo1", position=te1.output_position.getLinkPath(), template="Vec3d")
clip1.addObject('TriangleCollisionModel', group='1', selfCollision="0")
clip1.addObject('LineCollisionModel', group='1', selfCollision="0")
clip1.addObject('PointCollisionModel', group='1', selfCollision="0")
clip2=root.addChild("clip2")
clip2.addObject('MeshObjLoader', name='loader2', filename='./clip/Laz45.750p2newcoll.obj', rotation= [0, 0, -15])
clip2.addObject('MeshTopology', src='@loader2')
te2=clip2.addObject('TransformEngine', name ='te2', input_position='@loader2.position')
engines_clip2.append(te2)
clip2.addObject('MechanicalObject', name="mo2", position=te2.output_position.getLinkPath(), template="Vec3d")
clip2.addObject('TriangleCollisionModel', group='1', selfCollision="0")
clip2.addObject('LineCollisionModel', group='1', selfCollision="0")
clip2.addObject('PointCollisionModel', group='1', selfCollision="0")

root.addObject(ClosureController(engines_clip1=engines_clip1, engines_clip2=engines_clip2, init_clip1= [0, 0, +15], init_clip2= [0, 0, -15], target_clip1=[0, 0, +3], target_clip2=[0, 0, -3], velocity=0.06))

USE_GUI = True
''
def main():
import SofaRuntime
import Sofa.Gui

# aggiunta plugins
# plugins = ["MultiThreading", "SofaCUDA","SofaGeneralEngine", "Sofa.GL.Component.Rendering3D","Sofa.Component.Collision.Geometry", "Sofa.Component.Topology.Container.Dynamic", "Sofa.Component.ODESolver.Backward", "Sofa.Component.LinearSolver.Iterative", 'Sofa.Component.LinearSolver.Direct',
#           "Sofa.Component.IO.Mesh", "Sofa.Component.Collision.Response.Contact", 'Sofa.Component.Collision.Detection.Intersection', "Sofa.Component.SolidMechanics.FEM.HyperElastic", "Sofa.Component.Collision.Detection.Algorithm", "Sofa.Component.AnimationLoop", "Sofa.Component.Constraint.Lagrangian.Solver", "Sofa.Component.Constraint.Lagrangian.Correction"]
plugins=["Sofa.Component.Collision.Detection.Algorithm", 'Sofa.Component.Constraint.Lagrangian.Correction', "Sofa.Component.Collision.Detection.Intersection", "Sofa.Component.Collision.Geometry", "Sofa.Component.Collision.Response.Contact", "Sofa.Component.IO.Mesh", "Sofa.Component.LinearSolver.Direct","Sofa.Component.LinearSolver.Iterative","Sofa.Component.Mapping.Linear",
         "Sofa.Component.Mass","Sofa.Component.ODESolver.Backward", "Sofa.Component.AnimationLoop","Sofa.Component.MechanicalLoad","Sofa.Component.SolidMechanics.FEM.Elastic", "Sofa.Component.StateContainer", "Sofa.Component.Topology.Container.Constant", "Sofa.GL.Component.Rendering3D", "SofaCUDA", "SofaGeneralEngine" ]
for p in plugins:
    SofaRuntime.importPlugin(p)

root = Sofa.Core.Node("root") 
createScene(root)
Sofa.Simulation.init(root) 

if not USE_GUI:
    for iteration in range(10):
        Sofa.Simulation.animate(root, root.dt.value)
else:
    Sofa.Gui.GUIManager.Init("myscene", "qglviewer")
    Sofa.Gui.GUIManager.createGUI(root, __file__)
    Sofa.Gui.GUIManager.SetDimension(1080, 1080)
    Sofa.Gui.GUIManager.MainLoop(root)
    Sofa.Gui.GUIManager.closeGUI()

if name == 'main':
main()

Good behaviour:

--> second simulation (GPU)
import Sofa
import Sofa.Simulation
import Sofa.Simulation, Sofa.Gui, Sofa.Core, SofaRuntime
import numpy as np
import os, time
from Sofa.constants import *

class ClosureController(Sofa.Core.Controller):
def init(self, *args, **kwargs):
super().init(*args,**kwargs)
self.engines_clip1 = kwargs["engines_clip1"]
self.engines_clip2 = kwargs["engines_clip2"]
self.init_clip1 = kwargs["init_clip1"]
self.init_clip2 = kwargs["init_clip2"]
self.target_clip1 = kwargs["target_clip1"]
self.target_clip2 = kwargs["target_clip2"]
self.velocity = kwargs["velocity"]

    self.listening = True
    self.name = "ClosureController"
    print("Python::__init__::" + str(self.name.value))
    self.axe=2
    self.accumulated_time  = 0
    
def rotate(self,engines,velocity, direction):
    for engine in engines:
        engine.rotation[self.axe] += direction* velocity
    
def  onAnimateEndEvent(self, event): # begin event, sito
    dt = 0.02  # Il passo temporale che hai impostato
    
    # Aggiungi dt al tempo accumulato
    self.accumulated_time += dt
    print(f" Time: {self.accumulated_time:.2f} seconds")

    self.rotate(self.engines_clip1, self.velocity, -1) 
    self.rotate(self.engines_clip2, self.velocity, 1)   
    print("Rotating clips")

    # Salva la rotazione corrente dopo la rotazione
    self.current_rotation_clip1 = self.init_clip1[self.axe] + self.engines_clip1[0].rotation[self.axe]  
    self.current_rotation_clip2 = self.init_clip2[self.axe] + self.engines_clip2[0].rotation[self.axe] 

    # Ora possiamo fare il confronto con il target
    print(f"Current rotation clip1: {self.current_rotation_clip1}, Target: {self.target_clip1[self.axe]}")
    print(f"Current rotation clip2: {self.current_rotation_clip2}, Target: {self.target_clip2[self.axe]}")

    # Confronta la rotazione corrente con il target
    if abs(self.current_rotation_clip1) > abs(self.target_clip1[self.axe]):
        print("Rotation in progress for clip1")
    else:
        print("Target reached or no movement for clip1")
       
    if abs(self.current_rotation_clip2) > abs(self.target_clip2[self.axe]):
        print("Rotation in progress for clip2")
    else:
        self.velocity=0
        self.root.animate= False 
    print(f"Accumulated time: {self.accumulated_time:.2f} seconds")

def createScene(root):
root.name = 'root'
root.gravity = [0.0, -9810, 0.0]
root.addObject('VisualStyle', displayFlags='showVisualModels showForceFields showBehaviorModels')
root.addObject('FreeMotionAnimationLoop')
root.addObject('GenericConstraintSolver', maxIt= 1000, tolerance=1e-4)
root.addObject('CollisionPipeline', name="CollisionPipeline", verbose='0')
root.addObject('BruteForceBroadPhase')
root.addObject('BVHNarrowPhase')

root.addObject('NewProximityIntersection', alarmDistance = 0.3, contactDistance = 0.1)
root.addObject('CollisionResponse', name = "ContactResponse", response = "FrictionContactConstraint" )

# cylinder
cylinder=root.addChild("Cylinder")
cylinder.addObject("EulerImplicitSolver", name="cg_odesolver", rayleighStiffness= "0.1", rayleighMass="0.1") 
# cylinder.addObject("CGLinearSolver", name="linear solver",  iterations="25",tolerance="1e-09", threshold="1e-09")
cylinder.addObject('SparseLDLSolver', template="CompressedRowSparseMatrixd")
cylinder.addObject("MeshGmshLoader", name="meshLoader", filename="./cilindri/TETRA/tetra1.0.msh", translation=[0, -8, 0])

cylinder.addObject("MechanicalObject", name="dofs", src="@meshLoader", template="CudaVec3f")
cylinder.addObject("TetrahedronSetTopologyContainer", name="topo", src="@meshLoader")
cylinder.addObject("TetrahedronSetGeometryAlgorithms", template="CudaVec3f", name="GeomAlgo")
cylinder.addObject("DiagonalMass", template="CudaVec3f,CudaVec3f",name="DiagonalMass", totalMass="7e-5" )
cylinder.addObject("TetrahedronFEMForceField", template="CudaVec3f",name="FEM", method="large", poissonRatio="0.45", youngModulus="216", computeGlobalMatrix="0" )

cylinder.addObject('BoxROI', name='box_ind1', box="-2.5 -10.5 -10 2.5 -5.5 -9.8",  drawBoxes="1")
cylinder.addObject('BoxROI', name='box_ind2', box="-2.5 -10.5 10 2.5 -5.5 9.8",  drawBoxes="1") 
cylinder.addObject("FixedConstraint", name="FixedProjectiveConstraint", template="CudaVec3f", indices= '@box_ind1.indices')
cylinder.addObject("FixedConstraint", name="FixedProjectiveConstraint", template="CudaVec3f", indices= '@box_ind2.indices')
cylinder.addObject("GenericConstraintCorrection")
# Collision cylinder
surf=cylinder.addChild("Surf")
Loader=surf.addObject("MeshSTLLoader", filename="./cilindri/SHELL/shell1.0.stl", translation=[0, -8, 0])
surf.addObject("MechanicalObject", name="collision", position="@[-1].position", template="CudaVec3f" )
surf.addObject("TriangleSetTopologyContainer", name="topology",src=Loader.getLinkPath())
surf.addObject("TriangleCollisionModel", name="CollisionModel")
surf.addObject("LineCollisionModel", name="CollisionModel")
# surf.addObject("PointCollisionModel", name="CollisionModel")
surf.addObject("BarycentricMapping", name="collision mapping", input="@../dofs", output="@collision")

# # clip
engines_clip1=[]
engines_clip2=[]

clip1=root.addChild("clip1")
clip1.addObject('MeshObjLoader', name='loader1', filename='./clip/Laz45.750p1newcoll.obj',rotation= [0, 0, +15])
clip1.addObject('MeshTopology', src='@loader1')
te1=clip1.addObject('TransformEngine', name ='te1', input_position='@loader1.position')
engines_clip1.append(te1)
clip1.addObject('MechanicalObject', name="mo1", position=te1.output_position.getLinkPath(), template="Vec3d")
clip1.addObject('TriangleCollisionModel', group='1', selfCollision="0")
clip1.addObject('LineCollisionModel', group='1', selfCollision="0")
clip1.addObject('PointCollisionModel', group='1', selfCollision="0")
clip2=root.addChild("clip2")
clip2.addObject('MeshObjLoader', name='loader2', filename='./clip/Laz45.750p2newcoll.obj', rotation= [0, 0, -15])
clip2.addObject('MeshTopology', src='@loader2')
te2=clip2.addObject('TransformEngine', name ='te2', input_position='@loader2.position')
engines_clip2.append(te2)
clip2.addObject('MechanicalObject', name="mo2", position=te2.output_position.getLinkPath(), template="Vec3d")
clip2.addObject('TriangleCollisionModel', group='1', selfCollision="0")
clip2.addObject('LineCollisionModel', group='1', selfCollision="0")
clip2.addObject('PointCollisionModel', group='1', selfCollision="0")

root.addObject(ClosureController(engines_clip1=engines_clip1, engines_clip2=engines_clip2, init_clip1= [0, 0, +15], init_clip2= [0, 0, -15], target_clip1=[0, 0, +3], target_clip2=[0, 0, -3], velocity=0.06))

USE_GUI = True
''
def main():
import SofaRuntime
import Sofa.Gui

# aggiunta plugins
# plugins = ["MultiThreading", "SofaCUDA","SofaGeneralEngine", "Sofa.GL.Component.Rendering3D","Sofa.Component.Collision.Geometry", "Sofa.Component.Topology.Container.Dynamic", "Sofa.Component.ODESolver.Backward", "Sofa.Component.LinearSolver.Iterative", 'Sofa.Component.LinearSolver.Direct',
#           "Sofa.Component.IO.Mesh", "Sofa.Component.Collision.Response.Contact", 'Sofa.Component.Collision.Detection.Intersection', "Sofa.Component.SolidMechanics.FEM.HyperElastic", "Sofa.Component.Collision.Detection.Algorithm", "Sofa.Component.AnimationLoop", "Sofa.Component.Constraint.Lagrangian.Solver", "Sofa.Component.Constraint.Lagrangian.Correction"]
plugins=["Sofa.Component.Collision.Detection.Algorithm", 'Sofa.Component.Constraint.Lagrangian.Correction', "Sofa.Component.Collision.Detection.Intersection", "Sofa.Component.Collision.Geometry", "Sofa.Component.Collision.Response.Contact", "Sofa.Component.IO.Mesh", "Sofa.Component.LinearSolver.Direct","Sofa.Component.LinearSolver.Iterative","Sofa.Component.Mapping.Linear",
         "Sofa.Component.Mass","Sofa.Component.ODESolver.Backward", "Sofa.Component.AnimationLoop","Sofa.Component.MechanicalLoad","Sofa.Component.SolidMechanics.FEM.Elastic", "Sofa.Component.StateContainer", "Sofa.Component.Topology.Container.Constant", "Sofa.GL.Component.Rendering3D", "SofaCUDA", "SofaGeneralEngine" ]
for p in plugins:
    SofaRuntime.importPlugin(p)

root = Sofa.Core.Node("root") 
createScene(root)
Sofa.Simulation.init(root) 

if not USE_GUI:
    for iteration in range(10):
        Sofa.Simulation.animate(root, root.dt.value)
else:
    Sofa.Gui.GUIManager.Init("myscene", "qglviewer")
    Sofa.Gui.GUIManager.createGUI(root, __file__)
    Sofa.Gui.GUIManager.SetDimension(1080, 1080)
    Sofa.Gui.GUIManager.MainLoop(root)
    Sofa.Gui.GUIManager.closeGUI()

if name == 'main':
main()

Behaviour obtained with second simulation:

Thank you for the help

(SOFA version used: 23.06)

anyCli0 · 2025-03-06T10:34:31Z

anyCli0
Mar 6, 2025
Author

I don't know if it can help you but looking at the terminal i saw that i don't have errors but I saw this warning:

CudaBarycentricMapping: mapT with 11 entries per DOF and 10 blocs.
[WARNING] [TetrahedronFEMForceField(FEM)] buildDampingMatrix not implemented: for compatibility reason, the deprecated API (addBToMatrix) will be used. This compatibility will disapear in the future, and will cause issues in simulations. Please update the code of TetrahedronFEMForceField to ensure right behavior: the function addBToMatrix has been replaced by buildDampingMatrix
[WARNING] [IntersectorMap] Element Intersector TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TriangleCollisionModel<CudaVectorTypes<Vec<3,float>,Vec<3,float>,float>> NOT FOUND within : LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TLine<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TLine<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TLine<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TPoint<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-SphereCollisionModel<StdRigidTypes<3,double>> TLine<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TSphere<StdRigidTypes<3,double>>
LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-SphereCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TLine<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TSphere<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TPoint<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TPoint<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
SphereCollisionModel<StdRigidTypes<3,double>>-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TSphere<StdRigidTypes<3,double>>-TPoint<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
SphereCollisionModel<StdRigidTypes<3,double>>-SphereCollisionModel<StdRigidTypes<3,double>> TSphere<StdRigidTypes<3,double>>-TSphere<StdRigidTypes<3,double>>
SphereCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TSphere<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TPoint<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
SphereCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-SphereCollisionModel<StdRigidTypes<3,double>> TSphere<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TSphere<StdRigidTypes<3,double>>
SphereCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-SphereCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TSphere<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TSphere<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TTriangle<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TLine<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TTriangle<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TPoint<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-SphereCollisionModel<StdRigidTypes<3,double>> TTriangle<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TSphere<StdRigidTypes<3,double>>
TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-SphereCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TTriangle<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TSphere<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> TTriangle<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TTriangle<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
CubeCollisionModel-CubeCollisionModel Cube-Cube
RayCollisionModel-LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> nullptr
RayCollisionModel-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> nullptr
RayCollisionModel-SphereCollisionModel<StdRigidTypes<3,double>> Ray-TSphere<StdRigidTypes<3,double>>
RayCollisionModel-SphereCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> Ray-TSphere<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
RayCollisionModel-TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> Ray-TTriangle<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>
RayCollisionModel-TetrahedronCollisionModel Ray-Tetrahedron
TetrahedronCollisionModel-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> Tetrahedron-TPoint<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

SOFA CUDA , collision implementation #5275

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Replies: 3 comments 1 reply

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

SOFA CUDA , collision implementation #5275

Uh oh!

Uh oh!

anyCli0 Feb 17, 2025

Replies: 3 comments · 1 reply

Uh oh!

zhanglepy Feb 25, 2025

Uh oh!

alxbilger Mar 4, 2025 Maintainer

Uh oh!

Uh oh!

anyCli0 Mar 6, 2025 Author

Uh oh!

anyCli0 Mar 6, 2025 Author

anyCli0
Feb 17, 2025

Replies: 3 comments 1 reply

zhanglepy
Feb 25, 2025

alxbilger
Mar 4, 2025
Maintainer

anyCli0 Mar 6, 2025
Author

anyCli0
Mar 6, 2025
Author