Scalable multigrid methods for immersed finite element methods and immersed isogeometric analysis

TL;DR

This paper develops a scalable multigrid preconditioner for immersed finite element and isogeometric analysis methods, effectively addressing ill-conditioning issues and enabling efficient large-scale parallel computations.

Contribution

It introduces a geometric multigrid preconditioner that is mesh- and cut-element-independent, suitable for higher-order discretizations and various basis functions.

Findings

Provides mesh-independent convergence rates.

Enables linear-scaling parallel solutions.

Effective for higher-order and locally refined discretizations.

Abstract

Ill-conditioning of the system matrix is a well-known complication in immersed finite element methods and trimmed isogeometric analysis. Elements with small intersections with the physical domain yield problematic eigenvalues in the system matrix, which generally degrades efficiency and robustness of iterative solvers. In this contribution we investigate the spectral properties of immersed finite element systems treated by Schwarz-type methods, to establish the suitability of these as smoothers in a multigrid method. Based on this investigation we develop a geometric multigrid preconditioner for immersed finite element methods, which provides mesh-independent and cut-element-independent convergence rates. This preconditioning technique is applicable to higher-order discretizations, and enables solving large-scale immersed systems in parallel, at a computational cost that scales linearly with the number of degrees of freedom. The performance of the preconditioner is demonstrated for conventional Lagrange basis functions and for isogeometric discretizations with both uniform B-splines and locally refined approximations based on truncated hierarchical B-splines.

Figures

Equations

