<figure><img src="image_1.png" /></figure>

HAL
open science

# Efficient parallelization strategy for real-time FE simulations

Ziqiu Zeng, Hadrien Courtecuisse

► To cite this version:

Ziqiu Zeng, Hadrien Courtecuisse. Efficient parallelization strategy for real-time FE simulations. 2023. [hal-04122582](https://hal.science/hal-04122582)

HAL Id: hal-04122582

https://hal.science/hal-04122582

Preprint submitted on 8 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# EFFICIENT PARALLELIZATION STRATEGY FOR REAL-TIME FE SIMULATIONS*

ZIQIU ZENG † AND HADRIEN COURTEUISSE ‡‡

**Abstract.** This paper introduces an efficient and generic framework for finite-element simulations under an implicit time integration scheme. Being compatible with generic constitutive models, a fast matrix assembly method exploits the fact that system matrices are created in a deterministic way as long as the mesh topology remains constant. Using the sparsity pattern of the assembled system brings about significant optimizations on the assembly stage. As a result, developed techniques of GPU-based parallelization can be directly applied with the assembled system. Moreover, an asynchronous Cholesky precondition scheme is used to improve the convergence of the system solver. On this basis, a GPU-based Cholesky preconditioner is developed, significantly reducing the data transfer between the CPU/GPU during the solving stage. We evaluate the performance of our method with different mesh elements and hyperelastic models and compare it with typical approaches on the CPU and the GPU.

**Key words.** Real-time simulation, parallel algorithms, finite-element method

**AMS subject classifications.**

**1. Introduction.** Medical simulations have received strong interest in providing unlimited access to learn and rehearse complex interventions in a safe environment, without the ethical issues associated. Medical simulations have become more and more realistic, offering the possibility to simulate complex interactions in real-time such as contacts and friction between deformable structures. The general trend is currently focused on the possibility to bring medical simulations closer to the Operating Room (OR) for planning interventions, or even directly in the OR with visual assistance, registration, and augmented reality (AR). For this purpose, simulations must meet the antagonistic requirements of accuracy and fast computation time at the same time. Indeed, advanced finite-elements (FE) simulations are necessary to predict the complex behavior of organs during surgery and provide relevant information for surgeons in real-time.

The material behavior of tissues is generally admitted as non-linear. Hyperelastic FE models are nowadays compatible with real-time computations [30]. However, to provide large-scale simulations of detailed meshes, parallelization strategies must be employed to maintain the computational expense sufficiently low to account for user interactions. In this context, *general-purpose computing on graphics processing units (GPGPU)* has been widely studied because it provides access to massively parallel architecture with very low-cost memory transfers compared to distributed machines.

Low-level parallelization strategies are necessary to exploit the computational power of the GPUs efficiently. When applied to FEM, standard approaches aim to accumulate all the elements' contributions in parallel. Nevertheless, since several elements share the nodes of meshes, additional operations are necessary to handle concurrent memory access. Many solutions have been proposed to address this issue. For instance, [1] proposed a GPU-based matrix-free approach, enabling high-speed FE simulations with tetrahedral co-rotational elements. However, such a technique of parallelizing FE models is very invasive in the code. Specific solutions are necessary because each material law leads to different arithmetic operations needed to compute per-element matrices. Moreover, since memory

<footnote>†ICUBE, UNIVERSITY OF STRASBOURG, FRANCE</footnote>

<footnote>‡CNRS, FRANCE</footnote>

<footnote>*Submitted to the editors on May 2 2022.</footnote>

<footnote>**Funding:** This work was supported by French National Research Agency (ANR) within the project SPERRY ANR-18-CE33-0007.</footnote>
consumption also depends on the dimension of elements, the ratio between memory access and arithmetic operations brings additional concerns to manage the little cache memory available on the GPUs. As a result, the algorithm presented in [1] can hardly be applied in other constitutive models. A generic parallelization strategy for different models remains a challenge.

On the other hand, as a popular technique used in system solvers, preconditioning boosts convergence and improves performance. [13] introduced an asynchronous preconditioning scheme where a Cholesky preconditioner is factorized in a parallel thread. This strategy significantly accelerates the convergence since the preconditioner gives a close approximation to the system matrix. The overhead of factorization is removed from the main simulation loop, making the solving stage very efficient. However, applying the preconditioner (consisting of solving triangular systems) remains a CPU-based operation, leading to considerable data transfer between the CPU/GPU. Parallelizing the triangular systems on the GPU remains a challenge, as the forward/backward substitutions lead to data dependency all over the solving stage.

This paper introduces a framework for the system solver of FE simulations. Based on the currently fastest solving strategy in the *Sofa* framework, our main contributions are:

1. An efficient matrix assembly method compatible with generic constitutive models is proposed. The method requires no specific implementation of the material law on the GPU, but allows efficient solver with typical GPU-based implementation.

2. A fully GPU-based solving strategy, including the application of the asynchronous preconditioner, is proposed. The transfer between the CPU and the GPU is minimal due to an efficient GPU-based Cholesky solver.

These improvements enable efficient GPU-based parallelization for generic constitutive models.

The rest of this paper is organized as follows. After reviewing the related works in section 2, section 3 presents the relevant deformable models and preconditioning techniques used in this work. Section 4 is dedicated to the fast matrix assembling operation, and section 5 describes the parallelization strategies. In section 6 the method is extended to handle collisions and impose boundary conditions. Finally, the method is evaluated in section 7 using different FE models.

**2. Related works.** The technical level of computer-based training systems is increasing. Early works in this field proposed simplified models such as mass-spring systems [28]. Such discrete methods are simple to implement and fast, but material properties are difficult to parameterize. For this reason, they have been progressively replaced by Finite Element (FE) models. FE models provide a better understanding of the mechanisms involved in physiological or pathological cases, mainly because the soft-tissue behavior is directly explained through constitutive relations. With the rapid growth of computational power, FE models have become compatible with real-time and interactivity. First limited to linear elastic models [10], it was later extended to large displacements with the co-rotational formulation [18]. FE models are now used for the simulation of hyperelastic or viscoelastic materials in real-time [30], with advanced and complex interactions [14] between multiple structures. On the other hand, meshless methods, Position-Based Dynamics (PBD), and Neural Networks are other strategies to model soft tissues in real-time. A detailed review of this topic goes far beyond the scope of this article, but a survey can be found in [40].

**2.1. Time discretization.** In the context of interactive simulations, an important choice is the time integration scheme. Indeed, explicit methods have been widely used for medical simulations [26]. In this case, the solution only involves the (diagonalized)
mass matrix leading to very fast, simple to implement, and parallelizable solutions [12]. Unfortunately, user interactions may introduce sudden and stiff contacts at arbitrary location/frequency, which raises stability issues.

On the opposite, implicit methods are unconditionally stable, i.e., stable (but not necessarily accurate) for any time step and arbitrary stiff materials [3]. Implicit schemes provide better control of the residual vector and hence that the external and internal forces are balanced at the end of the time steps. Although these advantages come at the cost of solving a set of linear equations at each time step, implicit integration schemes offer a reasonable trade-off between robustness, stability, convergence, and computation time, particularly when combined with a GPU implementation.

**2.2. Solving the set of nonlinear equations.** The nonlinear problem obtained by implicit schemes is usually solved using an iterative Newton Raphson method. Each iteration of the Newton method consists of solving a linear problem whose solution reduces the error between internal and external forces. Since the mechanical matrices depend on the current position (and potentially velocities) of FE meshes, the linear problems must then be recomputed for each simulation step and ideally for each Newton iteration. Therefore, the performances of the simulation are directly linked to the efficiency of the solver, which explains why earlier studies mainly focused on sparse linear algebra. Two families of algorithms are proposed in the literature: direct and iterative methods.

Direct solvers provide the exact solution by computing a factorization (for instance, the Cholesky factorization [5]) or a decomposition (QR decomposition), or eventually, the actual inverse of the system matrix [10] (though not recommended for large matrices). As proposed in [20], the nested dissection ordering has been widely used in direct solvers by exploring the parallelism of subproblems while reducing the fill-in of the matrix. Partitioning and ordering are usually implemented through external tools such as METIS [27].

The solving phase can then be performed with so-called *forward/backward substitution*, using the two triangular systems (in the case of the Cholesky factorization). Efficient libraries exist both on the CPU (Pardiso, MUMPS, Taucs) and GPU (cuSPARSE, MAGMA, AmgX). The solving stage can be improved by partitioning and reordering the system [23]. Despite the stability of direct solvers, the complete factorization or decomposition of large matrices is usually too time-consuming to be recomputed at each time step, and it is very difficult to parallelize. Specific optimizations, inspired by the co-rotational model, have been proposed to incrementally update the sparse Cholesky factorization [22] but this approach does not extend to other material laws or element types.

In the interactive context, iterative methods are usually preferred because they limit the number of iterations to compute an approximated solution and better control the time spent during the solving process. The most popular method is the Conjugate Gradient (CG) algorithm [36], because of the fast convergence and its simple implementation. Parallel implementations both on CPU [34, 24] and GPU [7, 11, 1] were proposed. Nevertheless, the convergence of iterative methods can be significantly impacted for ill-conditioned problems, i.e., when the ratio of the largest and smallest eigenvalues is large.

**2.3. Matrix assembly and parallelized solver.** The main issue to improve the CG is to gain speedup on sparse matrix-vector multiplication (*SpMV*) operations. As it is presented in [6], to accelerate the *SpMV* operations, many methods are explored to implement them on throughput-oriented processors such as GPU. Several methods rely on the fact that CG iterations can be performed without explicitly assembling the system matrix [33, 31]. Matrix-free methods significantly reduce the memory bandwidth and are proven to be fast and stable. However, as a price of speedup, it lacks generality. As
an example, the method introduced in [1] is designed for the co-rotational formulation and relies on specific cache optimization to compute rotation matrices directly on the GPU. However, the specific cache optimizations proposed for the rotation matrices do not extend to other types of material, such as hyperelastic laws.

Explicit assembly of global matrices is necessary for direct solvers to compute the factorization or decomposition of the system. The assembly step is usually less critical than the solving process itself, but it may become the bottleneck when combined with efficient solvers. There are several ways to construct sparse matrices; the most popular method is first to collect triplets (the row/column index and the value); then compress the triplets in a sparse format. A very efficient implementation is provided in the Eigen library. Recently [25] proposed a row by row assembling method for isogeometric linear elasticity problems. To accelerate the assembling step and minimize memory transfers, several approaches proposed to assemble the matrix directly on the GPU [17, 39, 19]. However, specific GPU-based implementation of the assembling procedure is needed for each particular model.

**2.4. Preconditioner.** Another intense area of research aims to improve the performance of the CG algorithm with the use of preconditioners to speed up its convergence. There are several typical preconditioners: diagonal matrix is simple to build but has limited effect [4]; in contrast, precise ones such as incomplete Cholesky factorization are complex and costly to make but can significantly reduce the condition number [21].

For a typical synchronous preconditioner, the construction of the preconditioner has to be performed before the solving stage of each time integration, leading to additional computation costs. Some of the recent works aim to find a balance between the cost of applying the preconditioner and the effect of convergence boost, such as efficient preconditioners using the result of incomplete factorization [2] and inner Gauss-Seidel preconditioners [37].

On the other hand, the asynchronous preconditioners proposed in [13] exploit the continuity of the time line in physically-based simulations. Relying on the assumption that mechanical matrices undergo relatively small changes between consecutive time steps, the asynchronous preconditioning scheme processes the matrix factorization in a dedicated thread parallel to the main simulation loop and applies the factorization result as a preconditioner after a short delay. It enables access to a very efficient preconditioner with almost no overhead in the simulation loop. As a combination of a direct and iterative solver, the method requires explicitly assembling the matrix at a low frequency in the simulation loop to factorize the system in the dedicated thread. For both synchronous and asynchronous preconditioning schemes, applying the preconditioner requires processing the *forward/backward substitution*, leading to solving sparse triangular systems (STS). Parallelizing the solution of STS remains challenging in many applications. There are many works dedicated to improving the performance of STS solvers on the CPU [9] and on the GPU [35, 38, 29]. In [14], a GPU-based asynchronous preconditioner was designed to solve the STS with multiple right-hand sides (RHS) in the contact problem. However, the method cannot efficiently exploit parallelization when dealing with a single RHS. Therefore, despite the asynchronous preconditioning scheme being introduced with a GPU-based CG implementation of the co-rotational model, applying the preconditioner was performed on the CPU, requiring data transfers between CPU/GPU for each iteration of the preconditioned CG.

**2.5. Implementation.** It is important to note that even each model can be efficiently parallelized on the GPU with specific implementation, it will be hard to be developed and maintained in a large and generic framework such as *SOFA framework*.
A strong motivation of the current work lies in the fact that once the matrices are assembled, the solver can be parallelized independently from the FE models that generated the matrix (i.e., the material law or the type of elements of the FE mesh). For this purpose, a generic data structure is proposed to fill non-null values of mechanical matrices. The method exploits the fact that contributions are added to the system matrix in a deterministic way that only depends on the topology. Efficient GPU-based parallelization operations such as SpMV can be implemented with the assembled system. A specific parallelization strategy are also proposed to apply the preconditioner on the GPU, allowing this to parallelize the entire solving process independently from the constitutive laws or the type of elements used in the simulation.

**3. Background.** The current method is based on a general background for de-formable simulations using the implicit integration.

**3.1. FE models and constitutive law.** In order to underline the importance of the model, the method is tested with a co-rotational formulation [32] and with the hyperelastic models [30]. The method also applies for various types of elements such as beams, triangles, tetrahedra or hexahedra.

Using the co-rotational formulation, local stiffness matrices can be precomputed for each element e with the synthetic formulation:

$$ (3.1) \qquad \mathbf{K}_e = \int_{V_e} (\mathbf{C}_e \mathbf{D}_e \mathbf{C}_e^T dV_e) $$

where $\mathbf{D}_e$ corresponds to the stress-strain matrix, $V_e$ is the volume of the element and $\mathbf{C}_e$ is the strain-displacement matrix. The method is parametrized with the Young modulus $E$ and the Poisson's ratio $\nu$. The hyperelastic model is implemented with the Multiplicative Jacobian energy decomposition (MJED) method [30]. The method consists in decoupling the invariants of the right Cauchy-deformation tensor $C = \nabla\Phi^T\nabla\Phi$ from the expression of the deformation tensor $I_1 = \text{tr}C$, $I_2 = ((\text{tr}C)^2 - \text{tr}C^2)/2$ and the Jacobian $J = \det\nabla\Phi$, where $\Phi$ is the deformation function between the rest and deformed configurations. The method allows faster stiffness matrix assembly for a large variety of isotropic and anisotropic materials. In these formulations, the function $\mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})$ provides internal forces of the deformable body, given the nodal position $\mathbf{q}$ and velocities $\dot{\mathbf{q}}$.

**3.2. Time integration and implicit scheme.** For any time step $t$, the general way to describe the physical behavior of a deformable objective problem can be expressed using Newton's second law:

$$ (3.2) \qquad \mathbf{M}\ddot{\mathbf{q}} = \mathbf{p} - \mathcal{F}(\mathbf{q}, \dot{\mathbf{q}}) $$

Where $\mathbf{M}$ is the mass matrix, $\ddot{\mathbf{q}}$ the vector of the derivative of the velocity, $\mathbf{p}$ the external forces and $\mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})$ the function representing the internal forces.

