Efficient Simulation of Volumetric Deformable Objects in Unity3D: GPU-Accelerated Position-Based Dynamics

Abstract

This paper proposes an efficient approach for simulating volumetric deformable objects using the Position-Based Dynamics (PBD) method. Volumetric bodies generated by TetGen are used to represent three-dimensional objects, which accurately capture complex shapes and volumes. However, when a large number of constraints are applied to the system to solve using serialized algorithms on central processing units (CPU), the computational cost can become a bottleneck of the simulation. To address this issue, the proposed implementation algorithm takes advantage of graphic processing unit (GPU) acceleration and parallel processing to improve the efficiency of the simulation. We propose two specific contributions: firstly, the use of the PBD method with volume constraint for tetrahedral elements to simulate volumetric deformable objects realistically; secondly, an efficient GPU-accelerated algorithm for implementing the PBD method that significantly improves computational efficiency. We also applied the node-centric and constraint-centric algorithms to solve the stretch constraint in the GPU-based algorithm. The implementation was performed using Unity3D. The compute shader feature of Unity3D was utilized to perform thousands of parallel computations in a single pass, making it possible to simulate large and complex objects in real-time. The performance of the simulation can be accelerated by using GPU-based methods with stretch and bending constraints, which provides significant speedup factors compared to using only the CPU for deformable objects such as Bunny, Armadillo, and Dragon. The constraint-centric and node-centric GPU approaches provide speedup factors of up to 8.9x and 8x, respectively, while the GPU-based methods with all types of constraints exhibit a slight decrease but still operate at real-time speeds. Overall, this approach enables the simulation of complex and irregular shapes with plausible and realistic results, while also achieving speed, robustness, and flexibility. Additionally, the proposed approach can be applied to general simulation and other game engines that support GPU-based acceleration.

1. Introduction

The concept of the metaverse has grown in popularity recently since it represents a shared virtual environment where people can interact and engage in fully immersive experiences with one another as well as digital objects [1]. The Metaverse is a rapidly evolving concept that merges the physical and virtual worlds through the use of augmented reality (AR) and virtual reality (VR) technologies [2]. It provides a three-dimensional (3D) virtual universe where people can interact with each other and online resources, bridging the connectivity of social media with immersive virtual environments. The visual and interactive elements that bring the virtual world to life are provided by computer graphics, which are a crucial part of the metaverse [3]. As the metaverse develops and grows, computer graphics will continue to play a crucial role in constructing immersive and captivating virtual environments that inspire curiosity and offer limitless opportunities for exploration and future potential research.

Graphical simulation is an essential element of computer graphics since it allows for the creation of realistic and immersive virtual environments. Therefore, it is also a critical tool for industries that rely on virtual prototyping and testing before building products, such as automotive, aerospace, manufacturing, and medical industries [4,5,6,7,8]. Simulations help save time and costs, and also reduce the risk of errors and accidents, thus improving the safety of the products. In industries such as animation and the visual industry, simulation is used to create realistic motion and deformation of characters and objects. Similarly, in film and advertising, graphic simulations are essential for creating high-quality visual effects and animations that are difficult to achieve real-time performance with traditional techniques [9,10]. However, these simulations are often computationally expensive, requiring a significant amount of computing resources and time to produce realistic and accurate results [11].

Deformable objects, such as fluids, cloth, and soft bodies, present several challenges due to their complex physical properties. These objects can exhibit non-linear and time-varying behavior, making them difficult to model accurately using traditional methods. Additionally, deformable objects can interact with other objects in the scene, such as rigid bodies and obstacles, leading to even more complex interactions. It is important to consider the material properties of the object, such as elasticity, plasticity, viscosity, and collision response, as they have a significant impact on how deformable objects behave and interact similarly to their real-world counterparts. Deformable objects also often involve a large number of particles and elements that represent the object, which requires significant computation resources [12,13].

Over time, numerous approaches have been developed by researchers to create realistic models and simulations of soft bodies that are suitable for use in interactive applications. The finite element method (FEM) and mass-spring system are two widely used approaches for simulating deformable objects in computer graphics and simulation [14,15]. FEM is a numerical method that breaks down a deformable object into a mesh of smaller elements, each of which is governed by a set of equations describing its behavior. This allows for a more accurate and detailed simulation of the object’s physical properties, such as elasticity and plasticity. On the other hand, the mass-spring system is a simple and intuitive method that models a deformable object as a collection of interconnected particles, or springs, that are subject to forces such as tension, compression, and damping. While less computationally expensive than FEM, the mass-spring system can still produce realistic and visually appealing simulations, particularly for objects with simple geometries. On the other hand, position-based dynamics (PBD) is a relatively new approach to simulating deformable objects that has gained popularity in recent years [16]. PBD is a physics-based method that relies on the imposition of positional constraints on the particles or elements in a simulation rather than solving the underlying physical equations of motion. This makes PBD easier to implement and more computationally efficient than traditional methods. Additionally, PBD offers stability advantages over other approaches, such as avoiding unwanted numerical instabilities that can lead to a loss of visual realism. These advantages have made PBD a promising area of research in computer graphics and simulation, and it has been successfully applied to a range of problems, including cloth simulation, fluid simulation, and soft-body dynamics [17]. As the use of multicore processors becomes more widespread, it is expected that the performance and capabilities of modern computing systems will continue to improve, enabling new and innovative applications and solutions [18,19,20].

