CRESSim–MPM: A Material Point Method Library for Surgical Soft Body Simulation with Cutting and Suturing

Abstract

A number of recent studies have focused on developing surgical simulation platforms to train machine learning (ML) agents or models with synthetic data for surgical assistance. While existing platforms excel at tasks such as rigid body manipulation and soft body deformation, they struggle to simulate more complex soft body behaviors like cutting and suturing. A key challenge lies in modeling soft body fracture and splitting using the finite-element method (FEM), which is the predominant approach in current platforms. Additionally, the two-way suture needle/thread contact inside a soft body is further complicated when using FEM. In this work, we use the material point method (MPM) for such challenging simulations and propose new rigid geometries and soft-rigid contact methods specifically designed for them. We introduce CRESSim-MPM, a GPU-accelerated MPM library that integrates multiple MPM solvers and incorporates surgical geometries for cutting and suturing, serving as a specialized physics engine for surgical applications. It is further integrated into Unity, requiring minimal modifications to existing projects for soft body simulation. We demonstrate the simulator’s capabilities in real-time simulation of cutting and suturing on soft tissue and provide an initial performance evaluation of different MPM solvers when simulating varying numbers of particles. The source code is available at https://github.com/yafei-ou/CRESSim-MPM.

I INTRODUCTION

Realistic and efficient simulators are playing an increasingly important role in robotics research, not only because they serve as a playground for validating robot control and automation methods, but also because they provide as much synthetic data as possible that can be used for training and testing machine learning (ML) models [1, 2, 3, 4, 5, 6], reducing the need of real-world dataset.

While there are numerous high-performance and real-time simulators for other robotics applications, such as autonomous driving (e.g., CARLA) and general robot manipulation (e.g., Isaac Sim), there are limited resources in terms of surgical robotics simulation. Although simulating rigid surgical instruments is usually not challenging, the simulation of surgical scenes and operations requires specific simulation capabilities less present in other robotics fields.

The nature of surgery requires simulating diverse objects, involving rigid bodies (e.g., surgical instruments and bones), soft bodies (e.g., soft tissue), and fluids (e.g., blood and smoke). Additionally, a unique range of manipulation is needed, including cutting, suturing, burning (cauterization), fluid suction, and irrigation, among others. Simulating and rendering each of these objects and the manipulation types can be a unique research field, making it particularly challenging to develop simulators that meet a wide range of these requirements for surgical robotics research. As a large group of surgical robotics researchers focus on increasing the level of autonomy in surgeries using ML, the lack of efficient and realistic surgical simulators remains a significant barrier.

Figure 1: Simulation of soft tissue cutting (upper) and suturing (lower).

Simulating soft tissue and the corresponding surgical interactions such as cutting and suturing is an extensively explored area. A large amount of previous work utilizes the finite element method (FEM) for simulating soft bodies, including the widely used SOFA Framework [7]. While FEM is accurate for biomedical modeling, its expensive computational cost makes it difficult to achieve efficient and parallel simulation often needed by robotics researchers. Interactions such as cutting and suturing further require real-time tetrahedron re-meshing algorithms and complex contact constraints when using FEM, making it more computationally expensive. FEM also suffers from inaccuracies and instabilities under large deformations required in surgical simulation due to the distorted finite elements. Additionally, while GPUs can accelerate certain FEM computations, achieving full parallel efficiency remains difficult despite ongoing efforts [8].

The material point method (MPM) [9] is a promising yet underexplored approach for surgical simulation. The technique, together with its variations and improvements, has shown wide application in diverse fields including solid mechanics, chemical physics, and computer animation. MPM utilizes material points (particles) and a background grid for simulating continuum materials, including both solids and fluids. Despite being a hybrid approach using both Lagrangian particles and an Eulerian grid, it is often viewed as a meshfree method [10] and largely mitigates the disadvantages of the mesh-based FEM. Although it can be more computationally expensive than FEM if implemented on the CPU, its efficiency can be boosted significantly by GPU acceleration. It performs well under large deformations and can simulate fracture without re-meshing. As will be shown in this work, it also enables straightforward rigid coupling for simulating interactions such as cutting and suturing.