A backward Euler method is used to integrate the time step. The implicit scheme can be expressed as follows, where $h$ is the length of time interval $[t, t + h]$:

$$ (3.3) \qquad \dot{\mathbf{q}}_{t+h} = \dot{\mathbf{q}}_t + h\ddot{\mathbf{q}}_{t+h} \qquad \mathbf{q}_{t+h} = \mathbf{q}_t + h\dot{\mathbf{q}}_{t+h} $$

As $\mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})$ is a non-linear function, a first-order Taylor expansion is performed to linearize the problem [3]. This linearization corresponds to the first iteration of Newton-Raphson method. The incomplete approximation may cause numerical errors of the dynamic behavior but they lean towards to decrease at equilibrium.
The internal forces are expanded as following:

$$ (3.4) \qquad \mathcal{F}(\mathbf{q}_{t+h}, \dot{\mathbf{q}}_{t+h}) = \mathbf{f}_t + \frac{\partial \mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})}{\partial \mathbf{q}} h \dot{\mathbf{q}}_{t+h} + \frac{\partial \mathcal{F}(\mathbf{q}, \dot{\mathbf{q}})}{\partial \dot{\mathbf{q}}} h \ddot{\mathbf{q}}_{t+h} $$

with $\mathbf{f}_t = \mathcal{F}(\mathbf{q}_t, \dot{\mathbf{q}}_t)$.

During a time integration, the force function is considered as constant and the partial derivative terms could be expressed as matrices: $\frac{\partial \mathcal{F}}{\partial \ddot{\mathbf{q}}}$ the damping matrix **B** and $\frac{\partial \mathcal{F}}{\partial \dot{\mathbf{q}}}$ the stiffness matrix **K**.

By integrating the equations (3.2), (3.3) and (3.4) we obtain the dynamic equation:

$$ (3.5) \qquad (\mathbf{M} + h\mathbf{B} + h^2\mathbf{K}) \ddot{\mathbf{q}}_{t+h} = \mathbf{p}_t - (\mathbf{f}_t + \mathbf{B}\dot{\mathbf{q}}_t + h\mathbf{K}\dot{\mathbf{q}}_t) $$

With Rayleigh damping [3], the damping matrix can be expressed as a combination of matrices of mass and stiffness **B** = $\alpha\mathbf{M} + \beta\mathbf{K}$ with $\alpha$ and $\beta$ the proportional Rayleigh damping coefficients. By replacing **B** in the dynamic equation (3.5), it gives:

$$ (3.6) \qquad \underbrace{[(1 + h\alpha)\mathbf{M} + h(\beta + h)\mathbf{K}]}_{\mathbf{A}} \ddot{\mathbf{q}}_{t+h} = \underbrace{\mathbf{p}_t - \mathbf{f}_t + h\mathbf{K}\dot{\mathbf{q}}_t}_{\mathbf{b}} $$

Equation (3.6) provides a linear problem $\mathbf{A}\mathbf{x} = \mathbf{b}$ to solve. The left-hand side is a global system matrix **A** and the right-hand side a vector **b**. Both of them are constructed by some elements: **M** the mass matrix, **K** the stiffness matrix as well as scalar parameters (the time interval and the Rayleigh damping coefficients). The linear system must be solved at each time step as **K** depends on the position of FE models. Since **A** is large and sparse, the general-propose compressed sparse row (CSR) format is used to store the matrix information in three arrays.

**3.3. Asynchronous preconditioner.** Let $\mathbf{A}_t$ be the matrix built in a specific time $t$. Following [13], a preconditioner **P** can be built from an asynchronous $\mathbf{LDL}^T$ factorization:

$$ (3.7) \qquad \mathbf{P} = \mathbf{A}_t = \mathbf{LDL}^T $$

Where **D** is a diagonal matrix and **L** a sparse lower triangular matrix. The factorized matrices will be available after the factorization is done, normally several time steps after time $t$, and used as a preconditioner with the assumption that **P** remains a relatively good approximation to the current matrix $\mathbf{A}_{t+nh}$. In practice, the method is very efficient because the $\mathbf{LDL}^T$ factorization requires only few simulation steps (usually $n < 5$).

The application of the preconditioner consists mainly of solving the two triangular systems obtained after the factorization. The method is very efficient because in practice only 2 to 5 preconditioned CG iterations are necessary to converge (with threshold of $10^{-9}$). However, despite the triangular matrices are sparse and the solution can easily be implemented with a *gauss elimination* on the CPU, this step is difficult to parallelize on a GPU due to numerous data dependencies. A GPU-based solver with a CPU-based preconditioner leads to considerable data transfer between the different processors, making the solving process inefficient. Improving this process is an important issue that will be addressed in Section 5.

**4. Matrix assembly strategy.** Generic constitutive models can benefit from typical GPU-based matrix operations. But the matrix assembly usually leads to an overhead cost, which is not negligible. To address the issue, we propose a new assembly approach to meet the requirements of both efficiency and generality. The fast assembly method relies on the fact that the same assembly procedure is called in each time integration. When
building the system matrix in sequential order, the invariant topology structure brings an important property that the sequence of filling elements into the matrix and the sparsity pattern of **A** are definitive. A specific mapping from the filling element sequence to the final matrix pattern could be built. As sorting the initial filling sequence to the final sparse format is the most time-consuming stage in the matrix assembly, replacing it with the deterministic mapping brings a significant speedup.

<figure><img src="image_3.png" /><figcaption>Fig. 1: General workflow of the matrix assembly procedure. **Data Mapping** corresponds to the additional structure used to compress the unstructured contributions (see section 4.2) in CSR format. It is computed and sent on the GPU only once until no modifications of the fill ordering are detected during the collection phase.</figcaption></figure>

An overview of the general workflow of the assembly procedure is shown in the figure 1. The matrix assembly consists of the following steps:

1. **Collect data**: Collect the data of mass and stiffness for each element and fill the data in a triplet format (row index, column index, and value).

2. **Build matrix pattern**: Sort the collected triplet data by order of row and column. This step is necessary if and only if modifications of the structure have been detected during the collection phase.

3. **Compress**: Build the system matrix in CSR format.

The current matrix assembly method relies on the assumption that the topology remains invariant. Topological modifications are not addressed in this paper but the method remains generic since the topology modification only occurs in specific cases, such as cutting operations. Applying the asynchronous preconditioning in such case of sudden changes can be addressed with specific correction on the preconditioner [14]. Collisions and interactions may also change the fill ordering of the matrices. This specific issue is discussed in section 6.

**4.1. Collect data.** Matrices **M** and **K** in equations (3.6) are obtained by summing the local contributions of each element into the global matrices. Values are stored in a set of triplets at first, which is a structure containing 3 variables: row index, column index, and value. As the triplet vector corresponds to the original process of filling elements into the matrices, the sequence of row/column indices is unsorted and uncompressed<fnref n="1" />, but definitive in each time integration. Nevertheless, insertion of contributions into the triplet list will be called many times per second; it must therefore be optimized as much

<footnote n="1">¹i.e. the pair row/column may appear several times when filling matrices</footnote>
as possible. The pseudo-code of the **add** function is given in the algorithm 1, and exposed to the FE models in order to insert their contributions.

1 **Function add(col, row, val)**
2   **if** keepStruct and id < prevVal.size() and prevCol[id] = col and
      prevRow[id] = row **then**
3   **|** prevVal[id] = val;
4 **else**
5   **|** keepStruct = **false**;
6   **|** prevRow[id] = row;
7   **|** prevCol[id] = col;
8   **|** prevVal[id] = val;
9 **end**
10 **id**=id+1;

**Algorithm 1:** Procedure used to add value in the matrix. The boolean *keepStruct* is used to detect any modification in the filling order; *id* gives the next writing address in the uncompressed arrays (*prevRow*, *prevCol* and *prevVal* corresponding to the triplets list added in the previous time steps).

For each inserted value, the test performed line 2 checks the consistency of the pattern with respect to the previously built matrices. This test is a necessary overhead to detect changes in the structure. However, if the structure is not modified, only the value *val* is stored (line 7), allowing this way to take advantage of the cache of the CPU and minimize write operations.

**4.2. Build matrix pattern.** Let **X** be a generic matrix to be assembled (such as **M** and **K**). A method inspired by the Eigen's library is implemented to build the final CSR format for **X**. The method consists of computing twice the transpose of the matrix to sort the values. To store the temporary matrices, we introduce a temporary format called **uncompressed structure** which is similar to the CSR format: Like the CSR, an arranged row pointer encodes the index in the arrays of column index and values that are unsorted and uncompressed (duplicate indices exist). We summarize the states of the assembly in different stages in Table 1.

1. Firstly, the temporary transposed matrix $\mathbf{X}^T$ is built in the **uncompressed structure**. The computation of the transposed matrix requires beforehand to count the number of values per line, allowing this to allocate the necessary memory. Then, data can be moved to their correct location in the allocated structure. With the pre-defined matrix structure, the sequence of row index can be arranged with a time complexity of $O(2n)$, but inside each row, the sequence of column index remains unsorted.

2. Similar to the previous step, **X** is built in the **uncompressed structure** by transposing $(\mathbf{X}^T)^T$. The second transpose gives the initial matrix **X** with a sequence of values sorted both by rows and columns while the structure remains uncompressed.

3. Finally, the elements in the same position are merged, transferring the **uncompressed structure** into the CSR format.

One of the main differences with the Eigen's implementation is that the values of the transposed matrices are not directly stored in memory, making a strategy of **fast assembly** possible: Relying on the hypothesis that the mesh topology remains unchanged, the filling order, as well as the matrix pattern (row pointer and column index arrays in the CSR), could be reused. As long as the filling order remains unchanged at the collecting stage,
<table><caption>Table 1: State of storage format at different stages in matrix assembly process</caption><thead><tr><th>Matrix</th><th>X</th><th>X<sup>T</sup></th><th>X</th><th>X</th></tr></thead><tbody><tr><th>Format</th><td>triplet set</td><td>uncompressed<br>structure</td><td>uncompressed<br>structure</td><td>CSR</td></tr><tr><th>Row</th><td>unsorted<br>uncompressed</td><td>sorted<br>compressed</td><td>sorted<br>compressed</td><td>sorted<br>compressed</td></tr><tr><th>Column & Values</th><td>unsorted<br>uncompressed</td><td>unsorted<br>uncompressed</td><td>sorted<br>uncompressed</td><td>sorted<br>compressed</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 2: build matrix pattern and mapping vector from triplets set</figcaption></figure>