In the original approach, PBD consists of a set of constraints that need to be solved using common methods such as Jacobi and Gauss–Seidel solvers [21]. The process of solving constraints using the Gauss–Seidel method is more favorable since a constraint can be solved sequentially, one after the other. However, when a large number of constraints are applied to the system, the computational cost, which is CPU-based, increases and quickly becomes a bottleneck for the simulation.

In contrast, several existing solutions to improve simulation performance involve the use of meshless or mesh-free approaches that primarily consider only the vertices of the mesh to represent deformable objects [22]. The computational cost of the simulation can be decreased with this approach, but accuracy and realism may be degraded in the process. To balance accuracy and efficiency, additional methods can be applied, such as adaptive mesh refinement, to produce a simulation that is more reliable and optimized [23]. Therefore, the internal structure of the deformable object has also been studied to enhance the accuracy, stability, and simulation performance in real-time [24]. In the GPU-based method, the method’s computational efficiency was enhanced by utilizing fast summation techniques on both structured grids and unstructured meshes [25,26].

Volumetric bodies, such as tetrahedral models, are used in simulations to represent the deformation of an object in three dimensions, accurately capturing the complex shapes and volumes of objects in detail [27]. Tetrahedral models are particularly useful in cases where the geometry of an object is complex or irregular. Tetrahedral models consist of a set of interconnected tetrahedral elements, each of which is defined by a set of vertices and edges. These tetrahedra can be used to create a volumetric mesh that can be deformed under the influence of external forces or constraints. Therefore, the tetrahedral element can be enforced by volume constraints, which makes it well-suited for simulations using the PBD concept to obtain speed, robustness, flexibility, and a plausible result.

Therefore, in this paper, we implemented a method for deformable object simulation that uses volumetric bodies generated by TetGen to represent 3D objects. Since the PBD method is a stability method for large time-step conditions, we proposed an implementation algorithm that leverages GPU acceleration and parallel processing to improve the efficiency of the simulation. The proposed design and implementation of the simulation are performed in Unity3D, which is a widely used game engine that provides a powerful and flexible environment for developing real-time simulations. Our implementation uses the compute shader feature of Unity3D to utilize the parallel processing capabilities of modern GPUs. The compute shader allows us to perform thousands of parallel computations in a single pass, making it possible to simulate large and complex objects in real-time. Not only can the proposed approach be applied to Unity3D, but it can also be used with other game engines that support GPU-based acceleration. Although the process for other game engines may differ slightly, the fundamental concept of compute shaders is likely to be similar to other GPGPU libraries. The specific contributions of this paper are:

• The volume constraint for a tetrahedral element based on the PBD method is employed to achieve a plausible and realistic simulation of volumetric deformable objects.
• We propose an efficient algorithm for the implementation of the PBD method using GPU acceleration and parallel processing, which significantly improves the computational efficiency of the simulation.

The rest of the paper is organized as follows: in Section 2, we provide a review of related work in the field of deformable object simulation, focusing on previous approaches for simulating volumetric objects and their limitations. In Section 3, we provide a detailed overview of the position-based dynamics method, including its basic principles and advantages. In Section 4, we describe the design and implementation of our method for simulating volumetric deformable objects using PBD and GPU acceleration. In Section 5, we present the results of our experiments and discuss the limitations of our approach. Finally, we conclude the paper in Section 6, highlighting the potential applications of our method in various domains and suggesting areas for future research.

2. Related Work

A series of methods have been developed over time to address the complex challenge of simulating deformable objects in computer graphics. Zhang et al. [28] presented a comprehensive review of state-of-the-art deformable models for soft tissue deformation in interactive surgical simulations, providing an extensive analysis of the various modeling approaches, internal forces, numerical time integrations, and material properties. For the purpose of simulating deformable objects, a variety of techniques are available, each with specific pros and cons.

The most widely used methods in computer graphics for simulating dynamic systems are force-based methods. Force-based approaches, such as the finite element method (FEM) and mass-spring systems, rely on computing the forces acting on a system and then integrating those forces to update the positions and velocities of the objects being simulated [29]. While these methods can produce highly accurate simulations, they can be computationally expensive, especially for large or complex systems. PBD has been shown to be highly efficient and capable of handling large deformations and complex geometries, but it may not be as accurate in all cases and has its limitations [30]. Constraints, such as collisions, are also an important consideration, as different methods may handle them differently, and some methods may not be suitable for certain types of constraints. With the increasing computational power of modern GPUs, PBD can be efficiently implemented on parallel hardware to achieve high performance. A comprehensive survey of PBD was presented by Bender et al. [16].

Early research by Muller et al. [21] introduced the basic concept of PBD and presented different methods to solve the constraints in the PBD style. This approach allows for the creation of realistic and stable simulations of deformable objects, including cloth, fluids, and soft bodies. Since then, PBD has been widely adopted in the computer graphics industry and has been used in a variety of applications, including character animation, game physics, and virtual reality simulations. Over the years, researchers and developers have continued to refine and improve the PBD technique. For example, Macklin et al. [31] presented a unified dynamic framework for simulating multiple deformable objects in real-time. The same author also presented XPBD, which is a straightforward modification to the PBD method that enables it to simulate any type of elastic or dissipative energy potential in a precise and efficient manner using an implicit approach [32].