Leveraging these strengths, we now present CRESSim-MPM, a GPU-accelerated MPM library tailored for surgical simulation, enabling both soft tissue cutting and suturing within one solver. Our contributions include:

• A CUDA-based MPM library supporting multiple MPM solvers with C-style APIs for efficient integration with external applications and other physics solvers.

• Custom geometries and contact models designed to simulate surgical interactions between MPM soft bodies and scalpels, suture needles, and suture threads.

• A Unity package that allows easy use of the library and enables two-way interactions between MPM-based soft bodies and Unity’s built-in rigid body physics.

Based on these, we demonstrate:

• A preliminary evaluation of the simulation speed using different MPM solvers and number of particles.

• Real-time soft tissue cutting and suturing simulation within Unity, as shown in Fig. 1.

• Teleoperation using the Master Tool Manipulator (MTM) from the da Vinci Research Kit (dVRK) [11].

To the best of our knowledge, this is also the first work to simulate soft tissue suturing using the MPM. Our method couples 1D suturing geometries (i.e., an arc-shaped needle with segmented suture thread) with soft tissue dynamics, enabling two-way coupling between MPM soft bodies and 1D rigid bodies. This formulation also lays the groundwork for incorporating other solvers, such as position-based dynamics (PBD), for future extensions. The use of MPM further allows for easy future integration of fluids, such as blood and smoke in surgeries, as MPM naturally supports fluid simulation [10].

CRESSim-MPM is part of the CRESSim project [12, 6].

II RELATED WORK

II-A Surgical Simulators and Surgical Robotics

Surgical simulator development, particularly in the area of software-based (virtual) simulators, has mainly been influenced by two streams. The first one originates from the need for surgical training and skills assessment within a simulated environment, with the users being surgeons, residents, and medical researchers. This category primarily includes virtual simulators and hybrid simulators that are capable of simulating a wide range of surgical training tasks, such as cutting, suturing, and cauterization. Examples of these are the commercial ones, such as da Vinci SimNow and the open-source ones such as iMSTK. However, they are less used by robotics researchers due to their proprietary nature, the lack of real-time performance, or difficulty in integrating with existing software, such as the Robot Operating System (ROS).

The second stream arises from the recent trend in surgical robotics research, particularly for surgical automation and autonomy. Studies including AMBF [13], ATAR [14], and V-Rep Simulator for the dVRK [15] provide diverse applications for surgical robotics research, such as for benchmarking, synthesizing data, and trial-and-error learning. Other works such as ORBIT-Surgical [16], Surgical Gym [17], FF-SRL [18], LapGym [19], AMBF-RL [20], SurRoL [21], dVRL [22], and UnityFlexML [1] has specific focus on machine learning applications. These platforms typically rely on an existing physics engine, such as PhysX (used in Nvidia Omniverse), FleX, SOFA, and Bullet. Therefore, their simulation capability largely depends on the physics solvers. For instance, Bullet is primarily a rigid body engine and only provides limited FEM soft body support. Both PhysX 5 and SOFA use FEM for soft body simulation. While PhysX 5 supports FEM, it has yet to achieve realistic interactions such as cutting and suturing. SOFA does support cutting but incurs high computational costs. These challenges may stem from inherent limitations of FEM itself, as discussed earlier. A recent study [23] integrates basic MPM support into SurRoL using MLS-MPM but does not address cutting and suturing contact methods. In contrast, our work focuses on developing a comprehensive engine-level MPM library specifically designed for simulating surgical soft-body manipulations.

II-B Real-Time Soft Body Simulation

The commonly used soft body simulation methods for real-time applications include FEM, MPM, and PBD. As discussed in Section I, FEM has several disadvantages when used for large-scale simulations in robotic research, including computational expenses, the need for real-time re-meshing for topological change, and inaccuracies and instabilities under large deformation.