we propose to build a mapping C from the initial triplet set to the CSR format, making it very efficient to build the value array. Hence, the main operation is to merge the duplicated values in the triplet array that is reordered with the deterministic mapping C. This operation must be performed at each time step, but it can be easily parallelized both on CPU and GPU since the address of values in the CSR format are known and unique. More importantly, parallelization can be performed without impacting the code of the constitutive models that generate the matrices (i.e., our method is efficient for any generic constitutive model). The deterministic mapping C can be reused as long as no modifications in the filling order are detected at the previous stage. On the other hand, in case modifications are detected at the collecting stage, we process the complete rebuilding
of the matrix pattern (called **full assembly**). In this mode, the method provides similar performances as the default implementation of Eigen's library.

In order to get the global matrix **A** and the vector **b** in equation (3.6), one may note that both are generated from the sum of the same matrices **M** and **K** with various coefficients:

$$ (4.1) \qquad \mathbf{A} = (1 + h\alpha)\mathbf{M} + h(h + \beta)\mathbf{K} $$

$$ (4.2) \qquad \mathbf{b} = \mathbf{p}_t - \mathbf{f}_t - h\mathbf{K}\dot{\mathbf{q}}_t $$

Coefficients only depend on the time step and the Rayleigh damping constant during the entire simulation. Another consequence of our method is that the computation of the right-hand side and the left-hand side terms can be merged in a single procedure, allowing this way to exact a large amount of data that are well suited for GPU architectures, and benefits from cache optimization since the mapping **C** is accessed twice.

Once the vector of values is compressed, the CSR format can be used directly inside a parallel Conjugate Gradient (either on the CPU or the GPU<fnref n="2" />). For this purpose, the product of the sparse matrix with a vector (*SpMV*) needs to be parallelized, which is trivial with the assembled system. Many efficient CPU and GPU-based implementations exist for such a typical operation. In this paper, we use the *SpMV* implementation in the CUSPARSE library developed by NVIDIA. It can significantly improve the time spent in the CG iterations, providing a significant speedup to the entire simulation without modifying the code that generates the matrices.

**5. System solution.** The system solution can be efficiently processed with typical sparse matrix operations on GPU with the assembled matrix. Furthermore, as explained in 3.3, the solver can be boosted with a preconditioner. Compared to the matrix-free method, another important consequence of the fast assembly method lies in the possibility of directly using the assembled matrix to build a preconditioner (which most of the time requires the values of the system explicitly). For instance, the diagonal extraction for the Jacobi preconditioner or the lower triangular system for the SSOR preconditioner is exceptionally facilitated.

The asynchronous preconditioner method [13] needs to access the explicit values of the assembled system matrix at some specific time steps, which does not add any additional overhead with the proposed assembled solution. The factorization of the matrix being performed asynchronously, the preconditioner can be entirely computed on the CPU without blocking the main simulation thread. However, the application of the preconditioner at each iteration of the preconditioned CG implies solving sparse triangular systems. The data dependence between lines makes it difficult to be computed in parallel. Let **L** be the lower triangular system of the Cholesky factorization. We recall the main obstacle for solving a general lower triangular system **L**s = **r** is that the solution **r**j of a given row j depends on all previous solutions **s**i:

$$ (5.1) \qquad \mathbf{s}_j = \mathbf{r}_j - \sum_{i=0}^{i<j} (\mathbf{s}_i \mathbf{L}_{j,i}) $$

Consequently, the primary operations of the Conjugate Gradient algorithm are processed on the GPU, while the application of the preconditioner remains on the CPU. This hybrid solving strategy generates a huge amount of data transfer between the processors: in

<footnote n="2">²Note that the vector of values is already available on the GPU if the compression is performed on this architecture. The row index and column pointer need to be transferred only if the mapping is modified.</footnote>
each CG iteration, the processors need to send the residual vector from GPU to CPU to apply the preconditioner and then send the result vector back to GPU. In addition to their cost, data transfers impose numerous synchronizations between CPU/GPU, reducing the efficiency of the preconditioner. We aim to implement a GPU-based preconditioner that is at least as efficient as the CPU-based version to address the issue.

**5.1. GPU-based Cholesky preconditioner.** We propose a GPU-based Cholesky preconditioner, which is inspired by the solver in [14] for sparse triangular systems with multiple right-hand sides. Since the method in [14] was originally designed for multiple RHS. When applied in the problem with a single RHS, the level of parallelism for multiple RHS will be left unused. We propose to fill in this dimension with parallelism imported by the domain decomposition technique using the nested dissection algorithm. The nested dissection algorithm is used to reduce the filling of the matrix pattern, recursively dividing the mesh into two parts with nearly the same number of vertices while keeping the divider part at a small scale [20]. Consequently, **L** is reordered and partitioned into sub-domains with the indices given by the nested dissection algorithm. The reordering algorithm partitions the graph as follows:

$$ (5.2) \qquad \underbrace{\begin{bmatrix} \hat{\mathbf{L}}_a & \hat{\mathbf{L}}_b & \hat{\mathbf{L}}_c \\ \mathbf{V}_a & \mathbf{V}_b & \mathbf{V}_c \end{bmatrix}}_{\hat{\mathbf{L}}_{(a,b,c)}} \begin{bmatrix} \mathbf{s}_a \\ \mathbf{s}_b \\ \mathbf{s}_c \end{bmatrix} = \begin{bmatrix} \mathbf{r}_a \\ \mathbf{r}_b \\ \mathbf{r}_c \end{bmatrix} $$

where the *diagonal domains a* and **b** can be solved independently and the reordering algorithm guarantees that the *separator c* (which requires the solution of **a** and **b**) is as small as possible. The partition and reordering are processed recursively on *diagonal domains a* and **b** until the block size is small enough.

The specific parallelization level is assigned to each subdomain identified in the lower triangular system **L**s = **r**. The base rule is that the blocks with higher levels (left edge) require the solutions of low-level blocks. The upper triangular system problem **L**T**s** = **r** can be solved with the same method but with an opposite sequence of computation priority. The higher level a block has, the less dependence it has.

Within each level, the parallelization strategy presented in [14] can be used to solve each *diagonal block* ($\hat{\mathbf{L}}_a$, $\hat{\mathbf{L}}_b$) and the *separator* ($\mathbf{V}_a + \mathbf{V}_b + \tilde{\mathbf{L}}_c$) by the sequence of rows. It corresponds to the *Row Major* as illustrated in the figure 3 where $t * t$ threads are used to accumulate the contribution so that $t$ rows are processed in parallel ($t = 16$ in the current implementation). Due to the high dependencies, the diagonal part is treated separately as a dense matrix in shared memory. A parallel reduction is then used to sum the contribution for each row, and finally, the $t * t$ diagonal block is solved as a dense problem.

The opposite *Column Major* is also feasible by pre-accumulating the column's contributions. Instead of solving the combined block ($\mathbf{V}_a + \mathbf{V}_b + \tilde{\mathbf{L}}_c$) in a single kernel, the accumulation process of $\mathbf{V}_a$ and $\mathbf{V}_b$ is moved into the kernel of $\hat{\mathbf{L}}_a$ and $\hat{\mathbf{L}}_b$ respectively. Since it requires only the solution of $\hat{\mathbf{L}}_a$ (or $\hat{\mathbf{L}}_b$), the accumulation of block $\mathbf{V}_a$ (or $\mathbf{V}_b$) can be processed in the same kernel. The part $\tilde{\mathbf{L}}_c$ can be solved as a diagonal block after the accumulation of $\mathbf{V}_a$ and $\mathbf{V}_b$. Similarly, the diagonal and accumulation parts are treated with $t * t$ threads, and each $t$ column is processed simultaneously. This pre-accumulation leads to data writing conflicts since several columns may contribute to the same line simultaneously. The atomic add function defined in CUDA can automatically manage the data conflict.
<figure><img src="image_5.png" /><figcaption>Fig. 3: The solving stage for each subdomain is realized by GPU kernels, where contributions are accumulated in parallel. For the lower triangular system, the solution can be processed by column sequence (left), which pre-accumulates the data in higher levels, allowing sharing of the computation cost. On the other hand, when solving the upper triangular system, computation cost could be shared in lower levels, so the solution needs to be processed oppositely by row sequence (top-right).</figcaption></figure>

As illustrated in the figure 3, in order to share the computation cost in lower levels, the lower solver is implemented with *Column Major*, and the upper system is solved with *Row Major*. Our level-based parallelization strategy is similar to the approach in [38], with several main differences:

1. Our solver uses the block-row parallelization strategy in [14] to efficiently exploit the parallelism architecture of the GPU (see Figure 3).

2. Our solver is optimized for the problems in FE simulations (e.g. we keep using the analysis result of parallelization level until the matrix pattern is changed).