To improve performance, many approaches have been explored and implemented using parallel processing to partition the process of numerical solving. Fratarcangeli et al. [33] presented a massively parallel implementation of PBD for fast and interactive animation of deformable bodies. The implementation divides the set of constraints into independent partitions using a graph algorithm, solving each partition in parallel on the GPU for a significant performance improvement compared to its serial counterpart. Similarly, a highly parallelized version of PBD was demonstrated by Diziol et al. [34] that utilized a single GPU device. Lee et al. [24] introduced the internal shape-preserving constraint algorithm, which is designed to preserve the shape of a 3D deformable object and reduce the number of internal structures when integrated into PBD. The study also compared the performance of this approach to the traditional tetrahedral model. Rumman et al. [30] extended the application of PBD to simulate soft body characters by utilizing the tetrahedral volume constraint on a volumetric mesh and then reconstructing the surface information for rendering purposes.

A new algorithm called MSCOEA was also developed to solve multi-objective optimization problems while balancing uniformity and convergence. It was found to be effective, and it showed improved performance compared to other algorithms [35]. Xie et al. [36] presented the CFD-DEM coupling method to investigate the self-propulsion performance of an ice-strengthened Panamax bulk carrier in a brash ice channel. The results showed that this technique can provide a technical tool for hull line development of ice-strengthened polar ships and ship navigation guidance in brash ice channels. Recent studies have demonstrated the potential of integrating machine learning techniques with PBD to improve collision detection, material property estimation, and animation generation, among other capabilities [37,38]. This research could be relevant to studies on GPU-based, position-based dynamics algorithms for simulating deformable objects.

3. Position-Based Dynamic Algorithm Overview

3.1. The Principle of PBD

In deformable bodies, an object is represented by a set of $N$ nodes (vertices) and $M$ constraints. A node $i\in \left[1,\dots ,N\right]$ consist of position $x_i\in {ℝ}^3$, velocity $v_i\in {ℝ}^3$, and mass $m_i\in ℝ$. However, the inverted mass, which is denoted by $w_i\in ℝ$, is used instead to simplify the integration step. Similarly, an equality or inequality constraint $j\in \left[1,\dots ,M\right]$ is made using a cardinality $n_j$ and has different functions depending on the type of constraint ($C_j:{ℝ}^{3{\mathrm{n}}_{\mathrm{j}}}\to ℝ)$, and also uses coefficients to control the stiffness denoted by $k_j\in \left[0\dots 1\right]$. The result of the constraint solving is an estimation of the new position on the node, which satisfies all constraints conditions and is denoted by $p_i\in {ℝ}^3$. In other words, $p_i$ is the projected position, the position that satisfies all constraint functions $C_j$. The general algorithm of PBD is given in Algorithm 1.

Algorithm 1 Position-based dynamics algorithm

Required:

Position: $x$   Velocity: $v$   Inverse mass: $w$   External force: $f_{ext}$   Time step: $\Delta t$ 1.  for all nodes $\mathit{i}$ do: 2.  initialize ${\mathit{x}}_{\mathit{i}}={\mathit{x}}_{\mathit{i}}^{\mathbf{0}}$, ${\mathit{v}}_{\mathit{i}}={\mathit{v}}_{\mathit{i}}^{\mathbf{0}}$, ${\mathit{w}}_{\mathit{i}}=\mathbf{1}/{\mathit{w}}_{\mathit{i}}$ 3.  end for 4.  simulation loop 5.  for all nodes $i$ do: 6.  ${\mathit{v}}_{\mathit{i}}\leftarrow {\mathit{v}}_{\mathit{i}}+\Delta t{w}_i{f}_{ext}\left(x_i\right)$ 7.  $p_i\leftarrow x_i+\Delta t{v}_i$ 8.  end for 9.  for all nodes $i$ do Generate Collision Constraint $\left(x_i\to p_i\right)$ 10.  loop solver iteration times: 11.  projectConstraint($C_1,\dots , C_{M+M_{Coll}},p_1,\dots ,p_N)$ 12.  end loop 13.  for all nodes $i$ do: 14.  $v_i\leftarrow (p_i-x_i)/\Delta t$ 15.  $x_i\leftarrow p_i$ 16.  end for 17.  end loop

In lines 1–3, the variables are initialized with their initial conditions. Since the algorithm performs the explicit Euler integration step, the velocities and position are updated as shown in lines 5–8. The collision constraints are generated in line 9 for the current time step, which will be solved in the next stage. The core idea of PBD is to solve the constraints iteratively is shown in lines 10–12. The iterative solver is performed to iteratively correct the predicted position to satisfy both $M$ internal and $M_{Coll}$ external constraints. Finally, the positions and velocities can be updated based on the corrected position $p_i$.

In the PBD approach, the goal of the solver in lines 10–12 is to correct the predicted position $p$. There are several approaches to solving the system of constraints. However, in this paper, we utilized a non-linear Gauss–Seidel for the CPU-based method, as each constraint is linearized individually before its projection. This makes the solver more stable than a global approach, where the linearization is kept fixed for the entire global solution of a Newton iteration. We applied the Jacobi solver style to solve the constraint problem in the GPU-based method since each constraint can be processed in parallel.

In the process of solving the constraint, we want to find a $\Delta p$ that corrects the position. The constraint function can be written such that $C\left(\mathrm{p}+\Delta \mathrm{p}\right)=0$. We can approximate the constraint function as follows:

$$C\left(\mathrm{p}+\Delta \mathrm{p}\right)\approx C\left(\mathrm{p}\right)+{\nabla }_{p}C\left(\mathrm{p}\right)·\Delta p=0$$

(1)