PBD [24] and extended PBD (XPBD) [25] are efficient and numerically stable alternatives, which solve only positional constraints without relying on force integration, making the simulation highly stable. However, this also means less physically accurate behavior, and the realism of deformations depends on the formulation of constraints rather than the underlying physics. Although this is acceptable for computer animation, the application in biomedical and robotics research poses additional challenges, as careful parameter tuning is needed to match the simulation with real-world tissue behavior, as in [26]. Additionally, forming positional constraints can sometimes be challenging for certain types of materials, such as viscoelastic or elastoplastic models, which are common in biomedical applications.

MPM overcomes many disadvantages of FEM while preserving accuracy and versatility for simulating a wide range of material behaviors. MPM is considered within the particle-in-cell (PIC) methods family, where particles carry physical quantities and a background grid is used for computations. The initial PIC method [27] was developed for fluid mechanics, with later improvements by the fluid implicit particle method (FLIP) [28]. MPM extends the methods to simulating solids and multiphase materials, including elastic and hyperelastic materials, granular materials (e.g., snow–used by Disney in their film Frozen [29]–and sand), and elastoplastic materials such as biological tissues. In addition, since MPM does not explicitly use a mesh, topological changes are not a concern, making it easier to simulate cutting and fracture. The MPM and the PIC methods are related and have numerous variations such as generalized interpolation material point (GIMP) [30], convected particle domain interpolation (CPDI) [31], affine particle-in-cell (APIC) [32], moving least squares MPM (MLS-MPM) [33], and more recently position-based MPM (PB-MPM) [34]. However, MPM is less investigated and adopted for surgical simulation, possibly due to the lack of a high-performance library specifically built for surgical applications.

III MATERIAL POINT METHOD

III-A Standard MPM

MPM is an approach without explicit consideration of a computational mesh as in FEM. Instead, a background Eulerian (i.e., fixed) grid is used for force and momentum computations, whereas a continuum body is discretized into Lagrangian material points (i.e., moving particles) to store local material data, such as mass, velocity, stress tensor, and deformation gradient.

Figure 2: Overview of the standard MPM formulation, adapted from [10].

The standard formulation of the MPM follows 4 steps: (1) Particle to grid (P2G): Particle properties are transferred and distributed to the background grid using an interpolation shape function. (2) Grid update: the momentum of each grid node is updated using numerical integration. (3) Grid to particle (G2P): The updated node momenta are interpolated back to the particles. Particle properties such as the deformation gradient and stress are also updated. (4) Reset grid: clear grid node data to prepare for the next step. This is illustrated in Fig. 2. As an example, consider a Neo-Hookean hyperelastic model using a PIC-style particle velocity update. At time $t$, particle data is transferred to the grid using

$$m_{I}^{t}=\sum_{i}\phi(\mathbf{x}_{i}^{t})m_{i},\;(m\mathbf{v})_{I}^{t}=\sum_{i}\phi(\mathbf{x}_{i}^{t})(m\mathbf{v})_{i}^{t}.$$ (1)

Here, we use the subscript $i$ to denote the particle indices and $I$ for grid indices. $\mathbf{x}$ is the position, $m$ is the mass, $\mathbf{v}$ is the velocity, and $\phi$ is the shape function. As the shape function uses a basis function $N(x)$ with finite support, only grid nodes within this region are affected by a specific particle. With a quadratic B-spline basis, a particle can transfer data to 3 neighboring nodes per dimension. Force is transferred from the particle to the grid using

$$\mathbf{f}_{I}^{int,\,t}=-\sum_{i}V_{i}^{t}\boldsymbol{\sigma}_{i}^{t}\nabla\phi(\mathbf{x}_{i}^{t}),\;\mathbf{f}_{I}^{ext,\,t}=\sum_{i}\phi(\mathbf{x}_{i}^{t})m_{i}\mathbf{b}(\mathbf{x}_{i}^{t}),$$ (2)