3. Our solver benefits from the pre-accumulation technique which allows to share the computation cost in lower levels, making the solver more efficient (see Figure 3).

**5.2. Data Transfer between processors.** We evaluate the performance of our new GPU-based preconditioner in the following section 7. As illustrated in the table 5, our method is faster than the CPU-based implementation in various examples. Replacing the CPU-based preconditioner with our GPU-based implementation brings speedup for the solver and addresses the data transfer issue between the processors. Consequently, our new preconditioner makes it possible to execute a fully GPU-based preconditioned CG, requiring only one scalar to be transferred at each iteration from the GPU to the CPU in order to check the convergence (see Figure 4).

**6. Contact and interactions.** One can hardly talk about medical simulations without considering interactions between objects. The simulation of interacting deformable structure is an extremely large topic; a detailed review can be found in [8].
<figure><img src="image_6.png" /><figcaption>Fig. 4: Workflow of the asynchronous preconditioning scheme. The **collect** phase is performed on the CPU with any generic implementation of FE models. The **Op** process corresponds to the necessary operation to perform one CG iteration. All of them are either *Spmv* or linear algebra operations on vectors that can be easily parallelized on the GPU. The application of preconditioner **Apply Precond** is also performed on the GPU, resulting in a preconditioned CG fully implemented on GPU. Only a scalar needs to be copied to the CPU in each iteration to check the convergence state.</figcaption></figure>

**6.1. Projective Constraints.** The more straightforward solution to fix a set of nodes and impose boundary conditions is to use projective constraints. Such constraints consist of reducing the matrix dimension to remove the fixed points from the equations of motion. This can be implemented by clearing the corresponding rows and columns of the matrix and setting the value 1 on the diagonal. Clearing a row in the CSR format is trivial, but in order to clear the columns, either a filter could be added in the **add** procedure of algorithm 1 to skip the undesired values, or a search can be performed on each line afterward in order to erase non-null values on the column if it exists. Despite the fact that both strategies are highly time-consuming, they represent the only available alternatives when using the Eigen library.

With the proposed approach, the undesired values can be identified during the construction of the mapping *C*. The collection phase is not modified, and all the values associated with fixed points are added in the incoming triplets arrays with no overhead, while indices of fixed rows/columns are stored separately in specific vectors. During the computation of the matrix pattern, after the first transpose $X^T$, the rows (corresponding to the column of the final matrix) associated with fixed points can be skipped. Likewise, the undesired lines in the final matrix $(X^T)^T$ can be skipped and replaced by a single value 1 on the diagonal. Finally, the mapping *C* is built to only assemble the desired values of the projected matrix with no additional overhead in the simulation.

**6.2. Lagrangian multipliers.** Augmented Lagrangian Multipliers is an efficient solution to deal with constraints accurately and robustly. The size of the linear systems is increased with specific constraint equations, resulting in a Karush-Kuhn-Tucker (KKT) system:

$$ (6.1) \qquad \begin{cases} \mathbf{A}_1\mathbf{x}_1 - \mathbf{H}_1^\mathrm{T}\boldsymbol{\lambda} = \mathbf{b}_1 \\ \mathbf{A}_2\mathbf{x}_2 + \mathbf{H}_2^\mathrm{T}\boldsymbol{\lambda} = \mathbf{b}_2 \\ \mathbf{H}_1\mathbf{x}_1 - \mathbf{H}_2\mathbf{x}_2 = \boldsymbol{\Delta}\boldsymbol{\delta} \end{cases} $$
with subscript 1 and 2 representing two interacting objects, **H** are the linearized constraint equations, **λ** the associated Lagrangian Multipliers (contact forces) and $\Delta\delta$ the difference between interpenetration of the end and the beginning of the time step.

Contact constraints can be solved in Linear/Nonlinear Complementary Problem formulations [15], forming an LCP (linear) to simulate frictionless contact or an NLCP (nonlinear) in case of friction contact [16]. The solving process can be performed in several steps:

1. *Free motion* The motions are computed without considering the interactions between objects. It requires solving the linear systems $\mathbf{A}_1\mathbf{x}_1^{\text{free}} = \mathbf{b}_1$ and $\mathbf{A}_2\mathbf{x}_2^{\text{free}} = \mathbf{b}_2$.

2. *Constraints resolution* The constraints are defined and the compliance matrix $\mathbf{W} = \mathbf{H}_1\mathbf{A}_1^{-1}\mathbf{H}_1^T + \mathbf{H}_2\mathbf{A}_2^{-1}\mathbf{H}_2^T$, is built to solve the contact problem $\mathbf{W}\mathbf{\lambda} = \delta - \mathbf{H}_1\mathbf{x}_1^{\text{free}} - \mathbf{H}_2\mathbf{x}_2^{\text{free}}$ with the projected Gauss-Seidel.

3. *Motion correction* The motion is corrected solving equations: $\mathbf{x}_1 = \mathbf{x}_1^{\text{free}} - \mathbf{A}_1^{-1}\mathbf{H}_1^T\mathbf{\lambda}$ and $\mathbf{x}_2 = \mathbf{x}_2^{\text{free}} - \mathbf{A}_2^{-1}\mathbf{H}_2^T\mathbf{\lambda}$.

A significant advantage of using the Lagrange multipliers is that collision events never modify the system matrix. Indeed, the *Free motion* and the *Motion Correction* involve the solution of the same linear system as described in section 4 allowing this way to direct benefits for the Fast Assembling method. The computation of the *compliance matrix* is a time-consuming step that can be significantly improved using the asynchronous preconditioner and the multiple right-hand side solver proposed in [14]. Again, our fast assembly technique significantly improves the construction of the preconditioner resulting in a global speedup of constraint-based simulations.

**7. Results.** The simulation tests are conducted in the open-source SOFA framework with a CPU Intel® core i9-9900k at 3.60GHz and a GeForce RTX 2070 8 Gb.

**7.1. Matrix Assembly.** Our matrix assembly strategy aims to reach a compromise between the computation cost and the versatility of the code by assembling the matrix **A** with low cost. This section compares the matrix building time between the current assembly method and the standard assembly method implemented in the Eigen library. The simulation tests for the assembly stage are executed with a group of deformable mesh representing the shape of a raptor with various mesh resolutions (see table 2).

<table><caption>Table 2: Number of nodes and tetrahedral elements of the meshes.</caption><thead><tr><th>Example</th><th>Raptor 1</th><th>Raptor 2</th><th>Raptor 3</th></tr></thead><tbody><tr><td>Nodes</td><td>2996</td><td>4104</td><td>5992</td></tr><tr><td>Tetra</td><td>8418</td><td>12580</td><td>19409</td></tr></tbody></table>

The figure 5 shows the performances of the assembling stage, including the accumulation of triplets and the compression to the CSR format but excluding the computation of the mapping C. With the exception of the first time step where the mapping is actually computed, it corresponds to the standard performances obtained during the entire simulation with the various assembly methods. Compared with the standard method using Eigen library, the current method on CPU reduces by 72% time cost of building on average. This cost reduction rises to 81% for the fast assembly method on the GPU. The compression on the GPU provides a speedup of between 2.7× to 3× with respect to the parallel implementation of the compression on the CPU using 8 threads.

If topological modifications are performed or if the filling order is modified, the matrix
<figure><img src="image_7.png" /><figcaption>Fig. 5: Computational costs for the matrix assembly in **fast_assembly** mode (operations excluding the re-computation of the compression mapping C) and **full_assembly** mode (when rebuilding the matrix pattern) introduced in Section 4.2. We compare the performance between the Eigen's library implementation (**Eigen assembly**) and our fast assembly strategy (**FA**) with the compression performed on the CPU (with 8 CPU threads) or on the GPU.</figcaption></figure>

pattern needs to be rebuilt. In this case, the building cost, including the computation of the pattern, is measured in the figure 5. The time cost of the current method on CPU when the matrix pattern is rebuilt is slightly slower than the Eigen implementation, but it remains in the same order. The overhead is due to the additional computation of the index vector mapping C providing the position of the triplets in the CSR format. However, the cost is balanced because the mapping can be reused for the next time steps. Indeed, reusing the mapping for only two consecutive time steps already provides an acceleration compared to the Eigen implementation. Since the computation of the mapping is performed on the CPU, the GPU-based compression suffers a slowdown due to data transfers between the CPU and the GPU.

