Sofa cgal tetrahedron push with repulsive force by FadiWahbeh · Pull Request #4 · InfinyTech3D/CGALPlugin

FadiWahbeh · 2026-07-09T12:41:37Z

Based on #2

PS : don't merge or take in consideration PR #3 because it's not used and won't be used. Also, GitHub didn't let me add the .msh objects.

PR 02 — CGAL Detection → SOFA ForceField

Goal

This PR turns MeshRepulsionEngine into a real ForceField. It uses the PR 01 CGAL detection to generate penalty forces that push the slave out of the master, and lets the SOFA solver integrate those forces together with mass, FEM and damping. CGAL is used only for geometric detection and contact direction; SOFA integrates the motion.

The component now:

is a ForceField<DataTypes> (templated, Vec3Types), implementing addForce, addDForce, addKToMatrix and getPotentialEnergy;
builds the master CGAL surface mesh, its AABB tree and its Side_of_triangle_mesh inside/outside tester once (cached), since the master is static;
detects four kinds of contact each step, each one individually switchable:
- Inside: slave vertex strictly inside the master;
- Proximity: outside vertex still within safetyDistance of the surface;
- EdgeCrossing: endpoints of a slave edge crossing the master surface;
- FaceCrossing: vertices of a slave triangle intersecting the master;
merges contacts per vertex, keeping the most severe one (mergeContact);
for each contact, applies a penalty force F = stiffness * penetration * n along the master-oriented normal (orientedNormalFromClosestFace), with:
- crossingDepth giving edge/face crossings a force even when the vertex is outside the safety band,
- maxForce clamping the per-vertex magnitude,
- contactDamping opposing the vertex velocity;
optionally spreads each contact force over a local topological neighborhood (neighborhoodRings, neighborWeight) to better preserve the slave shape (applyDistributedForce);
provides the consistent stiffness term to the implicit solver (addDForce, addKToMatrix);
prints diagnostics and writes a CSV.

Data / API

New Data exposed by the component:

stiffness, safetyDistance, crossingDepth, maxForce, contactDamping
neighborhoodRings, neighborWeight
pushInsideVertices, pushEdgeCrossings, pushFaceCrossings, pushProximityVertices
rebuildMasterEveryStep
logEveryNSteps, writeCsv, csvFilePath

Read-only diagnostics:

lastStep, lastInsideContacts, lastEdgeContacts, lastFaceContacts,
lastProximityContacts, lastMaxForce, lastMeanForce

Performance

Two fixes were required to keep the scene interactive:

the Side_of_triangle_mesh tester and the master AABB tree are built once, not per contact (they were previously rebuilt for every contact, which dropped the scene to ~0.5-1 FPS);
on a tetrahedral slave, only surface vertices are tested (m_isSurfaceVertex mask); interior vertices follow the surface through the FEM.

Modify Test Scene

scenes/MeshRepulsion_PR02.scn

The slave now carries a solver, a mass, and the repulsion ForceField:

<EulerImplicitSolver/>
<CGLinearSolver iterations="40" tolerance="1e-6" threshold="1e-6"/>
<MechanicalObject/>
<UniformMass totalMass="0.5"/>
<MeshRepulsionEngine
    masterState="@../MasterMesh/dofs"
    masterTopology="@../MasterMesh/topology"
    slaveTopology="@topology"
    stiffness="400" safetyDistance="2.0" crossingDepth="3.0"
    maxForce="800" contactDamping="20"
    neighborhoodRings="2" neighborWeight="0.18"
    pushInsideVertices="true" pushEdgeCrossings="true"
    pushFaceCrossings="true" pushProximityVertices="true"/>

Test Scenes

Three patient scenes are provided. Each one loads its own liver (master) and diaphragm (slave) and runs the same repulsion component.

The slave is a tetrahedral volume driven by a TetrahedronFEMForceField: CGAL reads the boundary surface triangles to detect contacts, while the volumetric FEM integrates the resulting forces and preserves the organ shape. Surface and volume share the same DOFs, so a force applied by CGAL on a surface vertex propagates through the volume.

Click to expand XML configuration — Patient 01