where $\mathbf{f}_{I}^{int}$ is the internal force and $\mathbf{f}_{I}^{ext}$ is the external force. $V_{i}$ is the particle volume, $\boldsymbol{\sigma}_{i}$ is the stress tensor, $\mathbf{b}$ is the body force. $\mathbf{b}=\mathbf{g}$ when only gravity is applied. $\nabla\phi$ is the gradient of the shape function.

Next, the grid is updated for each node with

$$(m\mathbf{v})_{I}^{t+\Delta t}=(m\mathbf{v})_{I}^{t}+\mathbf{f}_{I}^{t}\Delta t,\;\mathbf{v}_{I}^{t+\Delta t}=(m\mathbf{v})_{I}^{t+\Delta t}/m_{I}^{t},$$ (3)

where $\mathbf{f}_{I}=\mathbf{f}_{I}^{int}+\mathbf{f}_{I}^{ext}$ is the total force. Boundary conditions (e.g., fixed nodes) can be applied by altering the momentum.

Updated nodal velocities are then transferred back to particles for integration, with

$$\mathbf{v}_{i}^{t+\Delta t}=\sum_{I}\phi(\mathbf{x}_{i}^{t})\mathbf{v}_{I}^{t+\Delta t},\;\mathbf{x}_{i}^{t+\Delta t}=\mathbf{x}_{i}^{t}+\mathbf{v}_{i}^{t+\Delta t}\Delta t.$$ (4)

Finally, particle velocity gradient $\mathbf{L}$, deformation gradient $\mathbf{F}$, volume, and stress are updated using

$$\mathbf{L}_{i}^{t+\Delta t} =\sum_{I}\nabla\phi(\mathbf{x}_{i}^{t})\mathbf{v}_{I}^{t+\Delta t},$$ (5)

$$\mathbf{F}_{i}^{t+\Delta t} =(\mathbf{I}+\mathbf{L}_{i}^{t+\Delta t}\Delta t)\mathbf{F}_{i}^{t},$$ (6)

$$V_{i}^{t+\Delta t} =JV_{i}^{t},\;\text{with }J=\det(\mathbf{F}_{i}^{t+\Delta t}),$$ (7)

$$\boldsymbol{\sigma}_{i}^{t+\Delta t} =\frac{1}{J}\left[\mu(\mathbf{F}\mathbf{F}^{T}-\mathbf{I})+\lambda\ln(J)\mathbf{I}\right].$$ (8)

$\mathbf{I}$ denotes the identity matrix. Equation (8) is specific to the Neo-Hookean model for calculating the Cauchy stress, where $\mu$ and $\lambda$ are the Lamé parameters. Some implementations use the Kirchhoff stress instead.

One benefit of MPM is that contact between multiple soft bodies is inherently handled, as each body is represented by a cluster of particles interacting through the shared background grid. This eliminates the need for explicit collision detection.

Since more advanced algorithms are now commonly used, standard MPM is implemented in this work only as a reference with limited feature support.

III-B Moving Least Squares MPM

The moving least squares MPM (MLS-MPM) [33] uses APIC-style transfers and unifies it with force computation by using a moving least squares shape function. This leads to modifications of the velocity and force transfers in the standard MPM formulation. For the internal force,

$$\mathbf{f}_{I}^{int,\,t}=-\sum_{i}V_{i}^{t}M_{i}^{-1}\boldsymbol{\sigma}_{i}^{t}N_{I}(\mathbf{x}_{i}^{t})(\mathbf{x}_{I}-\mathbf{x}_{i}^{t}),$$ (9)

where $N_{I}(x)$ is the same shape basis function used in standard MPM, centered at node $I$. $M_{i}=\frac{1}{4}{\Delta x}^{2}$ for quadratic basis and grid cell size $\Delta x$. APIC affine momentum $\mathbf{C}_{i}^{t}$ and force transfer is then fused in P2G, with