<table><caption>Table 3: Configurations of different scenario examples.</caption><thead><tr><th>Example</th><th>Liver</th><th>Cloth</th><th>Cube</th><th>Raptor</th></tr></thead><tbody><tr><th>Model</th><td colspan="3">Co-rotational</td><td>Hyperelastic</td></tr><tr><th>Type</th><td>Tetrahera</td><td>Triangle</td><td>Hexahedra</td><td>Tetrahedra</td></tr><tr><th>Nodes</th><td>2660</td><td>4900</td><td>8000</td><td>2996</td></tr><tr><th>Nb.element</th><td>12328</td><td>9522</td><td>6859</td><td>8418</td></tr><tr><th>Mesh</th><td></td><td></td><td></td><td></td></tr></tbody></table>
**7.2. Performances with the CG solver.** The performances of the global simulation are now compared in a complete simulation of a deformable body, including the time for the computation of the FE model, the assembling step and the solving process. Performances of the fast assembly method combined with a Conjugate Gradient solver (**CG GPU fast assembly**) is measured and compared with both a CPU-based matrix-free implementation of the Conjugate Gradient (**CG CPU matrix-free**) and with the method introduced in [1] which includes a matrix-free GPU-based Conjugate Gradient (**CG GPU matrix-free**) for the tetrahedral co-rotational model.

In order to verify the generality of the proposed solution, the specific GPU-based implementation introduced in [1] has been extended for other types of elements (triangles and hexahedron), requiring the development of specific code for each model on the GPU. In addition, the fast assembly method is also tested for hyperelastic material laws. However, since developing an efficient GPU parallelization is not trivial, the method is only compared with CPU-based matrix-free solvers. The scenarios are illustrated in Table 3.

For the scenarios in Figure 6, the run time increases linearly along with the number of iterations. The fast assembly method combined with the GPU-based CG is up to 16× faster than the standard CPU method implemented in SOFA and reaches the same computation cost level as the GPU-optimized method. The Fast Assembly method suffers a slowdown compared to the GPU matrix-free method with fewer iterations, but this case is inverted when the iteration increases. This is due to the fact that the fast assembly method takes time to build the matrix, but this overhead is compensated at each CG iteration since the parallel implementation of the *SpMV* operation is faster with the assembled matrix.

It's important to note that although performances are comparable to the GPU-based matrix-free implementation, the code of the co-rotational model is written for the CPU where optimizations are simply obtained by calling the **add** function of the algorithm 1, which is completely transparent for the code and enforces the compatibility with the rest of the models implemented in the SOFA framework. In addition, for computers without GPU-compatible hardware, the *SpMV* operation can also be parallelized on the CPU. The method **CG CPU fast assembly** uses 8 threads to perform the matrix-vector product, which leads to a speedup of up to 4.13× compared to the sequential method (see the figure 6a).

In the figure 7a, the method is directly tested with the Mooney-Rivlin material using the implementation of the MJED [30] provided in SOFA, without any modification of the code. The main difference with the co-rotational formulation lies in the fact that the computation of the hyperelastic formulation is significantly slower, and thus the time spent in the assembling and solving processes is smaller. Therefore, the benefits of the CPU parallelization with 8 threads (**CG CPU fast assembly**) is balanced by the overhead of assembling the matrix compared to the matrix-free version (**CG CPU matrix-free**). However, the GPU-based internal parallelization of the assembling and solving process provides a speedup between 1.31× and 2.05×. This represents the fastest method for nonlinear materials available in SOFA because no specific GPU-based parallelization of the MJED method is available. The method is also tested with the St Venant-Kirchhoff model using the MJED implementation. Compared to the Mooney-Rivlin material, the model is less complex so that the computation of the hyperelastic formulation is less costly. In the figure 7b, the fast assembly method gains a speedup between 1.60× and 3.34× compared to the matrix-free method, which is the fastest current implementation for the nonlinear model in SOFA.
<figure><img src="image_8.png" /><figcaption>Fig. 6: Computation time of a single time step for different examples modeled with the co-rotational formulation for various (fixed) number of CG iterations. The figures share the same legend on the top. The details of different examples can be found in Table 3. We note that the y axis is logarithmic in these figures.</figcaption></figure>

**7.3. Performances with the preconditioned CG solver.** The performances of different sparse LDL$^T$ solvers (including both the lower and upper triangular systems) are reported in the table 4. The proposed GPU-based parallelization relying on the nested dissection method (GPU ND) introduced in section 5 is 20.3 – 24.0× faster than the GPU-based implementation provided in NVIDIA’s CUSPARSE library. The main reasons lie in the fact that the CUSPARSE method requires performing the analysis of the data dependencies before actually solving the problem, and the parallelization strategies are optimal for much larger problems than the ones used in the context of real-time simulations. Such speedup compared to the golden-standard implementation (CUSPARSE library) is reported as maximally 5.8× in [35] and 19.5× in [38].

The method is also compared with the GPU-based implementation proposed in [14].
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>(a) Mooney-Rivlin model</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) St-Venant-Kichhoff model</figcaption></figure>

Fig. 7: Computation time of a single simulation step with hyperelastic models implemented in SOFA for various (fixed) number of CG iterations per time step. The figures share the same legend on the top. The scenarios simulate a tetrahedron raptor mesh (see Table 3).

As reported in this previous work, the GPU-based $LDL^T$ solver is 3x slower than a sequential CPU implementation, whereas the (GPU ND) provides a speedup of 1.4-2x, enabling the possibility to solve the problem directly on the GPU.

<table><caption>Table 4: Computation time (in ms) of various STS solvers</caption><thead><tr><th>Mesh</th><th>Method</th><th>LDL solver</th><th>Lower</th><th>Upper</th></tr></thead><tbody><tr><td rowspan="4">Raptor 1</td><td>CUSPARSE</td><td>13.46</td><td>3.33</td><td>2.48</td></tr><tr><td>[14] GPU</td><td>3.63</td><td>1.88</td><td>1.66</td></tr><tr><td>CPU</td><td>1.13</td><td>0.52</td><td>0.58</td></tr><tr><td>GPU ND</td><td>0.56</td><td>0.29</td><td>0.24</td></tr><tr><td rowspan="4">Raptor 2</td><td>CUSPARSE</td><td>22.77</td><td>5.63</td><td>3.91</td></tr><tr><td>[14] GPU</td><td>6.36</td><td>3.18</td><td>2.78</td></tr><tr><td>CPU</td><td>1.97</td><td>0.90</td><td>1.04</td></tr><tr><td>GPU ND</td><td>1.12</td><td>0.60</td><td>0.49</td></tr><tr><td rowspan="4">Raptor 3</td><td>CUSPARSE</td><td>44.96</td><td>11.02</td><td>6.97</td></tr><tr><td>[14] GPU</td><td>10.07</td><td>5.33</td><td>4.71</td></tr><tr><td>CPU</td><td>3.94</td><td>1.77</td><td>2.14</td></tr><tr><td>GPU ND</td><td>2.15</td><td>1.07</td><td>1.05</td></tr></tbody></table>