where the derivatives of the function C with respect to the $N$ component of the position are written as a vector ${\nabla }_{p}C\left(\mathrm{p}\right)$. A scalar value $\lambda$ is chosen to limit the direction of the change $\Delta p$ to be consistent with the gradient ${\nabla }_{p}C$, thereby ensuring their alignment as follows:

$$\Delta p=\lambda {\nabla }_{p}C\left(\mathrm{p}\right)$$

(2)

By substituting Equation (2) in Equation (1), $\lambda$ can be solved as follows:

$$\lambda =-\frac{C\left(\mathrm{p}\right)}{{\left|{\nabla }_{p}C\left(q\right)\right|}^2}$$

(3)

We also obtain the displacement of each node $i$ as follows:

$$\Delta p_i=-\frac{{\mathrm{w}}_{\mathrm{i}}C\left(\mathrm{p}\right)}{{{\displaystyle \sum }}_j{w}_j{\left|{\nabla }_{p}C\left(p\right)\right|}^2}{\nabla }_{p_i}C\left(\mathrm{p}\right)$$

(4)

However, different constraints can be used to represent the object with a variety of materials.

3.2. Stretch Constraint

The stretch constraint is a commonly used constraint in deformable body simulations which ensures that the distance between two nodes remains within a certain threshold during the simulation step [12]. The stretch constraint restores and resists the deflected nodes to their original distance d, as shown in Figure 1.

We then obtain a stretch constraint function $C_s\left(p_1,p_2\right)=\left|p_1-p_2\right|-d$, where $d$ is the rest distance of the stretch constraint. The derivatives with respect to $p_1$ and $p_2$ are ${\nabla }_{p_1}C\left(p_1,p_2\right)=\frac{p_1-p_2}{\left|p_1-p_2\right|}$ and ${\nabla }_{p_2}C\left(p_1,p_2\right)=-\frac{p_1-p_2}{\left|p_1-p_2\right|}$, respectively. The final corrections for $p_1$ and $p_2$ are shown in Equations (4) and (5).

$$\Delta p_1=-\frac{w_1}{w_1+w_2}\left(\left|p_1-p_2\right|-d\right)\frac{p_1-p_2}{\left|p_1-p_2\right|}$$

(5)

$$\Delta p_2=\frac{w_2}{w_1+w_2}\left(\left|p_1-p_2\right|-d\right)\frac{p_1-p_2}{\left|p_1-p_2\right|}$$

(6)

The coefficient to control the stiffness of the stretch constraint is denoted by $k_{strect}\in \left[0\dots 1\right]$. We then use $k_{strect}$ to multiply the corrections $\Delta p$.

3.3. Bending Constraint

In a practical simulation, the bending constraint is used to resist each pair of adjacent triangles. There are two types of bending constraints: dihedral and isometric bending. However, we only discuss the dihedral bending constraint in this section due to its simplicity. A bending constraint must have two adjacent triangles $\left(p_1,p_2,p_3\right)$ and $\left(p_1,p_2,p_4\right)$ and consist of four nodes, as illustrated in Figure 2. Since we have two pairs of triangles, two normals $n_1$ and $n_2$ create the initial angle ${\phi }_0$. Similar to the stretch constraint, we have $k_{bend}\in \left[0\dots 1\right]$ to define the bending stiffness of the object.

The function of the bending constraint can be written as follows:

$$\begin{gathered} C_b\left(p_1,p_2,p_3,p_4\right)=\arccos \left(n_1·n_2\right)-{\phi }_0 \\ =\arccos \left(\frac{\left(p_2-p_1\right)\times \left(p_3-p_1\right)}{\left|p_2-p_1\right|\times \left|p_3-p_1\right|}·\frac{\left(p_2-p_1\right)\times \left(p_4-p_1\right)}{\left|p_2-p_1\right|\times \left|p_4-p_1\right|}\right)-{\phi }_0 \end{gathered}$$

(7)

With $\frac{d}{dx}\arccos \left(x\right)=-\frac{1}{\sqrt{1-d^2}}$, the derivative with respect to the points are:

$${\nabla }_{{\mathrm{p}}_4}C=-\frac{1}{\sqrt{1-d^2}}\left({\left(\frac{\partial n_2}{\partial p_4}\right)}^T{n}_1\right)$$

(8)

$${\nabla }_{{\mathrm{p}}_3}C=-\frac{1}{\sqrt{1-d^2}}\left({\left(\frac{\partial n_1}{\partial p_3}\right)}^T{n}_2\right)$$

(9)

$${\nabla }_{{\mathrm{p}}_2}C=-\frac{1}{\sqrt{1-d^2}}\left({\left(\frac{\partial n_1}{\partial p_2}\right)}^T{n}_2+{\left(\frac{\partial n_2}{\partial p_2}\right)}^T{n}_1\right)$$

(10)

$${\nabla }_{{\mathrm{p}}_1}C=-{\nabla }_{{\mathrm{p}}_2}C-{\nabla }_{{\mathrm{p}}_3}C-{\nabla }_{{\mathrm{p}}_4}C$$

(11)

By using the gradient of the normalized cross-product, we obtain the following gradients:

$$q_4=\frac{p_2\times n_1+\left(n_2+p_2\right)d}{\left|p_2\times p_4\right|}$$

(12)

$$q_3=\frac{p_2\times n_2+\left(n_1+p_2\right)d}{\left|p_2\times p_3\right|}$$

(13)

$$q_2=-\frac{p_3\times n_2+\left(n_1+p_3\right)d}{\left|p_2\times p_3\right|}-\frac{p_4\times n_1+\left(n_2+p_4\right)d}{\left|p_2\times p_4\right|}$$

