How to calculate contact area

[KumaAirmoeba]

Title

How to calculate contact area

---

Description

Hello. I'm trying to calculate the contact area between a soft cube and a floor using SOFA v24.12.00 with Python3.
Although the cube visibly touches the floor during the simulation, the contact area is always reported as "no contact".

I'm using ContactListener and a custom Python controller to compute the area based on contact triangles.

---

Scene Overview

• A soft cube (tetrahedral mesh) is dropped onto a flat floor
• Collision models are added to both the cube and the floor
• ContactListener is used to detect contact
• A Python controller reads getContactPoints() and calculates area from triangle indices

---

ContactAreaController Code

class` ContactAreaController(Sofa.Core.Controller):
    def __init__(self, node):
        super().__init__()
        self.node = node
        self.cube = node.getRoot().getChild("cube")
        self.collisionNode = self.cube.getChild("modelCollis")

        self.positions = self.collisionNode.getObject("collisMO")
        self.topology = self.collisionNode.getObject("container_quadriped")
        self.contactListener = node.getObject("listener")
        self.frame_counter = 0

    def onAnimateBeginEvent(self, event):
        self.frame_counter += 1
        if self.frame_counter % 10 != 0:
            return

        contact_points = self.contactListener.getContactPoints()
        print(f"[Frame {self.frame_counter}] Number of contact points: {len(contact_points)}")

        if not contact_points:
            print(f"[Frame {self.frame_counter}] No contact")
            return

        contact_tri_ids = set()
        for i, contact in enumerate(contact_points):
            obj1, obj2 = contact.object1.getPathName(), contact.object2.getPathName()
            elem1, elem2 = contact.elem1, contact.elem2

            if "modelCollis" in obj1 and elem1 >= 0:
                contact_tri_ids.add(elem1)
            elif "modelCollis" in obj2 and elem2 >= 0:
                contact_tri_ids.add(elem2)

        if not contact_tri_ids:
            print(f"[Frame {self.frame_counter}] No valid contact triangles")
            return

        total_area = 0.0
        positions = self.positions.position.array()
        triangles = self.topology.triangles.array()

        print(f"[Frame {self.frame_counter}] Number of contact triangles: {len(contact_tri_ids)}")
        for tid in contact_tri_ids:
            if tid < 0 or tid >= len(triangles):
                print(f"  tri[{tid}] is an invalid index")
                continue
            tri = triangles[tid]
            p0, p1, p2 = positions[tri[0]], positions[tri[1]], positions[tri[2]]
            area = 0.5 * np.linalg.norm(np.cross(p1 - p0, p2 - p0))
            total_area += area
            print(f"  tri[{tid}]: vertices = {tri}, area = {area:.6f} m^2")

        print(f"[Frame {self.frame_counter}] Total contact area: {total_area:.6f} m^2")

---

What I Have Checked

• ContactListener is connected between /cube/modelCollis/triangle1 and /plane/triangle2
• Both objects use TriangleCollisionModel
• The cube visually touches the floor
• The position of collisMO indicates the cube is near the floor (Y ≈ -30)
• Gravity is applied, the cube falls onto the floor
• contactDistance is set to 2.0

---

Problem

Even though the cube clearly touches the floor, getContactPoints() always returns an empty list, and the script prints "No contact".

---

Question

Why is the contact not detected?
Is there anything wrong with how I'm using ContactListener, TriangleCollisionModel, or the collision setup?

---

Environment

• SOFA version: 24.12.00
• Operating System: [e.g. Ubuntu 22.04 or Windows 11]
• Python version: 3.10
• Mesh files: cube.msh, cube.stl, plane.vtk

Thank you in advance for your help!

[hugtalbot]

Hi @KumaAirmoeba

Thanks a lot for this detailed post.
I think the problem resides in using the onAnimateBeginEvent function in your controller. This function is called at the start of the time step while the collision detection did not take place yet. You should rather use the onAnimateEndEvent so that all collision info can be recovered.

Let us know if it helps.

[KumaAirmoeba]

Thank you for your reply.
I changed the name of the function as you suggested, but the result was the same and the contact was still not detected. Is there a problem with another part of the code rather than this function?

[bakpaul]

Hello, is it possible for you to share your scene ? IF so I could try it out on my end to see what is going wrong.

[KumaAirmoeba]

Thank you for your reply, Paul, and I apologize for the delayed response.
The code is as follows:

import Sofa
import numpy as np
import Sofa.Core
import SofaRuntime

# Controller to compute the contact area between softBody and the plane
class ContactAreaController(Sofa.Core.Controller):
    def __init__(self, node):
        super().__init__()
        self.node = node
        self.softBody = node.getRoot().getChild("softBody")
        self.listener = node.getObject("ContactListener")
        self.topology = self.softBody.getChild("modelCollis").getObject('container_quadriped')
        self.positions = self.softBody.getChild("modelCollis").getObject('collisMO')
        self.frame_counter = 0

    def onAnimateEndEvent(self, event):
        self.frame_counter += 1

        # Run every 10 frames to reduce verbosity
        if self.frame_counter % 10 != 0:
            return

        # Get contact information
        contacts = self.listener.getContactPoints()
        print(f"[Frame {self.frame_counter}] Contact points detected: {len(contacts)}")

        contact_tri_ids = set()
        for contact in contacts:
            tid1 = contact[0]
            tid2 = contact[2]
            if tid1 != -1:
                contact_tri_ids.add(tid1)
            if tid2 != -1:
                contact_tri_ids.add(tid2)

        # Compute contact area from triangle indices
        total_area = 0.0
        positions = self.positions.position.array()
        triangles = self.topology.triangles.array()

        for tid in contact_tri_ids:
            tri = triangles[tid]
            p0, p1, p2 = positions[tri[0]], positions[tri[1]], positions[tri[2]]
            area = 0.5 * np.linalg.norm(np.cross(p1 - p0, p2 - p0))
            total_area += area

        print(f"[Frame {self.frame_counter}] Contact area: {total_area:.6f} m^2, Contact triangles: {len(contact_tri_ids)}")


# Scene creation function
def createScene(rootNode):
    rootNode.gravity = [0, -9.81, 0]
    rootNode.dt = 0.01

    # Load required plugins
    rootNode.addObject('RequiredPlugin', name='SoftRobots')
    rootNode.addObject('RequiredPlugin', name='SofaPython3')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Visual')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Collision.Response.Contact')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Collision.Detection.Algorithm')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Collision.Detection.Intersection')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Collision.Geometry')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.StateContainer')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Mass')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.SolidMechanics.FEM.Elastic')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.ODESolver.Backward')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.LinearSolver.Direct')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.IO.Mesh')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Topology.Container.Constant')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Mapping.Linear')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.SolidMechanics.FEM.HyperElastic')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.AnimationLoop')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Constraint.Lagrangian.Correction')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Constraint.Lagrangian.Solver')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Constraint.Projective')
    rootNode.addObject('RequiredPlugin', name='Sofa.Component.Topology.Container.Dynamic')
    rootNode.addObject('RequiredPlugin', name='Sofa.GL.Component.Rendering3D')

    # Scene visual and simulation control setup
    rootNode.addObject('VisualStyle', displayFlags='showForceFields showBehaviorModels')
    rootNode.addObject('DefaultAnimationLoop')
    rootNode.addObject('FreeMotionAnimationLoop')
    rootNode.addObject('GenericConstraintSolver', tolerance=1e-7, maxIterations=3000)
    rootNode.addObject('CollisionPipeline', verbose=1)
    rootNode.addObject('BruteForceBroadPhase', name='N2')
    rootNode.addObject('BVHNarrowPhase')

    # Collision response configuration (rule-based and default)
    rootNode.addObject('RuleBasedContactManager', name='Response',
                       response='FrictionContactConstraint',
                       responseParams='mu=0.5',
                       rules=["1 2"])
    rootNode.addObject('CollisionResponse', response='FrictionContactConstraint', responseParams='mu=0.7')
    rootNode.addObject('LocalMinDistance', name='Proximity',
                       alarmDistance='20', contactDistance='15', angleCone='0.1')
    rootNode.addObject('OglSceneFrame', style='Arrows', alignment='TopRight')

    # Create deformable soft body cube
    softBody = rootNode.addChild('softBody')
    softBody.addObject('EulerImplicitSolver', name='odesolver', firstOrder=False)
    softBody.addObject('SparseLDLSolver', name='linearSolver', template='CompressedRowSparseMatrixMat3x3d')
    softBody.addObject('MeshGmshLoader', name='loader', filename='mesh/cube.msh')
    softBody.addObject('MeshTopology', src='@loader', name='container')
    softBody.addObject('MechanicalObject', name='tetlas', template='Vec3', showObject=True, showObjectScale=1)
    softBody.addObject('TetrahedronFEMForceField', template='Vec3', name='FEM',
                       method='large', poissonRatio=0.3, youngModulus=2000)
    softBody.addObject('UniformMass', totalMass=0.1)
    softBody.addObject('GenericConstraintCorrection')

    # Collision model for soft body
    collis = softBody.addChild('modelCollis')
    collis.addObject('MeshSTLLoader', name='loader', filename='mesh/cube.stl')
    collis.addObject('TriangleSetTopologyContainer', position='@loader.position', triangles='@loader.triangles', name='container_quadriped')
    collis.addObject('MechanicalObject', name='collisMO', template='Vec3')
    model1 = collis.addObject('TriangleCollisionModel', group=1, name="triangle1")
    collis.addObject('LineCollisionModel', group=1)
    collis.addObject('PointCollisionModel', group=1)
    collis.addObject('BarycentricMapping',
                     input='@../tetlas',
                     output='@collisMO',
                     mapForces=True,
                     mapMasses=False)

    # Create floor (static plane)
    plane = rootNode.addChild('plane')
    plane.addObject('MeshVTKLoader', name='loader', filename='mesh/plane.vtk')
    plane.addObject('MeshTopology', src='@loader', name='container')
    plane.addObject('MechanicalObject', name='dofs', translation='0 -30 0')
    plane.addObject('FixedConstraint', indices='0')  # Fix the plane in space
    model2 = plane.addObject('TriangleCollisionModel', simulated='0', moving='0', group=2, name="triangle2")
    plane.addObject('LineCollisionModel', simulated='0', moving='0', group=2)
    plane.addObject('PointCollisionModel', simulated='0', moving='0', group=2)
    plane.addObject('OglModel', name='Visual', src='@loader',
                    color='30 20 150', rotation='0 0 0',
                    translation='0 -30 0', scale='1')

    # Add a ContactListener to detect contact between soft body and plane
    listener = rootNode.addObject(
        "ContactListener",
        collisionModel1="@/softBody/modelCollis/triangle1",
        #collisionModel2="@/plane/triangle2",  # Optional if using RuleBasedContactManager
    )

    # Add custom Python controller to compute contact area
    rootNode.addObject(ContactAreaController(node=rootNode))

    return rootNode

Using this code, the contact area is now being printed correctly.
However, the number of contact points and the number of contact triangles seems incorrect — I believe it should be much higher. For example:

[Frame 90] Contact points detected: 3 [Frame 90] Contact area: 63.831415 m^2, Contact triangles: 2
How can I improve this?

[bakpaul]

HEllo, I am trully sorry for never answering back. Did you finally manage to make it work ? If yes, can you tell me what finally did the trick ? If not, where are you up to now ?