<?xml version="1.0"?>
<Node name="root" gravity="0 0 0" dt="0.003">
    <RequiredPlugin pluginName="Sofa.Component.AnimationLoop"/>
    <RequiredPlugin pluginName="Sofa.Component.IO.Mesh"/>
    <RequiredPlugin pluginName="Sofa.Component.Mapping.Linear"/>
    <RequiredPlugin pluginName="Sofa.Component.Setting"/>
    <RequiredPlugin pluginName="Sofa.Component.StateContainer"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Container.Constant"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Container.Dynamic"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Mapping"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Utility"/>
    <RequiredPlugin pluginName="Sofa.Component.Visual"/>
    <RequiredPlugin pluginName="Sofa.GL.Component.Rendering3D"/>
    <RequiredPlugin pluginName="Sofa.Component.LinearSolver.Iterative"/>
    <RequiredPlugin pluginName="Sofa.Component.Mass"/>
    <RequiredPlugin pluginName="Sofa.Component.ODESolver.Backward"/>
    <RequiredPlugin pluginName="Sofa.Component.SolidMechanics.FEM.Elastic"/>
    <RequiredPlugin pluginName="Sofa.Component.SolidMechanics.Spring"/>
    <RequiredPlugin pluginName="CGALPlugin"/>

    <BackgroundSetting color="0 0.16862745098 0.21176470588"/>
    <VisualStyle displayFlags="showVisualModels"/>
    <DefaultAnimationLoop/>

    <Node name="MasterLiver">
        <MeshGmshLoader name="loader" filename="C:/Users/Fadi WAHBEH/Desktop/STL_Intersection_Objects/Patient_01/liver_coarse_patient_01.msh"/>
        <MechanicalObject name="dofs" template="Vec3d" position="@loader.position" showObject="false"/>
        <TriangleSetTopologyContainer name="topology" src="@loader"/>
        <TriangleSetTopologyModifier name="topologyModifier"/>
        <TriangleSetGeometryAlgorithms name="geometryAlgorithms" template="Vec3d"/>

        <Node name="Visual">
            <OglModel name="visu" src="@../loader" color="1 0.05 0.05 1.0"/>
            <IdentityMapping input="@../dofs" output="@visu"/>
        </Node>
    </Node>

    <Node name="SlaveDiaphragmVolume">
        <EulerImplicitSolver name="odesolver" rayleighStiffness="0.25" rayleighMass="0.45"/>
        <CGLinearSolver name="linearsolver" iterations="60" tolerance="1.0e-6" threshold="1.0e-6"/>

        <MeshGmshLoader name="volumeLoader" filename="C:/Users/Fadi WAHBEH/Desktop/STL_Intersection_Objects/Patient_01/diaphragm_patient_01.msh"/>
        <MechanicalObject name="dofs" template="Vec3d" position="@volumeLoader.position" showObject="false"/>
        <TetrahedronSetTopologyContainer name="topology" src="@volumeLoader"/>
        <TetrahedronSetTopologyModifier name="topologyModifier"/>
        <TetrahedronSetGeometryAlgorithms name="geometryAlgorithms" template="Vec3d"/>

        <UniformMass totalMass="1.2"/>
        <TetrahedronFEMForceField name="tetraFEM" template="Vec3d" method="large" youngModulus="5000" poissonRatio="0.45"/>

        <Node name="ContactSurface">
            <MechanicalObject name="surfaceDofs" template="Vec3d" position="@../volumeLoader.position" showObject="false"/>
            <TriangleSetTopologyContainer name="surfaceTopology" src="@../volumeLoader"/>
            <TriangleSetTopologyModifier name="surfaceTopologyModifier"/>
            <TriangleSetGeometryAlgorithms name="surfaceGeometryAlgorithms" template="Vec3d"/>
            <BarycentricMapping name="surfaceToVolumeMapping" input="@../dofs" output="@surfaceDofs" mapForces="true" mapMasses="false"/>

            <RepulsiveSpringForceField name="diaphragmRepulsiveSprings" template="Vec3d"
                object1="@surfaceDofs"
                object2="@../../MasterLiver/dofs"
                stiffness="850"
                damping="75"
                lengths=""
                enabled=""
                springsIndices1=""
                springsIndices2=""
                drawMode="0"
                showArrowSize="0.01"/>

            <MeshRepulsionEngine name="diaphragmCgalSpringPairs" template="Vec3d"
                masterState="@../../MasterLiver/dofs"
                masterTopology="@../../MasterLiver/topology"
                slaveState="@surfaceDofs"
                slaveTopology="@surfaceTopology"
                repulsiveSpringField="@diaphragmRepulsiveSprings"
                safetyDistance="1.2"
                crossingDepth="2.2"
                pushInsideVertices="true"
                pushEdgeCrossings="true"
                pushFaceCrossings="false"
                rebuildMasterEveryStep="false"
                computeEveryNSteps="4"
                maxPairs="2500"
                maxEdgesCheckedPerUpdate="8000"
                maxFacesCheckedPerUpdate="0"
                logEveryNSteps="20"/>

            <Node name="Visual">
                <OglModel name="visu" src="@../../volumeLoader" color="0.15 0.55 1.0 0.35"/>
                <IdentityMapping input="@../surfaceDofs" output="@visu"/>
            </Node>
        </Node>
    </Node>