$$(m\mathbf{v})_{I}^{t}=\sum_{i}N_{I}(\mathbf{x}_{i}^{t})\left(m_{i}\mathbf{C}_{i}^{t}-\Delta tV_{i}^{t}M_{i}^{-1}\boldsymbol{\sigma}_{i}^{t}\right)(\mathbf{x}_{I}-\mathbf{x}_{i}^{t}),$$ (10)

where $\mathbf{C}_{i}^{t}$ can be updated in the APIC way during G2P using

$$\mathbf{B}_{i}^{t+\Delta t}=\sum_{i}N_{I}(\mathbf{x}_{i}^{t})v_{i}^{t+\Delta t}(\mathbf{x}_{I}-\mathbf{x}_{i}^{t})^{T},\;\mathbf{C}_{i}^{t+\Delta t}=\mathbf{B}_{i}^{t}\mathbf{D}_{i}^{t}.$$ (11)

For quadratic basis, $\mathbf{D}_{i}^{t}\equiv\frac{1}{4}{\Delta x}^{2}\mathbf{I}$, resulting in a scaling on $\mathbf{B}_{i}^{t}$. The deformation gradient is then updated with

$$\mathbf{F}_{i}^{t+\Delta t}=(\mathbf{I}+\mathbf{C}_{i}^{t+\Delta t}\Delta t)\mathbf{F}_{i}^{t}.$$ (12)

Algorithm 1 shows the structure of APIC-style MLS-MPM.

Algorithm 1 APIC-Style MLS-MPM Simulation Step

1:  for each particle $i$ do

2:     P2G mass transfer using the first of (1).

3:     P2G fused momentum-force transfer using (10).

4:  end for

5:  for each grid node $I$ do

6:     Update nodal velocity using (3) with external force.

7:     Modify nodal velocity based on contact and boundary conditions.

8:  end for

9:  for each particle $i$ do

10:     Affine momentum transfer using (11).

11:     Update particle positions using (4).

12:     Update deformation gradient using (12).

13:     Compute stress using the constitutive law.

14:  end for

15:  Reset grid: Clear grid momentum $(m\mathbf{v})_{I}$.

16:  $t\leftarrow t+\Delta t$

III-C Position-Based MPM

PB-MPM, proposed in 2024 [34], is a new method that formulates MPM as a constraint-solving problem similar to PBD, allowing larger integration steps and better numerical stability. The modifications needed to convert any MPM algorithm to PB-MPM are minimal: at each simulation step, multiple P2G-G2P cycles are conducted to iteratively solve a candidate particle velocity gradient. In P2G, instead of calculating the stress, a material constraint is solved for all particles independently. For example, for an elastic material with MLS-MPM formulation, the co-rotational constraint-solving step includes

$$\mathbf{F}_{i}^{t+\Delta t,*} \leftarrow(\mathbf{I}+\mathbf{C}_{i}^{t+\Delta t}\Delta t)\mathbf{F}_{i}^{t},$$ (13)

$$\mathbf{R},\mathbf{U} =\text{PolarDecomp}(\mathbf{F}_{i}^{t+\Delta t,*}),$$ (14)

$$\mathbf{F}_{i}^{t+\Delta t,*} \leftarrow\beta\mathbf{R}+(1-\beta)\frac{\mathbf{F}_{i}^{t+\Delta t,*}}{\det(\mathbf{F}_{i}^{t+\Delta t,*})},$$ (15)

$$\mathbf{C}_{i}^{t+\Delta t} \leftarrow\frac{1}{\Delta t}\left(\mathbf{F}_{i}^{t+\Delta t,*}(\mathbf{F}_{i}^{t})^{-1}-\mathbf{I}\right),$$ (16)

which formulates a Jacobi-style constraint projection solver.

PB-MPM integrates the benefits of both PBD and MPM but also inherits certain drawbacks from PBD, such as reliance on non-physical parameters and stiffness that depends on the number of iterations. Therefore, its application should be considered based on specific use cases.

III-D Computational Speed Comparison

