SOFAPython3_DynamicLiverModel #4862

VeLoXXHH · 2024-07-29T15:50:23Z

VeLoXXHH
Jul 29, 2024

Computational Dynamic Liver Model with SOFA Framework During my internship I have to create an average liver model that can react to external forces in a realistic way. My idea has been divided into several operational steps:

selection of a database of N target livers and 1 template liver
rigid registration of the template for each target followed by a non-rigid registration using the "Thin Plate Spline" method
calculation of the "Principal Component Analysis" to reduce the size of the livers to their eigenvalues and eigenvectors
removal of outliers
with the remaining livers, use SOFA framework in order to calculate a stiffness matrix for each liver
perform an average stiffness matrix that will regulate the behavior of my average liver
I am having problems writing the code on SOFA Framework, in particular:
- build a liver that behaves inside the volume as TetrahedralCorotationalFEMForceField and behaves on the surface as TriangularBendingSprings i.e. as the Glisson capsule.
- be able to apply 3 generic forces: 1 compression 1 traction 1 cut binding the entire posterior part of the liver.

import Sofa
import numpy as np
from scipy import sparse
from scipy import linalg
import numpy as np
from matplotlib import pyplot as plt


USE_GUI = True
exportCSV = True
showImage = False

def main():
    import SofaRuntime
    import Sofa.Gui

    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()