</Node>

Click to expand XML configuration — Patient 08

<?xml version="1.0"?>
<Node name="root" gravity="0 0 0" dt="0.003">
    <RequiredPlugin pluginName="Sofa.Component.AnimationLoop"/>
    <RequiredPlugin pluginName="Sofa.Component.IO.Mesh"/>
    <RequiredPlugin pluginName="Sofa.Component.Mapping.Linear"/>
    <RequiredPlugin pluginName="Sofa.Component.Setting"/>
    <RequiredPlugin pluginName="Sofa.Component.StateContainer"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Container.Constant"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Container.Dynamic"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Mapping"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Utility"/>
    <RequiredPlugin pluginName="Sofa.Component.Visual"/>
    <RequiredPlugin pluginName="Sofa.GL.Component.Rendering3D"/>
    <RequiredPlugin pluginName="Sofa.Component.LinearSolver.Iterative"/>
    <RequiredPlugin pluginName="Sofa.Component.Mass"/>
    <RequiredPlugin pluginName="Sofa.Component.ODESolver.Backward"/>
    <RequiredPlugin pluginName="Sofa.Component.SolidMechanics.FEM.Elastic"/>
    <RequiredPlugin pluginName="Sofa.Component.SolidMechanics.Spring"/>
    <RequiredPlugin pluginName="CGALPlugin"/>

    <BackgroundSetting color="0 0.16862745098 0.21176470588"/>
    <VisualStyle displayFlags="showVisualModels"/>
    <DefaultAnimationLoop/>

    <Node name="MasterLiver">
        <MeshGmshLoader name="loader" filename="C:/Users/Fadi WAHBEH/Desktop/STL_Intersection_Objects/Patient_08/liver_coarse_patient_08.msh"/>
        <MechanicalObject name="dofs" template="Vec3d" position="@loader.position" showObject="false"/>
        <TriangleSetTopologyContainer name="topology" src="@loader"/>
        <TriangleSetTopologyModifier name="topologyModifier"/>
        <TriangleSetGeometryAlgorithms name="geometryAlgorithms" template="Vec3d"/>

        <Node name="Visual">
            <OglModel name="visu" src="@../loader" color="1 0.05 0.05 1.0"/>
            <IdentityMapping input="@../dofs" output="@visu"/>
        </Node>
    </Node>

    <Node name="SlaveDiaphragmVolume">
        <EulerImplicitSolver name="odesolver" rayleighStiffness="0.25" rayleighMass="0.45"/>
        <CGLinearSolver name="linearsolver" iterations="60" tolerance="1.0e-6" threshold="1.0e-6"/>

        <MeshGmshLoader name="volumeLoader" filename="C:/Users/Fadi WAHBEH/Desktop/STL_Intersection_Objects/Patient_08/diaphragm_patient_08.msh"/>
        <MechanicalObject name="dofs" template="Vec3d" position="@volumeLoader.position" showObject="false"/>
        <TetrahedronSetTopologyContainer name="topology" src="@volumeLoader"/>
        <TetrahedronSetTopologyModifier name="topologyModifier"/>
        <TetrahedronSetGeometryAlgorithms name="geometryAlgorithms" template="Vec3d"/>

        <UniformMass totalMass="1.2"/>
        <TetrahedronFEMForceField name="tetraFEM" template="Vec3d" method="large" youngModulus="5000" poissonRatio="0.45"/>

        <Node name="ContactSurface">
            <MechanicalObject name="surfaceDofs" template="Vec3d" position="@../volumeLoader.position" showObject="false"/>
            <TriangleSetTopologyContainer name="surfaceTopology" src="@../volumeLoader"/>
            <TriangleSetTopologyModifier name="surfaceTopologyModifier"/>
            <TriangleSetGeometryAlgorithms name="surfaceGeometryAlgorithms" template="Vec3d"/>
            <BarycentricMapping name="surfaceToVolumeMapping" input="@../dofs" output="@surfaceDofs" mapForces="true" mapMasses="false"/>

            <RepulsiveSpringForceField name="diaphragmRepulsiveSprings" template="Vec3d"
                object1="@surfaceDofs"
                object2="@../../MasterLiver/dofs"
                stiffness="850"
                damping="75"
                lengths=""
                enabled=""
                springsIndices1=""
                springsIndices2=""
                drawMode="0"
                showArrowSize="0.01"/>

            <MeshRepulsionEngine name="diaphragmCgalSpringPairs" template="Vec3d"
                masterState="@../../MasterLiver/dofs"
                masterTopology="@../../MasterLiver/topology"
                slaveState="@surfaceDofs"
                slaveTopology="@surfaceTopology"
                repulsiveSpringField="@diaphragmRepulsiveSprings"
                safetyDistance="1.2"
                crossingDepth="2.2"
                pushInsideVertices="true"
                pushEdgeCrossings="true"
                pushFaceCrossings="true"
                rebuildMasterEveryStep="true"
                computeEveryNSteps="4"
                maxPairs="2500"
                maxEdgesCheckedPerUpdate="8000"
                maxFacesCheckedPerUpdate="0"
                logEveryNSteps="20"/>

            <Node name="Visual">
                <OglModel name="visu" src="@../../volumeLoader" color="0.15 0.55 1.0 0.35"/>
                <IdentityMapping input="@../surfaceDofs" output="@visu"/>
            </Node>
        </Node>
    </Node>
