# Efficient mesh deformation based on randomized RBF solvers

W. Bader<fnref n="1" />,<fnref n="2" />, A. Parret-Fréaud<fnref n="2" />, S. Da Veiga<fnref n="3" /> and Y. Mesri<fnref n="1" />

## Abstract

Mesh deformation methods [8] have been widely used for the past decades in various fields such as fluid-structure interaction, aerodynamic shape optimization, unsteady and aeroelastic computational fluid dynamics. Such methods are particularly interesting in order to update meshes during a simulation without the need to perform an (often expensive) full regeneration of the mesh, e.g. when facing moving boundaries or geometry update during a structural optimization loop.

Among the numerous existing methods, radial basis functions interpolation (RBF) [1] is particularly suitable for unstructured mesh applications due to its simplicity and the high quality of the resulting mesh. One key aspect of RBF-based mesh deformation is the resolution of a dense linear system, which tends to be computationally expensive and high memory demanding when dealing with large-scale meshes [2, 3], thus being a major drawback of the method. This could be mitigated using an iterative solver instead of a direct one during the resolution step, thus saving the memory needed to store the factorization. However, some radial basis functions lead to ill-conditioned systems, requiring the use of an efficient preconditioner which tends to complexify the problem.

In this work, we aim to speed-up the resolution of this linear system using alternative randomization techniques coming from probabilistic linear algebra to solve the as-

<footnote n="1">¹ CEMEF – Centre for material forming, MINES ParisTech, PSL, Sophia-Antipolis, France.</footnote>

<footnote n="2">² Safran Tech, Digital Sciences Technologies Department, Rue des Jeunes Bois, Châteaufort, 78114 Magny-Les-Hameaux, France.</footnote>

<footnote n="3">³ Ensai, CREST - UMR 9194, Rennes, France.</footnote>
sociated dense linear system. Indeed, such methods have been studied for two decades
and are being increasingly popular in various fields, including numerical linear algebra
and optimization [5]. Their key aspect is to reduce the complexity of solving large
scale linear systems by exploiting the spectral properties of the underlying operator.

In this study, we propose an alternative approach for dealing with the input matrix
by generating an approximate "sketch" of the initial problem. This sketch is easier
to solve compared to working with the original matrix directly, albeit at the expense
of reduced precision. Our focus lies specifically on matrices arising from RBF-based
mesh deformation procedures, which typically exhibit a rapid spectral decay and tend
to be numerically low rank. Leveraging these characteristics, we explore the potential
of probabilistic linear algebra techniques in this domain.

To embed the rows of the linear system into a lower-dimensional space while pre-
serving their geometric properties, we employ a dimension reduction map. By doing
so, we maintain the underlying geometry of the original space, thereby ensuring that
the approximated sketch exhibits similar behavior in terms of singular values and sin-
gular vectors as the original matrix. Our chosen method for constructing this map
involves utilizing highly structured random matrices, commonly referred to as ran-
domized linear embeddings or random projections. Subsequently, we confine the ma-
trix to the approximated subspace and compute a standard factorization of the reduced
matrix.

The proposed approach will be discussed on the basis of 2D and 3D applications.

**Keywords:** Mesh Deformation, RBF, Randomized Linear Algebra, Subspace embeddings, Low rank approximation, Large scale Systems

# 1 Introduction

In this paper, our primary objective is to significantly enhance the efficiency of solving the linear system involved in RBF-based mesh deformation approaches. We achieve this by employing alternative randomization techniques from probabilistic linear algebra, which have shown great promise in the context of RBF-based mesh deformation procedures. This is particularly relevant because a significant number of matrices encountered in such applications exhibit rapid spectral decay and tend to have low numerical rank.

To address this challenge, we propose a novel method that involves embedding the
rows of the linear system into a lower-dimensional space using a dimension reduc-
tion map, while ensuring the preservation of their geometric properties. By preserving
the initial space's geometry, we can demonstrate that the approximated "sketch" ex-
hibits similar behavior in terms of singular values and singular vectors compared to
the original matrix.

This paper is organised as follow. Section 2 shortly recalls the basis of mesh de-
formation based on RBF and derive the associated (dense) linear system. Section 3
introduce the application of randomized linear algebra on the previously obtained sys-
tem, focusing on the randomized singular value decomposition (RSVD) algorithm. Section 4 investigates the performance of the proposed method on the basis of a mesh moving application. Finally, Section 5 concludes this paper.

## 2 Mesh deformation based on Radial Basis Functions