CRESSim-MPM currently implements the aforementioned three types of MPM solvers in 3D for users to select from. More solver algorithms can be added in a modular manner.

We provide an initial performance evaluation when simulating a cubic object falling under gravity with different numbers of particles for the implemented solvers on Nvidia RTX 4090 (compute capability 8.9), with CUDA Toolkit version 12.4.1. For all solvers, each step advances the simulation by $0.02$ second. For MLS-MPM and standard MPM, each step contains 10 integration sub-steps with $\Delta t=0.002$. For PB-MPM, each step uses $\Delta t=0.02$ with 10 or 20 iterations. The grid size is fixed in all situations with $175,616$ nodes. The results are shown in Fig. 3. While comparing with existing implementations is reasonable, this work does not aim to outperform existing code in computational speed, as our main focus is on including multiple MPM solver options within one library. Further code optimization and performance benchmarking are planned as future work.

Figure 3: Simulation speed with different numbers of particles.

IV RIGID GEOMETRIES AND CONTACT MODELS FOR SURGICAL SIMULATION

IV-A Soft-Rigid Coupling

Figure 4: Grid level collision with SDF. Nodal velocities are modified if they are inside a rigid body.

Two-way coupling between MPM soft bodies and externally-solved rigid bodies can be achieved at the grid level [29, 33]. The most straightforward approach is to use a signed distance field (SDF) for each rigid body and modify the predicted nodal velocity based on the contact model if a grid node is inside the rigid geometry. This is similar to a predictor-corrector method [10], as shown in Fig. 4. We set the nodal velocity to zero in the normal direction, and apply frictional correction to the tangential component, unless the node is separating from the rigid body. Mathematically, the relative velocity between a node and a rigid body is $\mathbf{v}_{rel}=\mathbf{v}_{I}-\mathbf{v}_{r}$, with $\mathbf{v}_{r}$ being the rigid body’s velocity. Given a rigid surface normal $\mathbf{n}$, $\mathbf{v}_{rel}$ is decomposed into the tangent component $\mathbf{v}_{tg}$ and the normal part $v_{n}$:

$$\mathbf{v}_{rel}=\mathbf{v}_{tg}+\mathbf{n}v_{n}.$$ (17)

If the node is separating from the rigid body, the velocity is not changed, i.e., $\mathbf{v}^{\prime}_{rel}=\mathbf{v}_{rel}$. Otherwise,

$$\mathbf{v}^{\prime}_{rel}=c_{d}\max\left(1-\frac{\mu_{k}v_{n}}{\lVert\mathbf{v}_{tg}\rVert},0\right)\mathbf{v}_{tg},$$ (18)

where $\mu_{k}$ is the kinetic friction coefficient and $c_{d}$ applies a linear drag force on the velocity, which is useful for sticky effects such as tissue sticking to the scalpel and needle.

The velocity change $\mathbf{v}^{\prime}_{rel}-\mathbf{v}_{rel}$ can then be recorded and used as a momentum impulse applied back to the rigid body, yielding a two-way force coupling. Additional particle correction is performed to push out the remaining particles inside the rigid geometry, similar to [29]. The SDF-based approach works well for MPM-rigid collisions with regular geometries, as shown in Fig. 5.

Figure 5: Soft body (blue) coupling with rigid bodies (white sphere and cube) of different masses.

IV-B Geometries and Contact Method for Cutting

For cutting, authors of [33] propose using a colored distance field (CDF). The CDF stores points that contain an unsigned distance and a color information, which encodes a set of nearby surfaces and the side it is located at. If a particle and a node are located on different sides of a surface, the node is considered incompatible, and no P2G and G2P transfer occurs. This forms a discontinuity that blocks energy transfer at the cutting surface.

While this approach allows cutting with a random surface geometry, it requires the additional storage of nearby geometries and the use of rigid particles. This work simplifies the approach and uses a regular SDF representation for all rigid geometries. For a cutting geometry, the sign of the SDF value represents the side of the point relative to the nearest surface. This eliminates the need to store additional data for cutting geometries and particles and unifies the grid-level contact resolution by simply using the absolute value of the SDF for contact. However, as the absolute value is always positive, a fattened collision region is needed during contact resolution, as opposed to using CDF.