</Node>

Click to expand XML configuration — Patient 16

<?xml version="1.0"?>
<Node name="root" gravity="0 0 0" dt="0.003">
    <RequiredPlugin pluginName="Sofa.Component.AnimationLoop"/>
    <RequiredPlugin pluginName="Sofa.Component.IO.Mesh"/>
    <RequiredPlugin pluginName="Sofa.Component.Mapping.Linear"/>
    <RequiredPlugin pluginName="Sofa.Component.Setting"/>
    <RequiredPlugin pluginName="Sofa.Component.StateContainer"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Container.Constant"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Container.Dynamic"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Mapping"/>
    <RequiredPlugin pluginName="Sofa.Component.Topology.Utility"/>
    <RequiredPlugin pluginName="Sofa.Component.Visual"/>
    <RequiredPlugin pluginName="Sofa.GL.Component.Rendering3D"/>
    <RequiredPlugin pluginName="Sofa.Component.LinearSolver.Iterative"/>
    <RequiredPlugin pluginName="Sofa.Component.Mass"/>
    <RequiredPlugin pluginName="Sofa.Component.ODESolver.Backward"/>
    <RequiredPlugin pluginName="Sofa.Component.SolidMechanics.FEM.Elastic"/>
    <RequiredPlugin pluginName="Sofa.Component.SolidMechanics.Spring"/>
    <RequiredPlugin pluginName="CGALPlugin"/>

    <BackgroundSetting color="0 0.16862745098 0.21176470588"/>
    <VisualStyle displayFlags="showVisualModels"/>
    <DefaultAnimationLoop/>

    <Node name="MasterLiver">
        <MeshGmshLoader name="loader" filename="C:/Users/Fadi WAHBEH/Desktop/STL_Intersection_Objects/Patient_16/liver_coarse_patient_16.msh"/>
        <MechanicalObject name="dofs" template="Vec3d" position="@loader.position" showObject="false"/>
        <TriangleSetTopologyContainer name="topology" src="@loader"/>
        <TriangleSetTopologyModifier name="topologyModifier"/>
        <TriangleSetGeometryAlgorithms name="geometryAlgorithms" template="Vec3d"/>

        <Node name="Visual">
            <OglModel name="visu" src="@../loader" color="1 0.05 0.05 1.0"/>
            <IdentityMapping input="@../dofs" output="@visu"/>
        </Node>
    </Node>

    <Node name="SlaveDiaphragmVolume">
        <EulerImplicitSolver name="odesolver" rayleighStiffness="0.25" rayleighMass="0.45"/>
        <CGLinearSolver name="linearsolver" iterations="60" tolerance="1.0e-6" threshold="1.0e-6"/>

        <MeshGmshLoader name="volumeLoader" filename="C:/Users/Fadi WAHBEH/Desktop/STL_Intersection_Objects/Patient_16/diaphragm_patient_16.msh"/>
        <MechanicalObject name="dofs" template="Vec3d" position="@volumeLoader.position" showObject="false"/>
        <TetrahedronSetTopologyContainer name="topology" src="@volumeLoader"/>
        <TetrahedronSetTopologyModifier name="topologyModifier"/>
        <TetrahedronSetGeometryAlgorithms name="geometryAlgorithms" template="Vec3d"/>

        <UniformMass totalMass="1.2"/>
        <TetrahedronFEMForceField name="tetraFEM" template="Vec3d" method="large" youngModulus="5000" poissonRatio="0.45"/>

        <Node name="ContactSurface">
            <MechanicalObject name="surfaceDofs" template="Vec3d" position="@../volumeLoader.position" showObject="false"/>
            <TriangleSetTopologyContainer name="surfaceTopology" src="@../volumeLoader"/>
            <TriangleSetTopologyModifier name="surfaceTopologyModifier"/>
            <TriangleSetGeometryAlgorithms name="surfaceGeometryAlgorithms" template="Vec3d"/>
            <BarycentricMapping name="surfaceToVolumeMapping" input="@../dofs" output="@surfaceDofs" mapForces="true" mapMasses="false"/>

            <RepulsiveSpringForceField name="diaphragmRepulsiveSprings" template="Vec3d"
                object1="@surfaceDofs"
                object2="@../../MasterLiver/dofs"
                stiffness="850"
                damping="75"
                lengths=""
                enabled=""
                springsIndices1=""
                springsIndices2=""
                drawMode="0"
                showArrowSize="0.01"/>

            <MeshRepulsionEngine name="diaphragmCgalSpringPairs" template="Vec3d"
                masterState="@../../MasterLiver/dofs"
                masterTopology="@../../MasterLiver/topology"
                slaveState="@surfaceDofs"
                slaveTopology="@surfaceTopology"
                repulsiveSpringField="@diaphragmRepulsiveSprings"
                safetyDistance="1.2"
                crossingDepth="2.2"
                pushInsideVertices="true"
                pushEdgeCrossings="true"
                pushFaceCrossings="false"
                rebuildMasterEveryStep="false"
                computeEveryNSteps="4"
                maxPairs="2500"
                maxEdgesCheckedPerUpdate="8000"
                maxFacesCheckedPerUpdate="0"
                logEveryNSteps="20"/>

            <Node name="Visual">
                <OglModel name="visu" src="@../../volumeLoader" color="0.15 0.55 1.0 0.35"/>
                <IdentityMapping input="@../surfaceDofs" output="@visu"/>
            </Node>
        </Node>
    </Node>