Let's consider a mesh $\Omega_h$ with an interface $\Gamma_h^c$ containing $n_c$ nodes subjected to a given displacement $d_c$. In the following, we denote by $x$ the coordinates of any (internal) node of the mesh, and by $(x_c^j)_j$ the coordinates of the interface control nodes subjected to the given displacement. The main principle of RBF-based mesh deformation is to infer the displacement of any internal nodes $x$ given the known displacement of the control nodes on the interface $\Gamma_h^c$ by means of the following interpolation function $s$:

$$s(x) = \sum_{j=1}^{n_c} \alpha_j \phi(||x - x_c^j||) + p(x), \quad (1)$$

where $\phi$ denotes a radial basis function to be chosen, $(\alpha_j)_j$ the associated weights and $p$ a polynomial of degree depending on the choice of $\phi$. The weights $\alpha_j$ and the polynomial $p$ are computed such that the function $s$ returns the exact displacement at control points:

$$s(x_c^j) = d_c^j, \quad 1 \le j \le n_c, \quad (2)$$

under the additional constraint:

$$\sum_{j=1}^{n_c} \alpha_j q(x_c^j) = 0, \quad (3)$$

for all polynomials $q$ with a degree less or equal to the degree of $p$.

A common choice, which holds in the following, is to choose $s$ to be conditionally positive definite of order $m \le 2$ together with a linear polynomial $p$. Examples of such a function $s$ is given in Table (1).

<table><caption>Table 1: Examples of classical radial basis functions.</caption><thead><tr><th>Name</th><th>s(x)</th></tr></thead><tbody><tr><td>CP C<sup>0</sup></td><td>(1 - x)<sup>2</sup></td></tr><tr><td>CP C<sup>2</sup></td><td>(1 - x)<sup>4</sup>(4x + 1)</td></tr><tr><td>Gaussian</td><td>10<sup>-x<sup>2</sup></sup></td></tr></tbody></table>

Finally, Equations (2) and (3) lead to the following linear system of unknowns $\alpha_i$:

$$Ax = b, \quad \text{with} \quad A = \begin{bmatrix} M & P \\ P^T & 0 \end{bmatrix}, \quad x = \begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \quad b = \begin{bmatrix} d_c \\ 0 \end{bmatrix}, \quad (4)$$

where $M \in \mathbb{R}^{n_c \times n_c}$ is the interpolation matrix of coefficients $M_{i,j} = \phi(||x_c^i - x_c^j||)$ and $P \in \mathbb{R}^{n_c \times 4}$ is given by $P_{j,:} = [1, x_c^i, y_c^i, z_c^i]$.
### 3 Application of randomized linear algebra to RBF linear system

As stated before, applying both direct or iterative methods coming from deterministic linear algebra to the linear system (4) tends to be computationally expensive and memory demanding on large scale configurations. In this section, we present the application of one of the classical randomized linear algebra algorithms to the previous linear system.

Instead of working with the input matrix directly, one of the basic ideas of randomized linear algebra is to work with a "sketch" that approximates the initial problem. As the majority of RBF matrices exhibit a rapid spectral decay and tend to be numerically low rank, it is possible to embed the rows of the linear system into a lower-dimensional space using a dimension reduction map while conserving their geometrical properties. By preserving the geometry of the initial space, we can prove that the approximated "sketch" has similar behavior in terms of singular values and singular vectors compared to the original matrix.

The proposed method is conducted in two steps. First, we use highly structured random matrices to build such map which will be referred to as randomized linear embeddings, also known as random projection. Then, we restrict the matrix to the approximated subspace and compute a standard factorization of the reduced matrix.

More specifically, we start by randomly projecting the $n_c \times n_c$ RBF matrix $A$ using the randomized map to form the matrix $Y = A\Omega$, where $\Omega^*$ is a mixing random embedding such that the columns $(y_i)_{i \in [1,l]}$ of the matrix $Y$ form a random sampling from the range of $A$, where $l \ll n_c$.

For the selection of a random embedding matrix, it is recommended by the authors [5] to use a Gaussian random matrix. In this approach, each entry of the matrix is sampled from a Gaussian distribution. Another viable method that can be employed is the Subsampled Randomized Hadamard Transform (SRHT), denoted as $\Omega = \sqrt{\frac{n_c}{l}} DHR$. This transformation encompasses the following components:

* $H$: The Walsh Hadamard matrix with dimensions $n_c \times n_c$, which is deterministically defined.

* $D$: A diagonal $n_c \times n_c$ matrix with entries independently and uniformly distributed over $1, +1$.

* $R$: A $n_c \times l$ matrix obtained by randomly selecting $l$ columns from the $n_c \times n_c$ identity matrix.