$$\left\{\begin{array}[]{c l}\mbox{div}\left(\boldsymbol{\sigma}\right)+\boldsymbol{f}=\boldsymbol{0}&\mbox{in }\Omega,\\ \boldsymbol{u}=\boldsymbol{g}^{D}&\mbox{on }\Gamma^{D},\\ \boldsymbol{\sigma}\boldsymbol{n}=\boldsymbol{g}^{N}&\mbox{on }\Gamma^{N},\end{array}\right.$$

$$\left\{\begin{array}[]{l}\mbox{Find }\boldsymbol{u}_{h}\in\mathcal{V}_{h}\mbox{ such that for all }\boldsymbol{v}_{h}\in\mathcal{V}_{h}\mbox{:}\\[5.0pt] \qquad a_{h}\left(\boldsymbol{v}_{h},\boldsymbol{u}_{h}\right)=b_{h}\left(v_{h}\right),\end{array}\right.$$

$$a_h\left({\mathbf{v}}_h,{\mathbf{u}}_h\right)$$

$$b_h\left({\mathbf{v}}_h\right)$$

$$\mathbf{Ax}=\mathbf{b},$$

$$\kappa \left(\mathbf{A}\right)=\mathcal{O}\left({\eta }^{-\left(2p+1-2/d\right)}\right),$$

$$\eta =\underset{K_i\in {\mathcal{T}}_h|K_i\cap \Omega \ne \varnothing }{\min }\frac{|K_i\cap \Omega |}{|K_i|}.$$

$$\kappa \left(\mathbf{A}\right)=\mathcal{O}\left(h^{-2}\right).$$

$$\mathbf{r}=\mathbf{b}-\mathbf{A}{\mathbf{x}}^{\prime }.$$

$$\bar{\mathbf{x}}={\mathbf{A}}^{-1}\mathbf{r}.$$

$${\mathbf{\Phi }}_{\ell -1}={\mathbf{R}}_{\ell }{\mathbf{\Phi }}_{\ell }.$$

$${\tilde{\mathbf{x}}}_{\ell -1}^T{\mathbf{\Phi }}_{\ell -1}={\tilde{\mathbf{x}}}_{\ell -1}^T\left({\mathbf{R}}_{\ell }{\mathbf{\Phi }}_{\ell }\right)={\left({\mathbf{R}}_{\ell }^T{\tilde{\mathbf{x}}}_{\ell -1}\right)}^T{\mathbf{\Phi }}_{\ell }$$

$$0\le {\lambda }_{\min }\left(\gamma {\mathbf{M}}^{-1}\mathbf{A}\right)\le {\lambda }_{\max }\left(\gamma {\mathbf{M}}^{-1}\mathbf{A}\right)\le 2,$$

$$\rho \left(\mathbf{I}-\gamma {\mathbf{M}}^{-1}\mathbf{A}\right)<1,$$

$$\left(\mathbf{I}-{\mathbf{M}}^{-\mathrm{T}}\mathbf{A}\right)\left(\mathbf{I}-{\mathbf{M}}^{-1}\mathbf{A}\right){\mathbf{y}}_{{\lambda }_i}=\left(1-{\lambda }_i\right){\mathbf{y}}_{{\lambda }_i},$$

$$\left(\mathbf{I}-{\mathbf{M}}_{\ell }^{-\mathrm{T}}{\mathbf{A}}_{\ell }\right)\left(\mathbf{I}-{\mathbf{R}}_{\ell }^{\mathrm{T}}{\mathbf{A}}_{\ell -1}^{-1}{\mathbf{R}}_{\ell }{\mathbf{A}}_{\ell }\right)\left(\mathbf{I}-{\mathbf{M}}_{\ell }^{-1}{\mathbf{A}}_{\ell }\right){\mathbf{y}}_{{\lambda }_i}=\left(1-{\lambda }_i\right){\mathbf{y}}_{{\lambda }_i},$$

$$\mbox{V-cycle}\big{(}\ell=2,\mathbf{r}_{\ell}=\mathbf{A}\mathbf{y}_{\lambda_{i}}\big{)}=\lambda_{i}\mathbf{y}_{\lambda_{i}}.$$

$${\mathbf{M}}^{-1}=\sum _{j=1}^N{\mathbf{P}}_j\underset{{\mathbf{A}}_j^{-1}}{\underbrace{{\left({\mathbf{P}}_j^{\mathrm{T}}\mathbf{A}{\mathbf{P}}_j\right)}^{-1}}}{\mathbf{P}}_j^{\mathrm{T}},$$

$$\left\{{\varphi }_k:{\mathrm{supp}}_{\Omega }\left({\varphi }_k\right)\subseteq {\mathrm{supp}}_{\Omega }\left({\varphi }_j\right)\right\},$$

$${\lambda }_{\max }\left({\mathbf{M}}^{-1}\mathbf{A}\right)$$

$$=\underset{\mathbf{z}}{\max }\frac{{\mathbf{z}}^{\mathrm{T}}\mathbf{Az}}{\underset{{\sum }_{j=1}^N{\mathbf{P}}_j{\mathbf{z}}_j=\mathbf{z}}{\min }{\sum }_{j=1}^N{\mathbf{z}}_j^{\mathrm{T}}{\mathbf{P}}_j^{\mathrm{T}}\mathbf{A}{\mathbf{P}}_j{\mathbf{z}}_j}$$

$$=\underset{\{{\mathbf{z}}_j{\}}_{j=1}^N}{\max }\frac{\left({\sum }_{j=1}^N{\mathbf{P}}_j^{\mathrm{T}}{\mathbf{z}}_j^{\mathrm{T}}\right)\mathbf{A}\left({\sum }_{j=1}^N{\mathbf{P}}_j{\mathbf{z}}_j\right)}{{\sum }_{j=1}^N{\mathbf{z}}_j^{\mathrm{T}}{\mathbf{P}}_j^{\mathrm{T}}\mathbf{A}{\mathbf{P}}_j{\mathbf{z}}_j}$$

$$\le \underset{K_i}{\max }\underset{\{{\mathbf{z}}_j{\}}_{j=1}^{N_{K_i}}}{\max }\frac{\left({\sum }_{j=1}^{N_{K_i}}{\mathbf{P}}_j^{\mathrm{T}}{\mathbf{z}}_j^{\mathrm{T}}\right){\mathbf{A}}^{K_i}\left({\sum }_{j=1}^{N_{K_i}}{\mathbf{P}}_j{\mathbf{z}}_j\right)}{{\sum }_{j=1}^{N_{K_i}}{\mathbf{z}}_j^{\mathrm{T}}{\mathbf{P}}_j^{\mathrm{T}}{\mathbf{A}}^{K_i}{\mathbf{P}}_j{\mathbf{z}}_j}$$

$$\le \underset{K_i}{\max }\underset{\{{\mathbf{z}}_j{\}}_{j=1}^{N_{K_i}}}{\max }\frac{N_{K_i}{\sum }_{j=1}^{N_{K_i}}{\mathbf{z}}_j^{\mathrm{T}}{\mathbf{P}}_j^{\mathrm{T}}{\mathbf{A}}^{K_i}{\mathbf{P}}_j{\mathbf{z}}_j}{{\sum }_{j=1}^{N_{K_i}}{\mathbf{z}}_j^{\mathrm{T}}{\mathbf{P}}_j^{\mathrm{T}}{\mathbf{A}}^{K_i}{\mathbf{P}}_j{\mathbf{z}}_j}=\underset{K_i}{\max }{N}_{K_i},$$

$$\delta {\tilde{\mathbf{y}}}^j$$

$${\tilde{\mathbf{y}}}^j$$

$${\tilde{\mathbf{r}}}^j$$

$$\delta {\tilde{\mathbf{z}}}^c$$

$${\tilde{\mathbf{y}}}^c$$

$${\tilde{\mathbf{r}}}^c$$

$${\psi }_1(x,y,z)=16\left(1-\sum _{i=0}^6{e}^{-{\rho }_i(x,y,z)}\right)-x^4-y^4-z^4,$$

$$\begin{array}[]{l c c c c c}\rho_{0}(x,y,z)=&x^{2}&+&y^{2}&+&(z-2)^{2},\\ \rho_{1}(x,y,z)=&(x-2)^{2}&+&(y-2)^{2}&+&\left(\frac{z+2}{2}\right)^{2},\\ \rho_{2}(x,y,z)=&(x-2)^{2}&+&(y+2)^{2}&+&\left(\frac{z+2}{2}\right)^{2},\\ \rho_{3}(x,y,z)=&(x+2)^{2}&+&(y+2)^{2}&+&\left(\frac{z+2}{2}\right)^{2},\\ \rho_{4}(x,y,z)=&(x+2)^{2}&+&(y-2)^{2}&+&\left(\frac{z+2}{2}\right)^{2},\\ \rho_{5}(x,y,z)=&x^{2}&+&\left(\frac{y}{2}\right)^{2}&+&\left(\frac{z+2}{2}\right)^{2},\\ \rho_{6}(x,y,z)=&\left(\frac{x}{2}\right)^{2}&+&y^{2}&+&\left(\frac{z+2}{2}\right)^{2}.\end{array}$$

$${\psi }_2(x,y,z)=z+1.$$

$${\mathbf{g}}^N(x,y,z)=-\mathbf{n}{e}^{-\frac{1}{4}\left((x-2{)}^2+(y-2{)}^2+(z-2{)}^2\right)},$$

$${\psi }_1(x,y,z)=1-{\left(\frac{z}{1.7}\right)}^2-{\left(\frac{r(x,y)}{1+z/17}-0.7\right)}^2,$$

Full text

††journal: Computational Mechanics

1 Introduction

Immersed methods are useful tools to avoid laborious and computationally expensive procedures for the generation of body-fitted finite element discretizations or analysis-suitable NURBS geometries in isogeometric analysis, specifically for problems on complex, moving, or implicitly defined geometries. Immersed finite element techniques, such as the finite cell method [1, 2, 3], CutFEM [4, 5] and immersogeometric analysis [6, 7], have been successfully applied to a broad range of problems. Noteworthy applications include isogeometric analysis on trimmed CAD objects, e.g., [8, 9, 10, 11, 12, 13], fluid-structure interaction with large displacements, e.g., [14, 15, 16, 17, 18, 19, 20, 21], scan based analysis [22, 23, 24, 25, 26, 27] and topology optimization, e.g., [28, 29, 30, 31, 32, 33, 34].

An essential aspect of finite element methods and isogeometric analysis is the computation of the solution to a system of equations. This is specifically challenging for systems derived from immersed methods, since such methods generally yield severely ill-conditioned system matrices [35]. For this reason, many researchers resort to direct solvers, e.g., [2, 10, 9, 11, 3, 12, 24], the efficiency of which is not affected by the conditioning of the system matrix. Nevertheless, the computational cost of direct solvers, both in terms of memory and floating point operations, scales poorly with the size of the system. With iterative solvers, the scaling between the computational cost and the size of the system is generally better, making these more suitable for large systems of equations [36]. However, the efficiency and reliability of iterative solution methods depends on the conditioning of the system. Without dedicated treatments, the severe ill-conditioning of linear systems derived from immersed finite element methods generally forestalls convergence of iterative solution procedures. Multiple resolutions for these conditioning problems have been proposed, the most prominent of which are the ghost penalty, e.g., [37, 4, 5], constraining, extending, or aggregation of basis functions, e.g., [38, 39, 40, 41, 42, 43, 13, 44, 45, 46], and preconditioning, which is discussed in detail below.

Several dedicated preconditioners have been developed for immersed finite element methods. It is demonstrated in [47] that a diagonal preconditioner in combination with the constraining of very small basis functions results in an effective treatment for systems with linear bases. With certain restrictions to the cut-element-geometry, [48] derives that a scalable preconditioner for linear bases is obtained by combining diagonal scaling of basis functions on cut elements with standard multigrid techniques for the remaining basis functions. In [49] the scaling of a BDDC (Balancing Domain Decomposition by Constraints) is tailored to cut elements, and this is demonstrated to be effective with linear basis functions. An algebraic preconditioning technique is presented in [35], which results in an effective treatment for smooth function spaces. References [50] and [51] establish that additive Schwarz preconditioners can effectively resolve the conditioning problems of immersed finite element methods with higher-order discretizations, for both isogeometric and $hp$-finite element function spaces and for both symmetric positive definite (SPD) and non-SPD problems. Furthermore, the numerical investigation in [50] conveys that the conditioning of immersed systems treated by an additive Schwarz preconditioner is very similar to that of mesh-fitting systems. In particular, the condition number of an additive Schwarz preconditioned immersed system exhibits the same mesh-size dependence as mesh-fitting approaches [52, 53], which opens the doors to the application of established concepts of multigrid preconditioning. It should be mentioned that similar conditioning problems as in immersed methods occur in XFEM and GFEM. Dedicated preconditioners have been developed for these problems as well, a survey of which can be found in [50].

Multigrid methods effectively resolve the mesh-size dependence of the conditioning of linear systems and its effect on the convergence of iterative solution methods. In particular, the use of overlapping Schwarz smoothers can lead to multigrid methods with provably mesh-independent convergence rates [54, 55, 56]. There exists a rich literature on multigrid techniques, and interested readers are directed to the reference works [57, 58, 59, 60]. Multigrid methods have not been studied extensively in the context of immersed finite element methods, but detailed studies regarding closely related aspects are available. In isogeometric analysis, multigrid is an established concept, e.g., [61, 62, 63, 64, 65, 66, 67, 68]. In regard of the present manuscript [69, 70, 71, 72] are particularly noteworthy, as these all employ smoothers that are based on Schwarz-type techniques. Another interesting contribution is [73], which presents a multigrid technique for locally refined function spaces with truncated hierarchical B-splines, that are also employed in this manuscript. Further noteworthy references are multigrid preconditioners for XFEM [74, 75], discontinuous Galerkin (DG) [76], unfitted interface problems [77], and an algebraic multigrid (AMG) preconditioner that is applied to immersed systems which have been treated by an aggregation procedure [78].

The main objective of this contribution is to develop a geometric multigrid preconditioning technique that is applicable to higher-order immersed finite element methods with conventional, isogeometric, and locally refined basis functions. This preconditioner enables iterative solution methods with a convergence rate that is unaffected by either the cut elements or the grid size. The intrinsic dependence on mesh regularity in geometric multigrid approaches is non-restrictive for immersed finite element methods, as these methods generally employ structured grids. Based on the observations regarding the mesh-size dependence of the additive Schwarz preconditioner developed in [50], this contribution further investigates the spectral properties of immersed systems treated with Schwarz-type preconditioners, in order to establish the suitability of these as smoothers in a multigrid method. This results in a preconditioning technique with the desired properties, the performance of which is demonstrated on a range of test cases with multi-million degrees of freedom. This numerical investigation includes the application to locally refined bases with truncated hierarchical B-splines, and involves a detailed study of aspects that are specific for immersed finite elements.

This contribution only considers SPD problems. Multigrid methods are an established concept in fluid mechanics as well [60], however, and similar Schwarz-type methods have successfully been applied to flow problems with both immersed finite element methods [50] and mesh-fitting multigrid solvers [71], i.e., Vanka-smoothers [79]. Therefore, it is anticipated that the presented preconditioning technique extends matatis mutandis to non-SPD and mixed formulations.

In Section 2 of this contribution the employed immersed finite element method is presented, including the quadrature on cut elements and the applied discretization spaces. Section 3 investigates spectral properties of immersed finite element methods and presents the developed geometric multigrid preconditioner. In Section 4 this preconditioning technique is assessed on a range of test cases and conclusions are drawn in Section 5.

2 Immersed finite element formulation

We consider problems on a two-dimensional or three-dimensional domain $\Omega\subset\mathbb{R}^{d}$ ($d\in\{2,3\}$), that is referred to as the physical domain. The physical domain is encapsulated by a fictitious extension, to obtain an embedding domain of simple shape, as illustrated in Figure 1. Because of its simple shape, it is trivial to generate a structured, e.g., tensor product, mesh for the embedding domain. These structured meshes render immersed finite elements ideally suitable for geometric multigrid approaches. The set of elements that intersects the physical domain is referred to as the background grid, and can serve as a substructure to construct different types of basis functions. We denote this mesh by $\mathcal{T}_{h}$, where $h$ refers to the grid size, and the basis functions that are supported on the physical domain span the approximation space $\mathcal{V}_{h}$. A non-trivial aspect of immersed finite element methods is the integration over cut elements. Herein we employ the procedure as outlined in [23], which is illustrated in Figure 1.

The numerical examples in this contribution consider problems in linear elasticity:

[TABLE]

with Cauchy stress tensor $\boldsymbol{\sigma}=\boldsymbol{\sigma}(\boldsymbol{u})=\lambda\mbox{div}\left(\boldsymbol{u}\right)\mathbf{I}+2\mu\nabla^{s}\boldsymbol{u}$, Lamé parameters $\lambda$ and $\mu$, $\nabla^{s}$ denoting the symmetric gradient operator, $\mathbf{n}$ the exterior unit normal vector and $\Gamma^{D}\cup\Gamma^{N}=\partial\Omega$ complementary parts of the boundary on which Dirichlet and Neumann conditions are prescribed. The formulations in this section are restricted to pure Dirichlet or Neumann boundary conditions, but can easily be modified to mixed boundary conditions with a prescribed normal displacement and tangential traction or Robin-type conditions. We consider approximations of (1) based on a symmetric and coercive variational form:

[TABLE]

with the bilinear and linear operators defined as:

[TABLE]

The weak formulation in (2) and (3) imposes the Dirichlet conditions by the penalty method with parameters $\beta_{h}^{\lambda}$ and $\beta_{h}^{\mu}$ that are chosen inversely proportional to the grid size. The analysis of the conditioning of immersed finite elements in [35] conveys that cut-element-specific conditioning problems are not essentially affected by the type of weak enforcement of the boundary conditions, provided that with a Nitsche-type enforcement of Dirichlet conditions [81] coercivity of the variational form is retained [82]. Therefore all computations in this contribution apply the penalty method, which bypasses the computation of local stabilization parameters. It should be noted that, due to the penalty parameters, the operators $a_{h}\left(\cdot,\cdot\right)$ and $b_{h}\left(\cdot\right)$ are mesh-size dependent. This is relevant in regard of multigrid techniques, in which the same problem is considered on different meshes. The coarse problems in this contribution inherit the weak formulation from the finest grid, without modification of the penalty parameters, see Remark 3.1 for details.

An advantageous property of immersed methods is that different approximation spaces can relatively easily be employed, by virtue of the regularity of the underlying mesh. We herein consider uniform discretizations with both traditional Lagrange basis functions and B-splines [83], and locally refined discretizations with truncated hierarchical B-splines (THB-splines) [84]. The construction of THB-splines is illustrated in Figure 2. THB-splines are posed on a hierarchy of meshes, each of which has a certain number of active B-splines. The truncated basis is then obtained by truncating active B-splines with respect to the active B-splines on the finer grids in which they are nested. We refer the reader to [84, 85] for details regarding the construction of THB-splines. Truncating the basis functions reduces the computational cost, because the truncated basis functions have smaller supports than the standard (non-truncated) hierarchical basis functions. This results in a sparser system matrix. In the context of the multigrid preconditioner developed herein, THB-splines have a particularly important advantage concerning the computational cost. As will be elaborated in Section 3.4, the block selection for the Schwarz-type smoother in the developed preconditioner yields smaller blocks with THB-splines than it would with a non-truncated basis, and results in a nearly diagonal treatment of untrimmed basis functions. Furthermore, it is demonstrated in [73] that mesh-fitting systems with truncated bases are generally better conditioned than systems with non-truncated bases, and that multigrid methods with a diagonal smoother are effective for systems with THB-splines.

3 Multigrid methods for immersed finite element methods

This section presents the developed geometric multigrid preconditioner for immersed finite element methods. Section 3.1 discusses the conditioning effects of both cut elements and the mesh size of the background grid. In Section 3.2 the employed preconditioning algorithm is introduced. Section 3.3 considers the prolongation and restriction operators, specifically in the context of locally refined grids, and effective smoothers for immersed methods are discussed in Section 3.4.

3.1 Conditioning aspects of immersed finite element methods

Introducing the basis $\{\boldsymbol{\phi}_{i}\}_{i=1}^{n}$ for the approximation space $\mathcal{V}_{h}$, the variational formulation in (2) leads to the linear system:

[TABLE]

with symmetric positive definite (SPD) matrix $\mathbf{A}$, with $A_{ij}=a_{h}\left(\boldsymbol{\phi}_{i},\boldsymbol{\phi}_{j}\right)$, solution vector $\mathbf{x}$, such that $\mathbf{u}_{h}=\sum_{i=1}^{n}x_{i}\boldsymbol{\phi}_{i}$, and right hand side vector $\mathbf{b}$, with $b_{i}=b_{h}\left(\boldsymbol{\phi}_{i}\right)$. Immersed finite element methods generally lead to systems of equations that are fundamentally difficult to solve. Without dedicated treatment of the cut-element-specific ill-conditioning of systems derived from immersed finite element methods, in general the convergence of iterative solvers is severely retarded [35, 50, 51]. An important indicator of the feasibility of iterative solution procedures for a linear system as in (4) is the condition number, $\kappa\left(\mathbf{A}\right)$. For the symmetric positive definite (SPD) systems considered in this contribution, the condition number is equal to the quotient of the largest to the smallest eigenvalue. It should be noted that the convergence of iterative solution methods is not merely dependent on the condition number, but rather depends on the entire spectrum of the system matrix, i.e., on the distribution of the complete set of eigenvalues. Systems with well-clustered spectra generally lead to faster convergence than systems with eigenvalues that are spread out. For a detailed discussion about the aspects affecting the performance of iterative solvers, the reader is directed to [52, 86].

To elucidate the general characteristics of spectra emanating from immersed finite element methods, Figure 3 displays the spectrum of an immersed approximation of the Laplace operator on the star-shaped domain from Figure 1, as well as 4 characteristic eigenmodes. The immersed approximation space is composed of quadratic Lagrange basis functions on a mesh with $h=\frac{1}{8}$. Dirichlet boundary conditions are imposed by means of a penalty method with parameter $\frac{2}{h}$. Because symmetric positive definite (SPD) systems only posses real and positive eigenvalues, their spectra can be conveniently represented by plotting the eigenvalues $\lambda_{i}\in\mathbb{R}^{+}$ vs. their index $i\in\mathbb{N}$, as in Figure 3(a). Figure 3(b) displays a very small eigenmode that is only supported on a small cut element, which is exemplary for eigenmodes that are common in immersed finite element methods. As described in detail in [35], basis functions on small cut elements can become almost linearly dependent, which yields very small eigenvalues and is not repaired by Jacobi preconditioning. It can be observed that this eigenmode is very similar to the largest eigenmode of the system, plotted in Figure 3(c), as both eigenmodes consist of essentially the same basis functions. However, while the almost linearly dependent basis functions are subtracted from each other such that these cancel out in Figure 3(b), the opposite happens in Figure 3(c). Note that with the quadratic Lagrange basis, there are 4 almost linearly dependent basis functions that are only supported on the small cut elements, such that the largest eigenvalue is bounded from below by approximately $4$. Based on an analysis of the typical small eigenmodes in immersed finite element methods, an estimate of the condition number of SPD systems for second-order PDEs is derived in [35]:

[TABLE]

with $p$ the polynomial degree of the approximation space, $d$ the number of dimensions and $\eta$ the smallest volume fraction, defined as the smallest relative intersection of an element with the physical domain:

[TABLE]

Since cut elements can be arbitrarily small, systems derived from immersed finite element formulations can be arbitrarily ill-conditioned. As a result, iterative solvers are generally ineffective in case that no dedicated treatment for the cut-element-induced conditioning problem is applied, see e.g., the introduction in Section 1. Figures 3(d) and 3(e) portray characteristic eigenfunctions that are not exclusive to immersed finite element methods. The smooth eigenmode in Figure 3(d) is very similar to the usual smallest eigenmode in mesh-fitting finite elements, and is in close correspondence with the smallest analytical eigenmode of the considered PDE. Figure 3(e) plots the oscillatory eigenfunction with the highest frequency that can be captured by the grid. This eigenmode is similar to the usual largest eigenmode in mesh-fitting systems. The ratio between the largest and smallest eigenvalue in mesh-fitting systems, and with that the condition number, is therefore mesh-size dependent. It can be shown that for second-order PDEs it holds that [53]:

[TABLE]

This deteriorates the conditioning for very fine meshes, which also retards the convergence of iterative solvers. In particular, smooth eigenmodes as in Figure 3(d) – which yield small eigenvalues of $\mathcal{O}(h^{2})$ with Jacobi preconditioning – converge slowly. For fixed point iteration methods a convergence rate of $1-\mathcal{O}(h^{2})$ can be derived [59], and also plain Krylov subspace methods generally require $\mathcal{O}\left(h^{-1}\right)$ iterations [52].

Multigrid methods are often applied to resolve the mesh-dependence of the condition number according to (7) and the corresponding slow convergence, and provide a conditioning and convergence rate that is independent of the mesh size. Due to the smoothness of small eigenmodes in mesh-fitting systems, these modes can be effectively approximated on a coarser grid. The concept of multigrid methods is to combine iterations on a fine grid – by which the oscillatory modes with large eigenvalues as in Figure 3(e) converge quickly – with coarse grid corrections in which the approximation of the solution is amended by the solution on a coarser grid. These coarse grid corrections enhance the convergence of the smooth eigenmodes as in Figure 3(d), such that a convergence behavior is obtained that is independent of the grid size. The reader is directed to [57, 58, 59, 60] for reference works on multigrid techniques. In [50] an additive Schwarz preconditioner is presented that is tailored to immersed finite element methods, and resolves the conditioning problems related to cut elements. Furthermore, it is observed that the systems treated with this preconditioner behave similarly with respect to the grid size as mesh-fitting systems, both in terms of the condition number and in terms of the number of iterations with Krylov subspace methods. This results in a computational cost for solving immersed systems that – while scaling better with the size of the system than direct solvers – is suboptimal, in the sense that it is not yet linear with the number of degrees of freedom. In the following sections we therefore incorporate aspects of this preconditioner in a multigrid framework, to obtain a solution method that is robust to cut elements and independent of the size of the system.

3.2 Multigrid V-cycle algorithm

We present the geometric multigrid method in the context of the correction scheme (CS), but the analyses, algorithms and results extend mutatis mutandis to the full approximation scheme (FAS). We consider an algebraic system of the form (4). Let $\mathbf{x}^{\prime}$ denote an approximation to the solution of (4) and let $\mathbf{r}$ denote the corresponding residual according to:

[TABLE]

To obtain the solution to (4), the approximation $\mathbf{x}^{\prime}$ must be corrected as $\mathbf{x}^{\prime}+\bar{\mathbf{x}}$ with:

[TABLE]

Multigrid methods are based on the notion that if $\mathbf{A}$ derives from an elliptic operator, then the inverse $\mathbf{A}^{-1}$ in (9) can be efficiently approximated by a combination of fixed point iterations and a coarse grid correction. The fixed point iterations efficiently reduce the oscillatory components of the error, and are therefore commonly referred to as smoothing [59]. The coarse grid correction again involves the (approximate) inverse of a matrix, analogous to (9), but corresponding to a coarser grid. This inverse can likewise be approximated by smoothing operations and a coarse grid correction on an even coarser grid, see Figure 4. In multigrid methods, the smoothing operations and coarse grid correction are therefore applied recursively, until a grid is reached that is sufficiently coarse to enable direct inversion at a negligible computational expense. Our analyses are based on the notion that one entire multigrid cycle generates an approximation to the inverse of matrix $\mathbf{A}$. The result of the cycle is the vector $\tilde{\mathbf{x}}$, that approximates the vector $\bar{\mathbf{x}}$ in (9).

In this contribution we consider the multigrid V-cycle as outlined in Algorithm 1. The input parameters of this algorithm are the number of (remaining) recursive multigrid levels $\ell$, see Figure 4, and the residual of the linear system $\mathbf{r}_{\ell}$ at the current level. It is emphasized that the number of multigrid levels follows a different convention than the number of hierarchical refinement levels. The coarsest multigrid level in the V-cycle is $\ell=1$, while the coarsest hierarchical refinement level is $m=0$, see Figure 2. The output of the algorithm is the approximate solution $\tilde{\mathbf{x}}_{\ell}$ to $\mathbf{A}_{\ell}^{-1}\mathbf{r}_{\ell}$. The algorithm first initializes the vector $\tilde{\mathbf{x}}_{\ell}$ in line 2, and then in line 4 conducts a pre-smoothing step111Note that in this contribution we only apply a single pre-smoothing and post-smoothing operation in each cycle, and do not consider the possibility of multiple smoothing operations., i.e., fixed point iteration, with $\gamma>0$ denoting a relaxation parameter for stability and $\mathbf{M}_{\ell}^{-1}$ an approximate factorization or inverse of $\mathbf{A}_{\ell}$, e.g., Jacobi or Gauss-Seidel. Subsequently, smooth components of the error – which were not effectively reduced by the aforementioned smoothing operation – are treated by applying a coarse grid correction with coarser level $\ell-1$ in lines 7 to 9. The recursive nature of multigrid is implemented in Algorithm 1, by applying the V-cycle also to obtain an approximate solution to the coarse grid correction problem. The direct solver in line 14 is only applied at the coarsest level $\ell=1$. Note that the residual in the input of the algorithm should be interpreted as in (8) only on the finest level, and is a restriction of a finer residual in the recursive applications of the algorithm in line 8. To enforce symmetry of the linear operator induced by the algorithm, the post-smoothing operation in line 12 is performed with the adjoint of $\mathbf{M}_{\ell}^{-1}$. This is relevant for Gauss-Seidel-type smoothers, and is realized by a reverse sweep. The approximation of the solution is returned in line 16. The coarse grid correction is treated in more detail in Section 3.3, and the smoothing operations are discussed in detail in Section 3.4.

The multigrid V-cycle in Algorithm 1 can itself be applied as a solver, by iteratively performing the cycle to reduce the residual in every step, i.e., the approximation $\mathbf{x}^{\prime}$ is simply updated by directly adding $\tilde{\mathbf{x}}$. This is generally true for multigrid cycles, e.g., also the W-cycle or FMG-cycle [59]. We consider the V-cycle, as this simple setup is already suitable to demonstrate the effectivity of the multigrid concept in immersed FEM. Additionally, when the pre-smoothing and post-smoothing operations are chosen such that these are adjoint, the V-cycle algorithm yields a symmetric positive definite (SPD) linear operator. This has the advantage that, instead of direct application of the V-cycle as a solver, it can also be employed as a preconditioner in a conjugate gradient (CG) algorithm. In the results presented in Section 4, this multigrid-preconditioned CG-solver is applied. The advantage of Krylov subspace solvers compared to multigrid as a standalone solver is that the convergence of Krylov methods is not purely governed by the smallest eigenmode in the system. Therefore, these are more robust to artifacts in the spectrum resulting from e.g., geometrical complexities such as the artificial coupling that is observed in Section 4.1 [87].

3.3 Restriction, prolongation and coarse grid correction

The restriction and prolongation operations to communicate between different grid sizes, such as those in Figure 4, are essential aspects of the coarse grid correction in lines 7 to 10 of Algorithm 1. Under the usual assumption that the grid size is doubled in the mesh coarsening, the grid lines of the coarser level $\ell-1$ in Figure 4(b) coincide with grid lines at the finer level $\ell$ in Figure 4(a). Therefore, the level $\ell-1$ space is nested in the level $\ell$ space and the basis functions on level $\ell-1$ can be represented identically by linear combinations of the basis functions on level $\ell$. Denoting by $\boldsymbol{\Phi}_{\ell}$ and $\boldsymbol{\Phi}_{\ell-1}$ vectors of basis functions on level $\ell$ and level $\ell-1$, respectively, there exists a matrix of coefficients $\mathbf{R}_{\ell}$ such that:

[TABLE]

The matrix $\mathbf{R}_{\ell}$ defines the restriction operator. This restriction operator is employed in line 7 to restrict the residual at level $\ell$ to the level $\ell-1$ coarse mesh. In line 8 of Algorithm 1, the solution to the level $\ell-1$ problem is approximated by recursive application of the V-cycle algorithm, except at the coarsest level $\ell=1$ where a direct solution is carried out. The (approximate) solution to the level $\ell-1$ problem constitutes the coarse grid correction, which is prolongated and added to $\tilde{\mathbf{x}}_{\ell}$ at level $\ell$ in line 9. Let us note, that by the nesting of the approximation spaces at the different levels, the prolongation from level $\ell-1$ to level $\ell$ corresponds to injection. By virtue of relation (10) between the basis functions, we have the identities:

[TABLE]

In terms of the coefficients, the prolongation thus corresponds to the adjoint of the restriction operator. In line 10 the level $\ell$ residual is updated after the coarse grid correction.

Remark 3.1

The coarse grid correction in Algorithm 1 is performed purely algebraically. Therefore, the coarse problem inherits the weak formulation from the fine problem, without adapting the mesh-dependent parameters in the operators in (3). In essence, the finer level basis functions are replaced by the coarser level basis functions, such that the system matrix and residual at the coarser level can simply be obtained as $\mathbf{A}_{\ell-1}=\mathbf{R}_{\ell}\mathbf{A}_{\ell}\mathbf{R}_{\ell}^{\rm T}$ and $\mathbf{r}_{\ell-1}=\mathbf{R}_{\ell}\mathbf{r}_{\ell}$. An alternative approach is to adapt the operators in weak form (2) to the coarser grid, as in e.g., [88, 89]. The implementation of the algebraic approach in Algorithm 1 is considerably simpler, however, and has been observed to adequately resolve the smooth eigenmodes in all considered cases. It should be noted that the extension of this simple algebraic coarsening approach to other grid-size dependencies is not verified in this contribution. Examples of other grid-size dependencies in the weak formulation are the stabilization parameter in Nitsche’s method, see [81, 82] and for spectral effects specifically [90], and the ghost penalty, see [37, 4].

We have so far restricted ourselves to uniform meshes. For hierarchically refined meshes, as illustrated in Figure 5(a), the coarsening procedure in the multigrid method requires reconsideration. To enable an adequate approximation of functions that are smooth with respect to the local mesh width on level $\ell$, we apply a type of hierarchical derefinement to obtain a nested coarse space on level $\ell-1$. In practice, the mesh on level $\ell-1$ is the coarsest mesh for which the mesh of level $\ell$ can be obtained by a single level of hierarchical refinements, i.e., one uniform subdivision of a certain set of elements. The construction of such non-uniform coarser spaces is summarized in Algorithm 2, and illustrated in Figure 5. The algorithm denotes components of the mesh at level $\ell$ by: $K_{i}^{m}\in\mathcal{T}_{h}^{m}\subset\mathcal{T}_{h}$. As introduced in Figure 2, $K_{i}^{m}$ denotes an active element at local refinement level $m$ with index $i$. The set of all active elements of local refinement level $m$ is denoted by $\mathcal{T}_{h}^{m}=\{K_{i}^{m}\}$, and $\mathcal{T}_{h}=\cup_{m=0}^{m=M}\mathcal{T}_{h}^{m}$ denotes the full mesh, with $M$ denoting the number of local refinement levels. Elements of the coarser mesh at level $\ell-1$ are denoted as $k_{i}^{m}\in\mathcal{T}_{2h}^{m}\subset\mathcal{T}_{2h}$. The subscript $2h$ for the coarse mesh is to be conceived of as a symbolic notation. Note that the elements in $\mathcal{T}_{2h}^{m+1}$ are of the same size as the elements in $\mathcal{T}_{h}^{m}$. In line 1 of the algorithm, the coarser mesh is initialized as a uniform mesh with only unrefined elements of local refinement level $m=0$, i.e., $\mathcal{T}_{2h}=\mathcal{T}_{2h}^{0}$. Line 3 initiates a loop over the refinement levels $0\leq m\leq M-1$. Within this loop, the algorithm loops over all the elements $k_{i}^{m}\in\mathcal{T}_{2h}^{m}$, that are currently $m$ times refined. If $k_{i}^{m}$ intersects an element of $\mathcal{T}_{h}$ that is refined more than $m$ times, element $k_{i}^{m}$ is refined in lines 7 to 9. Note that line 7 abuses notation to simplify the expression, and that the refining of element $k_{i}^{m}$ implies removing $k_{i}^{m}$ from $\mathcal{T}_{2h}^{m}$ and adding the refinements to $\mathcal{T}_{2h}^{m+1}$. Note that this hierarchical derefinement procedure differs from the approach in [73], which employs an existing hierarchy of refined grids for the coarse grid corrections.

3.4 Smoothers for immersed finite element methods

In lines 4 and 12 of Algorithm 1, smoothing operations, i.e., fixed point iterations, are performed to resolve the components of the error that cannot be adequately captured by the coarse grid. Effective application of the multigrid algorithm requires that the eigenvalues of $\gamma\mathbf{M}^{-1}\mathbf{A}$ that correspond to non-smooth eigenfunctions are close to $1$. Note that in the formulations from here on, the subscripts indicating the level in the multigrid solver are omitted to simplify the notation. All these formulations are independent of the level $\ell$, however. Stability of fixed point iterations requires:

[TABLE]

with $\lambda_{\rm min}(\cdot)$ and $\lambda_{\rm max}(\cdot)$ denoting the smallest and largest eigenvalues, such that the spectral radius of the fixed point iteration is bounded:

[TABLE]

with $\mathbf{I}$ denoting the identity matrix with the same size as system matrix $\mathbf{A}$. In the case that $\lambda_{\rm max}\left(\gamma\mathbf{M}^{-1}\mathbf{A}\right)>2$, the error component in the direction of the eigenvector corresponding to the largest eigenvalue will increase with smoothing, which may result in divergence. For this reason, smoothers such as Jacobi and additive Schwarz require a sufficiently small relaxation parameter $\gamma$. Smoothers such as Gauss-Seidel and multiplicative Schwarz are unconditionally stable and do not require relaxation, see e.g., [52, Theorem 4.10 & 14.9]. As already mentioned in the description of Algorithm 1 in Section 3.2, it should be noted that symmetry of the linear operator that is induced by the V-cycle insists that the post-smoothing operation is the adjoint of the pre-smoothing operation. Furthermore, it should be mentioned that the computational efficiency can potentially be improved by performing multiple smoothing operations in each cycle. Such enhancements are, however, not considered in this contribution.

Jacobi iterations are not suitable as a smoother for immersed finite elements, which is illustrated by the example in Figure 3. The smallest eigenmode in Figure 3(b) with a very small eigenvalue is barely affected by the smoothing, and cannot be captured on a coarser grid. Furthermore, the relatively large eigenmodes caused by almost linear dependencies as in Figure 3(c) impose a small relaxation parameter, which further impairs the conditioning. The following subsections examine the suitability of a Gauss-Seidel smoother and Schwarz-type smoothers based on the preconditioner for immersed finite elements developed in [50].

Gauss-Seidel

A typical spectrum and characteristic eigenmodes with standard Gauss-Seidel preconditioning are shown in Figure 6(b). Similar to Figure 3 for Jacobi preconditioning, these figures correspond to the Laplace operator on the geometry in Figure 1 with quadratic Lagrange basis functions and boundary conditions imposed by the penalty method. To obtain a symmetric preconditioner, a double fixed point iteration with adjoint Gauss-Seidel operations is applied in these figures:

[TABLE]

with $\left(\mathbf{I}-\mathbf{M}^{-1}\mathbf{A}\right)$ corresponding to the initial Gauss-Seidel sweep, $\left(\mathbf{I}-\mathbf{M}^{-{\rm T}}\mathbf{A}\right)$ corresponding to the reverse sweep, and $\mathbf{y}_{\lambda_{i}}$ denoting the eigenvector corresponding to the $i$th eigenvalue $\lambda_{i}$. Hence, Figure 6(b) presents the eigenmodes of the system $\left(\mathbf{M}^{-1}+\mathbf{M}^{-{\rm T}}-\mathbf{M}^{-{\rm T}}\mathbf{A}\mathbf{M}^{-1}\right)\mathbf{A}$. As will follow in (15), this is actually very similar to the V-cycle, except for the omission of the coarse grid correction. To reduce the computational cost, the Gauss-Seidel routine is implemented with a graph coloring algorithm222This algorithm divides the basis functions in sets, or colors, such that the supports of basis functions with the same color do not intersect. The union over the sets of all colors constitutes the full approximation space. Note that, for Lagrange basis functions on uniform grids, this requires $(p+1)^{d}$ different colors. As basis functions of the same color do not intersect, updating the approximation of the solution, $\tilde{\mathbf{x}}$, at an index corresponding to a certain color, does not affect the residual, $\mathbf{r}$, at the other indices corresponding to that color. This enables a Gauss-Seidel routine that sweeps over all indices of a certain color at once, instead of sequentially updating each index of the approximation of the solution and updating residual accordingly. [91]. In Figure 6(b) a very small eigenvalue is plotted, which is similar to the smallest eigenmode with Jacobi preconditioning in Figure 3(b). Note that the eigenvalues of these eigenfunctions differ by a factor of (approximately) $2$, which is an expected consequence of the double iteration with Gauss-Seidel versus the single iteration with Jacobi. In Figure 6(c) it is shown that also the spectrum with Gauss-Seidel preconditioning contains an eigenmode that is in close correspondence with the smallest analytical eigenmode of the PDE, and is similar to the usual smallest eigenmode with mesh-fitting techniques. Since the colors are selected such that basis functions of the same color do not intersect, the unit vectors corresponding to basis functions with the last color of the graph coloring algorithm yield an eigenspace with an eigenvalue of exactly $1$. This can be observed in the spectrum in Figure 6(a) and is illustrated in Figure 6(d).

Gauss-Seidel is not suitable as a smoother for immersed finite element methods, despite the unconditional stability by which – in contrast to Jacobi – it does not require more relaxation in an immersed setting than in a mesh-fitting setting. The problem that renders Gauss-Seidel ineffective for immersed finite element methods is the small eigenmode plotted in Figure 6(b). Similar to Jacobi, the eigenvalue of this mode is too small to be adequately treated by the smoothing operations, and it is also not resolved by the coarse grid correction. This is shown in Figure 7, which displays the spectrum of the same problem preconditioned by the V-cycle in Algorithm 1 with $\ell=2$ levels and a single Gauss-Seidel sweep as smoother. This eigenvalue problem can be formulated as:

[TABLE]

or as:

[TABLE]

While a comparison of the spectra in Figures 6(a) and 7(a), i.e., without and with the coarse grid correction, reveals that several small eigenmodes are resolved by the coarse grid correction, Figures 6(b) and 7(b) demonstrate that both spectra contain approximately the same small eigenmode. The eigenvalues of this mode are barely affected by the coarse grid correction, with eigenvalue $3.20\cdot 10^{-6}$ for the double fixed point iteration and eigenvalue $3.78\cdot 10^{-6}$ for the full V-cycle.

Additive Schwarz

The spectra and eigenfunctions with Jacobi and Gauss-Seidel in Figures 3 and 6(b) clearly demonstrate that these are not robust to cut elements. This is consistent with the analysis of the conditioning problems in [35], which points out that diagonal preconditioners do not adequately mitigate the almost linear dependencies that occur in immersed finite element methods. In [50] it is derived that almost linearly dependent basis functions can be effectively treated collectively when these are inverted in a block manner by a Schwarz-type method. Based on the additive Schwarz lemma, [50] shows that additive Schwarz preconditioning is actually a very natural approach to resolve the small eigenmodes caused by almost linear dependencies. Furthermore, it is demonstrated that in terms of the condition number and the number of iterations in an iterative solution method, immersed methods with additive Schwarz preconditioning behave similar to mesh-fitting techniques. This section therefore examines the spectrum of immersed methods with additive Schwarz preconditioning, to establish the suitability as a smoother in a multigrid method.

In additive Schwarz preconditioning, a set of $N$ index blocks is selected, which correspond to sets of basis functions. For each index block $j\leq N$, the system matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ is restricted to the indices in the block, denoted by $\mathbf{A}_{j}\in\mathbb{R}^{n_{j}\times n_{j}}$. The block matrices are then inverted and prolongated to a matrix of size $n\times n$. The additive Schwarz preconditioner is obtained by summing these matrices:

[TABLE]

with $\mathbf{P}_{j}\in\mathbb{R}^{n\times n_{j}}$ a matrix that prolongates a block vector $\mathbf{y}_{j}\in\mathbb{R}^{n_{j}}$ to a vector $\mathbf{P}_{j}\mathbf{y}_{j}=\mathbf{y}\in\mathbb{R}^{n}$ – with nonzero entries only at the indices in block $j$ – and the transpose of $\mathbf{P}_{j}$ a restriction operator that restricts a vector $\mathbf{z}\in\mathbb{R}^{n}$ to a block vector $\mathbf{P}_{j}^{\rm T}\mathbf{z}=\mathbf{z}_{j}\in\mathbb{R}^{n_{j}}$ – containing only the indices in block $j$.

An essential aspect of additive Schwarz preconditioners is the choice of the index blocks. It is pointed out in [50] that almost linearly dependent basis functions are required to be in an index block together. Furthermore, it is demonstrated that devising a block for each cut element with all basis functions supported on it is an effective strategy to satisfy this requirement for uniform grids. As demonstrated in [51], however, it is not trivial to generalize this concept to locally refined meshes. Therefore this contribution applies an alternative strategy to select the Schwarz blocks based on so-called encapsulating supports, which is inspired by the Schwarz-type smoother developed for divergence-conforming discretizations in [71]. Accordingly, for every basis function a block is devised, containing all the basis functions whose support completely lies inside the support of the basis function associated to the block, see Figure 8. Note that in this contribution the support of a basis function refers to the support within the physical domain. The block that is associated to function $\phi_{j}$ is defined as:

[TABLE]

with ${\rm supp}_{\Omega}(\phi_{k})$ denoting the support of basis function $\phi_{k}$ within physical domain $\Omega$. For vector-valued problems, separate blocks are devised for basis functions describing different vectorial components of the solution, similar to the blocks in [50]. By construction, each block therefore contains the basis function associated to it, and for untrimmed (truncated hierarchical) B-splines this approach yields an approximately diagonal preconditioner. Let us note here the importance of the truncation of the basis functions in the locally refined approximations. As non-truncated hierarchical bases yield a very large number of basis functions with overlapping supports, this correspondingly results in very large blocks in the additive Schwarz preconditioner, which leads to significant computational costs. Trimmed basis functions on small cut elements – which can be almost linearly dependent – are also assigned to blocks associated to other functions. This satisfies the requirement formulated in [50] that almost linearly dependent basis functions need to be in a block together, and therefore resolves the small eigenmodes that are characteristic for immersed finite element methods. This strategy to select the Schwarz blocks is directly applicable to both B-splines on uniform grids and truncated hierarchical B-splines on non-uniform grids, in contrast to the element-wise strategy in [50]. It is, however, not natural to directly apply this strategy to Lagrange basis functions, as these are not all supported on the same number of elements. For the Lagrange bases, blocks are therefore only devised for the nodal basis functions333In this contribution different Lagrange basis functions are indicated as nodal, edge, face and volume functions. With $d$ denoting the number of dimensions: nodal functions attain the value $1$ at a vertex of the grid and span $2^{d}$ elements, edge functions attain the value $1$ on an edge and span $2^{d-1}$ elements, face functions attain the value $1$ on a face and span $2^{d-2}$ elements, and volume or element internal functions attain the value $1$ inside an element and span $1$ element. Note that face functions are not considered in the two-dimensional example.. This implies that with Lagrange bases a block is devised for every cluster of $2^{d}$ elements. Note that this yields approximately the same number of blocks for uniform Lagrange bases as for uniform B-spline bases, but, in contrast, does not reduce to a purely diagonal treatment of untrimmed basis functions for Lagrange bases. It should be mentioned that the block selection described here is not the only possible and effective method to select blocks for immersed finite elements, and interested readers are directed to [69, 70, 72] for a study of suitable block selections in isogeometric analysis and to the reference works [56, 92] for considerations regarding the block selections with traditional finite element bases.

Efficiency and stability of additive Schwarz as a smoother requires adequate selection of the relaxation parameter $\gamma$. For the smoothing operations to efficiently reduce the error components in the directions of modes with eigenfunctions that can not be adequately captured on coarser grids, it is required that the relaxed eigenvalues $\gamma\lambda_{i}\left(\mathbf{M}^{-1}\mathbf{A}\right)$ corresponding to such eigenmodes are close to $1$. The requirement with regard to stability is formulated in (12), and states that $\lambda_{\rm min}\left(\mathbf{M}^{-1}\mathbf{A}\right)\geq 0$ and that $\gamma\lambda_{\rm max}\left(\mathbf{M}^{-1}\mathbf{A}\right)\leq 2$. The positivity of the eigenmodes follows from the symmetric positive definiteness of both $\mathbf{M}^{-1}$ and $\mathbf{A}$. Since the eigenvalues of $\mathbf{M}^{-1}\mathbf{A}$ coincide with the eigenvalues of $\mathbf{M}^{-\frac{1}{2}}\mathbf{A}\mathbf{M}^{-\frac{1}{2}}$, the eigenmodes can be bounded from above by:

[TABLE]

with $\mathbf{A}^{K_{i}}\in\mathbb{R}^{n\times n}$ the part of system matrix $\mathbf{A}$ that results from the integration over element $K_{i}$, and $N_{K_{i}}$ denoting the number of blocks containing basis functions that are supported on this element. Note that in (19b) the additive Schwarz lemma is applied, see e.g., [56, 92]. In (19d) the maximal quotient with system matrix $\mathbf{A}$ is replaced by the maximum over the quotients with the element contributions $\mathbf{A}^{K_{i}}$, and in (19e) the Cauchy-Schwarz inequality is applied. For Lagrange basis functions this bound is observed to be considerably sharp, as volume basis functions are contained in $N_{K_{i}}$ blocks and form eigenfunctions similar to the one in Figure 9(d). These functions are not captured in the coarse grid correction, such that efficient smoothing requires $\gamma=N_{K_{i}}^{-1}=4^{-1}$ for a two-dimensional Lagrange basis. This relaxation parameter yields $\gamma\lambda_{\rm max}\left(\mathbf{M}^{-1}\mathbf{A}\right)\leq 1<2$, such that also the stability condition is satisfied.

Figure 9(b) presents the spectrum and characteristic eigenmodes of the Laplace operator on the geometry in Figure 1 with a quadratic Lagrange basis that is preconditioned by a double fixed point iteration with additive Schwarz and the relaxation parameter $\gamma=\frac{1}{4}$. The double fixed point iteration is applied such that later on these results can easily be related to the results with the V-cycle with additive Schwarz smoothing in Figure 10. Figure 9(b) shows that the smallest eigenvalue in the system with additive Schwarz does not correspond to an eigenfunction on a small cut element. Instead, the smallest eigenmode is similar to the usual smallest eigenmode with mesh-fitting methods, which can be considered as the smoothest possible mode satisfying the boundary conditions. Figure 9(c) plots an example of an eigenfunction that almost entirely consists of nodal basis functions. The spectrum contains multiple of such modes, and these are important because these are the smallest eigenmodes that cannot be adequately captured by the coarse grid correction. Therefore, these nodal modes form the bottleneck for the condition number of immersed systems preconditioned by the V-cycle with additive Schwarz smoothing, as will follow in Figure 10. In fact, the eigenvalue of such modes can be clarified. Since nodal basis functions are in only $1$ index block, a fixed point iteration reduces the contribution of such functions by approximately a factor $1-\gamma=\frac{3}{4}$. The double fixed point iteration therefore results in an eigenvalue of approximately $1-(1-\gamma)^{2}=\frac{7}{16}\approx 0.438$. The largest eigenmodes in the spectrum have an eigenvalue of approximately $1$, and correspond to eigenfunctions consisting almost entirely out of basis functions that are only supported on 1 element and therefore are contained in 4 index blocks. On the interior this only involves volume basis functions, resulting in the volume mode in Figure 9(d). By virtue of the eigenvalue of approximately $1$, error components in the direction of the eigenvectors of these modes are effectively eliminated by the additive Schwarz smoother.

Based on the smallest eigenmodes in the spectrum with additive Schwarz, this technique is suitable as a smoother in a multigrid method for immersed finite elements. The results of this smoother in a two-level V-cycle are presented in Figure 10. It can be observed from the spectra in Figures 9(a) and 10(a) – i.e., respectively without the coarse grid correction and with the coarse grid correction – that the smooth eigenmodes are effectively resolved, such that a method is obtained that is robust to both cut elements and the grid size. It is noteworthy that the smallest eigenvalues with the multigrid method correspond to nodal modes, see Figure 10(b). The limited effectiveness of the multigrid procedure for these modes derives from the fact that these modes are relatively insensitive to the smoothing operations, see Figure 9(c), and that these modes cannot be adequately captured on a coarser grid. With multigrid as a standalone solver, these modes would yield a convergence rate between $0.5$ and $0.6$. This rate can be improved by applying the multigrid cycle as a preconditioner in a Krylov subspace solver.

Figure 11 presents the spectra of the same problem with finer grids, preconditioned by a double fixed point iteration with additive Schwarz and preconditioned by the full V-cycle with additive Schwarz smoothing. Comparing these, and also the spectra with a grid of $16\times 16$ elements in Figures 9(a) and 10(a), conveys that systems without the coarse grid correction closely follow the grid-size dependence of the conditioning formulated in (7), and that with the full V-cycle a conditioning is obtained that is independent of the grid size, with for this problem eigenvalues $0.4<\lambda_{\rm min}<0.5$ and $\lambda_{\rm max}=1$.

While smoothing with additive Schwarz results in a conditioning that is independent of the grid size, the method is affected by the small relaxation parameter. This will be even more severe in three dimensions and for B-spline bases, which are supported on more elements and therefore require smaller and degree-dependent relaxation parameters, i.e., for B-splines $N^{K_{i}}=(p+1)^{d}$ with $p$ the spline degree and $d$ the number of dimensions. Therefore, the next section considers multiplicative Schwarz as a smoother, which is unconditionally stable [52] and thereby circumvents the small relaxation parameter required for stability.

Multiplicative Schwarz

Multiplicative Schwarz can be considered as the block equivalent of Gauss-Seidel. While with additive Schwarz all the locally inverted block matrices are applied to the same residual – similar to the diagonal elements in Jacobi – with multiplicative Schwarz the residual is updated after each block – similar to the update of the residual after each diagonal element in Gauss-Seidel. Multiplicative Schwarz can be formulated by initializing the zero vector $\tilde{\mathbf{y}}^{0}=\mathbf{0}$, defining the initial residual $\tilde{\mathbf{r}}^{0}=\mathbf{r}$ and looping over the blocks $j\leq N$:

[TABLE]

The linear operator of multiplicative Schwarz is then defined as $\mathbf{M}^{-1}\mathbf{r}=\tilde{\mathbf{y}}^{N}$. Similar to Gauss-Seidel, multiplicative Schwarz does not require stabilization [52], which is the most important motivation to examine multiplicative Schwarz as a smoother in multigrid methods for immersed finite elements.

The linear operator induced by a fixed point iteration with multiplicative Schwarz is not symmetric. Therefore, it is important that in the post-smoothing the direction in the loop over the index blocks is reversed, to restore symmetry of the linear operator induced by the V-cycle. As can be observed in (20), the application of multiplicative Schwarz requires updating the residual at every step, which is computationally expensive and impedes parallelization. Therefore, similar to standard Gauss-Seidel, the index blocks are ordered by a graph coloring algorithm [91]. To this end, the blocks are divided in $C$ colors indicated by $c\leq C$. The basis functions in a block of a certain color do not intersect functions in a different block with the same color. Therefore, blocks of the same color are not affected by each other’s update of the residual. As a result, the residual only needs to be updated after executing all blocks of a certain color. This reduces the computational cost and enables parallelization of the routine. The multiplicative Schwarz procedure with graph coloring can be formulated by initializing the zero vector $\tilde{\mathbf{z}}^{0}=\mathbf{0}$, defining the initial residual again as $\tilde{\mathbf{r}}^{0}=\mathbf{r}$ and looping over the colors $c\leq C$:

[TABLE]

with $\mathcal{J}_{c}$ denoting the set of blocks with color $c$. The linear operator of multiplicative Schwarz with graph coloring is defined as $\mathbf{M}^{-1}\mathbf{r}=\tilde{\mathbf{z}}^{C}$. The required minimal number of colors can be deduced by considering the overlap between the supports of basis functions. Because the supports of identically-colored basis functions are not allowed to intersect, for scalar problems and uniform grids the required number of blocks equals the number of elements on which basis functions are supported. For Lagrange bases, uniform B-spline bases and truncated hierarchical B-spline bases, the employed number of colors therefore amounts to $2^{d}$, $(p+1)^{d}$, and $M(p+1)^{d}$, respectively. For hierarchical bases the number of colors required for a uniform grid is multiplied by the number of hierarchical levels and for vector-valued problems the number of colors is multiplied by the number of components of the vector.

Figure 12 presents typical spectra of immersed systems, preconditioned by a symmetric double fixed point iteration with multiplicative Schwarz and preconditioned by a two-level V-cycle with multiplicative Schwarz smoothing. These results again pertain to the Laplace operator on the domain in Figure 1 with quadratic Lagrange basis functions. The single-grid results without the coarse grid correction show a largest eigenmode of $1$ and a mesh-dependent smallest eigenmode of order $\mathcal{O}(h^{2})$, such that the condition number follows the relation for mesh-fitting systems in (7). Similar to the results with additive Schwarz in Figure 11, the results with the two-level V-cycle show a spectrum that is robust to both cut elements and the mesh size of the background grid. The conditioning with multiplicative Schwarz is much better and even nearly optimal, however, because it is not compromised by the small relaxation parameter that was required for the stability of additive Schwarz. Therefore, we conclude that multiplicative Schwarz is the most suitable smoothing procedure for multigrid methods for immersed systems. In view of its superior smoothing properties over additive Schwarz, in the numerical examples in Section 4 we restrict our considerations to multiplicative Schwarz.

Remark 3.2

The resulting preconditioning techniques for Lagrange bases and B-splines have much in common, but also small differences are to be noted. As mentioned previously, Lagrange bases require a considerably smaller number of colors in the implementation. Also, the solver for the Lagrange bases requires approximately $2.5$ times fewer iterations, as can be observed in the results in Section 4. While both bases require approximately the same number of Schwarz blocks on the same grid, an advantage of the B-splines is that untrimmed basis functions result in a diagonal treatment, while untrimmed Lagrange basis functions still result in blocks with size $(2p-1)^{d}$. Also the considerably smaller number of degrees of freedom with B-splines is an aspect to consider. The methods for both bases, however, can be further optimized in regard of the computational efficiency. Therefore, we would like to emphasize that this contribution is intended to demonstrate the feasibility of multigrid methods for immersed finite elements, and not as a qualitative comparison between the different bases. Furthermore, while the computation time is linear with the number of degrees of freedom for both bases, the computation time is highly dependent on the implementation. For that reason, the computation time to solve the systems in Section 4 is not presented.

Remark 3.3

Similar to the Schwarz-type methods for immersed finite elements presented in [50] and [51], it is possible that block matrices contain eigenvalues of the order of the machine precision. Directly inverting such block matrices is unstable, as due to round-off errors these inverses can be inaccurate and can even contain negative eigenvalues for positive definite systems. As described in detail in [50, Remark 3.3], in case a block matrix contains an eigenvalue that is a factor $10^{16}$ smaller than the largest diagonal entry of the matrix, the basis function that is dominant in the corresponding eigenvector is removed from the block. Because this only pertains to basis functions with extremely small contributions, this does not affect the convergence of the iterative solver or the accuracy of the solution.

Remark 3.4

In principle, the ill-conditioning effects of small cut elements only need to be resolved on the finest grid, as the coarser levels are only required to resolve smooth components of the error. Diagonal smoothers therefore suffice on the coarser levels, and Schwarz-type smoothing is only required on the finest level. Applying Schwarz-type smoothing on all levels, however, retains the natural recursive character of the multigrid algorithm. For conciseness, we have therefore opted not to include results with diagonal smoothing on the coarser levels in the numerical results in Section 4.

4 Numerical examples

In this section we assess the developed geometric multigrid preconditioning technique for immersed finite element methods on a range of numerical examples. In all these examples the pre-smoothing and post-smoothing operations consist of a single multiplicative Schwarz sweep as described in Section 3.4. First, Section 4.1 considers three-dimensional elasticity problems on uniform grids with both B-splines and Lagrange basis functions. The goal of these simulations is to demonstrate the performance of the preconditioner on increasingly complex geometries. An aspect that is not covered in this section is the treatment of local mesh refinements. Therefore the preconditioner is applied to a level set based topology optimization problem with truncated hierarchical B-splines in Section 4.2. All simulations are performed with quadratic discretizations. This is sufficient to demonstrate the difference between Lagrange basis functions and B-splines, and quadratic function spaces already result in severely ill-conditioned systems in case a dedicated treatment is not applied. The assembly of quadratic systems is considerably less expensive than systems of third degree or higher, however, enabling finer grids to be assessed.

4.1 Linear elasticity problems

This section assesses the developed preconditioning technique in a conjugate gradient solver and discusses the obtained results. All problems are discretized with both quadratic Lagrange basis functions and quadratic B-splines on uniform grids. We first consider the deformation of a tooth-shaped domain subject to a distributed boundary traction. Next, an apple-shaped geometry containing two geometrically singular points is presented, which is deformed by a gravitational load. The third example considers a complex geometry of a $\mu$CT-scanned trabecular bone specimen, that is compressed by a prescribed displacement at the boundary. In this third example an effect regarding the number of levels in the preconditioner is observed. This effect is further investigated in the specifically designed test case in the last example, which is posed on a geometry in the shape of a triple helix.

Tooth-shaped geometry subject to a distributed boundary traction

This first example considers a dimensionless problem posed on an immersed geometry with the shape of a tooth, see Figure 13. The embedding domain is the cube $(-2,2)^{3}$, and the tooth-shaped geometry is obtained by trimming with the level set function:

[TABLE]

with:

[TABLE]

In this level set $\rho_{0}(x,y,z)$ creates the dent in the surface of the tooth, $\rho_{1}(x,y,z)$, $\rho_{2}(x,y,z)$, $\rho_{3}(x,y,z)$, and $\rho_{4}(x,y,z)$ create the dents in the sides of the tooth and the roots are obtained by $\rho_{5}(x,y,z)$ and $\rho_{6}(x,y,z)$. The tips of the roots below $z=-1$ are trimmed by a second trimming operation with the level set function:

[TABLE]

This creates four surfaces on which homogeneous Dirichlet conditions are applied. On the rest of the boundary a normal traction is applied with the magnitude:

[TABLE]

which concentrates around the corner of the tooth. The Lamé parameters are set to $\lambda=\mu=10^{3}$, and the Dirichlet conditions are enforced by the penalty method with penalty parameters $\beta_{h}^{\lambda}=\beta_{h}^{\mu}=\frac{2}{h}$. The resulting displacements and stresses are shown in Figure 13.

To demonstrate the robustness of the preconditioning technique to the number of elements, we discretize the embedding domain with grids of $20^{3}$, $40^{3}$, and $80^{3}$ elements. This results in, respectively, $129\cdot 10^{3}$, $874\cdot 10^{3}$, and $6.43\cdot 10^{6}$ degrees of freedom (DOFs) with quadratic Lagrange basis functions and $21.6\cdot 10^{3}$, $129\cdot 10^{3}$, and $878\cdot 10^{3}$ DOFs with quadratic B-splines. The integration depth as defined in Section 2 is set to $2$ for the grid with $20^{3}$ elements, $1$ for the grid with $40^{3}$ elements, and [math] for the grid with $80^{3}$ elements. This implies that the grid with $20^{3}$ elements is first partitioned by $2$ consecutive bisectioning operations before it is triangulated, the grid with $40^{3}$ elements is first partitioned by $1$ bisectioning operation before it is triangulated, and cut elements in the grid with $80^{3}$ elements are directly triangulated. Note that this results in identical integrated geometries for all grid sizes. The multigrid preconditioner is applied with $2$ and $3$ levels as defined in Section 3 for the grid with $20^{3}$ elements, $2$, $3$ and $4$ levels for the grid with $40^{3}$ elements and $2$, $3$, $4$ and $5$ levels for the grid with $80^{3}$ elements. These numbers of levels are chosen such that for all grid sizes, the largest number of levels results in a direct solution in the preconditioner for a system derived from $4^{3}$ elements.

The convergence of the multigrid-preconditioned conjugate gradient solver is plotted in Figure 14. The results demonstrate a convergence behavior that is virtually independent of the number of elements. Also, the number of iterations is nearly independent of the number of levels in the preconditioner. This can be attributed to the simple compact shape of the geometry, which is resolved well even on the coarsest grid of $4^{3}$ elements, such that with all numbers of levels the coarse grid corrections provide an effective approximation of the smooth eigenmodes. Finally, it can be observed that the systems with B-splines require more iterations than the systems with Lagrange basis functions. In this regard it should be noted, however, that for untrimmed basis functions the algorithm to select the Schwarz blocks yields a diagonal smoother for B-splines and a block smoother for Lagrange basis functions, cf. the discussion on the computational cost of Lagrange and B-spline spaces in Remark 3.2 that also considers the difference in the number of colors and the number of DOFs.

Apple-shaped geometry subject to a gravitational load

This second example is designed to establish the suitability of the preconditioner for non-convex geometries with sharp reentrant corners where stress singularities are to be expected, such that the effectivity of multigrid methods is not generally evident [93]. A dimensionless problem is posed on the geometry with the shape of an apple, see Figure 15. The embedding domain is again the cube $(-2,2)^{3}$, and the shape of the apple is derived through a trimming operation with the level set function:

[TABLE]

with $r(x,y)^{2}=x^{2}+y^{2}$. A second level set function models a bite being taken out of the apple:

[TABLE]

Homogeneous Dirichlet conditions are imposed on the surface of the bite, and a volumetric load $\mathbf{f}=(0,0,-1)$ is applied to model gravity acting on the apple. The Lamé parameters are again set to $\lambda=\mu=10^{3}$ and the Dirichlet conditions are again enforced by the penalty method with parameters $\beta_{h}^{\lambda}=\beta_{h}^{\mu}=\frac{2}{h}$. Figure 15 displays the resulting displacements and stresses.

The problem is discretized with quadratic Lagrange basis functions and quadratic B-splines on a background grid with $80^{3}$ elements. This yields $4.87\cdot 10^{6}$ DOFs supported in the physical domain with the Lagrange basis and $667\cdot 10^{3}$ DOFs with the B-splines. The integration depth on the cut elements is set to [math], implying that the the cut elements are directly triangulated and not first partitioned with a bisectioning operation. The convergence of the multigrid-preconditioned conjugate gradient solver with $2$, $3$, $4$, and $5$ levels for both bases is shown in Figure 16. The obtained numbers of iterations do not show an effect of the cusps in the geometry, and are similar to those for the tooth-shaped geometry. Also, the convergence is again virtually independent of the number of levels in the multigrid cycle.

Trabecular bone specimen loaded in compression

This third three-dimensional test case considers the challenging geometry of a $\mu$CT-scanned trabecular bone specimen. This geometry was first presented in [23], and is displayed in Figure 17, together with the embedding domain of dimension $(0,1.28)^{3}$ mm3. A linear elastic material model is employed with Young’s modulus $E=10$ GPa and Poisson’s ratio $\nu=0.3$. The specimen is compressed with an average uniaxial compressive strain of $1\%$, by imposing a homogeneous Dirichlet condition at the top boundary, and at the bottom boundary prescribing a normal displacement of $0.0128$ mm while constraining the tangential displacement. These boundary conditions are weakly enforced by the penalty method with penalty parameters $\beta_{h}^{\lambda}=\beta_{h}^{\mu}=\frac{2}{h}$.

We consider different grids on the embedding domain with $32^{3}$, $64^{3}$, and $128^{3}$ elements. The linear systems derived from these grids contain, respectively, $182\cdot 10^{3}$, $1.03\cdot 10^{6}$, and $6.65\cdot 10^{6}$ DOFs with the quadratic Lagrange bases and $39.9\cdot 10^{3}$, $189\cdot 10^{3}$, and $1.05\cdot 10^{6}$ DOFs with the quadratic B-splines. The integration depth is set to $2$ for the grid with $32^{3}$ elements, $1$ for the grid with $64^{3}$ elements, and [math] for the grid with $128^{3}$ elements. The multigrid cycle by which the systems are preconditioned applies $2$ and $3$ levels for the discretizations with $32^{3}$ elements, $2$, $3$ and $4$ levels for the discretizations with $64^{3}$ elements and $2$, $3$, $4$ and $5$ levels for the discretizations with $128^{3}$ elements. With these numbers of levels, the preconditioner with the largest number of levels applies a direct solver to a system derived from $8^{3}$ elements for all three grid sizes, similar to the first example.

The convergence of the systems with the different numbers of levels in the preconditioner is shown in Figure 18. It can be observed that the convergence of the systems with $32^{3}$ elements and $2$ levels, $64^{3}$ elements with $2$ and $3$ levels, and $128^{3}$ elements with $2$, $3$, and $4$ levels is very similar to that in the previous examples. For this test case, slightly more iterations are required. This is conjecturally connected with the complicated multiscale geometry, by which smaller geometric features are less adequately represented on coarse meshes. The convergence with the largest number of levels in the preconditioner, i.e., the preconditioners in which a direct solver is applied to a system with $8^{3}$ elements, is significantly slower. This is caused by underresolution and corresponding nonphysical behavior in the systems with $8^{3}$ elements, which was also observed in [23]. Both with the Lagrange and the B-spline bases, these coarse systems contain basis functions with a disjoint support in the physical domain, i.e., basis functions that cover the gap between disconnected parts of the geometry as e.g., in the circle in Figure 17(a). Physically, disconnected parts of the geometry should be able to move freely with respect to each other. When a basis function is supported on both parts, however, these disconnected parts are artificially coupled in the numerical model. This results in very different stiffness properties between the systems with a mesh of $8^{3}$ elements and systems with meshes of $16^{3}$ elements and finer. Therefore, the coarse grid corrections with $8^{3}$ elements do not approximate the smooth eigenmodes of the finer meshes adequately, which retards the convergence of the multigrid-preconditioned iterative solver. A similar effect can be observed in nearly incompressible elasticity, where the convergence deteriorates when the coarse grids experience volumetric locking, as in e.g., [94]. To illustrate the artificial coupling, the eigenfunction with the smallest eigenvalue in the three-level preconditioned system with B-splines on a grid of $32^{3}$ elements is shown in Figure 19. This function was obtained by $100$ iterations in a power algorithm, and it is clearly observable that this smallest preconditioned eigenmode contains peaks in the displacement and stress fields at points where basis functions of the coarse grid are supported on disconnected parts of the geometry. Such an artificial-coupling mode does not occur if the system is preconditioned with two levels, i.e., with a direct solver applied to the system with $16^{3}$ elements. While the artificial-coupling effect only moderately increases the number of iterations in this test case, and the iterative solver still converges in an acceptable number of iterations, it can not be ruled out that it can potentially hinder the convergence in other test cases more severely. Therefore, this effect is investigated in more detail in the next test case, which presents a geometry that is specifically designed to study this effect.

Triple helix loaded in compression

This final three-dimensional example is specifically designed to investigate the effect observed in the trabecular bone specimen, by which disconnected parts of the domain are artificially coupled on coarse grids. We consider a dimensionless problem on the geometry in Figure 20, which consists of a triple helix that is connected by a half ring at the top and bottom boundary. The embedding domain again consists of the cube $(-2,2)^{3}$. The helixes are obtained by rotations of the level set function:

[TABLE]

with $r_{\rm inner}=0.25$, $R_{\rm outer}=1.5$, and $r(x,y)^{2}=x^{2}+y^{2}$. For the half ring at the top and bottom boundary the level set function:

[TABLE]

is applied. The Lamé parameters are again set to $\lambda=\mu=10^{3}$. At the top boundary homogeneous Dirichlet conditions are imposed, and at the bottom boundary a normal displacement of $0.04$ prescribed. It should be noted that, as opposed to the bottom boundary condition on the trabecular bone specimen, the tangential displacement at the bottom boundary is not constrained. The boundary conditions are again weakly imposed by the penalty method with parameters $\beta_{h}^{\lambda}=\beta_{h}^{\mu}=\frac{2}{h}$. The solution with B-splines on a grid with $128^{3}$ elements is shown in Figure 20.

The embedding domain is discretized with $128^{3}$ elements, yielding $6.96\cdot 10^{6}$ and $1.11\cdot 10^{6}$ DOFs in the quadratic Lagrange and B-spline bases, respectively. The integration depth is set to [math], and the preconditioner applies $2$, $3$, $4$, and $5$ levels in the V-cycle. The convergence plots are shown in Figure 21. It is clearly observable that the convergence is severely retarded with $5$ levels in the preconditioner. This could be anticipated, because the separation of the helixes is approximately equal to the mesh size on the coarsest level of $8^{3}$ elements. Hence, artificial coupling between the helixes will occur on the coarsest mesh, reducing the effectiveness of the coarse grid correction. To further illustrate this effect, Figure 22 displays the smallest eigenmode in a system with B-splines on $32^{2}$ elements that is preconditioned with $3$ levels in the V-cycle, which behaves similarly as it also contains $8^{3}$ elements at the coarsest level. This mode is the equivalent of the smallest mode with the trabecular bone geometry in Figure 19, as both these smallest eigenmodes clearly demonstrate the aforementioned artificial coupling effect. The preconditioner with $4$ levels – for which approximately $2$ elements of the coarsest level fit between separate helixes – also converges significantly slower. For the quadratic B-splines that span $3$ elements this can obviously be attributed to the artificial coupling effect. For the Lagrange basis functions that only span $2$ elements, it should be noted that basis functions covering the gap between a helix and a half ring can still artificially increase the stiffness of the connection between these. In the systems preconditioned by a V-cycle with $2$ and $3$ levels, the coarsest level contains $64^{3}$ and $32^{3}$ elements, respectively, such that the nonphysical coupling effect is not observed and the convergence is similar to that of the previous examples. We expect that the observed deterioration of the convergence behavior can be mitigated by a dedicated coarsening algorithm, which prevents disjoint supports in coarse basis functions by applying local refinements or XFEM-type enrichments on the coarse grid. In this regard the work presented in [75] is noteworthy, since it considers an algebraic multigrid approach in which a similar form of artificial coupling through a crack on the coarse level is precluded.

4.2 A level set based topology optimization problem with truncated hierarchical B-splines

In this example we apply the developed multigrid preconditioning technique to a dimensionless level set based topology optimization problem, which is inspired by the classical MBB beam [95]. This test case is of particular interest, because it demonstrates the robustness to evolving geometries, and it establishes the suitability to locally refined grids with truncated hierarchical B-splines. Because this involves a large number of computations on an evolving geometry, this numerical experiment is performed in two dimensions to reduce the computational cost. The design space and boundary conditions of the optimization problem are shown in Figure 23, and samples of the grids and geometries during the procedure are presented in Figure 24. The objective of the design problem is to minimize the strain energy of the structure, subject to a volume constraint that restricts the volume to 30% of the design domain. The level set function is discretized by linear basis functions which are smoothed by a linear filter [96]. To mitigate the dependence of the optimization results on the initial level set function, a hole seeding method is used that considers the co-evolution of a density field [97]. The parameters of the discretized level set and density fields are treated as optimization variables and are updated in the optimization process by the globally convergent method of moving asymptotes (GCMMA) [98]. This test case is intended to establish that elasticity problems on these geometries and grids can be robustly solved in an iterative manner by means of the developed multigrid preconditioner. This opens the doors to level set based topology optimization problems beyond the reach of direct solvers.

On the locally refined grids supplied by the topology optimization procedure, quadratic truncated hierarchical B-splines are constructed, which are trimmed with the provided level set functions. In steps 0 to 49 of the routine, a uniform grid of $60\times 20$ elements is applied. In steps 50 to 124, grids with a single level of hierarchical refinements are used, see Figure 24. The refinement region is determined based on the design at iterations 49 (i.e., the final design on the uniform grid), 74, and 99. In steps 125 to 149, a grid with two levels of hierarchical refinements is employed. This two-level refinement is based on the final design with a single level of refinements in step 124. The two-level locally refined mesh is updated in step 150 based on the design in the previous step. After this mesh update, the topology optimization algorithm reaches convergence and terminates. The final design is displayed in Figure 24(f).

For each level set and each mesh occurring in the topology-optimization procedure, we consider a linear elasticity problem with Lamé parameters $\lambda=\mu=10^{3}$. The Dirichlet conditions are weakly imposed by the penalty method, using the penalty parameters $\beta_{h}^{\lambda}=\beta_{h}^{\mu}=4\frac{2}{h}$, with $h=\frac{1}{20}$ denoting the unrefined element size and the factor $4$ to cover for the local refinements. The vertical support is applied over the width of one unrefined element, and the load at the left top is applied over the same width and has a magnitude of $20$. We apply the conjugate gradient solver, preconditioned by the multigrid cycle with the tailored multiplicative Schwarz smoother. Three levels are applied in the preconditioner, such that the direct solver is applied to a coarsest grid of $15\times 5$ elements, with a similar pattern of [math], $1$, or $2$ levels of local refinements as supplied by the optimization procedure for that step.

Figure 25 plots the number of DOFs and the number of CG iterations during the procedure. The number of DOFs clearly shows a sharp increase after $50$ and $125$ iterations, when the number of local refinement levels is increased. The number of iterations in Figure 25(b) demonstrates that the preconditioning technique is robust to cut elements, and is not sensitive to changes in the geometry and topology. It can be observed that the number of iterations moderately increases from steps $25$ and $50$ and between steps $50$ and $125$, as the complexity of the evolving physical domain increases. Furthermore, the number of iterations is reduced when more levels of local refinements are applied. These effects are similar to those observed in the three-dimensional test cases; since the extra levels of local refinements are also applied to the coarser grids, the coarse grid correction terms resolve the smooth eigenmodes of the finest grid more accurately, which clarifies the reduction in the number of iterations.

5 Conclusion

This contribution develops a geometric multigrid preconditioner that enables iterative solutions for higher-order immersed finite element methods at a computational cost that is linear with the number of degrees of freedom. This preconditioning technique is robust to the cut elements in immersed methods, and is applicable to traditional, isogeometric, and locally refined discretizations. This is an improvement with respect to state-of-the-art preconditioning techniques for immersed finite element methods and immersed isogeometric analysis, as these are either restricted to linear discretizations [47, 48, 49] or provide convergence rates that are dependent on the grid size [35, 50, 51].

A spectral analysis of immersed finite element methods reveals a spectrum that consists of a combination of i) generic modes that can also be observed in mesh-fitting methods and ii) modes that are characteristic to immersed finite elements and are only supported on small cut elements. Furthermore, immersed systems that are treated by a Schwarz-type preconditioner – with blocks selected such that the linear dependencies on small cut elements are resolved – possess essentially identical spectral properties as mesh-fitting methods. This confirms the observations with additive Schwarz preconditioners in [50], and opens the door to the application of the established and highly efficient framework of multigrid techniques in immersed formulations.