The method is tested in complete simulations of deformable bodies with various constitutive laws (see Table 5). The tests are conducted with the same mesh group of raptors and solved with the asynchronous preconditioned CG. On average, the asynchronous preconditioner is updated every 2 to 4 simulation steps, which lead between 5 to 20 iterations (#it) according to different cases. Therefore, the asynchronous preconditioner already provides a significant speedup with respect to the standard GPU-based CG. With
<table><caption>Table 5: Computation time (in ms) for various models: Corotational (Corot), Mooney-Rivlin (MR) and St-Venant-Kichhoff (SVK).</caption><thead><tr><th rowspan="2">#it</th><th rowspan="2"></th><th colspan="3">Method</th><th rowspan="2">Raptor 1</th><th rowspan="2">Raptor 2</th><th rowspan="2">Raptor 3</th></tr><tr><th>Assembly</th><th>CG</th><th>Precond</th></tr></thead><tbody><tr><td rowspan="3">Corot</td><td rowspan="3">15</td><td>Fast Assembly</td><td>CPU</td><td>CPU</td><td>25.94</td><td>45.85</td><td>85.86</td></tr><tr><td><strong>Fast Assembly</strong></td><td><strong>GPU</strong></td><td><strong>GPU</strong></td><td><strong>14.58</strong></td><td><strong>26.46</strong></td><td><strong>46.84</strong></td></tr><tr><td>Matrix Free</td><td>CPU</td><td>CPU</td><td>43.18</td><td>65.65</td><td>118.33</td></tr><tr><td></td><td></td><td>Matrix Free</td><td>GPU</td><td>CPU</td><td>31.00</td><td>48.52</td><td>80.50</td></tr><tr><td rowspan="3">MR</td><td rowspan="3">12</td><td>Fast Assembly</td><td>CPU</td><td>CPU</td><td>59.02</td><td>90.64</td><td>153.90</td></tr><tr><td><strong>Fast Assembly</strong></td><td><strong>GPU</strong></td><td><strong>GPU</strong></td><td><strong>48.73</strong></td><td><strong>77.84</strong></td><td><strong>124.53</strong></td></tr><tr><td>Matrix Free</td><td>CPU</td><td>CPU</td><td>56.76</td><td>89.76</td><td>152.11</td></tr><tr><td rowspan="3">SVK</td><td rowspan="3">8</td><td>Fast Assembly</td><td>CPU</td><td>CPU</td><td>23.16</td><td>37.96</td><td>67.03</td></tr><tr><td><strong>Fast Assembly</strong></td><td><strong>GPU</strong></td><td><strong>GPU</strong></td><td><strong>16.60</strong></td><td><strong>26.46</strong></td><td><strong>42.86</strong></td></tr><tr><td>Matrix Free</td><td>CPU</td><td>CPU</td><td>24.18</td><td>37.35</td><td>64.67</td></tr></tbody></table>

the preconditioner, the matrix operations needed during the CG iterations are performed either with the fast assembly method or with the matrix-free method, either on the CPU or the GPU when available. The preconditioner is explicitly built using the fast assembly method and applied on the CPU as done in [14] or on the GPU with the method introduced in section 5.

The method **fast assembly + CG GPU + preconditioner GPU** is the fastest method and provides a speedup of between 1.7× and 2.1× for the co-rotational model compared to the solution proposed in [14]. An important advantage of the current solution is that the preconditioned CG is applied entirely on the GPU, without any need for data transfers or synchronizations between the CPU/GPU during the solving stage. In addition, since the matrix is assembled every time step, no additional overhead is introduced when the factorization needs to be recomputed.

As no GPU-based matrix-free method is implemented for hyperelastic models in SOFA, the comparison is made with the CG performed on the CPU. Although the computation cost of the hyperelastic formulation is significantly higher, the result of the proposed GPU version also provides a speedup from 1.15× to 1.22× for the Mooney-Rivlin material. For the St Venant-Kirchhoff material, where the model is more straightforward than the Mooney-Rivlin material, this speedup is raised from 1.41× to 1.51×.

**7.4. Contact and interactions.** The fast assembly method is applied in various simulations with interactions such as simulating an object falling on rigid planes and needle insertion in soft bodies. Interactions are based on Lagrange multipliers. The compliance matrix W is built using the asynchronous preconditioner, and constraints are solved using the GPU-based solver introduced in [14]. However, these simulations still benefit from the fast assembly method for the *free motion* and *motion correction*.

According to different scenarios, the speedup for the whole simulation depends on the ratio between the computation time for the *constraint resolution* and the rest of the simulation. For example, the simulation of the colliding deformable liver (see Figure 8a) involves 246 constraints per step on average. The *constraint resolution* stage takes 56.7% of the all computation time of a single time step while the *free motion* represents only 22.5%.

Although the benefits of the fast assembly method are provided only when the fill
<figure><img src="image_12.png" /><figcaption>Fig. 8: Simulation including contacts and interactions.</figcaption></figure>

ordering of the matrix remains constant, complex interaction can still be simulated. Figure 8b shows a complex needle insertion in a soft body with complex interaction involving friction. Lagrangian multipliers are dynamically added to constrain the relative displacement of the needle and the cube when the needle penetrates the soft structure. Despite the dynamic nature of the scene, the computation of the mapping is only performed at the initial step. The needle insertion test into human organs is simulated in a heterogeneous scenario where the liver presents a co-rotational model while the skin covering the liver is modeled with Mooney-Rivlin material (see 8c). The diaphragm and the intestine are modeled with hexahedral elements, whereas the liver and the skin are composed of tetrahedra. Therefore, the method proves its compatibility with the contact problem and significant flexibility to different elements and materials. For these scenarios with complex interactions, similar speedup as shown previously in Table 5 is observed during the free motion.

**8. Conclusion.** This paper introduces a framework for real-time finite element simulations. Besides its efficiency, our method remains generic for different constitutive models. We propose a new matrix assembly strategy, which gains a significant speedup when the topology structure remains invariant and keeps the building cost on the same level as standard methods when the matrix pattern needs to be rebuilt. The fast matrix assembly gives a possibility for parallelizations in the solving stage without any specific parallel implementation of the constitutive model. Moreover, we replace the CPU-based preconditioner with a new GPU-based implementation in the solving stage. This improvement significantly reduces the data transfer between CPU/GPU and makes running a fully GPU-based CG solver possible. Finally, we evaluate our matrix assembly and parallelization strategy in various examples, including different element types and constitutive models. Our approach is also proven to be compatible with contact problems. We hope that our work will help researchers and engineers improve the performance of their works on FE simulations.

## REFERENCES

[1] J. ALLARD, H. COURTEUISSE, AND F. FAURE, *Implicit FEM and fluid coupling on GPU for interactive multiphysics simulation*, in ACM SIGGRAPH 2011 Talks, SIGGRAPH’11, New York, New York, USA, 2011, ACM Press, p. 1, https://doi.org/10.1145/2037826.2037895, http://dl.acm.org/citation.cfm?doid=2037826.2037895.

[2] H. ANZT, E. CHOW, AND J. DONGARRA, *Iterative sparse triangular solves for preconditioning*, Lecture Notes in Computer Science (including subseries Lecture Notes in
Artificial Intelligence and Lecture Notes in Bioinformatics), 9233 (2015), pp. 650-661, https://doi.org/10.1007/978-3-662-48096-0-50/FIGURES/5, https://link.springer.com/chapter/10.1007/978-3-662-48096-0-50.

[3] D. BARAFF, *Linear-time dynamics using Lagrange multipliers*, in Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH 1996, ACM, 1996, pp. 137-146, https://doi.org/10.1145/237170.237226.

[4] D. BARAFF AND A. WITKIN, *Large steps in cloth simulation*, Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH 1998, (1998), pp. 43-54, https://doi.org/10.1145/280814.280821, http://portal.acm.org/citation.cfm?doid=280814.280821.

[5] J. BARBIĆ AND D. JAMES, *Real-time subspace integration for St.Venant-Kirchhoff deformable models*, ACM Transactions on Graphics, 24 (2005), pp. 982-990, https://doi.org/10.1145/1073204.1073300.

[6] N. BELL AND M. GARLAND, *Implementing sparse matrix-vector multiplication on throughput-oriented processors*, in Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, SC '09, New York, New York, USA, 2009, ACM Press, p. 1, https://doi.org/10.1145/1654059.1654078, http://dl.acm.org/citation.cfm?doid=1654059.1654078.

[7] J. BOLZ, I. FARMER, E. GRINSPUN, AND P. SCHRÖDER, *Sparse matrix solvers on the GPU: Conjugate gradients and multigrid*, in ACM SIGGRAPH 2005 Courses, SIGGRAPH 2005, 2005, https://doi.org/10.1145/1198555.1198781.

[8] N. BOURAGO AND V. KUKUDZHANOV, *A review of contact algorithms*, Mechanics of solids, 40 (2005), pp. 35-71.

[9] A. M. BRADLEY, *A hybrid multithreaded direct sparse triangular solver*, 2016 Proceedings of the Seventh SIAM Workshop on Combinatorial Scientific Computing, (2016), pp. 13-22, https://doi.org/10.1137/1.9781611974690.CH2, https://epubs.siam.org/terms-privacy.

[10] M. BRO-NIELSEN AND S. COTIN, *Real-time volumetric deformable models for surgery simulation using finite elements and condensation*, Computer Graphics Forum, 15 (1996), pp. 57-66.

[11] L. BUATOIS, G. CAUMON, AND B. LÉVY, *Concurrent number cruncher: a GPU implementation of a general sparse linear solver*, International Journal of Parallel, Emergent and Distributed Systems, 24 (2009), pp. 205-223, https://doi.org/10.1080/17445760802337010.

[12] O. COMAS, Z. A. TAYLOR, J. ALLARD, S. OURSELIN, S. COTIN, AND J. PASSENGER, *Efficient Nonlinear FEM for Soft Tissue Modelling and Its GPU Implementation within the Open Source Framework SOFA*, in Biomedical Simulation, vol. 5104 LNCS, 2008, pp. 28-39, https://doi.org/10.1007/978-3-540-70521-5-4, http://link.springer.com/10.1007/978-3-540-70521-5-4.

[13] H. COURTEUCISSE, J. ALLARD, C. DURIEZ, AND S. COTIN, *Asynchronous preconditioners for efficient solving of non-linear deformations*, in VRIPHYS 2010 - 7th Workshop on Virtual Reality Interactions and Physical Simulations, nov 2010, pp. 59-68, https://doi.org/10.2312/PE/vriphys/vriphys10/059-068.

[14] H. COURTEUCISSE, J. ALLARD, P. KERFRIDEN, S. P. BORDAS, S. COTIN, AND C. DURIEZ, *Real-time simulation of contact and cutting of heterogeneous soft-tissues*, Medical Image Analysis, 18 (2014), pp. 394-410, https://doi.org/10.1016/j.media.2013.11.001.

[15] C. DURIEZ, C. ANDRIOT, AND A. KHEDDAR, *Signorini's contact model for deformable objects in haptic simulations*, IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 4 (2004), pp. 3232-3237, https://doi.org/10.1109/IROS.2004.1389915, http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1389915.

[16] C. DURIEZ, F. DUBOIS, A. KHEDDAR, AND C. ANDRIOT, *Realistic haptic rendering of interacting deformable objects in virtual environments*, IEEE Transactions on Visualization and Computer Graphics, 12 (2006), pp. 36-47, https://doi.org/10.1109/TVCG.2006.13, http://arxiv.org/abs/0804.0561, https://arxiv.org/abs/0804.0561.

[17] A. DZIEKONSKI, P. SYPEK, A. LAMECKI, AND M. MROZOWSKI, *Finite element matrix generation on a GPU*, in Progress in Electromagnetics Research, vol. 128, Electromagnetics Academy, 2012, pp. 249-265, https://doi.org/10.2528/PIER12040301.

[18] C. FELIPPA, *A systematic approach to the element-independent corotational dynamics of finite elements*, Tech. Report January, College Of Engineeringuniversity Of Colorado, 2000, http://www.colorado.edu/engineering/cas/Felippa.d/FelippaHome.d/Publications.d/Report.CU-CAS-00-03.pdf.

[19] M. FU, A. KUNTZ, R. J. WEBSTER, AND R. ALTEROVITZ, *Safe Motion Planning for Steerable Needles Using Cost Maps Automatically Extracted from Pulmonary Images*, IEEE International Conference on Intelligent Robots and Systems, (2018), pp. 4942-4949, https://doi.org/10.1109/IROS.2018.8593407.

[20] A. GEORGE, *Nested Dissection of a Regular Finite Element Mesh*, SIAM Journal on Numerical
Analysis, 10 (1973), pp. 345-363, https://doi.org/10.1137/0710032.

[21] M. HAUTH, O. ETZMUSS, AND W. STRASSER, *Analysis of numerical methods for the simulation of deformable models*, Visual Computer, 19 (2003), pp. 581–600, https://doi.org/10.1007/s00371-003-0206-2.

[22] F. HECHT, Y. J. LEE, J. R. SHEWCHUK, AND J. F. O'BRIEN, *Updated sparse Cholesky factors for corotational elastodynamics*, ACM Transactions on Graphics, 31 (2012), pp. 1–13, https://doi.org/10.1145/2231816.2231821, http://dl.acm.org/citation.cfm?doid=2231816.2231821.

[23] P. HERHOLZ AND M. ALEXA, *Factor once: Reusing Cholesky factorizations on sub-meshes*, SIGGRAPH Asia 2018 Technical Papers, SIGGRAPH Asia 2018, 37 (2018), pp. 1–9, https://doi.org/10.1145/3272127.3275107, https://dl.acm.org/doi/10.1145/3272127.3275107.

[24] E. HERMANN, B. RAFFIN, AND F. FAURE, *Interactive Physical Simulation on Multicore Architectures*, in *Time*, 2009, pp. 1–9, https://doi.org/10.2312/EGPGV/EGPGV09/001-008, http://tel.archives-ouvertes.fr/docs/00/36/01/31/PDF/egpgv09.pdf.

[25] R. R. HIEMSTRA, G. SANGALLI, M. TANI, F. CALABRÒ, AND T. J. HUGHES, *Fast formation and assembly of finite element matrices with application to isogeometric linear elasticity*, Computer Methods in Applied Mechanics and Engineering, 355 (2019), pp. 234–260, https://doi.org/10.1016/j.cma.2019.06.020.

[26] G. R. JOLDES, A. WITTEK, AND K. MILLER, *Suite of finite element algorithms for accurate computation of soft tissue deformation for surgical simulation*, Medical Image Analysis, 13 (2009), pp. 912–919, https://doi.org/10.1016/j.media.2008.12.001, http://dx.doi.org/10.1016/j.media.2008.12.001.

[27] G. KARYPIS AND V. KUMAR, *METIS*\* A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices, Manual, (1998), pp. 1–44, https://www.researchgate.net/publication/254424775_METIS-A_Software_Package_for_Partitioning_Unstructured_Graphs_Partitioning_Meshes_and_Computing_Fill-Reducing_Ordering_of_Sparse_Matrices.

[28] U. KÜHNAPFEL, H. K. ÇAKMAK, AND H. MAASS, *Endoscopic surgery training using virtual reality and deformable tissue simulation*, Computers and Graphics (Pergamon), 24 (2000), pp. 671–682, https://doi.org/10.1016/S0097-8493(00)00070-4, https://linkinghub.elsevier.com/retrieve/pii/S009784930000704.

[29] R. LI AND C. ZHANG, *Efficient parallel implementations of sparse triangular solves for gpu architectures*, Proceedings of the 2020 SIAM Conference on Parallel Processing for Scientific Computing, (2020), pp. 106–117, https://doi.org/10.1137/1.9781611976137.10.

[30] S. MARCHESSEAU, T. HEIMANN, S. CHATELIN, R. WILLINGER, AND H. DELINGETTE, *Multiplicative Jacobian Energy Decomposition Method for Fast Porous Visco-Hyperelastic Soft Tissue Model*, in Lecture notes in computer science, vol. 6361, Springer, 2010, pp. 235–242, https://doi.org/10.1007/978-3-642-15705-9_29, http://link.springer.com/10.1007/978-3-642-15705-9_29.

[31] J. MARTÍNEZ-FRUTOS, P. J. MARTÍNEZ-CASTEJÓN, AND D. HERRERO-PÉREZ, *Fine-grained GPU implementation of assembly-free iterative solver for finite element problems*, Computers and Structures, 157 (2015), pp. 9–18, https://doi.org/10.1016/j.compstruc.2015.05.010.

[32] F. MORIN, H. COURTEUISSE, I. REINERTSEN, F. LE LANN, O. PALOMBI, Y. PAYAN, AND M. CHABANAS, *Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation*, Medical Image Analysis, (2017), https://doi.org/10.1016/j.media.2017.06.003.

[33] E. MÜLLER, X. GUO, R. SCHEICHL, AND S. SHI, *Matrix-free gpu implementation of a preconditioned conjugate gradient solver for anisotropic elliptic pdes*, Computing and Visualization in Science, 16 (2013), pp. 41–58, https://doi.org/10.1007/s00791-014-0223-x, http://arxiv.org/abs/1302.7193, https://arxiv.org/abs/1302.7193.

[34] E. G. PARKER AND J. F. O'BRIEN, *Real-time deformation and fracture in a game environment*, in Computer Animation, Conference Proceedings, 2009, pp. 165–175, https://doi.org/10.1145/1599470.1599492.

[35] A. PICCIAU, G. E. INGGS, J. WICKERSON, E. C. KERRIGAN, AND G. A. CONSTANTINIDES, *Balancing locality and concurrency: Solving sparse triangular systems on gpus*, Proceedings - 23rd IEEE International Conference on High Performance Computing, HiPC 2016, (2017), pp. 183–192, https://doi.org/10.1109/HIPC.2016.030.

[36] Y. SAAD, *Iterative Methods for Sparse Linear Systems*, Society for Industrial and Applied Mathematics, society fo ed., 2003, https://doi.org/10.1137/1.9780898718003.

[37] S. THOMAS, I. YAMAZAKI, L. BERGER-VERGIAT, B. KELLEY, J. HU, P. MULLOWNEY, S. RAJAMANICKAM, AND K. K. KATARZYNA 'SWIRYDOWICZ, *Two-stage gauss-seidel preconditioners and smoothers for krylov solvers on a gpu cluster*, (2021), https://doi.org/10.48550/arxiv.2104.01196, https://arxiv.org/abs/2104.01196v2.
[38] I. YAMAZAKI, S. RAJAMANICKAM, AND N. ELLINGWOOD, Performance portable supernode-based sparse triangular solver for manycore architectures, ACM International Conference Proceeding Series, (2020), https://doi.org/10.1145/3404397.3404428.

[39] R. ZAYER, M. STEINBERGER, AND H. P. SEIDEL, Sparse matrix assembly on the GPU through multiplication patterns, in 2017 IEEE High Performance Extreme Computing Conference, HPEC 2017, 2017, https://doi.org/10.1109/HPEC.2017.8091057.

[40] J. ZHANG, Y. ZHONG, AND C. GU, Deformable models for surgical simulation: A survey, IEEE Reviews in Biomedical Engineering, 11 (2018), pp. 143–164, https://doi.org/10.1109/RBME.2017.2773521.