Robust and scalable domain decomposition solvers for unfitted finite element methods

TL;DR

This paper develops robust, scalable domain decomposition preconditioners tailored for unfitted finite element methods, effectively addressing ill-conditioning issues and enabling large-scale 3D simulations.

Contribution

It introduces customized BDDC preconditioners with stiffness weighting and improved coarse degrees of freedom for unfitted finite element discretizations.

Findings

Preconditioners significantly improve solver robustness.

Method scales well for large 3D problems.

Numerical experiments confirm efficiency and robustness.

Abstract

Unfitted finite element methods, e.g., extended finite element techniques or the so-called finite cell method, have a great potential for large scale simulations, since they avoid the generation of body-fitted meshes and the use of graph partitioning techniques, two main bottlenecks for problems with non-trivial geometries. However, the linear systems that arise from these discretizations can be much more ill-conditioned, due to the so-called small cut cell problem. The state-of-the-art approach is to rely on sparse direct methods, which have quadratic complexity and are thus not well suited for large scale simulations. In order to solve this situation, in this work we investigate the use of domain decomposition preconditioners (balancing domain decomposition by constraints) for unfitted methods. We observe that a straightforward application of these preconditioners to the unfitted case has a very poor behavior. As a result, we propose a {customization} of the classical BDDC methods based on the stiffness weighting operator and an improved definition of the coarse degrees of freedom in the definition of the preconditioner. These changes lead to a robust and algorithmically scalable solver able to deal with unfitted grids. A complete set of complex 3D numerical experiments show the good performance of the proposed preconditioners.

Tables

Table 1. Table 1 . Information of the computational meshes used in the examples.

	Id.	$n^{\mathrm{sd}}$	$n^{\mathrm{el}}$	$n^{\mathrm{sdel}}$	$n^{\mathrm{dof}}$	$n^{\mathrm{sddof}}$	$h$	$H$	$H/h$
Sphere	1	8	1064	512	1461	729	0.1250	1.0000	8.0
	2	32	7160	512	8577	729	0.0625	0.5000	8.0
	3	160	51968	512	57123	729	0.0312	0.2500	8.0
	4	1064	396040	512	415775	729	0.0156	0.1250	8.0
	5	7160	3090104	512	3167421	729	0.0078	0.0625	8.0
Popcorn flake	1	8	1920	512	2559	729	0.1125	0.9000	8.0
	2	60	12936	512	15181	729	0.0563	0.4500	8.0
	3	324	95256	512	103701	729	0.0281	0.2250	8.0
	4	1924	729768	512	762497	729	0.0141	0.1125	8.0
	5	12936	5708900	512	5837623	729	0.0070	0.0563	8.0
Block	1	8	2120	512	3280	729	0.2312	1.8500	8.0
	2	56	14336	512	18772	729	0.1156	0.9250	8.0
	3	304	103152	512	120164	729	0.0578	0.4625	8.0
	4	2120	776584	512	842652	729	0.0289	0.2312	8.0
	5	14336	6011832	512	6272228	729	0.0145	0.1156	8.0
Array	1	216	53904	512	72796	729	0.2271	1.8167	8.0
	2	1512	376696	512	451516	729	0.1135	0.9083	8.0
	3	8208	2747880	512	3064716	729	0.0568	0.4542	8.0
Spiral	1	2	435	512	838	729	0.1500	1.2000	8.0
	2	16	2277	512	3501	729	0.0750	0.6000	8.0
	3	96	14309	512	18511	729	0.0375	0.3000	8.0
	4	441	99991	512	115378	729	0.0188	0.1500	8.0
	5	2277	744170	512	802975	729	0.0094	0.0750	8.0
	6	14309	5733813	512	5963098	729	0.0047	0.0375	8.0

Table 2. Table 2 . Description of the BDDC preconditioners used in the numerical examples.

Identifier	Object definition	Selected objects	Weighting operator
Standard BDDC(cef)	As in Sect. 3.2	corners, edges, faces	As in Sect. 4.3
Standard BDDC(ce)	As in Sect. 3.2	corners, edges	As in Sect. 4.3
1st variant BDDC(cef)	As in Sect. 4.4	corners, edges, faces	As in Sect. 4.3
1st variant BDDC(ce)	As in Sect. 4.4	corners, edges	As in Sect. 4.3
2nd variant BDDC(cef)	As in Sect. 4.5	corners, edges, faces	As in Sect. 4.3
2nd variant BDDC(ce)	As in Sect. 4.5	corners, edges	As in Sect. 4.3

Equations

$$-\Delta u=f\text{in }\Omega ,u=g^{\mathrm{D}}\text{on }{\Gamma }_{\mathrm{D}},\nabla u\cdot n=g^{\mathrm{N}}\text{on }{\Gamma }_{\mathrm{N}},$$

$${\mathcal{A}}_e(u,v)\doteq {\int }_{e\cap \Omega }\nabla u\cdot \nabla v\mathrm{d}V+{\int }_{{\Gamma }_{\mathrm{D}}\cap e}\left(\beta uv-v\left(\mathbf{n}\cdot \nabla u\right)-u\left(\mathbf{n}\cdot \nabla v\right)\right)\mathrm{d}S,$$

$${\ell }_e(v)\doteq {\int }_{{\Gamma }_{\mathrm{D}}\cap e}\left(\beta v{g}^{\mathrm{D}}-\left(\mathbf{n}\cdot \nabla v\right)g^{\mathrm{D}}\right)\mathrm{d}S,$$

$${\beta }_e\ge C_e\doteq \underset{v\in \bar{\mathbb{V}}}{\sup }\frac{{\mathcal{B}}_e(v,v)}{{\mathcal{D}}_e(v,v)}$$

$${\mathcal{D}}_e(u,v)\doteq {\int }_{e\cap \Omega }\nabla u\cdot \nabla v\mathrm{d}V,\text{ and }{\mathcal{B}}_e(u,v)\doteq {\int }_{{\Gamma }_{\mathrm{D}}\cap e}\left(\mathbf{n}\cdot \nabla u\right)\left(\mathbf{n}\cdot \nabla v\right)\mathrm{d}S.$$

$${\mathcal{B}}_{e}x=\tilde{\lambda }{\mathcal{D}}_{e}x.$$

$$\bar{\mathcal{A}}(u,v)\doteq \sum _{e\in {\theta }_{\mathrm{act}}}{\mathcal{A}}_e(u,v),\ell (v)\doteq \sum _{e\in {\theta }_{\mathrm{act}}}{\ell }_e(v),foru,v\in \bar{\mathbb{V}}.$$

$$k_2(\bar{\mathcal{A}})\sim \underset{e\in {\theta }_{\mathrm{cut}}}{\min }{\eta }_e^{-(2p_e+1-2/d)},{\eta }_e\doteq \frac{|e\cap \Omega |}{|e|},$$

$${\alpha }_{\lambda }^{\omega }(v_{\omega })\doteq {\int }_{\lambda }{\mathcal{R}}_{\lambda }^{\omega }(v_{\omega })is\; identical\; for\; all\omega \in \mathrm{neigh}(\lambda ),$$

$$\mathcal{W}u(\xi )\doteq \frac{{\sum }_{\omega \in \mathrm{neigh}(\xi )}{u}_{\omega }(\xi )}{|\mathrm{neigh}(\xi )|},$$

$$\mathcal{M}\doteq {\mathcal{I}}_0{\mathcal{A}}_0^{-1}{\mathcal{I}}_0^T+\mathcal{Q}{\tilde{\mathcal{A}}}^{-1}{\mathcal{Q}}^T.$$

$$k_2(\mathcal{M}\bar{\mathcal{A}})\le \tilde{C}\left(1+{\log }^2\left(\frac{H}{h}\right)\right),$$

$${\lambda }_{\max }$$

$${\mathcal{R}}_{\lambda }^{\omega }(v)=0,on\partial \omega \setminus \lambda ,\forall \lambda \in \Lambda ,\forall \omega \in \mathrm{neigh}(\lambda ).$$

$$w\doteq \mathcal{Q}v-v=\sum _{\lambda \in \Lambda }\sum _{\omega \in \mathrm{neigh}(\lambda )}\mathcal{E}{\mathcal{R}}_{\lambda }^{\omega }(\mathcal{Q}v-v)=\sum _{\lambda \in \Lambda }\sum _{\omega \in \mathrm{neigh}(\lambda )}\sum _{{\omega }^{\prime }\mathrm{neigh}(\lambda )}\frac{1}{|\mathrm{neigh}(\lambda )|}{\mathcal{J}}_{\omega ,{\omega }^{\prime }}^{\lambda }(v).$$