Additionally, since a knife or scalpel can cut with only the edge side instead of the spine (blunt) side, it is inappropriate to have a geometry that cuts in all directions. Considering this, we apply a circular-shaped distance field to the spine side and use sticky contact (zero relative velocity) if a node is in this region. When the spine side touches an MPM soft body, it generates contact similar to a cylinder colliding with the soft body, as shown in Fig. 6. The SDF of a scalpel’s cross-section is illustrated in Fig. 7. The geometry is named slicer in our code implementation, which includes QuadSlicer and TriangleMeshSlicer.

(a)

(b)

Figure 6: Contact between a soft body and a slicer geometry on (a) the spine side and (b) the edge side.

Figure 7: SDF for a slicer geometry. The signs of the distance on both sides are different.

IV-C Geometries and Contact Method for Suturing

In this work, we assume that the suture thread can be discretized into connected straight-line segments. To achieve suturing, we introduce the arc and the line segment geometries, which can model an arc-shaped suture needle and a discrete thread segment. They are both 1D geometries with SDFs. The signed distance is calculated analytically from a given query point to its nearest projected point on the geometry, which is always positive. The surface normal $\mathbf{n}$ is the normalized vector from the projected point to the query point. Nodal velocities inside the SDF region are corrected differently than for regular geometries. The relative velocity is decomposed into three parts: one along the arc or line tangent $\mathbf{v}_{tg1}$, one that lies on the cross-section of the geometry $\mathbf{v}_{tg2}$, and one along the surface normal $\mathbf{n}v_{n}$, as shown in Fig. 8, and only the tangent component is kept. Mathematically, the relative velocity $\mathbf{v}_{rel}$ is decomposed to

$$\mathbf{v}_{rel}=\mathbf{v}_{tg1}+\mathbf{v}_{tg2}+\mathbf{n}v_{n},$$ (19)

where $\mathbf{v}_{tg1}$ is the true tangent velocity along the arc or line. Nodal velocity correction is applied by replacing $\mathbf{v}_{tg}$ with $\mathbf{v}_{tg1}$ in (18). This ensures that the rigid geometry can move inside a soft body if the motion is along the geometry (e.g., inserting a needle) while preventing lateral motion, as shown in Fig. 9. Geometries for suturing include Arc and ConnectedLineSegments, where the latter is a generalized geometry of a line segment.

Figure 8: Nodal velocity decomposition for an arc geometry. The velocity is projected on the tangent direction of the arc and then corrected with friction and drag force.

(a)

(b)

Figure 9: 1D geometries can move freely along their tangent directions in a soft body, but deform the soft body when moving in the lateral direction.

Various external methods, such as the kinematic chain of rigid links with joints [35], mass-spring model [36], and PBD [37], can be used to simulate the suture thread itself, as long as the coupling momentum from the soft body can be applied back to each small segment.

V IMPLEMENTATION

V-A Object Model and Simulation Lifecycle

The library follows an object model similar to the one in PhysX, such that integration with other PhysX-based software is straightforward, as shown in Fig. 10.

Figure 10: CRESSim-MPM object model.

Scene contains all simulation objects and provides the method advance(float dt) for starting to advance the simulation by a time step dt, and the method fetchResults() retrieves the results of the simulation, blocking the thread until the results are available. This is similar to PhysX’s simulation lifecycle which splits the simulation into two steps, allowing asynchronous tasks to run while the solver advances.

V-B Implementation Details and Usage

The library is implemented in CUDA C++ and is dependent only on the standard C++ library and the CUDA runtime. 3 types of Solver discussed in Section III are currently included. ParticleObject represents an object of material points simulated by MPM, and Shape is a rigid body. SimulationFactory is used to create all objects. C++ Users should only include the header simulation_factory.h. A backend C-style API engine is also provided for integration into existing software.