The SRHT presents the advantage of requiring only $O(n\log(n))$ storage space and $O(k\log(n))$ operations to be applied to a vector. This efficiency is achieved through the utilization of a fast subsampled trigonometric transform algorithm.

Due to the randomness, the columns of the test matrix $\Omega$ are likely to be linearly independent and no linear combination falls in the null space of $A$. This means that
the columns of $Y$ are also linearly independent and span a subspace of the range of $A$. Next, a QR factorization may be applied to $Y$ in order to produce an orthonormal basis $Q$. However, as the columns of $Y$ tend to be strongly aligned in practice (almost linearly dependent), it is recommended [5] to use a numerically stable orthogonalization procedure [7] as Householder reflectors, double Gram-Schmidt, rank-revealing QR or TSQR algorithm. Finally, we restrict the input matrix $A$ to the approximated subspace in order to obtain a reduced $l \times n_c$ matrix $B$ and use a standard factorization method such as the singular value decomposition (SVD) on this last matrix. Algorithm 1 shows the different steps to compute the randomized approximation to the SVD factors.

<caption>**Algorithm 1:** Randomized singular value decomposition (RSVD) [5]</caption>

**Input:** $A \in R^{n_c \times n_c}$, $k$ a factorization rank and $p$ an oversampling parameter
**Output:** $U \in R^{n_c \times k}$, $V \in R^{k \times n_c}$ orthonormal matrices and $\Sigma \in R^{k \times k}$ a diagonal matrix
such that $A \approx U\Sigma V^T$
1. Draw a random $n_c \times (k+p)$ test matrix $\Omega$;
2. $Y = A\Omega$;
3. $Y = QR$; // QR Factorization
4. $B = Q^T A$;
5. $[\hat{U}, \Sigma, V] = svd_econ(B)$;
6. $U = Q\hat{U}$;

For a given target rank $k$, in order to improve the quality of the corresponding approximation, it is advised to add an oversampling parameter $p \ge 2$ such that $l = k + p \le n_c$. Then, for $\Omega \in R^{n_c \times (k+p)}$ a standard normal test matrix, the expectation of the error of the randomized SVD satisfies the following inequality [6]:

$$E(||A - U\Sigma V^T||) \le (1 + 4\frac{\sqrt{n_c(k+p)}}{p-1})\sigma_{k+1} \quad (5)$$

Consequently, the random projection procedure computes a $(k+p)$-dimensional subspace that captures as much of the action of a matrix $A$ as the best $k$-dimensional subspace to a small polynomial term, hence the reason behind the oversampling operation ($p = 5$ or $p = 10$). Due to the numerical low rank nature of the RBF matrices, we are expecting to have an accurate approximation even for moderate values of $k$.

With $T_{mult}$ the cost of a matrix-vector multiplication, the number of flops $T_{basic}$ required by Algorithm 1 satisfies:

$$T_{basic} = (k + p)T_{mult} + O(n_c k^2), \quad (6)$$

which mean that this algorithm can benefit from matrices with a fast matrix-vector product (structured, sparse, ...).

Finally, once obtaining the approximated factors, we use them to get an approximated solution to the linear system (4) coming from the RBF method.
## 4 Numerical assessments

In order to evaluate the performance of randomized mesh deformation approach based on Equation and Algorithm 1, we consider a test case constituted by a NACA-0012 airfoil geometry (Figure 1 - a) subjected to a 45-degree angular rotation around its center. The mesh deformation is performed in 10 steps between the initial mesh, located at the center of the domain, and the final configuration. A Winland (CP C²) RBF function is used for the interpolation together with several random projectors and compression ratio values.

We evaluate the performance of the approach in terms of resulting mesh quality and approximation error with respect to a reference deterministic direct factorization method on the matrix A. The numerical experiment are computed using Python programming language with standard scientific computing libraries such as NumPy and SciPy.

In order to evaluate the quality of the deformed mesh, we propose to use the following indicator on each element E [4]:

$$Q(E) = C_0 \frac{V_E}{h^d}, \quad \text{with} \quad h = \left( \frac{2}{d(d+1)} \sum_{i<j} ||\partial E_i - \partial E_j||^2 \right)^{\frac{1}{2}}, \quad (7)$$

where d is the space dimension, $V_E$ the volume of the element E, $(\partial E_i)_i$ the edges of element E, h the average of lengths of edges and $C_0$ is chosen such that the quality measurement of an equilateral element is equal to 1. As for the accuracy of the approximated solution of the linear system, we use the relative residual error e with respect to the spectral norm:

$$e(x) = \frac{||Ax - b||}{||b||}. \quad (8)$$

Finally, we also estimate the expected error and its standard deviation over 50 trials in order to evaluate the reproducibility of the randomized solution.

Figure (1 - b) represents the final mesh using an exact solver, Figure (1 - c) (respectively (1 - d)) represents the deformed mesh for the same number of steps but using a randomized direct solver with Gaussian (respectively Hadamard) random projectors, while keeping only 10% of the RBF matrix columns (which corresponds to $l = 0.1n_c$). The randomized mesh deformation approach had a similar behavior in terms of mesh quality distribution to the exact solver with an expected relative error up to $10^{-2}$ as shown by the Figure (2).

By decreasing the compression ratio, we were able to improve the accuracy of the randomized solver with an estimated error up to $10^{-4}$ as shown in Table (2), but no significant improvement was noticed in terms of mesh quality distribution.
<figure><img src="image_1.png" /><figcaption>Figure 1: NACA-0012 airfoil geometry : (a) initial mesh, (b) deformed mesh with exact solver, (c) deformed mesh with randomized solver using Gaussian projector, (d) deformed mesh with randomized solver using Hadamard projector.</figcaption></figure>

<table><caption>Table 2: Relative residual for different sample size</caption><thead><tr><th rowspan="2">Sample size %</th><th colspan="2">mean</th><th colspan="2">std</th></tr><tr><th>x</th><th>y</th><th>x</th><th>y</th></tr></thead><tbody><tr><td>0.1</td><td>2.72 10<sup>-2</sup></td><td>2.54 10<sup>-2</sup></td><td>2.72 10<sup>-5</sup></td><td>2.69 10<sup>-5</sup></td></tr><tr><td>0.3</td><td>1.82 10<sup>-3</sup></td><td>1.68 10<sup>-3</sup></td><td>6.08 10<sup>-8</sup></td><td>5.41 10<sup>-8</sup></td></tr><tr><td>0.5</td><td>4.74 10<sup>-4</sup></td><td>4.45 10<sup>-4</sup></td><td>2.73 10<sup>-9</sup></td><td>2.31 10<sup>-9</sup></td></tr></tbody></table>

## 5 Conclusion

This work presents an application of randomized linear algebra algorithms in order to improve the efficiency of RBF-based mesh deformation procedures. By exploiting the specificity of the spectral properties of RBF matrices, these randomized approaches allows to drastically reduce the computational complexity associated with the RBF linear system, while preserving sufficient accuracy for this class of problems. Indeed, the numerical experiments conducted using two different flavors of randomized SVD algorithm exhibit good quality of the resulting deformed meshes in comparison with a deterministic solver, even with high compression ratios on the RBF operator.

Current works focus on other types of random solvers, including adaptive variants, for a wider range of applications.
<figure><img src="image_2.png" /><figcaption>Figure 2: Comparison between the mesh quality distribution *Q* of the deterministic and randomized approaches while keeping only 10% of the RBF matrix columns.</figcaption></figure>

References

[1] De Boer, A., Van der Schoot, M. S. and Bijl, H. (2007). Mesh deformation based on radial basis function interpolation. *Computers & structures*, 85(11-14), 784-795.

[2] El Gharbi, Y., Parret-Fréaud, A., Bovet, C. and Gosselet, P. (2021). Two-level substructuring and parallel mesh generation for domain decomposition methods. *Finite Elements in Analysis and Design*, 192, 103484.

[3] Daldoul, W., Hachem, E. and Mesri, Y. (2020). A ‘R-to-H’ Mesh Adaptation Approach for Moving Immersed Complex Geometries Using Parallel Computers. *International Journal of Computational Fluid Dynamics*, 7-8, 597–609.

[4] Mesri Y., Guillard H., Coupez T. (2012). Automatic coarsening of three dimensional anisotropic unstructured meshes for multigrid applications. *Appl Math Comput.*, 218(21):10500-10519.

[5] Martinsson, P. G. and Tropp, J. A. (2020). Randomized numerical linear algebra: Foundations and algorithms. *Acta Numerica*, 29, 403-572.

[6] Halko, N., Martinsson, P. G., Tropp, J. A. (2011). Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. *SIAM review*, 53(2), 217-288.

[7] Björck, Å. (1994). Numerics of gram-schmidt orthogonalization. *Linear Algebra and Its Applications*, 197, 297-316.

[8] Selim, M. M., & Koomullil, R. P. (2016). Mesh deformation approaches—a survey. *Journal of Physical Mathematics*, 7(2), 1-9.