$$\Vert {\mathcal{J}}_{\omega ,{\omega }^{\prime }}^{\lambda }(v){\Vert }_{\mathcal{A}}\lesssim \left(1+{\log }^{\rho }\left(\frac{H}{h}\right)\right)|v{|}_{H^1(\omega )}^2,$$

$${\mathcal{R}}_{\lambda }^{\omega }(v)=0,on\partial (\omega \cap \Omega )\setminus \lambda ,\forall \lambda \in \Lambda ,\forall \omega \in \mathrm{neigh}(\lambda ),$$

$$\Vert {\mathcal{J}}_{\omega ,{\omega }^{\prime }}^{\lambda }(v){\Vert }_{\mathcal{A}}$$

$$\le |{\mathcal{R}}_{\lambda }^{\omega }(v)-{\alpha }_{\lambda }^{\omega }(v){|}_{H^{\frac{1}{2}}(\partial \omega )}+|{\mathcal{R}}_{\lambda }^{{\omega }^{\prime }}(v)-{\alpha }_{\lambda }^{{\omega }^{\prime }}(v){|}_{H^{\frac{1}{2}}(\partial {\omega }^{\prime })},$$

$${\lambda }_{\max }=\underset{\tilde{v}\in \tilde{\mathbb{V}}}{\max }\frac{⟨\mathcal{A}\mathcal{Q}\tilde{v},\mathcal{Q}\tilde{v}⟩}{⟨\mathcal{A}\tilde{v},\tilde{v}⟩}$$

$$\mathcal{W}u(\xi )\doteq \frac{{\sum }_{\omega \in \mathrm{neigh}(\xi )}{A}_{\xi \xi }^{\omega }{u}_{\omega }(\xi )}{{\sum }_{\omega \in \mathrm{neigh}(\xi )}{A}_{\xi \xi }^{\omega }},$$

$$u(x,y,z)=\sin \left(5\pi {\left(x^2+y^2+z^2\right)}^{1/2}\right),(x,y,z)\in \Omega \subset {\mathbb{R}}^3.$$

Full text

This work has partially been funded by the European Research Council under the Proof of Concept Grant No. 737439 NuWaSim - On a Nuclear Waste Deep Repository Simulator and the EC - International Cooperation in Aeronautics with China (Horizon 2020) under the project: Efficient Manufacturing for Aerospace Components USing Additive Manufacturing, Net Shape HIP and Investment Casting (EMUSIC). SB gratefully acknowledges the support received from the Catalan Government through the ICREA Acadèmia Research Program.

${\dagger}$ Universitat Politècnica de Catalunya, Jordi Girona 1-3, Edifici C1, E-08034 Barcelona $\&$ Centre Internacional de Mètodes Numèrics en Enginyeria, Parc Mediterrani de la Tecnologia, Esteve Terrades 5, E-08860 Castelldefels, Spain E-mail: sbadia@cimne.upc.edu.

${\ddagger}$ Universitat Politècnica de Catalunya, Jordi Girona 1-3, Edifici C1, E-08034 Barcelona $\&$ Centre Internacional de Mètodes Numèrics en Enginyeria, Parc Mediterrani de la Tecnologia, Esteve Terrades 5, E-08860 Castelldefels, Spain E-mail: fverdugo@cimne.upc.edu.

Keywords: Unfitted finite elements, embedded boundary methods, linear solvers, parallel computing, domain decomposition

1. 1 Introduction

2. 2 Unfitted FE method

3. 2.1 Cell partition and subdomain partition

4. 2.2 Model problem and space discretization

5. 3 Domain decomposition for body-fitted FE meshes

6. 3.1 Sub-assembled problem

7. 3.2 Coarse DOFs

8. 3.3 Transfer operator

9. 3.4 Preconditioner and condition number bounds

10. 4 Domain decomposition solvers for unfitted FE meshes

11. 4.1 Breakdown example

12. 4.2 Enlarged coarse spaces and adaptive BDDC methods

13. 4.3 Stifness weighting operator

14. 4.4 New coarse objects: variant 1

15. 4.5 New coarse objects: variant 2

16. 5 Numerical experiments

17. 5.1 Setup

18. 5.2 “Popcorn flake” example.

19. 5.3 Multi-body example: sphere, popcorn flake, block, array of blocks and spiral.

20. 6 Conclusions

1. Introduction

The use of unfitted finite element methods (FEMs) is an appealing approach for different reasons. They are interesting in coupled problems that involve interfaces (e.g., fluid-structure interaction [7] or free surface flows), or situations in which one wants to avoid the generation of body-fitted meshes. These type of techniques have been named in different ways. Unfitted FEMs for capturing interfaces are usually denoted as eXtended FEM (XFEM) [10], whereas these techniques are usually denoted as embedded (or immersed) boundary methods, when the motivation is to simulate a problem using a (usually simple Cartesian) background mesh. Recently, different realizations of the method have been coined in different ways, depending on the way the numerical integration is performed, how Dirichlet boundary conditions are enforced on no-matching surfaces, or the type of stabilization, if any, being used. To mention some examples, the finite cell method combines an XFEM-type functional space, a numerical integration based on adaptive cell refinement and full sub-cell integration, to perform integration on cut cells, and a Nitsche weak imposition of boundary conditions [29]. CutFEM makes also use of Nitsche’s method, but includes additional “ghost penalty" stabilization terms to improve the the stability of the algorithm [13]. The huge success of isogeometrical analysis methods (spline-based discretization) and the severe limitations of this approach in complex 3D geometries will probably increase the interest of unfitted methods in the near future [19].

The main showstopper up to now for the succesful application of unfitted methods for realistic applications is the linear solver step. The condition number of the resulting linear system does not only depend on the characteristic size of the background mesh elements, but also on the characteristic size of the cut elements; cut elements can be arbitrarily small and can have arbitrarily high aspect ratios for small cut cell situations. Enforcing some reasonable thresholds, one can use robust sparse direct solvers for these problems. However, sparse direct methods are very expensive, due to their quadratic complexity, which is especially dramatic at large scales. Scalability is also hard to get, especially at very large scales. For large scale applications with body-fitted meshes in CSE, iterative solvers, usually Krylov solvers combined with preconditioners, is the natural way to go. Scalable preconditioners at large scales are usually based on algebraic multigrid (AMG) and domain decomposition solvers [31, 17, 4, 22]. Unfortunately, the lack of robust and scalable iterative solvers for unfitted FEM has limited the applicability of unfitted methods in real applications. It is well-known that some rudimentary methods, e.g., Jacobi preconditioners, can solve the issue of the ill-conditioning that comes from cut elements. However, these preconditioners are neither scalable nor optimal. Different serial linear solvers for unfitted FEMs have been recently proposed (see, e.g., [27, 11, 18, 15]). The methods in [27, 11] consider a segregation of nodes into healthy and ill nodes, and a domain decomposition solver for the ill nodes. The method is only applied to 2D problems in [27] and serial computations, and the domain decomposition solver is not scalable, due to the fact that no coarse correction is proposed. The method in [11] proposes to use an AMG method for the healthy nodes. A specific-purpose serial AMG solver is designed in [18] and a serial incomplete factorization solver is considered in [15].

Motivated by the lack of robust (with respect to the element cuts) and scalable parallel solvers for unfitted methods, we develop in this work a domain decomposition solver with these desired properties based on the balancing domain decomposition by constraints (BDDC) framework. First, we provide a short description of the mathematical analysis that proves the scalability of BDDC methods on body-fitted meshes. Next, we show an example that proves that the use of standard BDDC methods cannot be robust with respect to the element cuts. In fact, the method performs poorly when straighforwardly applied to unfitted meshes. Next, we propose (based on heuristic arguments and numerical experimentation, but motivated from the mathematical body of domain decomposition methods) a modified BDDC method. First, we consider stiffness weighting operator (already proposed in the original BDDC article [16]), relying on the diagonal of the (sub-assembled stiffness matrix). The use of this weighting operators proves to be essential for the robustness of the solver. Next, we consider enhanced coarse spaces, with a sub-partition of the edge constraints. Since no mathematical analysis is available for domain decomposition methods with stiffness weighting, we have performed a comprehensive set of numerical experiments on 3D geometries with different levels of complexity. The methods are not only robust and algorithmically weakly scalable, but what is more surprising, the number of iterations seems to be independent of the domain shape too. Furthermore, the cost of the modified weighted operator is identical to the standard one, and the algorithm has the same building blocks as the standard BDDC preconditioner. This is good, since one can use the extremely scalable implementation of these methods, e.g., in the FEMPAR library [5, 1] or in the PETSc library [36]. In any case, the cost of the new method can differ from the standard one in the number of coarse degrees of freedom (DOFs). We have observed that in practice, the CPU cost of the modified formulation is very close to the one of the standard preconditioner. In fact, as we increase the size of the coarse problem of the modified solver (which is robust for unfitted methods) tends to the one of the standard preconditioner (which performs very bad for unfitted meshes). We note that the preconditioner proposed herein could also be readily applied to XFEM-enriched interface problems and other discretization techniques (e.g., discontinuous Galerkin methods).

Let us describe the outline of this work. In Sect. 2 we show the unfitted FE method considered in this paper, in particular our choice to integrate in cut cells and the surface intersection algorithm. (In any case, the preconditioners proposed later on do not depend on these choices, and can be used in other situations). We also state the model problem, its discretization, and the mesh and subdomain partitions. The standard BDDC preconditioner and its building blocks for body-fitted meshes are presented in Sect. 3. In Sect. 4, we show a breakdown example of the preconditioner on unfitted meshes, an expensive solution that is provably robust, and a cheap solution based on the stiffness weighting and (optionally) two slightly larger coarse spaces. A complete set of complex 3D numerical experiments are included in Sect. 5 to show experimentally the good performance of the proposed preconditioner. Finally, some conclusions are drawn and future work is described in Sect. 6.

2. Unfitted FE method

2.1. Cell partition and subdomain partition

Let $\Omega\subset\mathbb{R}^{d}$ be an open bounded polygonal domain, with $d\in\{2,3\}$ the number of spatial dimensions. For the sake of simplicity and without loss of generality, we consider that the domain boundary is defined as the zero level-set of a given scalar function $\phi^{\mathrm{ls}}$, namely $\partial\Omega\doteq\{x\in\mathbb{R}^{d}:\phi^{\mathrm{ls}}(x)=0\}$. We note that the problem geometry could be described using 3D CAD data instead of level-set functions, by providing techniques to compute the intersection between cell edges and surfaces (see, e.g., [26]). In any case, the way the geometry is handled does not affect the domain decomposition solver presented below. Like in any other immersed boundary method, we build the computational mesh by introducing an artificial domain $\Omega_{\rm art}$ such that it has a simple geometry easy to mesh using Cartesian grids and it includes the physical domain $\Omega\subset\Omega_{\rm art}$ (see Fig. 1(a)).

Let us construct a partition of $\Omega_{\rm art}$ into cells, represented by $\theta_{\rm art}$. We are interested in $\theta_{\rm art}$ being a Cartesian mesh into hexahedra for $d=3$ or quadrilaterals for $d=2$, even though unstructured background meshes can also be considered. Cells in $\theta_{\rm art}$ can be classified as follows: a cell $e\in\theta_{\rm art}$ such that $e\subset\Omega$ is an internal cell; if $e\cap\Omega=\emptyset$, $e$ is an external cell; otherwise, $e$ is a cut cell (see Fig. 1(b)). The set of interior (resp. external and cut) cells is represented with $\theta_{\rm in}$ and its union $\Omega_{\rm in}\subset\Omega$ (resp. $(\theta_{\rm ext},\Omega_{\rm ext})$ and $(\theta_{\rm cut},\Omega_{\rm cut}))$. Furthermore, we define the set of active cells as $\theta_{\rm act}\doteq\theta_{\rm in}\cup\theta_{\rm cut}$ and its union $\Omega_{\rm act}$. Let us also consider a subdomain partition $\Theta$ of $\Omega_{\rm act}$, obtained by aggregation of active cells in $\theta_{\rm act}$, i.e., there is an element partition $\theta_{\rm act}^{\omega}\doteq\{e\in\theta_{\rm act}:e\subset\omega\}$ for any $\omega\in\Theta$. The interface of the subdomain partition is $\Gamma\doteq\cup_{\omega\in\Theta}\partial\omega\setminus\partial\Omega$.

To generate the subdomain partition in the context of immersed boundary methods is simple, since Cartesian meshes are allowed. For uniform Cartesian meshes, the elements can be easily aggregated into uniform subdomains as shown in Fig. 1(d), whereas for adaptive Cartesian meshes the elements can be efficiently aggregated using space-filling curves [2]. In both cases, the partitions can be generated efficiently in parallel without using graph-based partitioning techniques [20], which are otherwise required for complex unstructured body-fitted meshes. In fact, the elimination of both body-fitted mesh generation and graph partitioning is the main motivation of this work.

2.2. Model problem and space discretization

Let us consider as a model problem the Poisson equation with Dirichlet and Neumann boundary conditions: find $u\in H^{1}(\Omega)$ such that

[TABLE]

where $(\Gamma_{\mathrm{D}},\Gamma_{\mathrm{N}})$ is a partition of the domain boundary (the Dirichlet and Neumann boundaries, respectively), $f\in H^{-1}(\Omega)$, $g^{\mathrm{D}}\in H^{1/2}(\Gamma_{\mathrm{D}})$, and $g^{\mathrm{N}}\in H^{-1/2}(\Gamma_{\mathrm{N}})$.

For the space discretization, we consider $H^{1}$-conforming FE spaces on conforming meshes. For unfitted grids, it is not possible to include Dirichlet conditions in the approximation space in a strong manner. Thus, we consider Nitsche’s method [8, 28] to impose Dirichlet boundary conditions weakly on $\Gamma_{D}$. It provides a consistent numerical scheme with optimal converge rates (also for high-order elements) that is commonly used in the immersed boundary community [33]. We define the space $\bar{\mathbb{V}}\subset H^{1}(\Omega)$ as the global FE space related to the mesh $\theta_{\rm act}$. Further, we define the FE-wise operators:

[TABLE]

defined for a generic mesh element $e\in\theta_{\rm act}$. Vector $\boldsymbol{n}$ denotes the outwards normal to $\partial\Omega$. The bilinear form ${\mathcal{A}}_{e}(\cdot,\cdot)$ includes the usual form resulting for the integration by parts of (1) and the additional term associated with the weak imposition of Dirichlet boundary conditions with Nitsche’s method. The right-hand side operator $\ell_{e}(\cdot)$ includes additional terms related to Nitsche’s method. We will make abuse of notation, using the same symbol for a bilinear form, e.g., ${\mathcal{A}}:\mathbb{V}\rightarrow\mathbb{V}^{\prime}$, and its corresponding linear operator, i.e., $\langle{\mathcal{A}}u,v\rangle\doteq{\mathcal{A}}(u,v)$.

The coefficient $\beta>0$ is a mesh-dependent parameter that has to be large enough to ensure the coercivity of ${\mathcal{A}}_{e}(\cdot,\cdot)$. As shown in [15], coercivity is mathematically guaranteed by considering an element-wise constant coefficient $\beta$ such that

[TABLE]

for all the mesh elements $e\in\theta_{\rm cut}$ intersecting the boundary $\Gamma^{\mathrm{D}}$. In previous formula, $\beta_{e}$ is the value of $\beta$ restricted to element $e$, whereas $\mathcal{D}_{e}(\cdot,\cdot)$ and $\mathcal{B}_{e}(\cdot,\cdot)$ are the forms defined as

[TABLE]

Since $\bar{\mathbb{V}}$ is finite dimensional, and $\mathcal{D}_{e}(\cdot,\cdot)$ and $\mathcal{B}_{e}(\cdot,\cdot)$ are symmetric and bilinear forms, the value $C_{e}$ (i.e., the minimum admissible coefficient $\beta_{e}$) can be computed numerically as $C_{e}\doteq\tilde{\lambda}_{\mathrm{max}}$, being $\tilde{\lambda}_{\mathrm{max}}$ the largest eigenvalue of the generalized eigenvalue problem (see [15] for details):

[TABLE]

In this work, we chose $\beta_{e}\doteq 2\tilde{\lambda}_{\mathrm{max}}$ in all the examples in order to be on the safe side, as suggested in [15]. Let us note that both operators have the same kernel, the space of constants. Thus, the generalized eigenvalue problem (6) is well-posed (see [15] for more details).

The global FE operator ${\bar{\mathcal{A}}}:\bar{\mathbb{V}}\rightarrow\bar{\mathbb{V}}^{\prime}$. and right-hand side term $\ell\in\bar{\mathbb{V}}^{\prime}$ are stated as the sum of the element contributions, i.e.,

[TABLE]

Further, we define $b:\bar{\mathbb{V}}^{\prime}\rightarrow\bar{\mathbb{V}}$ as $b(v)\doteq f(v)+g^{\rm N}(v)+\ell(v)$, for $v\in\bar{\mathbb{V}}$.

The global problem can be stated in operator form as: find $u\in\bar{\mathbb{V}}$ such that ${\bar{\mathcal{A}}}u=b$ in $\bar{\mathbb{V}}^{\prime}$. After the definition of the set of DOFs $\bar{\mathcal{N}}$ and the corresponding FE basis (of shape functions) $\{\phi_{a}(\boldsymbol{x})\}_{a\in\bar{\mathcal{N}}}$ that span $\bar{\mathbb{V}}$, the previous problem leads to a linear system to be solved. At large scales, sparse linear systems are usually solved with Krylov iterative methods [32] due to the quadratic complexity of sparse direct methods. Here, we can consider a conjugate gradient (CG) iterative method since ${\bar{\mathcal{A}}}$ is symmetric positive definite. The convergence rate of the conjugate gradient method is very sensitive to the condition number of the operator, namely $k_{2}({\bar{\mathcal{A}}})\doteq|{\bar{\mathcal{A}}}|_{2}/|{\bar{\mathcal{A}}}^{-1}|_{2}$, where $|{\bar{\mathcal{A}}}|_{2}$ denotes the operator 2-norm. As shown in [15] (under over-optimistic assumptions), the condition number scales as

[TABLE]

where $p_{e}$ is the polynomial degree of the interpolation at element $e$. Thus, arbitrarily high condition numbers are expected in practice since the position of the interface cannot be controlled. In turn, the convergence rate of the iterative solver is expected to be very slow and highly sensitive to the position of the geometry unless a robust preconditioner is considered. Our goal is to develop a BDDC preconditioner able to deal with the bad conditioning properties associated with cut elements, be quasi-optimal with respect to $h$ (the characteristic element size of $\theta_{\rm act}$), and be weakly scalable (not dependent on the number of subdomains for a fixed subdomain problem size). Attaining these objectives, we will have at our disposal a solver for embedded boundary systems that will efficient exploit distributed memory resources for large scale problems.

3. Domain decomposition for body-fitted FE meshes

In this section, we briefly introduce the BDDC method as usually considered in the literature for conventional body-fitted meshes [16]. The description of the algorithm is concise and we mainly focus on the parts that are modified later in Sect. 4 to deal with unfitted meshes. We refer the reader to [25, 12] for a detailed exposition of BDDC methods and to [3, 4] for the practical implementation details. We have also included a short proof of the condition numbers for body-fitted meshes. It is important for the subsequent discussion to explain the breakdown of the standard algorithm for unifitted meshes and motivate robust modifications. In this section, $\Omega\equiv\Omega_{\rm act}$ and the Dirichlet boundary conditions are strongly enforced in the definition of the FE spaces.

3.1. Sub-assembled problem

Non-overlapping domain decomposition preconditioners rely on the definition of a sub-assembled FE problem, in which contributions between subdomains have not been assembled. In order to do so, at every subdomain $\omega\in\Theta$, we consider the FE space $\mathbb{V}_{\omega}$ associated to the element partition $\theta_{\rm act}^{\omega}$ with homogeneous Dirichlet boundary conditions on $\partial\omega\cap\partial\Omega_{\rm D}$. One can define the subdomain operator ${{\mathcal{A}}_{\omega}}(u,v)=\sum_{e\in\theta_{\rm act}^{\omega}}{\mathcal{A}}_{e}(u,v)$, for $u,v\in\mathbb{V}_{\omega}$.

Subdomain spaces lead to the sub-assembled space of functions $\mathbb{V}\doteq\Pi_{\omega\in\Theta}\mathbb{V}_{\omega}$. For any $u\in\mathbb{V}$, we define its restriction to a subdomain $\omega\in\Theta$ as $u_{\omega}$. Any function $u\in\mathbb{V}$ can be represented by its unique decomposition into subdomain functions as $\{u_{\omega}\in{{V}_{\omega}}\}_{\omega\in\Theta}$. We also define the sub-assembled operator $\mathcal{A}(u,v)\doteq\Pi_{\omega\in\Theta}{\mathcal{A}}_{\omega}(u_{\omega},v_{\omega})$.

With these definitions, $\bar{\mathbb{V}}$ can be understood as the subspace of functions in $\mathbb{V}$ that are continuous on the interface $\Gamma$, and $\bar{\mathcal{A}}$ is the Galerkin projection of $\mathcal{A}$ onto $\bar{\mathbb{V}}$. We note that $\theta_{\rm act}$ and the FE type defines $\bar{\mathbb{V}}$, whereas $\Theta$ is also required to define the local spaces $\{\mathbb{V}_{\omega}\}_{\omega\in\Theta}$ and the sub-assembled space $\mathbb{V}$, respectively. We consider Lagrangian FE spaces, where DOFs are associated to nodes (spatial points). We represent the set of nodes related to the FE space $\mathbb{V}$ with $\mathcal{N}$ (analogusly for $\mathcal{N}_{\omega}$).

3.2. Coarse DOFs

A key ingredient in domain decomposition preconditioners is to classify the set of nodes of the FE space $\bar{\mathbb{V}}$. The interface $\partial e$ of every FE in the mesh $\theta_{\rm act}$ can be decomposed into vertices, edges, and faces. By a simple classification of these entities, based on the set of subdomains that contain them, one can also split the interface $\Gamma$ into faces, edges, and vertices (at the subdomain level), that will be called geometrical objects. We represent the set of geometrical objects by $\Lambda$. In all cases, edges and faces are open sets in their corresponding dimension. By construction, faces belong to two subdomains and edges belong to more than two subdomains in three-dimensional problems. This classification of $\Omega$ into objects automatically leads to a partition of interface DOFs into DOF objects, due to the fact that every DOF in a FE does belong to only one geometrical entity. These definitions are heavily used in domain decomposition preconditioners (see, e.g., [35, p. 88]).

Next, we associate to some (or all) of these geometrical objects a coarse DOF. In BDDC methods, we usually take as coarse DOFs mean values on a subset of objects $\Lambda_{O}$. Typical choices of $\Lambda_{O}$ are $\Lambda_{O}\doteq\Lambda_{C}$, when only corners are considered, $\Lambda_{O}\doteq\Lambda_{C}\cup\Lambda_{E}$, when corners and edges are considered, or $\Lambda_{O}\doteq\Lambda$, when corners, edges, and faces are considered. These choices lead to three common variants of the BDDC method referred as BDDC(c), BDDC(ce) and BDDC(cef), respectively. This classification of DOFs into objects can be restricted to any subdomain $\omega\in\Theta$, leading to the set of subdomain objects $\Lambda_{O}(\omega)$.

With the classification of the interface nodes and the choice of the objects in $\Lambda_{O}$, we can define the coarse DOFs and the corresponding BDDC space. Given an object $\lambda\in\Lambda_{O}(\omega)$, let us define its restriction operator $\mathcal{R}_{\lambda}^{\omega}$ on a function $u\in\mathbb{V}_{\omega}$ as follows: $\mathcal{R}_{\lambda}^{\omega}(u)(\xi)=u(\xi)$ for a node $\xi\in\mathcal{N}_{\omega}$ that belongs to the geometrical object $\lambda$, and zero otherwise. We define the BDDC space $\tilde{\mathbb{V}}\subset{\mathbb{V}}$ as the subspace of functions $v\in{\mathbb{V}}$ such that the constraint

[TABLE]

where ${\rm neigh}(\lambda)$ stands for the set of subdomains that contain the object $\lambda$. (The integral on $\lambda$ is just the value at the vertex, when $\lambda$ is a vertex.) Thus, every $\lambda\in\Lambda_{O}$ defines a coarse DOF value (8) that is continuous among subdomains. Further, we can define the BDDC operator $\tilde{{\mathcal{A}}}$ as the Galerkin projection of ${{\mathcal{A}}}$ onto $\tilde{\mathbb{V}}$.

3.3. Transfer operator

The next step is to define a transfer operator from the sub-assembled space $\mathbb{V}$ to the continuous space $\bar{\mathbb{V}}$. The transfer operator is the composition of a weighting operator and a harmonic extension operator.

1. (1)

The weighting operator $\mathcal{W}$ takes a function $u\in\mathbb{V}$ and computes mean values on interface nodes, i.e.,

[TABLE]

at every node $\xi\in\mathcal{N}$ of the FE mesh $\theta_{\rm act}$, where ${\rm neigh}(\xi)$ stands for the set of subdomains that contain the node $\xi$. It leads to a continuous function $\mathcal{W}u\in\bar{\mathbb{V}}$. It is clear that this operator only modifies the DOFs on the interface. We denote this weighting as topological weighting, since it is defined from the topology of the mesh and the subdomain partition. 2. (2)

Next, let us define the bubble space $\mathbb{V}_{0}\doteq\{v\in\mathbb{V}:v=0\hbox{ on }\Gamma\}$ and the Galerkin projection ${{\mathcal{A}}}_{0}$ of ${{\mathcal{A}}}$ onto ${\mathbb{V}}_{0}$. We also define the trivial injection $\mathcal{I}_{0}$ from $\mathbb{V}_{0}$ to $\bar{\mathbb{V}}$. The harmonic extension reads as $\mathcal{E}v\doteq(I-\mathcal{I}_{0}{{\mathcal{A}}}_{0}^{-1}\mathcal{I}_{0}^{T}{\bar{\mathcal{A}}})v$ for a function $v\in\bar{\mathbb{V}}$, where $I$ is the identity operator in $\bar{\mathbb{V}}$. This operator corrects interior DOFs only. The computation of ${{\mathcal{A}}}_{0}^{-1}$ involves to solve local problems with homogeneous Dirichlet boundary conditions on $\Gamma$. The set of discrete global harmonic functions is $\mathbb{V}_{0}^{\perp}\doteq\{v\in\mathbb{V}:\ {\mathcal{A}}(v,w)=0,\ \forall w\in\mathbb{V}_{0}\}$.

The transfer operator ${\mathcal{Q}}:\mathbb{V}\rightarrow\bar{\mathbb{V}}$ is defined as ${\mathcal{Q}}\doteq\mathcal{E}{\mathcal{W}}$.

3.4. Preconditioner and condition number bounds

With all these ingredients, we are now in position to define the BDDC preconditioner. This preconditioner is an additive Schwarz preconditioner (see, e.g., [35, Ch. 2]), with corrections in $\mathbb{V}_{0}$ and the BDDC correction in $\tilde{\mathbb{V}}$ with the transfer ${\mathcal{Q}}$. As a result, the BDDC preconditioner reads as:

[TABLE]

For body-fitted meshes and second order elliptic problems, the BDDC preconditioned operator, namely ${\mathcal{M}}{\bar{\mathcal{A}}}$, has a condition number bounded as follows (see [25]).

Theorem 3.1**.**

The operator ${\mathcal{M}}{\bar{\mathcal{A}}}$, where ${\mathcal{M}}$ can be the preconditioner BDDC(c) or BDDC(ce) in 2D and BDDC(ce) and BDDC(cef) in 3D, has a condition number bounded by:

[TABLE]

where $\tilde{C}$ is a constant independent from the sizes $h$ and $H$ (the characteristic mesh and subdomain size). In the case of BDDC(ce), every face $\lambda_{F}\in\Lambda_{F}$ must have an edge $\lambda_{E}\in\Lambda_{E}$ on its boundary, i.e., $\lambda_{E}\subset\partial\lambda_{F}$.

Proof.

The BDDC method can be casted in the abstract additive Schwarz setting (see [12, Ch. 7] and [35, Ch. 2]). As a result, the condition number of the preconditioned operator is $k_{2}({\mathcal{M}}{\bar{\mathcal{A}}})=\frac{\lambda_{\rm max}}{\lambda_{\rm min}}$, with $\lambda_{\rm min}=1$ and

[TABLE]

where we have used the triangle inequality and the energy minimization property of discrete harmonic functions in the last bound (see [25]). In order to bound $\lambda_{\rm max}$, we need some further elaboration. We can decompose any discrete harmonic function $v\in\mathbb{V}_{0}^{\perp}$ as the sum of object functions $v=\sum_{\lambda\in\Lambda}\sum_{\omega\in{\rm neigh}(\lambda)}\mathcal{E}\mathcal{R}_{\lambda}^{\omega}(v)$. For body-fitted meshes, it is clear that

[TABLE]

Further, given $v\in\mathbb{V}_{0}^{\perp}$, $\lambda\in\Lambda$, and subdomains $\omega,\omega^{\prime}\in{\rm neigh}(\lambda)$, we define the jump function $\mathcal{J}_{\omega,\omega^{\prime}}^{\lambda}(v)\in\mathbb{V}_{0}^{\perp}$ that is equal to zero in all subdomains but $\omega$, where it take the value $\mathcal{J}_{\omega,\omega^{\prime}}^{\lambda}(v)(\xi)\doteq v_{\omega}(\xi)-v_{\omega^{\prime}}(\xi)$ if $\xi\in\lambda$ and zero otherwise. Thus, we get

[TABLE]

Using standard domain decomposition theory, one can prove the following bound (see, e.g., [25]):

[TABLE]

where $\rho=1$ for edges and 2 for faces. The bound for faces in BDDC(cef) differs from the fact that $\alpha_{\lambda}^{\omega}(v)\neq\alpha_{\lambda}^{\omega^{\prime}}(v)$ in general. It can be bounded similarly relying on the assumption that there exists an edge in $\partial\lambda$ that is in $\Lambda_{O}$ (see [35, p. 182] for more details).

∎

In a weak scaling test (i.e., refining the mesh and the partition at the same rate), $H/h$ is constant and so is the condition number bound regardless of the problem size. This property leads to an optimal weak scaling of the preconditioner (i.e., to increase $h$ and $H$ in the same proportion) when used within a conjugate gradient iteration. That is, the number of iterations needed to solve the problem up to a certain (relative) tolerance is proportional to the condition number $k_{2}({\mathcal{M}}{\bar{\mathcal{A}}})$, and as the condition number is bounded by a constant independent of the problem size, so are also the number of iterations. The coarse DOFs can be reduced, including only a subset of edges in $\Lambda_{O}$. We refer to [24] for more details and the extension to linear elasticity problems.

4. Domain decomposition solvers for unfitted FE meshes

The BDDC method presented in the previous section can be applied verbatim to unfitted meshes. In that case, the computational mesh is $\theta_{\rm act}$ and the extended domain $\Omega_{\rm act}$. It obviously requires specific strategies for integrating the problem matrices within cut elements and to impose Dirichlet conditions on unfitted boundaries, but the BDDC solver can be applied unaltered. The question is, however, whether the performance of BDDC is strongly deteriorated or not when moving from body-fitted to unfitted meshes.

The bound for the condition number in (11) cannot be extended to unfitted methods. One of the problems is that (13) for the physical subdomain $\omega\cap\Omega$, namely

[TABLE]

is not true for unfitted meshes, since a partition can have cut cells on the subdomain interface. As a result, the mathematical theory of substructuring domain decomposition methods, which is grounded on the object decomposition of functions and trace theorems, cannot be readily applied. In fact, it is easy to find examples that prove that the standard BDDC method breaks down for unfitted meshes.

4.1. Breakdown example

Assuming that (13) holds for the unfitted problem at hand (which cannot be enforced in practice), it is easy to check that the standard analysis of BDDC methods, relying on trace inequalities, holds. We observe that $w$ in (14) vanishes on corners for $v\in\tilde{\mathbb{V}}$, since corner values are identical on all subdomains. Thus, only edge and face objects must be considered. Let us consider edge terms. We consider $\lambda\in\Lambda_{E}$ and $\omega,\omega^{\prime}\in{\rm neigh}(\lambda)$. Using the standard analysis of BDDC methods, we note that the mean value of any $v\in\tilde{\mathbb{V}}$ is continuous between subdomains for BDDC(ce) and BDDC(cef), i.e., $\alpha_{\lambda}^{\omega}(v)=\alpha_{\lambda}^{\omega^{\prime}}(v)$. We have:

[TABLE]

where we have used the trace theorem in $H^{1}(\omega)$, the fact that the mean value of functions in $\tilde{\mathbb{V}}$ are continuous between subdomains, the fact that $\mathcal{J}_{\omega,\omega^{\prime}}^{\lambda}(v)$ is discrete harmonic, and (13). More specifically, (13) is used to say that $\mathcal{J}_{\omega,\omega^{\prime}}^{\lambda}(v)=0$ on $\partial\omega\setminus\lambda$. We proceed analogously for faces in BDDC(cef). However, if (13) does not hold, additional terms would appear on the boundary, e.g., the ones related to Nitsche’s method, in the first equality in (17), that cannot be bounded as above.111We note that Nitsche’s methods for unfitted FEMs can require arbitrarily large penalty methods also affected by the small cut cell problem. Even for Neumann boundary conditions on the unfitted boundary only, bounds cannot be independent of cuts, as described with the following example.

Let us consider the problem in Fig. 2 and the BDDC(c) preconditioner in 2D. The left side of the domain is at distance $\epsilon$ to the interface. The underlying PDE is a Poisson problem with Neumann boundary conditions on the left (moving) side and Dirichlet Boundary conditions on the rest of the boundary, defined on a 2D rectangular domain that is dicretized with an unfitted mesh of bilinear Lagrange elements and partitioned into eight uniform subdomains. The example is designed in such a way that, as $\varepsilon$ tends to zero, the area of subdomains $\omega^{1}$ and $\omega^{5}$ tend to zero (see Fig. 2(a)). Now, we consider a function $v\in\tilde{\mathbb{V}}$ that vanishes everywhere but in $\omega^{1}$, where it has value 1 on all DOFs that belong to the edge shared by $\omega^{1}$ and $\omega^{2}$ and zero elsewhere. Thus, the energy of the function before weighting is $\int_{\omega^{1}}|\nabla v|^{2}=C_{1}\epsilon$, where $C_{1}$ is a constant independent from $h$. On the other hand, after communication, the energy related to $\omega^{2}$ is $\int_{\omega^{2}}|\nabla v|^{2}\geq C_{2}$ for some $C_{2}$ independent of $\epsilon$ and $h$. As a result, the condition number is bounded by $\frac{C_{2}}{C_{1}\epsilon}$.222We note that it could not happen if (13) would hold (as in body-fitted methods), because the gradient would increase as $\frac{1}{\epsilon}$. As a result,

[TABLE]

tends to $\infty$ as $\varepsilon\to 0$, i.e., the condition number of the BDDC preconditioner cannot be uniform with respect to the surface intersection. (Similar examples can be designed for 3D problems and the other variants of the preconditioner.) It is clearly observed that the usual topological weighting operator leads to a solver that is not robust with respect to the position of the interface. It is also observed experimentally in Fig. 2 that the iteration count explodes as $\varepsilon$ tends to zero.

4.2. Enlarged coarse spaces and adaptive BDDC methods

In order to fix the BDDC preconditioner for unifitted FE meshes, one could consider as coarse corner DOFs all the nodes that violate (13) and are on the interface. Such approach would lead to a huge space. The authors note that reduced coarse spaces could be considered, especially when the cut cells are related to Neumann boundary conditions. In any case, even though improved solutions can be designed, the coarse space is still too large for practical computations. In the following section, we consider a different approach related to stifness weighting.

Recently, adaptive BDDC methods have been developed based on the idea to use specral solvers to define the coarse space in such a way that the bound for (18) can be controlled by the user (see, e.g., [30, 34, 23, 37, 21, 14]). The analysis of these methods does not rely on trace theorems (a difference with respect to classical BDDC methods). The analysis is algebraic (not relying on trace theorems, functional analysis tools, or FE inverse inequalities) and these algorithms could readily be used for unfitted FE methods to lead to provably robust solvers. However, such approach has two limitations: the additional cost related to the eigenvalue solvers is very large and it only pays the price in some cases, and it would lead to a severe load imbalance (only subdomains intersected by the surface would really require these eigenvalue solvers for many problems of interest). For these reasons, we do not pursue this line in this work.

4.3. Stifness weighting operator

In this section, we will present a modification of the BDDC preconditioner that has been experimentally observed to be robust with respect to small cuts in unfitted FE methods with a moderate increase of the coarse solver size. The use of stiffness matrix-based weighting has previously been used in other situations [16] with good results, e.g., to increase robustness of BDDC preconditioners in isogeometric analysis [9]. Unfortunately, there are no rigorous condition number bounds for this approach.

We can represent functions in $\bar{\mathbb{V}}$ as vectors in $\mathbb{R}^{|\mathcal{N}|}$, and define the symmetric positive definite matrix $A_{ab}\doteq\langle{\mathcal{A}}\phi_{b},\phi_{a}\rangle$ (analogously for the subassembled subdomain matrix $A_{ab}^{\omega}$). We introduce the so-called algebraic weighting operator $\mathcal{W}:\mathbb{V}\rightarrow\bar{\mathbb{V}}$ for a given function $u\in\mathbb{V}$ and a node $\xi$ of the mesh:

[TABLE]

Formula (19) provides the node values after the weighting, which uniquely defines a FE function in $\bar{\mathbb{V}}$. This type of weighting is referred to as the algebraic weighting since it is defined from the problem matrices, in contrast to the topological weighting (see Sect. 3.3) that is only based on the topology of the partition.

As we will show, the algebraic weighting leads to a much more robust BDDC method with respect to the relative position of domain boundary $\partial\Omega$ and the unfitted grid. We consider the 2D example in Fig. 2, for which the BDDC preconditioner with topological weighting is not robust. Fig. 2(b) shows the number of PGC iterations in function of the position of the interface. It is clearly observed that the usual topological weighting operator leads to a solver that is not robust with respect to the position of the interface. The iteration count explodes as $\varepsilon$ tends to zero. On the other hand, the algebraic weighting operator leads to an extremely robust solver. The number of iterations are nearly independent of the position of the interface, even when the value $\varepsilon/h$ reaches machine precision.

Using the algebraic weighting, we have successfully solved also complex 3D examples (see the numerical results in Sect. 5). However, we have observed that the number of iterations needed to solve the systems for meshes with cut elements are (about 2 times) higher that the iterations needed to solve the same problem without cutting the mesh elements. This increase in iterations is not dramatic, but further improvement is possible. To this end, we propose two further enhancements of the BDDC method based on the definition of the coarse objects that lead to almost identical iteration counts than the case in which all cut cells are filled and the standard BDDC method is used.

4.4. New coarse objects: variant 1

The current mathematical analysis of BDDC methods is for body-fitted meshes, and therefore, it implicitly assumes that the coarse objects are not intersected by the boundary $\partial\Omega$. However, when body-fitted meshes are allowed, the objects can be cut by $\partial\Omega$ (see Fig. 3(a)). In this situation, the main mathematical properties required to prove scalability of the BDDC method are lost. It can lead to situations in which the corner constraints are exterior corners of cells with close to zero $\eta_{e}$, possibly leading to close to singular local problems. Motivated by this fact, we propose a modification in the definition of the coarse objects, which eliminates intersections with the boundary $\partial\Omega$. Since fixing the elements of the kernel (the constant for the Laplacian or rigid-body modes for elasticity) for all the edges is enough in 3D (see, e.g., [24]), we will only consider the modification of edge coarse DOFs.

The new objects are computed in two main steps. The first one is to build the standard objects $\Lambda_{O}$ of the BDDC method for $\theta_{\rm act}$, as presented in Sect. 3.2. The second step is to identify the edges in $\Lambda_{O}$ cut by the boundary and split them into new edges/corners, leading to uncut edges (see Fig. 3). We consider a segregation of a coarse edge $\lambda\in\Lambda_{E}$ as follows. Given the set of FE edges $t\subset\lambda$, we segregate them into interior edges $t\subset\Omega$, exterior edges such that $t\cap\Omega=\emptyset$, and the remaining cut edges. All DOFs on cut edges are considered as coarse corner DOFs. The rest of FE edges are agregated into connected subedges (see Fig. 3(b)). This procedure is done efficiently by using the data structures commonly available in domain decomposition codes and immersed boundary methods. It requires to perform a loop on the edges objects (which are usually identified anyway in a domain decomposition code), and determine if those edges are intersected by the interface or not. Intersection points are usually already pre-computed in immersed boundary codes. As a result, no further significant overhead is introduced. Moreover, this operation has to be performed only in the subdomains cut by the boundary, which is typically a small set. The numerical examples (see Sect.5) show that by splitting edge objects in this way, optimal weak scaling is recovered if Neumann boundary conditions are imposed on the unfitted parts of the domain boundary $\partial\Omega$.

4.5. New coarse objects: variant 2

However, the case of weak imposition of Dirichlet conditions has been observed to be more challenging and requires further elaboration. In that case, we propose a second strategy to split the edges. The method is motivated by the fact that the domain decomposition analysis of non-spectral methods relies on trace theorems, and it requires constant weighting coefficients within each coarse object in order to derive condition number bounds (see Eq. (14) and, e.g., [35, p. 141], [25], or [12, p. 208]).333We note that the numerical analysis of more involved non-diagonal weighting operators (e.g., deluxe scaling, based on spectral methods) does not rely on trace theorems. We note that the coefficient $\beta$ required to ensure coercivity in Nitsche’s method might be very different between neighbor elements (even different orders of magnitude) in function of the location of the interface. This induces strong variations in the coefficients of the stiffness matrix, and therefore highly oscillatory weighting coefficients, when they are computed as in (19). As a result, the weighting coefficient might be non-constant within the BDDC objects cut by the interface (see Fig. 4(a)), which violates one of the basic assumptions of the domain decomposition analysis. Motivated by this fact, we try to build objects that have constant weighting coefficients. To this end, we follow two steps. First, we build the standard BDDC objects (Sect. 3.2). Then, we split all coarse objects classified as edges that have non-constant weighting coefficients. Precisely, we aggregate neighbor nodes on the edge that have the same weighting coefficient (note that many nodes on a coarse edge have the same coefficient, up to machine precision, for a structured meshes and linear elements), and then we set the remaining un-aggregated nodes as coarse corners (see Fig. 4(b)). Since edges are one dimensional curves, this aggregation can be easily performed. Moreover, this process is performed only for the small set of subdomains that are cut by the boundary. By modifying the coarse objects in this way, an optimal weak scaling is also recovered for weak Dirichlet boundary conditions, as showed in the following 3D numerical examples.

5. Numerical experiments

5.1. Setup

In the numerical examples below, we consider the Poisson problem with the discretization being used in Sect. 2.2 defined on five different complex 3D domains shown in Fig. 5. The chosen geometries are 1) a sphere, 2) a body that reminds to a popcorn flake, 3) a hollow block, 4) a 3-by-3-by-3 array of such blocks, and 5) a spiral. These geometries are used often in the literature to study the performance of unfitted FE methods (see, e.g., [13]). This variety of shapes is considered here to evaluate the versatility of the proposed solvers and their ability to deal with complex and diverse 3D domains. A detailed numerical experimentation is basic, since the solvers proposed are based on heuristic motivations, and numerical experimentation is the only way to prove their robustness.

For each geometry, weak scalability tests are executed in order to evaluate the performance of the solvers. To this end, a family of computational meshes are generated for each piece spanning from coarse meshes up to fine meshes with several millions of elements (see Table 1 for details). Note that the subdomain partitions are refined at the same rate than the meshes, leading to a constant ratio $H/h$, which should lead to an optimal weak scaling of the BDDC preconditioner. That is, the number of iterations of the underlying preconditioned conjugate gradient solver is expected to be (asymptotically) independent of the problem size. The main purpose of the numerical example below is to check that the proposed BDDC preconditioners for unfitted meshes have this optimal performance. In the remainder of this sub-section, we further detail the setup of the numerical examples.

We define the source term and boundary conditions of the Poisson equation such that the PDE has an exact solution, namely

[TABLE]

We use the exact solution to assess the quality of the numerical approximations provided by the solvers. In order to define the boundary conditions, we partition the boundary $\partial\Omega$ into two non-overlapping sets $\Gamma_{\mathrm{bb}}$ and $\Gamma_{\mathrm{cut}}$ (see Fig. 6). The set $\Gamma_{\mathrm{bb}}$ is the intersection of $\partial\Omega$ with the bounding box used to define the computational mesh, and $\Gamma_{\mathrm{cut}}$ is its complement, $\Gamma_{\mathrm{cut}}\doteq\partial\Omega-\Gamma_{\mathrm{bb}}$. After discretizing the bounding box with a Cartesian mesh, the set $\Gamma_{\mathrm{bb}}$ results to be aligned with the mesh faces. On the contrary, $\Gamma_{\mathrm{cut}}$ is a non-conforming boundary living inside the cut elements. We define strong Dirichlet conditions on the boundary $\Gamma_{\mathrm{bb}}$. Since $\Gamma_{\mathrm{bb}}$ is aligned with the mesh, strong Dirichlet boundary conditions are imposed by constraining the corresponding mesh nodes. On the other hand, we impose two different types of boundary conditions on the unfitted surface $\Gamma_{\mathrm{cut}}$, either weak Dirichlet boundary conditions using Nitsche’s method or Neumann boundary conditions. Note that for the popcorn flake and the spherical domain, the set $\Gamma_{\mathrm{bb}}$ is empty since all the boundary $\partial\Omega$ in non-conforming (see Fig. 6). In these cases, we strongly impose the value of the solution at an interior node of $\Omega$ in order to have a well posed problem when considering Neumann conditions on the entire boundary.

After the discretization process, the resulting systems of linear equations $Ax=b$ are solved with a BDDC-preconditioned conjugate gradient (PCG) solver. Convergence of the PCG solver is declared when the 2-norm of the algebraic residual $r\doteq b-Ax$ is smaller that the 2-norm of the right-hand-side $b$ times a tolerance, namely $|r|_{2}<\mathrm{tol}\cdot|b|_{2}$. In all cases, the selected tolerance is $\mathrm{tol}=10^{-9}$. With this value, the error committed in solving the linear system is negligible in front of the discretization error associated with the FE mesh. This is demonstrated by the fact that the optimal convergence rate of FEM is obtained by using this tolerance for all the geometries (see Fig. 7).

Six different versions of the BDDC preconditioner are studied in the examples, see Table 2. They differ in the definition of the coarse objects (either the Standard definition of section 3.2 or the extended definitions for unfitted meshes in Sects. 4.4 and 4.5) and the subset of objects used to define the coarse DOFs (either corners and edges, or corners edges and faces). In all the cases, the algebraic weighting operator of Sect. 4.3 is considered. We do not include the topological weighting operator (detailed in Sect. 3.3) since it leads to a solver that is not robust with respect to the relative position between the geometry and the mesh (see, e.g., the breakdown example in Sect. 4.1). We have experimentally observed that the topological weighting operator leads to unacceptably large iteration counts for all the 3D geometries studied here. In particular, more than 1000 iterations are needed to solve the “popcorn flake” with a medium-size mesh (mesh # 3 described in Table 1) for the case of Neumann boundary conditions on $\Gamma_{\mathrm{cut}}$. Thus, the breakdown of the method does not only appear in manufactured examples as the one in 4.1 but in all practical examples analyzed in this section.

Finally, in order to compare the performance of the proposed BDDC methods for unfitted grids, with respect to the optimal behavior of BDDC for conventional conforming meshes, we construct auxiliary meshes by replacing the cut elements by standard full elements (see Fig. 8). Our goal is to recover a similar performance for cut elements, that the one obtained with full elements.

5.2. “Popcorn flake” example.

For the sake of clarity, we first present the numerical results for a single geometry, namely the “popcorn flake” geometry in Fig. 5(b). The results for the other shapes are discussed in Sect. 5.3 below. Fig. 9 shows the result of the weak scaling test for this geometry, where Fig. 9(a) is for the case of Neumann boundary conditions on $\Gamma_{\mathrm{cut}}$, and Fig. 9(b) is for weak Dirichlet boundary conditions. It is observed that the standard BDDC method leads to a perfect weak scaling for the case of full elements (as expected) for the two types of boundary conditions. The performance of the standard BDDC method when dealing with cut elements (i.e., red lines in Fig. 9) is decent, in the sense that the number of iterations increases only mildly with the problem size, but it is clear that the absolute number of iterations needed to solve the problems is more than twice the iterations required for full elements (which is the target value). In contrast, both the 1st and 2nd variants of BDDC achieve a perfect weak scaling with iteration counts close to the ones associated with full elements for the case of Neumann boundary contitions (see blue and yellow lines in Fig. 9(a)).

However, for Dirichlet boundary conditions, only the 2nd variant provides a perfect weak scaling close to the one observed for full elements (see yellow line in Fig. 9(b)). In conclusion, the 1st variant of the proposed method is enough to achieve a perfect weak scaling for Neumann boundary conditions, whereas the 2nd variant is required for Dirichlet conditions imposed with Nitsche’s method. Similar results are obtained if the coarse space is defined with corners, edges, and faces (cef), or only with corners and edges (ce) (see Fig. 10). The main difference is that the iterations associated with corners and edges are slightly larger than for corners, edges, and faces, which is the expected behavior in BDDC methods. The iteration count is smaller in the later case since the associated coarse space is larger.

Fig. 11 shows the size of the coarse space of the new BDDC methods, namely variants 1 and 2, with respect to the size of the standard coarse space. Fortunately, the size of the coarse spaces associated with variants 1 and 2 tends to be equal to the size of the standard coarse space in BDDC as the mesh is refined (see Fig. 11). This behavior is explained by the fact that extra corners are only added on edges that are close to cut elements. Since the number of BDDC objects in the interior of $\Omega$ scales much faster that the objects cut by the boundary $\partial\Omega$, the number of extra corners added near the boundary become nearly negligible as the partition is refined. As a result, the finer is the mesh the more competitive are the proposed BDDC preconditioners for unfitted meshes. Thanks to this behavior and for the perfect weak scaling results in [6, 4], the proposed BDDC methods (variants 1 and 2) are suitable for large-scale problems requiring fine meshes and a large number of subdomains. For Neumann boundary conditions, variant 1 is the best choice since it leads to a smaller coarse space than variant 2 for approximately the same number of PCG iterations. For Dirichlet boundary conditions with Nitsche’s method, the 2nd variant is the preferred choice specially for fine meshes, since it leads to a better weak scaling.

All the results presented so far are obtained considering a fixed ratio $H/h=8$ between the characteristic sub-domain size $H$ and the characteristic mesh size $h$. We evaluate the effect of $H/h$ on the number of iterations of the three variants of the BDDC preconditioner for the popcorn flake geometry in Figures 12 and 13. As expected, the effect of $H/h$ (measure of the local problem size) on the number of iterations is very mild. An algorithmic strong scalability analysis, i.e., to fix the global problem and consider $P$, $8P$, and $64P$ processors, can be obtained by looking for a given value of DOFs at the number of iterations for the three cases under consideration. These figures have not been included for the sake of brevity.

It is also desirable that the corner detection mechanisms in variants 1 and 2 do not have a negative impact in the load balancing of the local operations involved in the BDDC preconditioner, namely the number of local Neumann problems required to compute the coarse basis functions, which is the main algorithmic source of load imbalance (the number of local Dirichlet problems and constrained Neumann problems in the application step are identical for all subdomains). The effect of the proposed method on the load balancing is studied in Fig. 14. Here, the extra extra corners added by 1st and 2nd variants of the new method can potentially affect the load balancing. However, the load imbalance in terms of the coarse DOFs per subdomain is not getting worse as the mesh is refined (see Fig 14). The results in Fig. 14 are for the 2nd variant of the BDDC(ce) method. Even better results are obtained for the 1st variant and the BDDC(cef) method. The plots are not included for the sake of brevity since they do not add extra information.

In a final phase of the experiment, we study the robustness of the methods with respect to the position of the cut geometry. To this end, we slightly move the location of the background mesh by an arbitrary value $\varepsilon$ in one of the three spatial directions. By translating the mesh in this way, the location of the geometry changes and creates different kinds of intersections with the mesh elements, leading to cut elements with different active volume fractions that can affect the conditioning of the problem in different ways. Our goal is to show that the proposed methods are not affected by this translation and the different configurations of cut elements. The values for $\varepsilon$ are chosen in the closed interval $[-h,+h]$, where $h$ is the element size of the background mesh. Fig. 15 shows the number of iterations of the PCG solver as a function of $\varepsilon$, for the mesh number 3 in Table 1. The results are for weak Dirichlet boundary conditions, which is the most challenging case. The results in Fig. 15 show that the 2nd variant is very robust with respect to the position of the interface. The maximum change in the number of iterations is only 1 iteration for all the 51 different values of $\varepsilon$ considered. The other methods have also a decent robustness with a maximum change of 4 iterations.

5.3. Multi-body example: sphere, popcorn flake, block, array of blocks and spiral.

The purpose of this final example is to show that the proposed BDDC preconditioners are able to deal with different complex geometries, and that the results obtained in previous example are qualitatively reproduced also in other settings. To this end, we study the BDDC preconditioners for all the geometries previously presented in Fig. 5. The results of the weak scalability test are presented in Fig. 16. The figure shows that performance of the methods for this set of geometries is similar to the one previously obtained for the “popcorn flake”. For Neumann boundary conditions a perfect weak scaling is obtained for both the 1st and 2nd variants, whereas for Dirichlet boundary conditions a perfect weak scaling is only obtained with the 2nd variant. It is remarkable that the number of iterations is almost identical for all the geometries for the finest meshes (see 1st and 2nd variants for Neumann conditions, and 2nd variant for Dirichlet conditions). This is true even for the array of blocks, which is the most complex geometry considered herein.

The size of the coarse spaces behaves similar to the previous example (see Fig. 17). The increment in the coarse space size associated with variants 1st and 2nd becomes smaller in front of the size of the coarse space in standard BDDC as the mesh is refined. This is again due to the fact that the number of objects in the bulk of the domain scales faster than the objects cut by the boundary. For that reason, the sphere is the shape that presents a smaller increment in the coarse space size. It minimizes the cut surface per unit of domain volume. Note, however, that the coarse space size of the variants 1 and 2 tends also to the coarse space size of standard BDDC for the other shapes (see, e.g., the spiral or the array of blocks).

Finally, we study the robustness of the methods with respect to the position of the unfitted boundary $\partial\Omega$. In all cases, we move the background mesh a value $\varepsilon\in[-h,+h]$ to see the influence of the boundary location. Fig. 18 shows the number of iterations of the PCG solver with respect to the value $\varepsilon$ for all the shapes. The results correspond to weak Dirichlet boundary conditions, which is the more challenging case. The results are similar to previous example. It is observed that the performance of the 2nd variant is very robust. It is almost independent of the position of the interface with a maximum change of 1 iteration in all the considered geometries.

In summary, the presented results show that the proposed variants of the BDDC method are able to provide excellent results for different and complex geometries discretized with unfitted FE meshes. As in previous example, the 1st variant is already enough to achieve a perfect weak scaling for all the shapes for the case of Neumann boundary conditions. The 2nd variant is the only one that leads to perfect weak scaling for Dirichlet conditions. It is remarkable that the number of PCG iterations is nearly independent from the geometry for the finest meshes. This shows that the methods are able to effectively get rid of the complexities introduced by the considered geometries and cut elements. The size of the coarse spaces of the new methods tends to the size for standard BDDC as the mesh is refined. Therefore, the efficiency of the new methods is comparable to the one of standard BDDC for fine meshes.

6. Conclusions

In this work, we have analyzed the main problems related to the use of BDDC methods to unfitted FE meshes. From a mathematical point of view, there are different assumptions in the analysis of these algorithms that do not hold anymore for unfitted meshes. As a result, the polylogarithmic condition number bounds in BDDC methods are not true in general. In fact, using a breakdown example, we have observed that the condition number of these methods cannot be robust with respect to the intersection between the surface and the background mesh. One can fulfill again the assumptions of the BDDC analysis by considering as corner constraint all the DOFs on the interface that belong to cut cells. Unfortunately, even thought the resulting algorithm is perfectly robust, such approach would lead to a dramatic increase in the coarse problem size, and it is not considered in practice.

Alternatively, and motivated by both numerical experimentation and the mathematical analysis of BDDC methods, we have considered a stiffness-based weighting and a subdivision of the BDDC edges. The first ingredient has turned to be an essential ingredient to get a robust BDDC solver for unfitted meshes. On the other hand, the enrichment of the coarse BDDC space is especially important in regions with weak imposition of Dirichlet boundary conditions. As usual, the use of stiffness weighting is not backed up with a full mathematical theory. Instead, a complete set of numerical experiments has been performed to show the robustness and algorithmic weak scalability of the proposed solvers.

Future work will include the extension to more complex (non-coercive) physical problems and to ghost penalty stabilization. The use of stabilized versions is expected to be easier than the problem considered herein, since the condition number bound in these cases can be bounded independently of the location of the intersection. Furthermore, we want to implement these solvers within the FEMPAR library [5, 1], in order to show weak scalability results for its highly scalable implementation [3, 4] of multilevel BDDC preconditioners. We also want to apply these methods to adaptive octree meshes and space-filling curves [2].