While we also provide a CPU version for each solver, they are not optimized for multi-core parallelism or single instruction, multiple data (SIMD). Users who require simulations that fully utilize the CPU should refer to alternative software, such as Taichi-MPM for MLS-MPM, and MPM-GIMP for standard MPM and GIMP.

V-C Unity Integration

Integration into Unity is straightforward using the C-style APIs. Only a few modifications are needed to allow the coupling between the MPM soft body and rigid bodies simulated in Unity’s default engine. For each rigid body in Unity, a corresponding Shape is created for two-way coupling with the soft body. During each step, the rigid body is simulated first with Unity, after which the Shape’s velocity is updated. The MPM engine is then advanced, recording a total momentum change of the contact grids with the rigid body. This is then used to apply a momentum impulse on the rigid body back in Unity. A Unity plugin is provided as part of CRESSim-MPM, as shown in Fig. 11a.

(a)

(b)

Figure 11: Usage of the library: (a) Unity integration; (b) Teleoperating the suturing scene using the dVRK MTM.

VI SURGICAL SIMULATION EXAMPLES

This section includes screenshots from two example scenes representing surgical cutting and suturing. While the visuals and rendering are not highly realistic, the examples are primarily designed to showcase the physics capabilities. For demonstration purposes, we use two different solvers, although either one works in both scenarios. Additional examples are provided in the supplementary video.

VI-A Cutting (PB-MPM)

This example demonstrates that we can simulate soft tissue cutting using a scalpel. The scalpel moves right to cut the tissue first and then down to split and reveal both sides of the cut tissue. The screenshots are shown in Fig. 12.

Figure 12: Soft tissue cutting with a scalpel moving right then down.

Figure 13: Suturing between two soft body blocks.

VI-B Suturing (MLS-MPM)

This example demonstrates the capability of simulating soft tissue suturing using an arc-shaped needle. For the sake of simplicity, we use a kinematic chain of rigid links with spherical joints to represent a suture thread, which is implemented using Unity’s default physics engine. This is a traditional approach that does not provide an adequate level of realism [35, 38]. Nevertheless, the focus of this work is not the thread simulation itself. Integration of more advanced approaches, such as PBD-based threads [37], is straightforward since contact force coupling can be achieved in a similar manner. Screenshots of this example are shown in Fig. 13. This example can be teleoperated using the MTM robot from the dVRK, as shown in Fig. 11b.

VII LIMITATIONS AND FUTURE WORK

We are aware of several areas where the library can be enhanced, and we’re actively working on improving these before an official release. These include (1) code optimization, particularly through increased use of fused multiply-add (FMA) operations, shared memory for rigid shapes, and improved thread block scheduling; (2) implementation of additional material models, such as elastoplastic materials and fluids, to simulate specific tissue types, blood, and smoke; (3) integration of other MPM and PIC solvers, including APIC and FLIP, as well as PBD solvers for fascia and suture thread simulation; and (4) improved rendering of particle objects.

There is also a lack of performance comparison with existing general-purpose MPM implementations, such as Taichi-MPM which does not support features such as suturing and only uses MLS-MPM. This is left as a future work after code optimization. We also recognize that relying on CUDA restricts broader compatibility with computing devices beyond Nvidia GPUs. Extending to other platforms can require a significant amount of future work, although porting this work to AMD devices should be straightforward using ROCm/HIP.

VIII CONCLUSION

In this work, we introduced CRESSim-MPM, a new MPM library designed for surgical simulation that implements multiple solvers and surgery-specific contact models for cutting and suturing. By using it as a separate physics module, integration with existing surgical simulation platforms is straightforward, and we further provide a Unity plugin to allow direct integration into Unity. With the capability of simulating cutting and suturing on soft tissue efficiently on the GPU, this work enables more realistic and efficient surgical simulations that can be potentially used in various applications, including soft tissue modeling, surgical skills training, and surgical robot learning.