The presented examples convey that the developed multigrid preconditioner results in a number of iterations that is i) independent of the mesh size, ii) only very slightly affected by the number of levels in the multigrid cycle, and iii) only marginally affected by the geometry of the problem. These observations indicate that, while the examples in this contribution already contain multi-million degrees of freedom, the preconditioned iterative solution method enables large-scale computations with immersed finite element methods. While the presented results were obtained with a sequential implementation, further upscaling of the number of degrees of freedom requires an efficient parallel implementation, for which the procedure is suitable.

The numerical results are obtained with a baseline multigrid algorithm and relatively straightforward selection of the Schwarz blocks. While the scaling of the computational cost is already optimal with respect to the system size, the framework of multigrid preconditioners with Schwarz-type smoothers allows for adjustments that can enhance the efficiency even further. First of all, as observed in [69, 70, 72], the efficiency of Schwarz-type smoothers in isogeometric methods is sensitive to the size and the overlap in the selection of the blocks, and different block selections have not been investigated in this contribution. Second, smoothing with multiplicative Schwarz is computationally more expensive than smoothing with additive Schwarz. While multiplicative Schwarz generally results in superior spectral properties and fewer iterations, a detailed comparison between these in the context of the overall computational cost has not been performed. Finally, the V-cycle with a single pre-smoothing and post-smoothing operation is the simplest symmetric multigrid cycle, and it is anticipated that improvements of the cycle design are possible.

Another recommendation pertains to the robustness with respect to the geometrical complexity. Different grid sizes can yield very different stiffness properties when coarse grid basis functions cover the gap between disconnected parts of the geometry. This reduces the effectivity of the multigrid method. Therefore, the robustness can be enhanced by a dedicated coarsening algorithm that identifies coarse grid basis functions with disjoint supports and resolves this with a local refinement or an XFEM-type enrichment.

This contribution only considers symmetric positive definite problems. Based on the investigation of the conditioning problems for different partial differential equations in [50], this is representative for the specific cut-element-related conditioning problems in immersed finite element methods and immersed isogeometric analysis. While symmetric positive definite problems cover a large variety of problems in computational mechanics, the developed preconditioning technique is not immediately applicable to nonsymmetric and mixed formulations in, in particular, flow problems. However, multigrid methods are commonly applied in mesh-fitting flow problems, see e.g., [60], and have been observed to be effective with very similar Schwarz (or Vanka [79]) blocks in [71]. Furthermore, it is demonstrated in [50] that Schwarz-type methods can effectively resolve the cut-element-specific conditioning problems in immersed flow problems. Therefore, it is anticipated that the developed multigrid preconditioner extends mutatis mutandis to problems that are not symmetric positive definite.

Acknowledgement

The research of F. de Prenter was funded by NWO under the Graduate Program Fluid & Solid Mechanics. The research of J.A. Evans, C. Messe, and K. Maute was funded by the Defense Advanced Research Projects Agency under Grant Number HR0011-17-2-0022. The simulations in this work were performed using the open source software package Nutils ([99], www.nutils.org).