def createScene(root):
    root.gravity=[0, -9.81, 0]
    root.dt=0.02

    root.addObject("RequiredPlugin", pluginName=[
        'Sofa.Component.Collision.Detection.Algorithm',
        'Sofa.Component.Collision.Detection.Intersection',
        'Sofa.Component.Collision.Geometry',
        'Sofa.Component.Collision.Response.Contact',
        'Sofa.Component.Constraint.Projective',
        'Sofa.Component.IO.Mesh',
        'Sofa.Component.LinearSolver.Iterative',
        'Sofa.Component.Mapping.Linear',
        'Sofa.Component.Mass',
        'Sofa.Component.ODESolver.Backward',
        'Sofa.Component.SolidMechanics.FEM.Elastic',
        'Sofa.Component.StateContainer',
        'Sofa.Component.Topology.Container.Dynamic',
        'Sofa.Component.Visual',
        'Sofa.GL.Component.Rendering3D',
        'Sofa.Component.SolidMechanics.Spring',
        'Sofa.Component.MechanicalLoad',
        'Sofa.Component.Topology.Container.Constant'
    ])

    root.addObject('DefaultAnimationLoop')

    root.addObject('VisualStyle', displayFlags="showCollisionModels")
    root.addObject('CollisionPipeline', name="CollisionPipeline")
    root.addObject('BruteForceBroadPhase', name="BroadPhase")
    root.addObject('BVHNarrowPhase', name="NarrowPhase")
    root.addObject('DefaultContactManager', name="CollisionResponse", response="PenalityContactForceField")
    root.addObject('DiscreteIntersection')

    root.addObject("MeshSTLLoader", name="LiverSurface", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_sofa\Liver_Target_2.stl")

    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    liver.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    liver.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_gmsh\Liver_Target_2.msh")  
    liver.addObject('TetrahedronSetTopologyContainer', name="tetra", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="TetraGeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    FEM = liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="TetraFEM", method="large", poissonRatio="0.48", youngModulus="3000", computeGlobalMatrix="0")
    

    liverCapsule = liver.addChild('GlissonCapsule')
    liverCapsule.addObject('MeshSTLLoader', name='CapsuleMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_segmentation\Liver_Target_2.stl")
    liverCapsule.addObject('MechanicalObject', name='dofsCapsule', src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetTopologyContainer', name="tri", src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetGeometryAlgorithms', template="Vec3d", name="TriGeomAlgo")
    liverCapsule.addObject('TriangularBendingSprings', template="Vec3d", name="BendingSprings", stiffness=10000, damping=0.2)
    liverCapsule.addObject('OglModel', name='visualCapsule', color='0.2 0.8 0.2 1.0', src='@CapsuleMesh')
    liverCapsule.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@dofsCapsule", output="@visualCapsule")
    

    visualization = liver.addChild('visualization')
    visualization.addObject('OglModel', name='VisualModelParenchyma', color='0.8 0.2 0.2 1.0', src='@../meshLoader')
    visualization.addObject('BarycentricMapping', name="VisualMappingParenchyma", input="@../dofs", output="@VisualModelParenchyma")
    
    liver.addObject(MatrixAccessController('MatrixAccessor', name='matrixAccessor', force_field=FEM))

    return root

class MatrixAccessController(Sofa.Core.Controller):


    def __init__(self, *args, **kwargs):
        Sofa.Core.Controller.__init__(self, *args, **kwargs)
        self.force_field = kwargs.get("force_field")

    def onAnimateEndEvent(self, event):
        stiffness_matrix = self.force_field.assembleKMatrix()

        print('====================================')
        print('Stiffness matrix')
        print('====================================')
        print('dtype: ' + str(stiffness_matrix.dtype))
        print('shape: ' + str(stiffness_matrix.shape))
        print('ndim: ' + str(stiffness_matrix.ndim))
        print('nnz: ' + str(stiffness_matrix.nnz))
        print('norm: ' + str(sparse.linalg.norm(stiffness_matrix)))

        if exportCSV:
            np.savetxt('stiffness.csv', stiffness_matrix.toarray(), delimiter=',')
        if showImage:
            plt.imshow(stiffness_matrix.toarray(), interpolation='nearest', cmap='gist_gray')
            plt.show(block=False)


if __name__ == '__main__':
    main()

hugtalbot · 2024-07-29T16:10:23Z

hugtalbot
Jul 29, 2024
Maintainer

Hi @VeLoXXHH

Thanks for posting
I seee the issue here in your Capsule description (see mainly the BarycentricMapping):

...

    liverCapsule = liver.addChild('GlissonCapsule')
    liverCapsule.addObject('MeshSTLLoader', name='CapsuleMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_segmentation\Liver_Target_2.stl")
    liverCapsule.addObject('MechanicalObject', name='dofsCapsule', src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetTopologyContainer', name="tri", src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetGeometryAlgorithms', template="Vec3d", name="TriGeomAlgo")
    liverCapsule.addObject('TriangularBendingSprings', template="Vec3d", name="BendingSprings", stiffness=10000, damping=0.2)
    liverCapsule.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofs", output="@dofsCapsule")

    capsuleVisu = liverCapsule.addChild('capsule-visualization')
    capsuleVisu.addObject('OglModel', name='visualCapsule', color='0.2 0.8 0.2 1.0', src='@CapsuleMesh')
    capsuleVisu.addObject('IdentityMapping')
...

0 replies

VeLoXXHH · 2024-07-30T12:58:02Z

VeLoXXHH
Jul 30, 2024
Author

thank you very much for your answer @hugtalbot . now the code works. I have one last question:
now I would like to constrain the entire rear part of the liver (the part not exposed to the surgeon during the operation, so to speak) but I don't know how to do it. I tried to use FixedConstraint but I don't know how to select all the indices of the nodes of the posterior part.
After that I would like to apply 3 different forces:

1 traction force
1 compression force
1 cutting force

I tried to apply the traction and compression forces with ConstantForceField but I don't think it's the best method, or I tried to use SurfacePressureForceField which gives me a more homogeneous result in terms of forces while for the cutting force I have no idea how to implement it. So I wanted to ask you how do I constrain only the rear nodes of the liver, what is the best method to apply a traction/compression force and how can I implement a cutting force.

Follow my current python code:

import Sofa
import numpy as np
from scipy import sparse
from scipy import linalg
import numpy as np
from matplotlib import pyplot as plt


USE_GUI = True
exportCSV = True
showImage = False

def main():
    import SofaRuntime
    import Sofa.Gui

    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()


def createScene(root):
    root.gravity=[0, 0, 0]
    root.dt=0.02

    root.addObject("RequiredPlugin", pluginName=[
        'Sofa.Component.Collision.Detection.Algorithm',
        'Sofa.Component.Collision.Detection.Intersection',
        'Sofa.Component.Collision.Geometry',
        'Sofa.Component.Collision.Response.Contact',
        'Sofa.Component.Constraint.Projective',
        'Sofa.Component.IO.Mesh',
        'Sofa.Component.LinearSolver.Iterative',
        'Sofa.Component.Mapping.Linear',
        'Sofa.Component.Mass',
        'Sofa.Component.ODESolver.Backward',
        'Sofa.Component.SolidMechanics.FEM.Elastic',
        'Sofa.Component.StateContainer',
        'Sofa.Component.Topology.Container.Dynamic',
        'Sofa.Component.Visual',
        'Sofa.GL.Component.Rendering3D',
        'Sofa.Component.SolidMechanics.Spring',
        'Sofa.Component.MechanicalLoad',
        'Sofa.Component.Topology.Container.Constant',
        'Sofa.Component.SolidMechanics.FEM'
    ])

    root.addObject('DefaultAnimationLoop')

    root.addObject('VisualStyle', displayFlags="showCollisionModels")
    root.addObject('CollisionPipeline', name="CollisionPipeline")
    root.addObject('BruteForceBroadPhase', name="BroadPhase")
    root.addObject('BVHNarrowPhase', name="NarrowPhase")
    root.addObject('DefaultContactManager', name="CollisionResponse", response="PenalityContactForceField")
    root.addObject('DiscreteIntersection')

    root.addObject("MeshSTLLoader", name="LiverSurface", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_sofa\Liver_Target_2.stl")

    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    liver.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    liver.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_gmsh\Liver_Target_2.msh")  
    liver.addObject('TetrahedronSetTopologyContainer', name="tetra", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="TetraGeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    FEM = liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="TetraFEM", method="large", poissonRatio="0.48", youngModulus="3000", computeGlobalMatrix="0")
    liver.addObject('FixedConstraint', name="FixedConstraint", indices="3 39 64")
    
    liverCapsule = liver.addChild('GlissonCapsule')
    liverCapsule.addObject('MeshSTLLoader', name='CapsuleMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_sofa\Liver_Target_2.stl")
    liverCapsule.addObject('MechanicalObject', name='dofsCapsule', src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetTopologyContainer', name="tri", src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetGeometryAlgorithms', template="Vec3d", name="TriGeomAlgo")
    liverCapsule.addObject('TriangularBendingSprings', template="Vec3d", name="BendingSprings", stiffness=10000, damping=0.2)
    liverCapsule.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofs", output="@dofsCapsule")
    liverCapsule.addObject('SurfacePressureForceField', name="tractionForce", template="Vec3d", pressure="-1000", triangleIndices = "@CapsuleMesh.triangles")

    visualization = liver.addChild('visualization')
    visualization.addObject('OglModel', name='VisualModelParenchyma', color='0.8 0.2 0.2 1.0', src='@../meshLoader')
    visualization.addObject('BarycentricMapping', name="VisualMappingParenchyma", input="@../dofs", output="@VisualModelParenchyma")
    capsuleVisu = liverCapsule.addChild('capsule-visualization')
    capsuleVisu.addObject('OglModel', name='visualCapsule', color='0.2 0.8 0.2 1.0', src='@../CapsuleMesh')
    capsuleVisu.addObject('IdentityMapping')
    
    liver.addObject(MatrixAccessController('MatrixAccessor', name='matrixAccessor', force_field=FEM))

    return root

class MatrixAccessController(Sofa.Core.Controller):


    def __init__(self, *args, **kwargs):
        Sofa.Core.Controller.__init__(self, *args, **kwargs)
        self.force_field = kwargs.get("force_field")

    def onAnimateEndEvent(self, event):
        stiffness_matrix = self.force_field.assembleKMatrix()

        print('====================================')
        print('Stiffness matrix')
        print('====================================')
        print('dtype: ' + str(stiffness_matrix.dtype))
        print('shape: ' + str(stiffness_matrix.shape))
        print('ndim: ' + str(stiffness_matrix.ndim))
        print('nnz: ' + str(stiffness_matrix.nnz))
        print('norm: ' + str(sparse.linalg.norm(stiffness_matrix)))

        if exportCSV:
            np.savetxt('stiffness.csv', stiffness_matrix.toarray(), delimiter=',')
        if showImage:
            plt.imshow(stiffness_matrix.toarray(), interpolation='nearest', cmap='gist_gray')
            plt.show(block=False)

if __name__ == '__main__':
    main()

1 reply

hugtalbot Jul 30, 2024
Maintainer

Hi @VeLoXXHH

To define a fixed boundary condition, you can indeed use the FixedProjectiveConstraint, see the video tutorial here.

Then, about applying forces:

1 traction or compression force can be done using either:
- the ConstantForceField applying a nodal constant force
- the SurfacePressureForceField applying a pressure over a surface (triangles)
cutting is something very specific, please see the tearing plugin developed by the consulting company InfinyTech3D

VeLoXXHH · 2024-08-01T11:47:38Z

VeLoXXHH
Aug 1, 2024
Author

thank you very much for your reply @hugtalbot , it was very helpful.
I have another question please.
I tried to apply 1 constant force and 1 pressure but starting the program and running the animation it is as if no force is applied and I can't understand why. I also tried to make a collision between the liver and a floor (which I would like to impose to infinite stiffness) and, during the animation, the liver crosses the floor as if it did not exist. My questions were therefore:

if I did something wrong in the implementation of the constant force and the implementation
how can I impose an infinite stiffness to the floor or in any case make it so that it does not model itself in contact with the liver?
how come the liver crosses the floor avoiding the collision?
thanks in advance for the answer, below are my codes:
CODE ABOUT CONSTANT FORCE FIELD:

import Sofa
import numpy as np
from scipy import sparse
from scipy import linalg
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

USE_GUI = True
exportCSV = True
showImage = False

def main():
    import SofaRuntime
    import Sofa.Gui

    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()


def createScene(root):
    root.gravity=[0, 0, 0]
    root.dt=0.02

    root.addObject("RequiredPlugin", pluginName=[
        'Sofa.Component.Collision.Detection.Algorithm',
        'Sofa.Component.Collision.Detection.Intersection',
        'Sofa.Component.Collision.Geometry',
        'Sofa.Component.Collision.Response.Contact',
        'Sofa.Component.Constraint.Projective',
        'Sofa.Component.IO.Mesh',
        'Sofa.Component.LinearSolver.Iterative',
        'Sofa.Component.Mapping.Linear',
        'Sofa.Component.Mass',
        'Sofa.Component.ODESolver.Backward',
        'Sofa.Component.SolidMechanics.FEM.Elastic',
        'Sofa.Component.StateContainer',
        'Sofa.Component.Topology.Container.Dynamic',
        'Sofa.Component.Visual',
        'Sofa.GL.Component.Rendering3D',
        'Sofa.Component.SolidMechanics.Spring',
        'Sofa.Component.MechanicalLoad',
        'Sofa.Component.Topology.Container.Constant',
        'Sofa.Component.SolidMechanics.FEM'
    ])

    root.addObject('DefaultAnimationLoop')

    root.addObject('VisualStyle', displayFlags="showCollisionModels")
    root.addObject('CollisionPipeline', name="CollisionPipeline")
    root.addObject('BruteForceBroadPhase', name="BroadPhase")
    root.addObject('BVHNarrowPhase', name="NarrowPhase")
    root.addObject('DefaultContactManager', name="CollisionResponse", response="PenalityContactForceField")
    root.addObject('DiscreteIntersection')

    root.addObject("MeshSTLLoader", name="LiverSurface", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.stl")
    root.addObject("MeshOBJLoader", name="ConstraintIDpoints", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\ConstraintIDPoints.obj") 
    root.addObject("MeshOBJLoader", name="ForceIDpoints", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\ForceIDPoints.obj") 

    #Liver Parenchyma with TetrahedralCorotationalFEMForceField
    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    liver.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    liver.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.msh")  
    liver.addObject('TetrahedronSetTopologyContainer', name="tetra", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="TetraGeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="TetraFEM", method="large", poissonRatio="0.48", youngModulus="3000")
    # #constraints
    liver.addObject('FixedConstraint', name="FixedConstraint", indices="@../ConstraintIDpoints.position")
    #forces
    liver.addObject('ConstantForceField', force = "0 0 20", indices="@../ForceIDpoints.position" )

    # Glisson capsule with TriangularBendingSprings
    liverCapsule = liver.addChild('GlissonCapsule')
    liverCapsule.addObject('MeshSTLLoader', name='CapsuleMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.stl")
    liverCapsule.addObject('MechanicalObject', name='dofsCapsule', src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetTopologyContainer', name="tri", src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetGeometryAlgorithms', template="Vec3d", name="TriGeomAlgo")
    liverCapsule.addObject('TriangularBendingSprings', template="Vec3d", name="BendingSprings", stiffness=10000, damping=0.2)
    liverCapsule.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofs", output="@dofsCapsule")
    
    #visualization Parenchyma
    liverVisu = liver.addChild('visualization')
    liverVisu.addObject('OglModel', name='VisualModelParenchyma', color='0.8 0.2 0.2 1.0', src='@../meshLoader')
    liverVisu.addObject('BarycentricMapping', name="VisualMappingParenchyma", input="@../dofs", output="@VisualModelParenchyma")

    #visualization Glisson's Capsule
    capsuleVisu = liverCapsule.addChild('capsule-visualization')
    capsuleVisu.addObject('OglModel', name='visualCapsule', color='0.2 0.8 0.2 1.0', src='@../CapsuleMesh')
    capsuleVisu.addObject('IdentityMapping')

    return root

if __name__ == '__main__':
    main()

CODE ABOUT PRESSURE:

import Sofa
import numpy as np
from scipy import sparse
from scipy import linalg
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

USE_GUI = True
exportCSV = True
showImage = False

def main():
    import SofaRuntime
    import Sofa.Gui

    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()


def createScene(root):
    root.gravity=[0, 0, 0]
    root.dt=0.02

    root.addObject("RequiredPlugin", pluginName=[
        'Sofa.Component.Collision.Detection.Algorithm',
        'Sofa.Component.Collision.Detection.Intersection',
        'Sofa.Component.Collision.Geometry',
        'Sofa.Component.Collision.Response.Contact',
        'Sofa.Component.Constraint.Projective',
        'Sofa.Component.IO.Mesh',
        'Sofa.Component.LinearSolver.Iterative',
        'Sofa.Component.Mapping.Linear',
        'Sofa.Component.Mass',
        'Sofa.Component.ODESolver.Backward',
        'Sofa.Component.SolidMechanics.FEM.Elastic',
        'Sofa.Component.StateContainer',
        'Sofa.Component.Topology.Container.Dynamic',
        'Sofa.Component.Visual',
        'Sofa.GL.Component.Rendering3D',
        'Sofa.Component.SolidMechanics.Spring',
        'Sofa.Component.MechanicalLoad',
        'Sofa.Component.Topology.Container.Constant',
        'Sofa.Component.SolidMechanics.FEM'
    ])

    root.addObject('DefaultAnimationLoop')

    root.addObject('VisualStyle', displayFlags="showCollisionModels")
    root.addObject('CollisionPipeline', name="CollisionPipeline")
    root.addObject('BruteForceBroadPhase', name="BroadPhase")
    root.addObject('BVHNarrowPhase', name="NarrowPhase")
    root.addObject('DefaultContactManager', name="CollisionResponse", response="PenalityContactForceField")
    root.addObject('DiscreteIntersection')

    root.addObject("MeshSTLLoader", name="LiverSurface", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.stl")
    root.addObject("MeshOBJLoader", name="ConstraintIDpoints", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Constraint2IDPoints.obj") 
    root.addObject("MeshOBJLoader", name="ForceIDpoints", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\ForceIDPoints.obj") 

    #Liver Parenchyma with TetrahedralCorotationalFEMForceField
    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    liver.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    liver.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.msh")  
    liver.addObject('TetrahedronSetTopologyContainer', name="tetra", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="TetraGeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="TetraFEM", method="large", poissonRatio="0.48", youngModulus="3000")
    # #constraints
    liver.addObject('FixedConstraint', name="FixedConstraint", indices="@../Constraint2IDpoints.position")

    # Glisson capsule with TriangularBendingSprings
    liverCapsule = liver.addChild('GlissonCapsule')
    liverCapsule.addObject('MeshSTLLoader', name='CapsuleMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.stl")
    liverCapsule.addObject('MechanicalObject', name='dofsCapsule', src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetTopologyContainer', name="tri", src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetGeometryAlgorithms', template="Vec3d", name="TriGeomAlgo")
    liverCapsule.addObject('TriangularBendingSprings', template="Vec3d", name="BendingSprings", stiffness=10000, damping=0.2)
    liverCapsule.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofs", output="@dofsCapsule")
    #pressure
    #liverCapsule.addObject('SurfacePressureForceField', name="compressionForce", template="Vec3d", pressure="-100", triangleIndices = "@CapsuleMesh.triangles")
    #liverCapsule.addObject('SurfacePressureForceField', name="tractionForce", template="Vec3d", pressure="1000", triangleIndices = "@CapsuleMesh.triangles")
    
    #visualization Parenchyma
    liverVisu = liver.addChild('visualization')
    liverVisu.addObject('OglModel', name='VisualModelParenchyma', color='0.8 0.2 0.2 1.0', src='@../meshLoader')
    liverVisu.addObject('BarycentricMapping', name="VisualMappingParenchyma", input="@../dofs", output="@VisualModelParenchyma")

    #visualization Glisson's Capsule
    capsuleVisu = liverCapsule.addChild('capsule-visualization')
    capsuleVisu.addObject('OglModel', name='visualCapsule', color='0.2 0.8 0.2 1.0', src='@../CapsuleMesh')
    capsuleVisu.addObject('IdentityMapping')

    return root

if __name__ == '__main__':
    main()

CODE ABOUT COLLISION:

import Sofa
import numpy as np
from scipy import sparse
from scipy import linalg
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

USE_GUI = True
exportCSV = True
showImage = False

def main():
    import SofaRuntime
    import Sofa.Gui

    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()


def createScene(root):
    root.gravity=[0, 0, -9.81]
    root.dt=0.02

    root.addObject("RequiredPlugin", pluginName=[
        'Sofa.Component.Collision.Detection.Algorithm',
        'Sofa.Component.Collision.Detection.Intersection',
        'Sofa.Component.Collision.Geometry',
        'Sofa.Component.Collision.Response.Contact',
        'Sofa.Component.Constraint.Projective',
        'Sofa.Component.IO.Mesh',
        'Sofa.Component.LinearSolver.Iterative',
        'Sofa.Component.Mapping.Linear',
        'Sofa.Component.Mass',
        'Sofa.Component.ODESolver.Backward',
        'Sofa.Component.SolidMechanics.FEM.Elastic',
        'Sofa.Component.StateContainer',
        'Sofa.Component.Topology.Container.Dynamic',
        'Sofa.Component.Visual',
        'Sofa.GL.Component.Rendering3D',
        'Sofa.Component.SolidMechanics.Spring',
        'Sofa.Component.MechanicalLoad',
        'Sofa.Component.Topology.Container.Constant',
        'Sofa.Component.SolidMechanics.FEM'
    ])

    root.addObject('DefaultAnimationLoop')

    root.addObject('VisualStyle', displayFlags="showCollisionModels")
    root.addObject('CollisionPipeline', name="CollisionPipeline")
    root.addObject('BruteForceBroadPhase', name="BroadPhase")
    root.addObject('BVHNarrowPhase', name="NarrowPhase")
    root.addObject('DefaultContactManager', name="CollisionResponse", response="PenalityContactForceField")
    root.addObject('DiscreteIntersection')

    root.addObject("MeshSTLLoader", name="LiverSurface", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.stl")
    root.addObject("MeshSTLLoader", name="Floor", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\Platform\plane.stl")

    #Liver Parenchyma with TetrahedralCorotationalFEMForceField
    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    liver.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    liver.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.msh", translation = "300 300 800")  
    liver.addObject('TetrahedronSetTopologyContainer', name="tetra", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="TetraGeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="TetraFEM", method="large", poissonRatio="0.48", youngModulus="3000")
    liver.addObject("TriangleCollisionModel", selfCollision="0")
    liver.addObject("LineCollisionModel", selfCollision="0")
    liver.addObject("PointCollisionModel", selfCollision="0")

    # Glisson capsule with TriangularBendingSprings
    liverCapsule = liver.addChild('GlissonCapsule')
    liverCapsule.addObject('MeshSTLLoader', name='CapsuleMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.stl", translation = "300 300 800")
    liverCapsule.addObject('MechanicalObject', name='dofsCapsule', src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetTopologyContainer', name="tri", src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetGeometryAlgorithms', template="Vec3d", name="TriGeomAlgo")
    liverCapsule.addObject('TriangularBendingSprings', template="Vec3d", name="BendingSprings", stiffness=10000, damping=0.2)
    liverCapsule.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofs", output="@dofsCapsule")
    
    #visualization Parenchyma
    liverVisu = liver.addChild('visualization')
    liverVisu.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.msh", translation = "@../meshLoader.transition")
    liverVisu.addObject('OglModel', name='VisualModelParenchyma', color='0.8 0.2 0.2 1.0', src='@../meshLoader')
    liverVisu.addObject('BarycentricMapping', name="VisualMappingParenchyma", input="@../dofs", output="@VisualModelParenchyma")

    #visualization Glisson's Capsule
    capsuleVisu = liverCapsule.addChild('capsule-visualization')
    capsuleVisu.addObject('OglModel', name='visualCapsule', color='0.2 0.8 0.2 1.0', src='@../CapsuleMesh')
    capsuleVisu.addObject('IdentityMapping')

    #floor 
    impactFloor = root.addChild('impactFloor')
    impactFloor.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    impactFloor.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    impactFloor.addObject('MeshGmshLoader', name='floorMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\Platform\plane.msh", rotation = "-90 0 0")
    impactFloor.addObject('TetrahedronSetTopologyContainer', name="tetraFloor", src="@floorMesh")
    impactFloor.addObject('MechanicalObject', name='dofsFloor', src="@floorMesh")
    impactFloor.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="TetraGeomAlgoFloor")
    impactFloor.addObject('UniformMass', name="Mass", totalMass="1.0")
    impactFloor.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="TetraFEM", method="large", poissonRatio="0.48", youngModulus="3000")
    impactFloor.addObject('FixedConstraint', name="FixedConstraint", fixAll = True)
    impactFloor.addObject("TriangleCollisionModel", simulated="false", moving="false")
    impactFloor.addObject("LineCollisionModel", simulated="false", moving="false")
    impactFloor.addObject("PointCollisionModel", simulated="false", moving="false")
    
    #visualization floor
    floorVisu = impactFloor.addChild('floor-visualization')
    floorVisu.addObject('MeshGmshLoader', name='floorMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\Platform\plane.msh", rotation = "@../floorMesh.rotation")
    floorVisu.addObject('OglModel', name='visualFloor', color='0.2 0.2 0.8 1.0', src='@../floorMesh')
    floorVisu.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofsFloor", output="@visualFloor")

    return root

if __name__ == '__main__':
    main()
    ```

0 replies

VeLoXXHH · 2024-08-12T15:13:03Z

VeLoXXHH
Aug 12, 2024
Author

Here the warning about the collision between liver and floor:

[WARNING] [IntersectorMap] Element Intersector PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>> NOT FOUND within :  LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr     
  LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
  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>>-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>>  nullptr
  TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
ouble>> NOT FOUND within :  LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr     
  LineCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
  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>>-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>>  nullptr
  TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
  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>>  nullptr
  TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
e>,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>>  nullptr
  TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
ble>,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>>  nullptr
  TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
  TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>-TriangleCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
  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-PointCollisionModel<StdVectorTypes<Vec<3,double>,Vec<3,double>,double>>  nullptr
  RayCollisionModel-SphereCollisionModel<StdRigidTypes<3,double>>  Ray-TSphere<StdRigidTypes<3,double>>
  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>>

2 replies

hugtalbot Aug 12, 2024
Maintainer

Hi @VeLoXXHH

Could you please try to replace DiscreteIntersection with LocalMinDistance and let me know if this solves your problem

VeLoXXHH Aug 12, 2024
Author

the warning disappears, but the collision does not occur yet.

VeLoXXHH · 2024-08-12T17:47:24Z

VeLoXXHH
Aug 12, 2024
Author

the warning disappears, but the collision does not occur yet. Il giorno lun 12 ago 2024 alle ore 17:17 Hugo ***@***.***> ha scritto:

…

Hi @VeLoXXHH <https://github.com/VeLoXXHH> Could you please try to replace DiscreteIntersection with LocalMinDistance and let me know if this solves your problem — Reply to this email directly, view it on GitHub <#4862 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BHL6D5M4O4IULWTXCLDOJ2TZRDG2NAVCNFSM6AAAAABLUTBGTWVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAMZRGU4DSMA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

1 reply

hugtalbot Aug 13, 2024
Maintainer

can you check whether collision really does not occur by activating the Collision > DetectionOuputs (or Behavior > Interactions) in the View panel of the runSofa GUI?
This should allow visualize any potential collision detected.

Let me know about what you see

VeLoXXHH · 2024-08-13T13:25:48Z

VeLoXXHH
Aug 13, 2024
Author

![Liver_Test_impact_](https://github.com/user-attachments/assets/071687dc-9b9e-4378-9746-485b4b524318) I saw that using the "interaction" option inside "view" in Sofa GUI in the initial screen with the snake you can see the contact iterations highlighted in red; unfortunately in my script screen there is no interaction detected between liver and floor (I share the screenshot). I also tried to change the "DefaultContactManager" with response="FrictionContactConstraint" and the program crashed when the two objects came into contact and I got the following message in the terminal: ########## SIG 11 - SIGSEGV: segfault ########## sofa::helper::BackTrace::dump sofa::helper::BackTrace::sig log2f log2f _C_specific_handler _chkstk RtlFindCharInUnicodeString KiUserExceptionDispatcher sofa::component::mapping::linear::BarycentricMapperMeshTopology<sofa::defaulttype::StdVectorTypes<sofa::type::Vec<3,double>,sofa::type::Vec<3,double>,double>,sofa::defaulttype::StdVectorTypes<sofa::type::Vec<3,double>,sofa::type::Vec<3,double>,double>

::

sofa::component::mapping::linear::BarycentricMapping<sofa::defaulttype::StdVectorTypes<sofa::type::Vec<3,double>,sofa::type::Vec<3,double>,double>,sofa::defaulttype::StdVectorTypes<sofa::type::Vec<3,double>,sofa::type::Vec<3,double>,double>

::apply

sofa::core::Mapping<sofa::defaulttype::StdVectorTypes<sofa::type::Vec<3,double>,sofa::type::Vec<3,double>,double>,sofa::defaulttype::StdRigidTypes<3,double>

::apply

sofa::component::collision::response::mapper::BarycentricContactMapper<sofa::component::collision::geometry::LineCollisionModel<sofa::defaulttype::StdVectorTypes<sofa::type::Vec<3,double>,sofa::type::Vec<3,double>,double>

,sofa::defaulttype::StdVectorTyp

sofa::component::collision::response::contact::FrictionContact<sofa::component::collision::geometry::TriangleCollisionModel<sofa::defaulttype::StdVectorTypes<sofa::type::Vec<3,double>,sofa::type::Vec<3,double>,double>

,sofa::component::collision::geometr

sofa::component::collision::response::contact::FrictionContact<sofa::component::collision::geometry::LineCollisionModel<sofa::defaulttype::StdVectorTypes<sofa::type::Vec<3,double>,sofa::type::Vec<3,double>,double>

,sofa::component::collision::geometry::S

sofa::component::collision::detection::algorithm::CollisionPipeline::doCollisionResponse sofa::simulation::PipelineImpl::computeCollisionResponse sofa::simulation::CollisionVisitor::processCollisionPipeline sofa::simulation::AnimateEndEvent::getEventTypeIndex sofa::simulation::CollisionAnimationLoop::preCollisionComputation sofa::simulation::CollisionVisitor::processNodeTopDown sofa::simulation::graph::DAGNode::executeVisitorTopDown sofa::simulation::graph::DAGNode::doExecuteVisitor sofa::simulation::Node::executeVisitor sofa::simulation::AnimateVisitor::processCollisionPipeline sofa::simulation::AnimateVisitor::processNodeTopDown sofa::simulation::graph::DAGNode::executeVisitorTopDown sofa::simulation::graph::DAGNode::doExecuteVisitor sofa::simulation::Node::executeVisitor sofa::simulation::DefaultAnimationLoop::step sofa::simulation::Simulation::animate sofa::gui::qt::RealGUI::step QMetaObject::activate QTimer::timerEvent QObject::event QApplicationPrivate::notify_helper QApplication::notify QCoreApplication::notifyInternal2 QEventDispatcherWin32::event QApplicationPrivate::notify_helper QApplication::notify QCoreApplication::notifyInternal2 QCoreApplicationPrivate::sendPostedEvents qt_plugin_query_metadata QEventDispatcherWin32::processEvents DispatchMessageW DispatchMessageW QEventDispatcherWin32::processEvents qt_plugin_query_metadata QEventLoop::exec QCoreApplication::exec sofa::gui::qt::RealGUI::mainLoop sofa::gui::common::GUIManager::MainLoop PyInit_Gui PyInit_Gui PyCFunction_Call PyEval_EvalFrameDefault PyEval_EvalFrameDefault PyEval_EvalCodeWithName PyEval_EvalCodeEx PyEval_EvalCode PyArena_New PyArena_New Py_wfopen PyUnicode_CompareWithASCIIString PyRun_SimpleFileExFlags PyRun_AnyFileExFlags Py_gitversion Py_RunMain Py_RunMain Py_Main Py_Main BaseThreadInitThunk RtlUserThreadStart Il giorno mar 13 ago 2024 alle ore 09:16 Hugo ***@***.***> ha scritto:

…

can you check whether collision really does not occur by activating the Collision > DetectionOuputs (or Behavior > Interactions) in the View panel of the runSofa GUI? This should allow visualize any potential collision detected. Screenshot_2024-08-13_09-14-26.png (view on web) <https://github.com/user-attachments/assets/a19ea8e5-fe39-486c-93d8-efee89191643> Let me know about what you see — Reply to this email directly, view it on GitHub <#4862 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BHL6D5KJAVXFKSGIGGT7UVTZRGXEJAVCNFSM6AAAAABLUTBGTWVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAMZSGIZDKNY> . You are receiving this because you were mentioned.Message ID: ***@***.***>

0 replies

hugtalbot · 2024-08-13T19:44:41Z

hugtalbot
Aug 13, 2024
Maintainer

Hi @VeLoXXHH

I see from you scene that you are not using a dedicated mesh for the collision detection. We usually prefer to use a surface mesh, representing the envelope of the object instead of the full 3D mesh (containing tetrahedra, and where the triangle orientations is uncertain). Could you try:

checking the normals (again, in the GUI the option Normals)
or using a dedicated surface mesh for the collision detection

3 replies

VeLoXXHH Aug 14, 2024
Author

I tried both to replace the volumetric mesh of the floor with a surface mesh like this:

import Sofa
import numpy as np
from scipy import sparse
from scipy import linalg
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

USE_GUI = True
exportCSV = True
showImage = False

def main():
    import SofaRuntime
    import Sofa.Gui

    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()


def createScene(root):
    root.gravity=[0, 0, -9.81]
    root.dt=0.02

    root.addObject("RequiredPlugin", pluginName=[
        'Sofa.Component.Collision.Detection.Algorithm',
        'Sofa.Component.Collision.Detection.Intersection',
        'Sofa.Component.Collision.Geometry',
        'Sofa.Component.Collision.Response.Contact',
        'Sofa.Component.Constraint.Projective',
        'Sofa.Component.IO.Mesh',
        'Sofa.Component.LinearSolver.Iterative',
        'Sofa.Component.Mapping.Linear',
        'Sofa.Component.Mass',
        'Sofa.Component.ODESolver.Backward',
        'Sofa.Component.SolidMechanics.FEM.Elastic',
        'Sofa.Component.StateContainer',
        'Sofa.Component.Topology.Container.Dynamic',
        'Sofa.Component.Visual',
        'Sofa.GL.Component.Rendering3D',
        'Sofa.Component.SolidMechanics.Spring',
        'Sofa.Component.MechanicalLoad',
        'Sofa.Component.Topology.Container.Constant',
        'Sofa.Component.SolidMechanics.FEM'
    ])

    root.addObject('DefaultAnimationLoop')

    root.addObject('VisualStyle', displayFlags="showCollisionModels")
    root.addObject('CollisionPipeline', name="CollisionPipeline")
    root.addObject('BruteForceBroadPhase', name="BroadPhase")
    root.addObject('BVHNarrowPhase', name="NarrowPhase")
    root.addObject('DefaultContactManager', name="CollisionResponse", response="FrictionContactConstraint")
    root.addObject('LocalMinDistance', name="Proximity", alarmDistance="1.5", contactDistance="1")

    root.addObject("MeshSTLLoader", name="LiverSurface", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_segmentation\Liver_Target_1.stl")
    root.addObject("MeshSTLLoader", name="Floor", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\Platform\plane.stl")

    #Liver Parenchyma with TetrahedralCorotationalFEMForceField
    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    liver.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    liver.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.msh", translation = "0 0 800")  
    liver.addObject('TetrahedronSetTopologyContainer', name="tetra", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="TetraGeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="TetraFEM", method="large", poissonRatio="0.48", youngModulus="3000")
    liver.addObject("TriangleCollisionModel", selfCollision="0")
    liver.addObject("LineCollisionModel", selfCollision="0")
    liver.addObject("PointCollisionModel", selfCollision="0")

    # Glisson capsule with TriangularBendingSprings
    liverCapsule = liver.addChild('GlissonCapsule')
    liverCapsule.addObject('MeshSTLLoader', name='CapsuleMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.stl", translation = "0 0 800")
    liverCapsule.addObject('MechanicalObject', name='dofsCapsule', src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetTopologyContainer', name="tri", src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetGeometryAlgorithms', template="Vec3d", name="TriGeomAlgo")
    liverCapsule.addObject('TriangularBendingSprings', template="Vec3d", name="BendingSprings", stiffness=10000, damping=0.2)
    liverCapsule.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofs", output="@dofsCapsule")

    #visualization Parenchyma
    liverVisu = liver.addChild('visualization')
    liverVisu.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.msh", translation = "@../meshLoader.transition")
    liverVisu.addObject('OglModel', name='VisualModelParenchyma', color='0.8 0.2 0.2 1.0', src='@../meshLoader')
    liverVisu.addObject('BarycentricMapping', name="VisualMappingParenchyma", input="@../dofs", output="@VisualModelParenchyma")

    #visualization Glisson's Capsule
    capsuleVisu = liverCapsule.addChild('capsule-visualization')
    capsuleVisu.addObject('OglModel', name='visualCapsule', color='0.2 0.8 0.2 1.0', src='@../CapsuleMesh')
    capsuleVisu.addObject('IdentityMapping')

    #floor 
    impactFloor = root.addChild('impactFloor')
    impactFloor.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    impactFloor.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    impactFloor.addObject('MeshOBJLoader', name='floorMesh', filename="mesh/floor3.obj", rotation = "-90 0 0")
    impactFloor.addObject("MeshTopology", src="@floorMesh")
    impactFloor.addObject('MechanicalObject', name='dofsFloor', src="@floorMesh", scale="8")
    impactFloor.addObject('FixedConstraint', name="FixedConstraint", fixAll = True)
    impactFloor.addObject("TriangleCollisionModel",name="FloorTriangle", simulated="false", moving="false")
    impactFloor.addObject("LineCollisionModel",name="FloorLine", simulated="false", moving="false")
    impactFloor.addObject("PointCollisionModel",name="FloorPoint", simulated="false", moving="false")
    
    #visualization floor
    floorVisu = impactFloor.addChild('floor-visualization')
    floorVisu.addObject('MeshOBJLoader', name='floorMesh', filename="mesh/floor3.obj", rotation = "@../floorMesh.rotation")
    floorVisu.addObject('OglModel', name='visualFloor', color='0.2 0.2 0.8 1.0', src='@../floorMesh')
    floorVisu.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofsFloor", output="@visualFloor")


    return root

if __name__ == '__main__':
    main()

or to replace both the floor and the liver with two surface meshes, but in both cases the two objects pass through each other without colliding.

hugtalbot Aug 15, 2024
Maintainer

Hi @VeLoXXHH

You do not use the surface mesh (e.g. "LiverSurface" in your simulation as collision model)
It should look like for the liver object:

...
    root.addObject("MeshSTLLoader", name="LiverSurface", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_segmentation\Liver_Target_1.stl")
    root.addObject("MeshSTLLoader", name="Floor", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\Platform\plane.stl")

    #Liver Parenchyma with TetrahedralCorotationalFEMForceField
    liver = root.addChild('Liver')
    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    liver.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    liver.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.msh", translation = "0 0 800")  
    liver.addObject('TetrahedronSetTopologyContainer', name="tetra", src="@meshLoader")
    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="TetraGeomAlgo")
    liver.addObject('DiagonalMass', name="Mass", massDensity="1.0")
    liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="TetraFEM", method="large", poissonRatio="0.48", youngModulus="3000")
    
   ############### COLLISION LIVER ##########################
    liverCollision = liver.addChild('liverCollision')
    liverCollision.addObject("MeshTopology", src="@../../LiverSurface")
    liverCollision.addObject("MechanicalObject", name="StoreContactForces")
    liverCollision.addObject("TriangleCollisionModel", selfCollision="0")
    liverCollision.addObject("LineCollisionModel", selfCollision="0")
    liverCollision.addObject("PointCollisionModel", selfCollision="0")
    liverCollision.addObject('BarycentricMapping', name="CollisionModelMapping", input="@../dofs", output="@StoreContactForces")
    ######################################################

    # Glisson capsule with TriangularBendingSprings
    liverCapsule = liver.addChild('GlissonCapsule')
    liverCapsule.addObject('MeshSTLLoader', name='CapsuleMesh', filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.stl", translation = "0 0 800")
    liverCapsule.addObject('MechanicalObject', name='dofsCapsule', src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetTopologyContainer', name="tri", src="@CapsuleMesh")
    liverCapsule.addObject('TriangleSetGeometryAlgorithms', template="Vec3d", name="TriGeomAlgo")
    liverCapsule.addObject('TriangularBendingSprings', template="Vec3d", name="BendingSprings", stiffness=10000, damping=0.2)
    liverCapsule.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofs", output="@dofsCapsule")

    #visualization Parenchyma
    liverVisu = liver.addChild('visualization')
    liverVisu.addObject("MeshGmshLoader", name="meshLoader", filename = r"C:\Users\monta\Desktop\Materie\Altro\Tirocinio\Liver_Database\3D_Slicer\liver_test\Liver_Test.msh", translation = "@../meshLoader.transition")
    liverVisu.addObject('OglModel', name='VisualModelParenchyma', color='0.8 0.2 0.2 1.0', src='@../meshLoader')
    liverVisu.addObject('BarycentricMapping', name="VisualMappingParenchyma", input="@../dofs", output="@VisualModelParenchyma")

    #visualization Glisson's Capsule
    capsuleVisu = liverCapsule.addChild('capsule-visualization')
    capsuleVisu.addObject('OglModel', name='visualCapsule', color='0.2 0.8 0.2 1.0', src='@../CapsuleMesh')
    capsuleVisu.addObject('IdentityMapping')

    #floor 
    impactFloor = root.addChild('impactFloor')
    impactFloor.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
    impactFloor.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
    impactFloor.addObject('MeshOBJLoader', name='floorMesh', filename="mesh/floor3.obj", rotation = "-90 0 0")
    impactFloor.addObject("MeshTopology", src="@floorMesh")
    impactFloor.addObject('MechanicalObject', name='dofsFloor', src="@floorMesh", scale="8")
    impactFloor.addObject('FixedConstraint', name="FixedConstraint", fixAll = True)
    

    ############### COLLISION FLOOR ##########################

    floorCollision = impactFloor.addChild('floorCollision')
    floorCollision.addObject("MeshTopology", src="@../../Floor")
    floorCollision.addObject("MechanicalObject", name="StoreContactForces")
    floorCollision.addObject("TriangleCollisionModel",name="FloorTriangle", simulated="false", moving="false")
    floorCollision.addObject("LineCollisionModel",name="FloorLine", simulated="false", moving="false")
    floorCollision.addObject("PointCollisionModel",name="FloorPoint", simulated="false", moving="false")
    floorCollision.addObject('BarycentricMapping', name="CollisionModelMapping", input="@../dofsFloor", output="@StoreContactForces")


    ######################################################


    #visualization floor
    floorVisu = impactFloor.addChild('floor-visualization')
    floorVisu.addObject('MeshOBJLoader', name='floorMesh', filename="mesh/floor3.obj", rotation = "@../floorMesh.rotation")
    floorVisu.addObject('OglModel', name='visualFloor', color='0.2 0.2 0.8 1.0', src='@../floorMesh')
    floorVisu.addObject('BarycentricMapping', name="VisualMappingCapsule", input="@../dofsFloor", output="@visualFloor")


    return root

if __name__ == '__main__':
    main()

hugtalbot Aug 15, 2024
Maintainer

@VeLoXXHH

You should definitely watch the SOFA basic tutorial before making your own simulation : see this part of the training (add a visual and a collision model)

@VeLoXXHH · 2024-10-21T06:16:19Z

github-actions[bot]
bot Oct 21, 2024

⚠️ ⚠️ ⚠️
@VeLoXXHH
Feedback has been given to you by the project reviewers, however we have not received a response from you. Without further news in the coming weeks, this discussion will be automatically closed in order to keep this forum clean and fresh 🌱 Thank you for your understanding

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

SOFAPython3_DynamicLiverModel #4862

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 8 comments 7 replies

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

Select a reply

SOFAPython3_DynamicLiverModel #4862

VeLoXXHH Jul 29, 2024

Replies: 8 comments · 7 replies

hugtalbot Jul 29, 2024 Maintainer

VeLoXXHH Jul 30, 2024 Author

hugtalbot Jul 30, 2024 Maintainer

VeLoXXHH Aug 1, 2024 Author

VeLoXXHH Aug 12, 2024 Author

hugtalbot Aug 12, 2024 Maintainer

VeLoXXHH Aug 12, 2024 Author

VeLoXXHH Aug 12, 2024 Author

hugtalbot Aug 13, 2024 Maintainer

VeLoXXHH Aug 13, 2024 Author

hugtalbot Aug 13, 2024 Maintainer

VeLoXXHH Aug 14, 2024 Author

hugtalbot Aug 15, 2024 Maintainer

hugtalbot Aug 15, 2024 Maintainer

github-actions[bot] bot Oct 21, 2024

VeLoXXHH
Jul 29, 2024

Replies: 8 comments 7 replies

hugtalbot
Jul 29, 2024
Maintainer

VeLoXXHH
Jul 30, 2024
Author

hugtalbot Jul 30, 2024
Maintainer

VeLoXXHH
Aug 1, 2024
Author

VeLoXXHH
Aug 12, 2024
Author

hugtalbot Aug 12, 2024
Maintainer

VeLoXXHH Aug 12, 2024
Author

VeLoXXHH
Aug 12, 2024
Author

hugtalbot Aug 13, 2024
Maintainer

VeLoXXHH
Aug 13, 2024
Author

hugtalbot
Aug 13, 2024
Maintainer

VeLoXXHH Aug 14, 2024
Author

hugtalbot Aug 15, 2024
Maintainer

hugtalbot Aug 15, 2024
Maintainer

github-actions[bot]
bot Oct 21, 2024