(14)

$$q_1=-q_2-q_3-q_4$$

(15)

The final correction of each node in the bending constraint can be written as follows:

$$\Delta p_i=-\frac{w_i\sqrt{1-d^2}\left(\arccos \left(d\right)-{\phi }_0\right)}{{{\displaystyle \sum }}_j{w}_j{\left|q_j\right|}^2}{\mathrm{q}}_{\mathrm{i}}$$

(16)

3.4. Volume Constraint

In most simulations of deformable objects composed of tetrahedron elements, it is straightforward to employ stretch and bending constraints. However, in cases where the object’s accuracy of deformation and volume conservation are required, volume constraints can be used to accurately capture the behavior of the deformable object. Volume constraints based on PBD are very useful for most cases such as volume maintenance, shape restoration, inverse volume deformation, etc. The tetrahedral model consists of a set of tetrahedra, and each tetrahedron is composed of 4 nodes $\left(p_1,p_2,p_3,p_4\right)$ and an initial volume $V_0$, as shown in Figure 3. The initial volume of a tetrahedron can be calculated using the parallelepiped concept and expressed as follows: $V_0=\frac{1}{6}\left[\left(x_2-x_1\right)\times \left(x_3-x_1\right)\right]·\left(x_4-x_1\right)$. Additionally, the stiffness of the volume constraint can be controlled using the parameter $k_{vol}\in \left[0\dots 1\right]$.

The function of the volume constraint can be written as follows:

$$C_v\left(p_1,p_2,p_3,p_4\right)=\frac{1}{6}\left[\left(p_2-p_1\right)\times \left(p_3-p_1\right)\right]·\left(p_4-p_1\right)-V_0$$

(17)

The gradient of each node in the tetrahedron can be written as follows:

$${\nabla }_{{\mathrm{p}}_4}C=\frac{1}{6}\left(p_1-p_2\right)\times \left(p_1-p_3\right)$$

(18)

$${\nabla }_{{\mathrm{p}}_3}C=\frac{1}{6}\left(p_1-p_2\right)\times \left(p_1-p_4\right)$$

(19)

$${\nabla }_{{\mathrm{p}}_2}C=\frac{1}{6}\left(p_1-p_3\right)\times \left(p_1-p_4\right)$$

(20)

$${\nabla }_{{\mathrm{p}}_1}C=-{\nabla }_{{\mathrm{p}}_2}C-{\nabla }_{{\mathrm{p}}_3}C-{\nabla }_{{\mathrm{p}}_4}C$$

(21)

Therefore, we obtain the correction of each node belonging to the tetrahedron, which can be written as follows:

$$\Delta p_i=-\frac{w_i{C}_v\left(p_1,p_2,p_3,p_4\right)}{{{\displaystyle \sum }}_j^4{w}_j{\left|{\nabla }_{{\mathrm{p}}_{\mathrm{j}}}C\right|}^2}{\nabla }_{{\mathrm{p}}_{\mathrm{i}}}C$$

(22)

4. Implementation Details

In this section, we present a technique for deformable object simulation using the original PBD method for the volumetric model generated by the TetGen library [39]. The implementation proposed in this paper is specifically designed for Unity3D, a well-known tool for game development and physics simulation. In addition, it is applicable to develop immersive content for superior VR/AR devices.

Although GPUs were originally designed for image rendering and animation, their architecture with many processing cores and high throughput has made them useful for accelerating other applications such as graphic modeling and simulation, general computing, and AI [40,41,42]. The parallel processing capabilities of GPUs can significantly increase the processing speed for tasks that can be divided into smaller sub-tasks. Therefore, we utilize compute shaders, which are programs that runs on the GPU outside of the normal rendering pipeline. Compute shaders acts as kernel programs with a defined number of workgroup blocks. These blocks are made up of threads that are organized into a one, two, or three-dimensional grid and executed in a non-sequential order, providing an advantage to the performance of extremely large data. Compute shaders work similarly to other general-purpose computing on graphics processing units (GPGPU) pipelines, in which a large amount of data can be processed simultaneously [43]. However, compute shaders are specifically designed to allow developers to write custom programs that can operate on these data in a highly parallelized manner. They can be used for a variety of tasks beyond graphics rendering, such as physics simulation and machine learning. This makes them a powerful tool for accelerating computationally intensive applications. Figure 4 shows the structure of a compute shader. Note that the compute shader is called in a Unity3D script (C# script).

In a single compute shader script in the Unity engine, the user can create multiple kernels. Therefore, we designed multiple kernels in order to perform Algorithm 1 purely on the GPU as a parallel algorithm.

Table 1. List of kernel programs for deformable object simulation using the PBD approach.

Kernel	Grid Size	Size of Thread Block
ExplicitEuler()	$N/1024$	$\left(1024, 1, 1\right)$
SolveStretchConstraint()	$M_{Stretch}/1024$	$\left(1024, 1, 1\right)$
SolveBendConstraint()	$M_{Bend}/1024$	$\left(1024, 1, 1\right)$
SolveVolumeConstraint()	$M_{Tet}/1024$	$\left(1024, 1, 1\right)$
AverageDelta()	$N/1024$	$\left(1024, 1, 1\right)$
UpdatePositions()	$N/1024$	$\left(1024, 1, 1\right)$

describes the kernel used to perform parallel PBD in Unity3D. In order to speed up the performance of the parallel process, we need to increase the number of threads per block. However, our machine can only utilize 1024 threads per block. In addition, each kernel is called in a sequential order in the C# script for every frame of the simulation, and we must provide the grid size information as well.

We also require the GPU buffer in order to perform the numerical calculations and render the deformable object based on the given vertices data.

Table 2. List of GPU buffers used to perform the PBD method in compute shader.

Buffer	Size	Description
Positions buffer	$N$	Contain the position of vertices
Velocities buffer	$N$	Contain the velocity of vertices
Projected Position buffer	$N$	Contain the predicted position of vertices
DeltaPos Buffer	$N$	Contain the corrected position from the constraint satisfaction process
Stretch Constraint buffer	$M_{Stretch}$	Contain the stretch constraint to be solved
Bend Constraint buffer	$M_{Bend}$	Contain the bending constraint to be solved
Volume Constraint buffer	$M_{tet}$	Contain the tetrahedron volume constraint to be solved

shows the list of the buffers used in this paper.

In order to perform the simulation purely on the GPU for both calculations and rendering, we propose a flowchart to simulate the volumetric mesh using the compute shader of Unity3D. In our proposed scheme, the GPU rendering is performed, which involves the rendering pipeline and requires a custom shader (vertex, fragment, and geometry) to fully render the deformable object efficiently. Figure 5 shows the flowchart of the proposed algorithm.

Note that in our proposed scheme, we modified some processes from the original concept of PBD such as the initialization of collision constraint, velocity damping, and collision constraint solving, as shown in Algorithm 1. However, it still serves the main purpose of the simulation using the PBD method. Since we used Jacobi solver for GPU-based algorithms, once all constraints are processed, the node’s total constraint delta is divided by $n_i$: the number of constraints affecting the node. The constraint averaging can be written as follows:

$$\Delta \widetilde{p_i}=\frac{1}{n_i}{\displaystyle \sum }\Delta p_i$$

(23)

The constraint averaging ensure convergence, but it converges significantly slower than the Gauss–Seidel solver. In some cases, the averaging is too aggressive, making it difficult to control the solver’s performance. To overcome this issue, we use a global user parameter $\omega$ to control the rate of successive over-relaxation (SOR), similar to Fratarcangeli et al. [44]. The constraint average equation can be written as follows:

$$\Delta \widetilde{p_i}=\frac{\omega }{n_i}{\displaystyle \sum }\Delta p_i$$

(24)

Note that the recommended range for the global parameter used to control the rate of SOR is $1\le \omega \le 2$ in order to achieve good results. In this paper, we use $\omega =1.5$ for the GPU-based PBD solver, as it can be simulated in any model stably and achieve good convergence.

In addition to the GPU-based constraint solving in parallel, there may be a problem of data racing in the stretch constraint solving section when multiple threads attempt to invoke an element of buffer at the same time. There are many possible solutions to solve the stretch constraint in parallel, such as graph coloring and sum reduction, which could be more effective, but they may require much memory and require additional initialization processes [45].

Therefore, we only consider node-centric and constraint-centric approaches. The node-centric algorithms are designed to efficiently update the positions and velocities of nodes based on their adjacent nodes with respect to the constraint list. In addition, the node-centric approach for stretch constraint solving functions as a gathering method where each node gathers information from its neighboring nodes and updates its position based on the stretch constraint equation of PBD. Algorithm 2 shows the pseudo-code to perform node-centric operations in a compute shader. In a constraint-centric algorithm for stretch constraint in parallel, each constraint is assigned to a processing unit. The processing units are responsible for updating the corresponding node positions. This approach allows for better parallelism than node-centric algorithms as each constraint can be computed independently. However, each processing unit might require an atomic add operation to concurrently update the position with guaranteeing that multiple threads can simultaneously access and modifying the same memory location. This prevents race conditions, where multiple threads can interfere with each other’s changes to the same memory location, potentially leading to unpredictable results. Algorithm 3 shows the pseudo-code to perform constraint-centric operations in a compute shader.

Algorithm 2 Node-Centric Algorithm

1.  for all nodes i in parallel do: 2.  for all constraint c affecting node i do: 3.  Accumulate the delta in Equation (2) 4.  end for 5.  end for

Algorithm 3 Constraint-Centric Algorithm

1.  for all constraint c in parallel do: 2.  Calculate the delta in Equation (2) 3.  Atomic add operation to update the delta 4.  end for

In the implementation, we apply both techniques to simultaneously handle stretch constraints. However, in the node-centric algorithm, we compute the list of constraints that affect each node during the initializing stage to make the thread directly accessible. However, we employ constraint-centric operations for all types of constraints as they are easier and directly update the positions.

5. Results and Discussion

We have implemented a simulation of the volumetric mesh based on the PBD method in Unity3D. The specifications used for our experiment on the desktop are presented in

Table 3. Experimental environment for simulation.

Component	Specification
Operation System	Windows 10 Pro 10.0.19044 Build 19044
CPU	Intel® Core™ i7-7700 3.60 GHz
RAM	32 GB
GPU	NVIDIA GeForce RTX 2070 SUPER
	8 GB V-RAM
Shader	HLSL Shader model 5.0 (Windows Graphics API)
IDE	Unity engine version 2020.3.8f1, Microsoft Studio Community 2019 version 16.11.2, TetGen 1.6.0

.

In our proposed approach, calculations are performed by a compute shader that is dispatched in the C# script, and the vertices’ positions are stored in the GPU buffer. The custom shader is required to shade the tetrahedral object based on the vertices’ position and normal vector. The shading stage is performed in a shader that includes vertex, fragment, and geometry shaders.

In order to obtain the volumetric object used in this paper, we first obtained 3D models from the Stanford repository, which are free and available to download. We then used TetGen to convert them into volumetric meshes that include the vertices and triangles on the surface and the list of Tetrahedra.

Table 4. The 3D models used in this simulation for performance comparisons.

3D Model	Vertices	Triangles	${\mathit{M}}_{\mathit{S}\mathit{t}\mathit{r}\mathit{e}\mathit{t}\mathit{c}\mathit{h}}$	${\mathit{M}}_{\mathit{b}\mathit{e}\mathit{n}\mathit{d}\mathit{i}\mathit{n}\mathit{g}}$	${\mathit{M}}_{\mathit{T}\mathit{e}\mathit{t}}$
Bunny	2990	5976	16,701	8964	10,724
Armadillo	6362	12,720	35,442	19,080	22,721
Dragon	31,800	10,582	42,072	15,419	36,690

shows the information on the 3D model used in this paper.

In order to compare how well the simulation performs, we measured the average fps over 400 frames for a single object in the scene. To ensure a fair comparison between simple and complex objects, we set the VSync count to “Don’t Sync” instead of the default setting of “Every V-Blank” in Unity3D. This default setting targets the animation framerate to match the monitor’s framerate. As shown in Table 3, the local workgroup’s grid size for each GPU kernel was set to 1024, which is the maximum size allowed for the one-dimension of workgroup size in HLSL’s compute shader. We only perform a single iteration of solving the constraint to compare the performance. Figure 6 shows a comparison of the performance between CPU-based and GPU-based PBD methods that utilize stretching and bending constraints on a tetrahedral model. Furthermore, the performance comparison between CPU-based and GPU-based PBD methods that use all types of constraints (stretch, bending, and volume constraint) on the tetrahedral model is illustrated in Figure 7.

The results when using only stretch and bending constraints show that using a GPU provides significant speedup factors compared to using only a CPU for all three models. The constraint-centric GPU approach provides a speedup factor of up to 2x for Bunny and up to 6.5x for Armadillo compared to CPU-only processing. Similarly, the node-centric GPU approach provides a speedup factor of up to 2.1x for Bunny and up to 8x for Armadillo compared to the CPU-only processing. When we enforce the volume constraint on the deformable object, the results demonstrate that using a GPU provides significant speedup factors compared to using only a CPU for all three models. The constraint-centric GPU approach provides a speedup factor of up to 3.7x for Bunny, 6.2x for Armadillo, and 8.9x for Dragon compared to CPU-only processing. Similarly, the node-centric GPU approach provides a speedup factor of up to 3.8x for Bunny, 6.2x for Armadillo, and 8.8x for Dragon compared to CPU-only processing. Additionally the results demonstrate that the utilization of all types of constraints exhibit only a slight decrease in performance compare to the utilization of stretch and bending constraint, however it still operates at real-time speeds.

The number of vertices for Dragon is larger than Bunny to reflect the result of CPU-based comparison. However, in this work, we didn’t limit the maximum FPS in any case. The results of the GPU-based comparison are slightly differen, as unity3D has the potential to restrict the number of fps. Note that the empty scene only obtains approximately 300 fps. In addition, the usage of multi-threading in the GPU makes the simulation of the deformable object to be fast enough, and the differences between GPUs are much smaller.

When comparing the performance results to previous works that used GPU-based PBD in Unity3D, Lee et al. [24], implemented the PBD method in Unity3D using a compute shader with distance, bending, and strain constraints. The real-time performance (≥30 fps) can handle approximately 6K vertices and 19K constraints, while our proposed approach can handle 31K vertices and approximately 94K constraints in real time. In addition, Va et al. [46] proposed a GPU-based algorithm for deformable object simulations using the OpenGL platform. The simulation using an 82k vertices and 93k triangles model still achieves real-time performance, with 154 fps using the mass-spring system, 149 fps for the distance constraint, and 150 fps for the volume preservation constraint method. Although this performance result is faster and more accurate than our GPU-based PBD method, our proposed method in unity3D is more popular and better suited for general problems, with a superior rendering method. In addition, the mass-spring method is more accurate than PBD; however, it requires a small integration time step to maintain a stable simulation. While the PBD is unconditionally stable, it is also proof that it can restore the deformed shape to its original, as shown in [47].

We performed additional analysis to compare the performance of the deformable object in terms of volume error under various conditions. Firstly, we investigated the volume error of the bunny model. Since the number of iterations significantly affects the accuracy of the constraint solving, we employed a range of iteration numbers to resolve the constraints for the stretch and bend constraints cases, as well as for the case involving all types of constraints, which encompasses the volume constraints for the tetrahedron. Figure 8 shows the volume error over time in the free-fall simulation of the bunny model. Note that the stiffness coefficient for every constraint is set to 1.0 in all scenarios. Additionally, the time step used in the experiments has an impact on the simulation results. In this paper, the time step used for all experiments is set to 0.005.

As a result of the experimental comparison, the maximum volume error of simulation using only stretch and bending constraints that use 5, 10, 15, 20, and 25 iteration numbers to solve the constraint obtains the maximum volume loss of 64.66%, 34.5%, 20.79%, 15.39%, and 12.31%, respectively. By enforcing the volume constraint on each tetrahedron element, the volume error is decreased, and the shape is restored to its original form. The maximum volume error using all types of constraint using 5, 10, 15, 20, and 25 iterations obtains 60.93%, 31.10%, 17.79%, 13.34%, and 11.47%, respectively. Snapshots of the simulation using the bunny model are illustrated in Figure 9.

An additional experimental test was also conducted on a different model. We used the armadillo model as a deformable object falling to the ground. We then computed the volume error over time. Figure 10 shows the volume error over time in a free-fall simulation of the armadillo model. After conducting an experimental comparison, the results showed that the use of only stretch and bending constraints, with iteration numbers of 5, 10, 15, 20, and 25, resulted in maximum volume losses of 86.95%, 64.45%, 51.11%, 43.90%, and 42.47%, respectively. However, by enforcing volume constraints on each tetrahedron element, the volume errors decreased, and the shape was restored to its original form. When using all types of constraints with the same iteration numbers, the maximum volume errors obtained were 79.12%, 53.03%, 37.83%, 32.19%, and 28.07%, respectively. Snapshots of the simulation we conducted using the armadillo model are illustrated in Figure 11.

The experimental results suggest that enforcing volume constraints on each tetrahedron element can significantly decrease the volume errors and restore the shape to its original form. Compared to using only stretch and bending constraints, which resulted in a higher maximum volume loss, the use of all types of constraint yielded lower maximum volume errors. Iteration numbers also played a crucial role, as increasing the number of iterations resulted in lower maximum volume errors for all types of constraints. Overall, these findings demonstrate the importance of considering volume constraints in simulations to ensure accurate and realistic results. In the volume error analysis, the number of iterations ranged from five to twenty-five, but in the performance comparison, the number of iterations was only one. However, the calculations using GPU-based method are fast enough that the given range of the iteration numbers makes no difference. In fact, we increased the iteration up to 150, as shown in Figure 12, which compares of the performance between the GPU-based method performing the deformable object using the PBD method with 150 iteration numbers. The result shows that even increase the iteration number to 150, the large deformable model still performs at real-time speed.

While we have used Unity3D as the case study in this paper, OpenGL also provides the same possibility to accelerate deformable object simulation using the PBD method [48]. Both GPU-based OpenGL and Unity offer hardware acceleration using compute shaders. GPU-based OpenGL is a low-level graphics API that provides direct access to the GPU, allowing developers to optimize rendering performance at a low level. This approach is more flexible and customizable, allowing developers to fine-tune their rendering pipeline to meet the specific needs of their application. However, it also requires more programming knowledge and can be more complex to implement. On the other hand, Unity3D is a high-level game engine that provides a more user-friendly and streamlined approach to 3D rendering. It offers a range of pre-built tools and features that simplify the development process, including a visual editor for building scenes, built-in physics simulation, and support for a range of platforms. However, it may not be as customizable as GPU-based OpenGL and may not offer the same level of optimization for specific hardware configurations. However, the choice between GPU-based OpenGL and Unity will depend on the specific needs of the application and the expertise of the development team. For applications that require a high degree of customization and optimization, GPU-based OpenGL may be the better choice. For those that prioritize ease of use and streamlined development, Unity may be the more suitable option.

6. Conclusions and Future Works

This paper proposes a novel methodology for designing and implementing a simulation of a deformable object that is based solely on volumetric mesh. To achieve this, we utilized a position-based dynamic approach and integrated a variety of constraints to enhance the object’s volume which was executed in the Unity3D environment. Furthermore, our simulation implementation was specifically designed to be executed on a GPU-based approach, leveraging the computational power of graphics processing units to accelerate the performance compared to traditional CPU-based methods. This implementation strategy allows for efficient and fast simulation of complex deformable objects with plausible results.

In this study, the performance of a GPU-based PBD approach for simulating deformable objects was compared to a CPU-based PBD approach. The results showed that the GPU-based approach provided significant speedup factors for all three models tested, with speedup factors ranging from 2x to 8.9x depending on the approach used and the model tested. Additionally, this study examined the effect of enforcing volume constraints on deformable objects. The results demonstrated that the use of volume constraints significantly decreased volume errors and restored the shape to its original form. The maximum volume errors were found to decrease as the number of iterations increased, emphasizing the importance of considering iteration numbers in simulations. Overall, these findings highlight the benefits of using a GPU-based approach and enforcing volume constraints in simulations to ensure accurate and realistic results. This study also provides valuable insights for researchers and practitioners working in the field of computer graphics and simulation. Future research could explore further optimization of GPU-based approaches and investigate the impact of other types of constraints on simulation performance and accuracy.

Despite these results, there are still some limitations to our approach that need to be acknowledged. One limitation is that our method is currently focused on volumetric mesh-based deformable objects, which may not be suitable for all types of simulations. Additionally, while our approach provides significant speedup compared to a CPU-based approach, it may still not be fast enough for certain real-time applications with more complex objects and higher-resolution meshes. Another limitation is that the accuracy of our method strongly relies on the number of iterations used to solve the constraints. While we have demonstrated that increasing the number of iterations can improve accuracy, this also comes at the cost of increased computation time. Furthermore, it should be noted that our discussion on the constraints used in deformable object simulation only focused on the dihedral bending constraint. However, it is worth mentioning that isometric bending constraints can also be utilized for this purpose.

Lastly, our current implementation does not take into account self-collisions between the deformable object and other objects in the environment, which could lead to unrealistic behavior. Incorporating collision detection and response would be an important area for future research.