</Node>

Observation on Patient 16

Patient 16 and 08 are the case that the slave rotates and re-enters the master from the opposite side before coming back out. This is an unanchored deformable body under one-sided contact forces, pushing one side out injects momentum that might makes the free body pivot, so the opposite side dips into the master, generates its own contacts, and gets pushed out in turn. Patients 01 is more stable. The solution to this issue is to try and fixe a region of the slave so it wouldn't pevot twords the master.

Build

"/c/Program Files/CMake/bin/cmake.exe" --build /c/projects/sofa-build --config Release --target CGALPlugin

Run

cd /c/projects/sofa-build/bin/Release ./runSofa.exe /c/projects/sofa-plugins/CGALPlugin/scenes/MeshRepulsion_PR02.scn

What This PR Proves

the CGAL detection can drive a SOFA-integrated force, not only a diagnostic;
the repulsion is solved together with the FEM, so the slave shape is preserved while the contact zone is pushed out;
edge/face crossings are handled even when no vertex is inside;
the correction is stable and the inside count converges toward 0.

Not Included Yet

.msh export;
automatic per-patient anchoring.

…shes

…t not tested

FadiWahbeh added pr: enhancement pr: status to review pr: new feature labels Jul 9, 2026

[Intersection Push] Adding Repulsive Force Field using tetrahedron me…

89860d5

…shes

epernod force-pushed the SOFA_CGAL_tetrahedron_push_with_Repulsive_Force branch from 7e795f5 to 89860d5 Compare July 9, 2026 13:48

epernod changed the base branch from master to dev_IT3D July 9, 2026 13:48

epernod added 2 commits July 9, 2026 17:13

Update scene to use a RestShapeSpringFF as boundary condition

803071b

[scenes] Add node to map collision obj mesh. Also add MeshExporter bu…

7ea8d2a

…t not tested

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Sofa cgal tetrahedron push with repulsive force#4

Sofa cgal tetrahedron push with repulsive force#4
FadiWahbeh wants to merge 3 commits into
dev_IT3Dfrom
SOFA_CGAL_tetrahedron_push_with_Repulsive_Force

FadiWahbeh commented Jul 9, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

FadiWahbeh commented Jul 9, 2026

PR 02 — CGAL Detection → SOFA ForceField

Goal

Data / API

Performance

Modify Test Scene

Test Scenes

Observation on Patient 16

Build

Run

What This PR Proves

Not Included Yet

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants