A minimal stabilization procedure for Isogeometric methods on trimmed geometries

TL;DR

This paper introduces a minimal stabilization technique for Isogeometric Analysis on trimmed geometries, addressing stability issues and ensuring accurate solutions for elliptic PDEs on complex CAD-derived domains.

Contribution

The paper presents a novel stabilization method for Isogeometric methods on trimmed geometries, improving stability and accuracy in solving elliptic PDEs.

Findings

Stability issues occur with Nitsche's method on trimmed domains.

The proposed stabilization guarantees well-posedness and optimal error estimates.

Numerical examples confirm theoretical stability and accuracy improvements.

Abstract

Trimming is a common operation in CAD, and, in its simplest formulation, consists in removing superfluous parts from a geometric entity described via splines (a spline patch). After trimming the geometric description of the patch remains unchanged, but the underlying mesh is unfitted with the physical object. We discuss the main problems arising when solving elliptic PDEs on a trimmed domain. First we prove that, even when Dirichlet boundary conditions are weakly enforced using Nitsche's method, the resulting method suffers lack of stability. Then, we develop novel stabilization techniques based on a modification of the variational formulation, which allow us to recover well-posedness and guarantee accuracy. Optimal a priori error estimates are proven, and numerical examples confirming the theoretical results are provided.

Equations

$$\Omega ={\Omega }_0\setminus \bigcup _{i=1}^N{\overline{\Omega }}_i,$$

$${\xi }_1=\cdots ={\xi }_{p+1}<{\xi }_{p+2}\le \cdots \le {\xi }_n<{\xi }_{n+1}=\cdots ={\xi }_{n+p+1}.$$

$${\tilde{I}}_j:=\mathrm{int}\bigcup \{\mathrm{supp}({\hat{B}}_{k,p}):\mathrm{supp}({\hat{B}}_{k,p})\cap I_j\ne \varnothing ,1\le k\le n\}=\left({\xi }_{i-p},{\xi }_{i+p+1}\right).$$

$${\hat{\mathcal{M}}}_0=\{Q_{\mathbf{j}}=I_{1,j_1}\times \cdots \times I_{d,j_d}:I_{\ell ,j_{\ell }}=({\zeta }_{\ell ,j_{\ell }},{\zeta }_{\ell ,j_{\ell }+1}):1\le j_{\ell }\le M_{\ell }-1\},$$

$${\mathcal{M}}_{0,h}:=\{K\subset \Omega :K=\mathbf{F}(Q),Q\in {\hat{\mathcal{M}}}_{0,h}\}.$$

$$V_h=\mathrm{span}\{B_{\mathbf{i},p}(\mathbf{x}):={\hat{B}}_{\mathbf{i},p}\circ {\mathbf{F}}^{-1}(\mathbf{x}):\mathbf{i}\in \mathbf{I}\},$$

$$\{\begin{array}{ll}-\Delta u=f& \text{in}\Omega ,\\ u=g_D& \text{on}{\Gamma }_D,\\ \frac{\partial u}{\partial n}=g_N& \text{on}{\Gamma }_N,\end{array}$$

$${\int }_{\Omega }\nabla u_h\cdot \nabla v_h$$

$$=$$

$$a_h(u_h,v_h):={\int }_{\Omega }\nabla u_h\cdot \nabla v_h-{\int }_{{\Gamma }_D}\frac{\partial u_h}{\partial n}{v}_h-{\int }_{{\Gamma }_D}{u}_h\frac{\partial v_h}{\partial n}+\beta {\int }_{{\Gamma }_D}{\mathsf{h}}^{-1}{u}_h{v}_h.$$

$${\left(u_h,v_h\right)}_{1,h,\Omega }:={\int }_{\Omega }\nabla u_h\cdot \nabla v_h+{\int }_{{\Gamma }_D}{\mathsf{h}}^{-1}{u}_h{v}_h,$$

$${\Vert u_h\Vert }_{1,h,\Omega }^2:={\Vert \nabla u_h\Vert }_{L^2(\Omega )}^2+{\Vert {\mathsf{h}}^{-\frac{1}{2}}{u}_h\Vert }_{L^2({\Gamma }_D)}^2.$$

$$\alpha {\Vert u_h\Vert }_{1,h,\Omega }^2\le a_h(u_h,u_h)\forall u_h\in {\tilde{V}}_h,$$

$$a_h(u_h,v_h)\le \gamma {\Vert u_h\Vert }_{1,h,\Omega }{\Vert v_h\Vert }_{1,h,\Omega }\forall u_h,v_h\in {\tilde{V}}_h.$$

$${\int }_{\Omega }\nabla u_h\cdot \nabla v_h-{\int }_{{\Gamma }_D}\frac{\partial u_h}{\partial n}{v}_h-{\int }_{{\Gamma }_D}{u}_h\frac{\partial v_h}{\partial n}+{\int }_{{\Gamma }_D}{\mathsf{h}}^{-1}{u}_h{v}_h={\lambda }_h(u_h,v_h{)}_{1,h,\Omega }\forall v_h\in {\tilde{V}}_h.$$

$${\Vert {\mathsf{h}}^{\frac{1}{2}}\frac{\partial v_h}{\partial n}\Vert }_{L^2({\Gamma }_K)}\le C{\Vert \nabla v_h\Vert }_{L^2(K\cap \Omega )},$$

$$\frac{\left|\hat{\Omega }\cap Q\right|}{\left|Q\right|}\ge \theta .$$

$$\mathcal{N}(K):=\{K^{\prime}\in\mathcal{M}_{h}:\operatorname{dist}\left(K,K^{\prime}\right)\leq C{\left.\kern-1.2pt\mathsf{h}\vphantom{\big{|}}\right|_{K}}\}\setminus\{K\},$$

$$R_h:{\tilde{V}}_h\to L^2({\Gamma }_D)$$

$${\overline{a}}_h(u_h,v_h)={\int }_{\Omega }f{v}_h+{\int }_{{\Gamma }_N}{g}_N{v}_h-{\int }_{{\Gamma }_D}{g}_D{R}_h(v_h)+\beta {\int }_{{\Gamma }_D}{\mathsf{h}}^{-1}{g}_D{v}_h\forall v_h\in {\tilde{V}}_h,$$

$${\overline{a}}_h(u_h,v_h):={\int }_{\Omega }\nabla u_h\cdot \nabla v_h-{\int }_{{\Gamma }_D}{R}_h(u_h)v_h-{\int }_{{\Gamma }_D}{u}_h{R}_h(v_h)+\beta {\int }_{{\Gamma }_D}{\mathsf{h}}^{-1}{u}_h{v}_h.$$

$${\Vert {\mathsf{h}}^{\frac{1}{2}}{R}_h(v_h)\Vert }_{L^2({\Gamma }_K)}\le C{\Vert \nabla v_h\Vert }_{L^2(\Omega \cap K^{\prime })}\forall v_h\in {\tilde{V}}_h,$$

$$\left|{\overline{a}}_h(u_h,v_h)\right|\le {\Vert \nabla u_h\Vert }_{L^2(\Omega )}{\Vert \nabla v_h\Vert }_{L^2(\Omega )}+{\Vert {\mathsf{h}}^{\frac{1}{2}}{R}_h(u_h)\Vert }_{L^2({\Gamma }_D)}{\Vert {\mathsf{h}}^{-\frac{1}{2}}{v}_h\Vert }_{L^2({\Gamma }_D)}$$

$$+{\Vert {\mathsf{h}}^{\frac{1}{2}}{R}_h(v_h)\Vert }_{L^2({\Gamma }_D)}{\Vert {\mathsf{h}}^{-\frac{1}{2}}{u}_h\Vert }_{L^2({\Gamma }_D)}+\beta {\Vert {\mathsf{h}}^{-\frac{1}{2}}{u}_h\Vert }_{L^2({\Gamma }_D)}{\Vert {\mathsf{h}}^{-\frac{1}{2}}{v}_h\Vert }_{L^2({\Gamma }_D)}$$

$$\le$$

$$+\beta {\Vert u_h\Vert }_{1,h,\Omega }{\Vert v_h\Vert }_{1,h,\Omega }\le C{\Vert u_h\Vert }_{1,h,\Omega }{\Vert v_h\Vert }_{1,h,\Omega },$$

$$a_h(u_h,u_h)\ge$$

$$+\beta {\Vert {\mathsf{h}}^{-\frac{1}{2}}{u}_h\Vert }_{L^2({\Gamma }_D)}^2\ge \left(1-\frac{C}{\delta }\right){\Vert \nabla u_h\Vert }_{L^2(\Omega )}^2+\left(\beta -\delta \right){\Vert {\mathsf{h}}^{-\frac{1}{2}}{u}_h\Vert }_{L^2({\Gamma }_D)}^2,$$

$$R_{K}(v_{h}):=\frac{\partial{\left.\kern-1.2ptv_{h}\vphantom{\big{|}}\right|_{K}}}{\partial n};$$

$$R_{K}(v_{h}):=\frac{\partial\left(\mathcal{E}\left({\left.\kern-1.2pt\widehat{v}_{h}\vphantom{\big{|}}\right|_{Q^{\prime}}}\right)\circ\mathbf{F}^{-1}\right)}{\partial n},$$

$$R_{K}(v_{h}):=\frac{\partial{\left.\kern-1.2ptv_{h}\vphantom{\big{|}}\right|_{K}}}{\partial n};$$

$$R_{K}(v_{h}):=\frac{\partial\left(\mathcal{E}\left(P({\left.\kern-1.2ptv_{h}\vphantom{\big{|}}\right|_{K^{\prime}}})\right)\right)}{\partial n},$$

$${\Vert {\mathsf{h}}^{\frac{1}{2}}{R}_h(v_h)\Vert }_{L^2({\Gamma }_K)}\le C{\Vert \nabla v_h\Vert }_{L^2(\Omega \cap K^{\prime })}\forall v_h\in {\tilde{V}}_h,$$

Full text

{AMS}

65N12, 65N15, 65N30, 65N85

1 Introduction

Complex models are processed within Computer Aided Design (CAD) tools where several geometric manipulations are possible. Geometries are described as collection of their boundary surfaces, often defined as tensor-product splines or NURBS, and during the design process these surfaces can be joined, intersected or simply superflous parts can be cut away. All these Boolean operations act on the original surfaces through a common procedure of trimming. When the superfluous surface areas are cut away, the visualization of the resulting surface changes, while its mathematical description does not. This description of the geometry is called “boundary representation” (B-rep) (see, e.g., [38, 45], or the recent review [31] and references therein) and is clearly not well suited for the simulation of PDEs.

Several efforts have been undertaken in the last years to improve the usability of CAD geometries in the solution of PDEs, especially thanks to the advent of Isogeometric Analysis (IGA) [15, 27] and its tremendous success (see, e.g., [1, 4, 6, 26, 42]). The geometric modelling community has also provided important inputs to this scientific challenge [41], and, in this respect, volumetric representations (V-rep) are a major contribution [33]. On the other hand, trimming remains a main tool for the design of complex models via Boolean operations, and basically all developments of IGA described above rely on strong requirements on the underlying geometric models and, in general, do not support trimmed geometric entities. The aim of this paper is to contribute to the design of isogeometric methods that robustly support trimming in the geometric description of the computational domain, and do not require the construction of a global re-parametrization (meshing).

Two main issues arising when dealing with trimmed geometries are the presence of elements unfitted with the boundary, making the research for efficient quadrature rules and the stable imposition of boundary conditions a challenge, and the existence of basis functions whose support has been cut, affecting the conditioning of the related linear system. Let us briefly review some of the most successful methods which have been proposed so far in this regard.

In connection with the construction of quadrature rules on trimmed elements let us mention, from the engineering side, the pioneering works [35, 48], the shell analysis on geometric models with B-reps [8, 37], and the finite cell method combined with IGA [17, 40, 43, 49].

Concerning the lack of stability and the ill-conditioning of the stiffness matrix we should mention stabilization based on polynomial extrapolation in the parametric domain [30, 32] and the so-called Cut-IGA method proposed in [20, 28]. The former traces its roots back to the groundbreaking work [25] where the authors, by modifying the discrete functional space, tackle both the stability and the conditioning issues. The theory was proved in [25] for finite elements with (extended) B-splines, which corresponds to consider below the map $\F$ equal to the identity (i.e. without the isogeometric map), and the polynomial extrapolation may lead to sub-optimal convergence properties, as it will be discussed in Section 6. The latter is a generalization of the Cut-FEM method [11, 12, 13, 14, 23] born in the framework of fictitious domain methods for finite elements. It is based on a modification of the weak formulation of the problem by the addition a penalisation term on the boundary. This technique aims to deal at the same time with the ill-conditioning and the lack of stability of the bilinear form. Note that the penalisation term acts on the jumps of the normal derivatives of all orders over cut elements’ boundaries, which can be very demanding if high-order B-splines are employed for the analysis (which is usually the case in IGA).

We now describe our simplified mathematical setting. Let $\Omega_{0}\subseteq\R^{d}$ (here $d=2,3$) be a domain described by a bijective spline map $\mathbf{F}:(0,1)^{d}\to\Omega_{0}$, i.e. a patch in the isogeometric terminology, and let $\Omega_{1},\ldots,\Omega_{N}$ be bounded domains of $\R^{d}$. We assume that $\bigcup_{i=1}^{N}\overline{\Omega}_{i}$ are to be cut away from $\Omega_{0}$, i.e.

[TABLE]

and we assume that the computational domain $\Omega$ is Lipschitz. After trimming the mathematical description of the domain remains unchanged, that is, the elements and basis functions fit the boundary of $\Omega_{0}$ instead of that of $\Omega$. In this paper, we focus on a simple Poisson problem, with weakly imposed boundary conditions, in the domain $\Omega$ described above. First of all, we discuss the difference between bad matrix conditioning and lack of stability. The former can be improved by modifying the chosen basis (preconditioning), while the latter needs to act on the bilinear form directly.

Regarding the lack of stability, we propose a stabilization technique, inspired by [24], that acts only on those cut elements that affect stability. For example, in the case of a Neumann condition on the trimmed boundary, no stabilization is needed. The stabilization is “minimal” in the sense that no additional parameters are introduced, in contrast with the CutFEM [13] and Finite Cell methods [17] for instance. Our stabilization is parameter free and its computation requires only local projections at the element level, and only for “bad” cut elements. We follow a “local approach” [31], i.e. we modify the analysis, rather than the geometry, in order to be able to face the challenges arising from trimming. Moreover, we remain faithful to the so-called isogeometric paradigm, in the sense that we just locally modify the weak formulation, while keeping the discrete functional space unaffected.

We present two different versions of the stabilization. The first one is based on polynomial extrapolation in the parametric domain, that is easier to implement from the numerical point of view, but suboptimal in some cases. The second one is a projection-based stabilization performed directly on the physical domain, which allows us to recover optimal a priori error estimates.

Concerning the conditioning issue, we do not have a sound solution to the problem, but we constructed tests to check the behaviour of the condition number of the stiffness matrix. Numerical evidences show that a rescaling of the stiffness matrix, coupled with our stabilization, greatly reduces the condition number, although it does not solve the conditioning issue in all configurations (see Section 8.4). A clear theoretical understanding of this issue is beyond the scope of this paper and, in this regard, the interested reader is referred to [34] where a $L^{p}$-stable basis is constructed in the context of B-splines and to [18, 19].

The document is organized as follows. Section 2 presents an overview of isogeometric analysis in trimmed domains. In Section 3 we set the model problem, and explain the main challenges we need to face, namely integration, conditioning and numerical stability of the associated linear system. After having explained in detail in Section 4 the causes for the lack of stability of the simple Nitsche’s formulation, in Section 5 we present our new stabilization technique. Two possible constructions of the stabilization operator are suggested and analysed in Section 6, and error estimates are provided in Section 7. Finally, we conclude by showing some numerical examples, obtained using the MATLAB library GeoPDEs [47], confirming the theoretical results.

2 Isogeometric analysis on trimmed domains

2.1 The univariate case

For a more detailed introduction to isogeometric analysis, we refer the interested reader to the review article [4]. Given two positive integers $p$ and $n$, we say that $\Xi:=\{\xi_{1},\dots,\xi_{n+p+1}\}$ is a p-open knot vector if

[TABLE]

We assume $\xi_{1}=0$ and $\xi_{n+p+1}=1$. We also introduce $Z:=\{\zeta_{1},\dots,\zeta_{M}\}$, the set of breakpoints, or knots without repetitions, which forms a partition of the unit interval $(0,1)$. Note that $\Xi=\{\underbrace{\zeta_{1},\dots,\zeta_{1}}_{m_{1}\;\text{times}},\underbrace{\zeta_{2},\dots,\zeta_{2}}_{m_{2}\;\text{times}},\dots,\underbrace{\zeta_{M},\dots,\zeta_{M}}_{m_{M}\;\text{times}}\},$ where $\sum_{i=1}^{M}m_{i}=n+p+1$. Moreover, we assume $m_{j}\leq p$ for every internal knot and we denote $I_{i}:=(\zeta_{i},\zeta_{i+1})$ and its measure $h_{i}:=\zeta_{i+1}-\zeta_{i}$, $i=1,\dots,M-1$.

We denote as $\widehat{B}_{i,p}:[0,1]\to\R$ the i-th B-spline, $1\leq i\leq n$, obtained using the Cox-de Boor formula, see for instance [4]. Moreover, let $S_{p}(\Xi)=\operatorname{span}\{\widehat{B}_{i,p}:1\leq i\leq n\}$ the vector space of univariate splines of degree $p$. $S_{p}(\Xi)$ can also be characterized as the space of piecewise polynomials of degree $p$ with $k_{j}:=p-m_{j}$ continuous derivatives at the breakpoints $\zeta_{j}$, $1\leq j\leq M$ (Curry-Schoenberg theorem).

Moreover, given an interval $I_{j}=\left(\zeta_{j},\zeta_{j+1}\right)=(\xi_{i},\xi_{i+1})$, we define its support extension $\widetilde{I}_{j}$ as

[TABLE]

2.2 The multivariate case

Let $d$ be the space dimension. Assume that $M_{\ell},n_{\ell}\in\N$, $p\in\N$, $\Xi_{\ell}=\{\xi_{\ell,1},\dots,\xi_{\ell,n_{\ell}+p+1}\}$ and $Z_{\ell}=\{\zeta_{\ell,1},\dots,\zeta_{\ell,M_{\ell}}\}$ are given for every $1\leq\ell\leq d$. We set the degree $\mathbf{p}:=(p,\dots,p)$ and $\mathbf{\Xi}:=\Xi_{1}\times\dots\times\Xi_{d}$. With no risk of ambiguity we will write $p$ in place of $\mathbf{p}$. Note that the breakpoints of $Z_{\ell}$ form a Cartesian grid in the parametric domain $\widehat{\Omega}_{0}=(0,1)^{d}$. We define the parametric Bézier mesh

[TABLE]

where each $Q_{\mathbf{j}}$ is called Bézier element, with $h_{Q_{\mathbf{j}}}:=\operatorname{diam}\left(Q_{\mathbf{j}}\right)$. Let $h:=\max\{h_{Q}:Q\in\widehat{\mathcal{M}}_{0}\}$, hence we denote $\widehat{\mathcal{M}}_{0,h}=\widehat{\mathcal{M}}_{0}$.

Throughout the manuscript we are going to rely on the following shape-regularity hypothesis, which allows us to assign $h_{Q}$ as the unique size of the element, without the necessity of dealing with the length of its edges separately. This hypothesis is implicitly used throughout the paper in any result involving the mesh size $h$. {assumption} The family of meshes $\{\widehat{\mathcal{M}}_{0,h}\}_{h}$ is assumed to be shape-regular, that is, the ratio between the smallest edge of $Q\in\widehat{\mathcal{M}}_{0,h}$ and its diameter $h_{Q}$ is uniformly bounded with respect to $Q$ and $h$.

Remark 2.1**.**

*The shape-regularity hypothesis implies that the mesh is locally-quasi uniform, i.e. the ratio of the sizes of two neighboring elements is uniformly bounded (see [3]). *

Let $\mathbf{I}:=\{\mathbf{i}=(i_{1},\dots,i_{d}):1\leq i_{\ell}\leq n_{\ell}\}$ be a set of multi-indices. For each $\mathbf{i}=(i_{1},\dots,i_{d})$, we define the set of multivariate B-splines

$\{\widehat{B}_{\mathbf{i},p}(\mathbf{\widehat{x}})=\widehat{B}_{i_{1},p}(\widehat{x}_{1})\dots\widehat{B}_{i_{d},p}(\widehat{x}_{d}):\mathbf{i}\in\mathbf{I}\}.$ Moreover, for an arbitrary Bézier element $Q_{\mathbf{j}}\in\widehat{\mathcal{M}}_{0,h}$, we define its support extension $\widetilde{Q}_{\mathbf{j}}=\widetilde{I}_{1,j_{1}}\times\dots\times\widetilde{I}_{d,j_{d}}$, where $\widetilde{I}_{l,j_{\ell}}$ is the univariate support extension of the univariate case defined above.

The multivariate spline space in $\widehat{\Omega}$ is defined as $S_{p}(\mathbf{\Xi})=\operatorname{span}\{\widehat{B}_{\mathbf{i},p}(\mathbf{\widehat{x}}):\mathbf{i}\in\mathbf{I}\},$ which can also be seen as the space of piecewise multivariate polynomials of coordinate degree $p$ and with regularity across the Bézier elements given by the knots multiplicities. Note that $S_{p}(\mathbf{\Xi})=\bigotimes_{\ell=1}^{d}S_{p}(\Xi_{\ell})$.

Remark 2.2**.**

*What has been said so far can be easily generalized to the case of Non-Uniform Rational B-Splines (NURBS) basis functions. See for instance [15]. *

Remark 2.3**.**

*The previous construction as well as what follows could be done in a more general setting considering different, but fixed, degrees in each Cartesian direction. In this case all the inequality constants appearing in the theoretical results would depend on the difference between the degrees, and would possibly explode if this difference is not kept bounded. Further generalizations to the anisotropic setting, either in terms of the mesh or in terms of the degree, would rely on the approximation theory in anisotropic Sobolev spaces (see, for instance, [5, 16]), These extensions are far from trivial, and out of the scope of this work. *

2.3 Parametrization, mesh and approximation space for trimming domains

Let $\Omega_{0}\subset\R^{d}$ be the original domain before trimming. We assume that there exists a map $\mathbf{F}\in\left(S_{p^{0}}(\mathbf{\Xi}^{0})\right)^{d}$ such that $\Omega_{0}=\mathbf{F}(\widehat{\Omega}_{0})$, for given degree $p^{0}$ and knot vector $\mathbf{\Xi}^{0}$. We define the (physical) Bézier mesh as the image of the elements in $\widehat{\mathcal{M}}_{0,h}$ through $\mathbf{F}$:

[TABLE]

We denote $h_{K}:=\operatorname{diam}\left(K\right)$ for each $K\in\mathcal{M}_{0,h}$. To prevent the existence of singularities in the parametrization we make the following assumption. {assumption} The parametrization $\mathbf{F}:\widehat{\Omega}_{0}\to\Omega_{0}$ is bi-Lipschitz. Moreover, {\left.\kern-1.2pt\mathbf{F}\vphantom{\big{|}}\right|_{\overline{Q}}}\in C^{\infty}(\overline{Q}) for every $Q\in\widehat{\mathcal{M}}_{0,h}$ and {\left.\kern-1.2pt\mathbf{F}^{-1}\vphantom{\big{|}}\right|_{\overline{K}}}\in C^{\infty}(\overline{K}) for every $K\in\mathcal{M}_{0,h}$.

Some consequences of Assumption 2.3 are the following.

$h_{Q}\approx h_{K}$, i.e. $\exists\ C_{1}>0,C_{2}>0$ such that $C_{1}h_{K}\leq h_{Q}\leq C_{2}h_{K}$; 2. 2.

$\exists\ C>0$ such that, $\forall\ Q\in\widehat{\mathcal{M}}_{0,h}$ such that $\mathbf{F}(Q)=K$, it holds

$\left\lVert D\mathbf{F}\right\rVert_{L^{\infty}(Q)}\leq C$ and $\left\lVert D\mathbf{F}^{-1}\right\rVert_{L^{\infty}(K)}\leq C$; 3. 3.

$\exists\ C_{1}>0,C_{2}>0$ such that $\quad C_{1}\leq\left\lvert\operatorname{det}(D\mathbf{F}(\widehat{\bf x}))\right\rvert\leq C_{2}$, for all $\widehat{\bf x}\in\widehat{\Omega}_{0}$.

Let $\widehat{V}_{h}=S_{p}(\mathbf{\Xi})$ be a refinement of $S_{p^{0}}(\mathbf{\Xi}^{0})$ and define

[TABLE]

where $\{\widehat{B}_{\mathbf{i},p}:\mathbf{i}\in\mathbf{I}\}$ is a basis for $\widehat{V}_{h}$. Note that throughout this document $C$ will denote generic constants that may change at each occurrence, but that are always independent of the local mesh size.

Let us clarify the interpolation strategy we are going to rely upon in the rest of this manuscript. Given a function $u\in H^{s}(\Omega)$, $s\geq 1$, we extend it using the Sobolev-Stein extension operator (see, for instance, Section 3.2 of [36]) $E:H^{s}(\Omega)\to H^{s}(\R^{d})$ and denote its restriction to the un-trimmed domain as \widetilde{u}:={\left.\kern-1.2ptE(u)\vphantom{\big{|}}\right|_{\Omega_{0}}}, for $u\in H^{s}(\Omega)$. The spline quasi-interpolant operator ([10]) associated to the uncut mesh $\mathcal{M}_{0,h}$ is $\Pi_{0}:H^{s}(\Omega_{0})\to V_{h}$. Hence, we are allowed to write $\left\lVert\Pi_{0}(u)\right\rVert_{H^{s}(\Omega)}\leq\left\lVert\Pi_{0}(\widetilde{u})\right\rVert_{H^{s}(\Omega_{0})}\leq C\left\lVert\widetilde{u}\right\rVert_{H^{s}(\Omega_{0})}\leq C\left\lVert u\right\rVert_{H^{s}(\Omega)}$ and, similarly, for every $0\leq t\leq s$, $\left\lVert u-\Pi_{0}(u)\right\rVert_{H^{t}(\Omega)}\leq\left\lVert\widetilde{u}-\Pi_{0}(\widetilde{u})\right\rVert_{H^{t}(\Omega_{0})}\leq Ch^{s-t}\left\lVert\widetilde{u}\right\rVert_{H^{s}(\Omega_{0})}\leq Ch^{s-t}\left\lVert u\right\rVert_{H^{s}(\Omega)}$.

3 The isogeometric formulation

At this point we suppose to trim $\Omega_{0}$ as explained in (1), for simplicity, with $N=1$ , i.e. the new domain is $\Omega=\Omega_{0}\setminus\overline{\Omega}_{1}$. We denote the trimming curve as $\Gamma_{trim}=\partial\Omega\cap\partial\Omega_{1}$. Let us consider the Poisson equation as model problem. Given $f\in L^{2}(\Omega)$, $g_{D}\in H^{\frac{1}{2}}(\Gamma_{D})$ and $g_{N}\in H^{-\frac{1}{2}}(\Gamma_{N})$, find $u:\Omega\to\R$ such that

[TABLE]

where $\Gamma_{D}\cup\Gamma_{N}=\Gamma=:\partial\Omega$ and $\mathring{\Gamma}_{D}\cap\mathring{\Gamma}_{N}=\emptyset$, and $\displaystyle{\frac{\partial u}{\partial n}}:=\nabla u\cdot\mathbf{n}$ is the normal derivative, with $\mathbf{n}$ the outward unit normal to $\Gamma$. Observe that, in general, $\Gamma_{trim}\cap\Gamma_{D}\neq\emptyset$.

Let us now develop and extend the notation introduced in Section 2, to adapt it to trimmed domains.

The approximation space on the trimmed domain $\Omega$ is \widetilde{V}_{h}:=\operatorname{span}\{{\left.\kern-1.2ptB_{\mathbf{i},\mathbf{p}}\vphantom{\big{|}}\right|_{\Omega}}:\mathbf{i}\in\mathbf{I}\}. The new parametric Bézier mesh is $\widehat{\mathcal{M}}_{h}=\{Q\in{\widehat{\mathcal{M}}}_{0,h}:Q\cap\widehat{\Omega}\neq\emptyset\},$ where $\widehat{\Omega}=\F^{-1}(\Omega)$. The physical mesh is $\mathcal{M}_{h}=\{\F(Q):Q\in\widehat{\mathcal{M}}_{h}\}.$ and the set of Bézier elements cut by the trimming curve is denoted as $\mathcal{G}_{h}=\{K\in\mathcal{M}_{h}:\overline{K}\cap\Gamma_{trim}\neq\emptyset\}.$

For every $K\in\mathcal{M}_{h}$, let $h_{K}:=\operatorname{diam}(K)$, $h_{\max}:=\max_{K\in\mathcal{M}_{h}}h_{K}$ and $h_{\min}:=\min_{K\in\mathcal{M}_{h}}h_{K}$. We define $\mathsf{h}:\Omega\to\R^{+}$ to be the piecewise constant mesh-size function of $\mathcal{M}_{h}$ given by {\left.\kern-1.2pt\mathsf{h}\vphantom{\big{|}}\right|_{K}}:=h_{K}.

First of all, we make an assumption on how the mesh is cut by the boundary. {assumption} There exists $C>0$ such that, $\forall\ h>0,\ \forall\ K\in\mathcal{M}_{h}$, it holds $\left\lvert\Gamma_{K}\right\rvert\leq Ch_{K}^{d-1}$, where $\Gamma_{K}:=\Gamma_{D}\cap\overline{K}\neq\emptyset$.

We denote as $\widehat{\Gamma}_{D}:=\F^{-1}\left(\Gamma_{D}\right)$ and as $\widehat{\mathbf{n}}$ its outward unit normal.

Since the point is to avoid a reparametrization and a remeshing of the trimmed domain, it is natural to see the analogy with fictitious domain methods, where the physical domain, with a possibly complicated topology, is immersed into a simpler, but unfitted, background mesh. Similarly to fictitious domain methods, we need to be able to impose essential boundary conditions when the mesh is not fitted with boundary. Following [21, 44], we decide to employ Nitsche’s method, which in its symmetric form reads as follows.

Find $u_{h}\in\widetilde{V}_{h}$ such that

[TABLE]

where $\beta>0$ is a penalization parameter.

We define

[TABLE]

Our main goal is to provide a minimal stabilization to make formulation (3) uniformly well-posed with respect to the mesh-size.

4 Stability

Firstly, we need to clarify what we actually mean by “stability” of the discrete variational problem (3). We introduce the following mesh-dependent scalar product

[TABLE]

which induces the discrete norm

[TABLE]

Definition 4.1**.**

Problem (3) is stable if there exist $\overline{\beta}>0$ and $\alpha>0$ such that for every $\beta\geq\overline{\beta}$, for every $h>0$ it holds that

[TABLE]

and for every fixed $\beta\geq\overline{\beta}$ there exists $\gamma>0$ such that for every $h>0$ it holds that

[TABLE]

Remark 4.2**.**

*The main point of Definition 4.1 is to find $\overline{\beta}$, $\alpha$ and $\gamma$ that do not depend on the trimming configuration. In the following lines we will show with a numerical example that the formulation (3) is not stable. *

Remark 4.3**.**

*In Definition 4.1 we have followed [44]. Note that the coercivity constant only depends on $\overline{\beta}$, while the continuity constant depends on the penalization parameter $\beta$, and in particular it grows with $\beta$. This dependence of the constant on $\beta$ also occurs in Theorem 7.1 and in Proposition 7.6. In practice, $\beta$ has to be chosen large enough (i.e. larger than $\overline{\beta}$), but as close as possible to $\overline{\beta}$, to avoid that the continuity constant deteriorates. *

Proposition 4.4**.**

*If for all $K\in\mathcal{G}_{h}$ we have $\Gamma_{K}=\emptyset$, then problem (3) is stable. *

Proof 4.5**.**

*We refer the reader to [44]. *

The following numerical experience shows that there exists a trimming configuration for which the formulation (3) is not stable according to Definition 4.1. In particular we show that for every fixed $\beta$ the continuity constant $\gamma$ may be arbitrarily large, for a given $h>0$. First, we notice that if $\gamma_{\beta}$ is the continuity constant corresponding to $\beta$, then $\gamma_{\beta}>\gamma_{1}$ for every $\beta>1$. So, we fix $\beta=1$ and show that $\gamma_{1}$ can be arbitrarily large.

Let us consider the following eigenvalue problem. Find $u_{h}\in\widetilde{V}_{h}\setminus\{0\}$ and $\lambda_{h}\in\mathbb{R}$ such that

[TABLE]

As the problem is symmetric, the continuity constant $\gamma_{1}$ equals the maximum eigenvalue of (5). Let us consider $\Omega_{0}=(0,1)^{2}$ and as trimmed domain $\Omega=(0,1)\times(0,0.757)$. We fix $h=2^{-5}$ as mesh size and $p=3$ as degree. We construct a sequence of discrete spaces $\left(\widetilde{V}_{h,\eps}\right)_{\eps}$ of degree $p$ and of class $C^{2}$ at the internal knots, starting from the uniform knot vectors $\Xi_{x}$, $\Xi_{y}$ and substituting in the latter the knot $0.75$ with $\overline{\xi}=0.757-\eps$, see Figure LABEL:domain_eps. Basically, the horizontal knot line $\{(x,y):y=0.75\}$ is replaced by $\{(x,y):y=\overline{\xi}\}$, which is such that the smaller $\eps>0$ is, the closer to the trimming curve it becomes.

In Figure LABEL:fig_eigs_bad we can see the dependence of the spectrum of (5) on the magnitude of $\varepsilon$. In particular, the magnitude of the largest generalized eigenvalue goes to infinity as $\varepsilon$ goes to [math], implying that the discrete formulation (3) is not stable, as the continuity constant can be made arbitrarily large by reducing $\eps$.

Going through the proof of stability, we clearly miss a discrete trace inequality which is uniform with respect to any mesh-trimming curve configuration, namely:

[TABLE]

where $C$ does not depend neither on $K$ nor on $K\cap\Omega$.

At this point, in order to be able to deal with a well-posed, hence stable, problem we want to find a way to improve this discrete trace inequality, where the constant does not depend on how the trimmed boundary intersects the mesh.

5 A new stabilization technique

The goal of this section is to present a new stabilization technique for the problem (3). Our construction is inspired by the work of J. Haslinger and Y. Renard in [24].

Let us partition the elements of the Bézier mesh into two disjoint sub-families.

Definition 5.1**.**

Let $\theta\in(0,1]$ and $Q\in\widehat{\mathcal{M}}_{h}$. We say that $Q$ is a good element if

[TABLE]

Otherwise, $Q$ is a bad element. Thanks to the regularity Assumption 2.3 on $\F$, this classification on the parametric elements induces naturally a classification on the physical elements. $\mathcal{M}_{h}^{g}$ stands for the collection of the good physical Bézier elements and $\mathcal{M}_{h}^{b}$ for the one of the bad physical elements. Note that $\mathcal{M}_{h}\setminus\mathcal{G}_{h}\subseteq\mathcal{M}_{h}^{g}$ and $\mathcal{M}_{h}^{b}\subseteq\mathcal{G}_{h}$. We denote the set of neighbors of $K$ as:

[TABLE]

*where $C$ does not depend on the mesh size. *

The following assumption is not restrictive since it holds true if the mesh is sufficiently refined. {assumption} We assume that for any $K\in\mathcal{M}_{h}^{b}$, $\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}\neq\emptyset$.

In what follows, we will use Assumption 5 to construct a stable representation of the normal flux of discrete functions. Let us assume that there exists an operator

[TABLE]

which approximates the normal derivative on $\Gamma_{D}$ in a sense that will be specified. We propose the following stabilized formulation of problem (3).

Find $u_{h}\in\widetilde{V}_{h}$ such that

[TABLE]

where

[TABLE]

Theorem 5.2**.**

Suppose the following stability property is satisfied: there exists a uniform $C>0$ such that for every $K\in\mathcal{G}_{h}$

[TABLE]

*where $K^{\prime}=K$ if $K\in\mathcal{M}_{h}^{g}$, otherwise $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$. Then problem (7) is stable in the sense of Definition 4.1 (modified accordingly). *

Proof 5.3**.**

For the continuity, let $u_{h},v_{h}\in\widetilde{V}_{h}$ and estimate

[TABLE]

where we employed first Cauchy-Schwarz inequality, the definition of the norm (4) and the stability property (8). Take $u_{h}\in\widetilde{V}_{h}$. Using Young inequality, with $\delta>0$, and, again, the stability property (8), we obtain:

[TABLE]

*from which we deduce the coercivity, provided $C<\delta<\beta$. *

Remark 5.4**.**

*In order for the solution of (7) to be a good approximation of $u$, it is clear that we will also need to quantify the error between $R_{h}(u_{h})$ and $\frac{\partial u}{\partial n}$. This fact will be addressed in the next section. *

6 Construction of the stabilization operator

The definition of the operator $R_{h}$ is not unique. As already observed, we seek for a stable approximation of the normal derivative on the trimmed part of the boundary, namely on $\Gamma_{K}$ for every $K\in\mathcal{G}_{h}$. Here, we propose two different constructions of such an operator.

• •

a stabilization in the parametric domain: for each $K\in\mathcal{M}_{h}^{b}$ we take the (unique) polynomial extension of the pull-back of the functions of $\widetilde{V}_{h}$ from $Q^{\prime}=\F^{-1}(K^{\prime})$ to $Q=\F^{-1}(K)$, where $K^{\prime}$ is a good neighbor;

• •

a stabilization in the physical domain: for each $K\in\mathcal{M}_{h}^{b}$, we first $L^{2}$-project the spline functions restricted to the good neighbor $K^{\prime}$ onto the polynomial space $\mathbb{Q}_{p}(K^{\prime})$, then we take their (unique) polynomial extension up to $K$.

Definition 6.1** (Stabilization in the parametric domain).**

We define the operator $R_{h}$ locally as {\left.\kern-1.2ptR_{h}(v_{h})\vphantom{\big{|}}\right|_{K}}:=R_{K}(v_{h}) $\;\forall\ K\in\mathcal{G}_{h}$, $\forall\ v_{h}\in\widetilde{V}_{h}$, where

• •

if $K\in\mathcal{M}_{h}^{g}$,

[TABLE]

• •

if $K\in\mathcal{M}_{h}^{b}$, $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$,

[TABLE]

where $\mathcal{E}:\mathbb{Q}_{p}(Q^{\prime})\to\mathbb{Q}_{p}(Q^{\prime}\cup Q)$ is the polynomial natural extension.

Definition 6.2** (Stabilization in the physical domain).**

An alternative stabilization operator can be defined by using the $L^{2}$-projection in the physical domain. We define the operator $R_{h}$ locally as {\left.\kern-1.2ptR_{h}(v_{h})\vphantom{\big{|}}\right|_{K}}:=R_{K}(v_{h}) $\forall\ K\in\mathcal{G}_{h}$, $\forall\ v_{h}\in\widetilde{V}_{h}$:

• •

if $K\in\mathcal{M}_{h}^{g}$,

[TABLE]

• •

if $K=\mathbf{F}(Q)\in\mathcal{M}_{h}^{b}$, $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$,

[TABLE]

where $P:L^{2}\left(K^{\prime}\right)\to\mathbb{Q}_{p}\left(K^{\prime}\right)$ is the $L^{2}$-orthogonal projection and

$\mathcal{E}:\mathbb{Q}_{p}(K^{\prime})\to\mathbb{Q}_{p}(K^{\prime}\cup K)$* is the polynomial natural extension.*

Remark 6.3**.**

*Note that in the trivial case where $\F=\operatorname{\bf{Id}}$, the $L^{2}$-projection $P$, restricted to $\widetilde{V}_{h}$, reduces to the identity operator and the two stabilizations coincide. *

6.1 Properties of the stabilization in the parametric domain

We are now up to verify if our choice of $R_{h}$ verifies the stability property (8). Its proof relies on a series of quite technical results that are reported in the Appendix.

Theorem 6.4**.**

The stability property (8) holds for $R_{h}$ defined as in Definition 6.1, i.e., there exists $C>0$ such that for every $K\in\mathcal{G}_{h}$

[TABLE]

*where $K^{\prime}=K$ if $K\in\mathcal{M}_{h}^{g}$, otherwise $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$. *

Proof 6.5**.**

Fixed $K\in\mathcal{G}_{h}$, it is enough to prove

[TABLE]

for $v_{h}\in\widetilde{V}_{h}$ such that {\left.\kern-1.2pt\overline{v}_{h}\vphantom{\big{|}}\right|_{K}}:=\mathcal{E}\left({\left.\kern-1.2pt\widehat{v}_{h}\vphantom{\big{|}}\right|_{Q^{\prime}}}\right)\circ\mathbf{F}^{-1}, where $\mathcal{E}:\mathbb{Q}_{p}(Q^{\prime})\to\mathbb{Q}_{p}(Q^{\prime}\cup Q)$ and $K=\mathbf{F}(Q)$, $K^{\prime}=\mathbf{F}(Q^{\prime})\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$. We can restrict ourselves to the case $K\in\mathcal{M}_{h}^{b}$ with good neighbor $K^{\prime}$. It holds:

[TABLE]

where we have used $\mathbf{F}^{-1}(\Gamma_{K})=\mathbf{F}^{-1}(\Gamma_{D})\cap\mathbf{F}^{-1}(K)=\widehat{\Gamma}_{D}\cap Q$, because $\mathbf{F}$ preserves boundaries (as homeomorphisms do). We then have:

[TABLE]

In the first inequality we have used Hölder inequality. Now, we employ Lemma A.5 and Assumption 3:

[TABLE]

At this point, notice that we can use Lemma A.7 because $\frac{\left\lvert\Omega\cap K^{\prime}\right\rvert}{\left\lvert K^{\prime}\right\rvert}\geq\theta$ implies $\frac{\left\lvert\widehat{\Omega}\cap Q^{\prime}\right\rvert}{\left\lvert Q^{\prime}\right\rvert}\geq C\theta_{\min}$, where $C$ depends just on $\mathbf{F}$, thanks to Assumption 2.3.

Let us continue with the inequalities:

[TABLE]

*Gathering together (9) and (10), we conclude the proof. *

In what follows, we analyse the approximation properties of the operator $R_{h}$, and provide estimates that will be used in Section 7 to deduce a complete error estimate.

Proposition 6.6**.**

Let $\frac{1}{2}<k\leq p$. There exists $C>0$ such that for every $K\in\mathcal{G}_{h}$

• •

if $\frac{1}{2}<k<p-\frac{1}{2}$, for every $v\in H^{k+1}(\Omega)$,:

[TABLE]

where $K^{\prime}=K$ if $K\in\mathcal{M}_{h}^{g}$, otherwise $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$;

• •

if $p-\frac{1}{2}\leq k\leq p$ and each internal knot line is not repeated, for every $v\in H^{k+1}(\Omega)$, for all $\eps>0$,

[TABLE]

where $K^{\prime}=K$ if $K\in\mathcal{M}_{h}^{g}$, otherwise $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$.

Proof 6.7**.**

First of all, let $v\in H^{k+1}(\Omega)$, with $\frac{1}{2}<k\leq p$. We take $K\in\mathcal{G}_{h}$. Let us distinguish two cases: either $K\in\mathcal{M}_{h}^{g}$ or $K\in\mathcal{M}_{h}^{b}$.

If $K\in\mathcal{M}^{g}$. We use Lemma A.1 and standard approximation results:

[TABLE]

If $K=\mathbf{F}(Q)\in\mathcal{M}_{h}^{b}$ and $K^{\prime}=\mathbf{F}(Q^{\prime})\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$ be its good neighbor. We easily obtain:

[TABLE]

The second term converges as expected because of the properties of spline quasi-interpolants [10]. We focus on the first one. Let $\widehat{q}=q\circ\F\in\mathbb{Q}_{p}(\mathbb{R}^{d})$ be a global polynomial. Note that, trivially, \mathcal{E}({\left.\kern-1.2pt\widehat{q}\vphantom{\big{|}}\right|_{Q^{\prime}}})={\left.\kern-1.2pt\widehat{q}\vphantom{\big{|}}\right|_{Q^{\prime}}}. By triangular inequality:

[TABLE]

Using Corollary A.3, we can bound the last term of (12) as follows:

[TABLE]

The first term of (12) can be bounded using the stability property of $R_{h}$, given in Theorem 6.4:

[TABLE]

Thus, combining (12), (13) and (14), we obtain:

[TABLE]

Again, the first term converges as expected by standard approximation results. Concerning the other term, there are some issues, related to the regularity of the parametrization. By the theory of bent Sobolev spaces (see [3]), we have $\widetilde{v}\in H^{k+1}(\Omega_{0})$, but, in general, {\left.\kern-1.2pt\widetilde{v}\circ\mathbf{F}\vphantom{\big{|}}\right|_{Q\cup Q^{\prime}}}\notin H^{k+1}(Q\cup Q^{\prime}), since it is bent by $\F$, a spline of degree $p$ and regularity $p-1$ (under the assumption that internal knot lines are not repeated). It holds, indeed, that {\left.\kern-1.2pt\widetilde{v}\circ\mathbf{F}\vphantom{\big{|}}\right|_{Q\cup Q^{\prime}}}\in H^{r+1}(Q\cup Q^{\prime}), where $r+1:=\min\{k+1,p+\frac{1}{2}-\eps\}$, hence $0\leq r\leq k$ and $0\leq r\leq p-\frac{1}{2}-\eps$. So, the following inequality follows:

[TABLE]

where $0\leq r\leq k$ and $0\leq r\leq p-\frac{1}{2}-\eps$, for any $\eps>0$. Hence, pushing forward to the physical domain:

[TABLE]

Hence, from (15) and (16), we deduce:

[TABLE]

We want to rewrite inequality (LABEL:ineq_incasinata) by distinguishing two cases.

• •

$\frac{1}{2}<k<p-\frac{1}{2}$. In this case,

[TABLE]

for any $0\leq r\leq k$. Hence,

[TABLE]

• •

If $p-\frac{1}{2}\leq k\leq p$, then

[TABLE]

for any $\eps>0$.

Remark 6.8**.**

*Note that if $\frac{1}{2}<k<p-\frac{1}{2}$, the estimate is optimal. In the case $p-\frac{1}{2}\leq k\leq p$ the estimate is sub-optimal, instead. As already mentioned during the proof, this is due to the fact that $u\in H^{k+1}(K\cup K^{\prime})$ does not imply $u\circ\mathbf{F}\in H^{k+1}(Q\cup Q^{\prime})$: if the knot line between $K$ and $K^{\prime}$ is not repeated, namely $\F\in C^{p-1}(Q\cup Q^{\prime})$, then it holds $u\circ\mathbf{F}\in H^{r+1}(Q\cup Q^{\prime})$ with $r+1:=\min\{k+1,p+\frac{1}{2}-\eps\}$. Moreover, if the parametrization is less regular than requested in the hypotheses of Proposition 6.6, then the sub-optimality may be even worse. More precisely, if $\F\in C^{s}\left(Q\cup Q^{\prime}\right)$, which is the case if the knot line is repeated $p-s$ times, then we have $r+1:=\min\{k+1,s+\frac{3}{2}-\eps\}$. We will see an example of this sub-optimal behaviour in the worst case scenario of $s=0$ in Section 8.3. *

Remark 6.9**.**

*Any method based on polynomial extrapolation of the B-splines in the parametric domain may also suffer of this sub-optimality depending on the regularity of the isogeometric map $\mathbf{F}$, because the theory of bent Sobolev spaces from [3] cannot be applied. In particular the method of extended B-splines which works very well in the parametric domain [25], may suffer a lack of accuracy in the isogeometric setting [30, 32]. *

6.2 Properties of the stabilization in the physical domain

Theorem 6.10**.**

The stability property (8) holds for $R_{h}$ defined as in Definition 6.2, i.e., there exists $C>0$ such that for every $K\in\mathcal{G}_{h}$

[TABLE]

*where $K^{\prime}=K$ if $K\in\mathcal{M}^{g}$, otherwise $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}^{g}$. *

Proof 6.11**.**

Let us start applying Hölder inequality and Lemma A.5:

[TABLE]

We finish with Lemma A.7, Assumption 3 and the stability of the $L^{2}$-orthogonal projection $P$, see for instance [7],

[TABLE]

*We conclude by summing over $K\in\mathcal{G}_{h}$. *

Proposition 6.12**.**

Let $\frac{1}{2}<k\leq p$. There exists $C>0$ such that for every $K\in\mathcal{G}_{h}$

[TABLE]

*where $K^{\prime}=K$ if $K\in\mathcal{M}_{h}^{g}$, otherwise $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$. *

Proof 6.13**.**

We can focus on the case $K\in\mathcal{M}_{h}^{b}$. Let $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$ and $q\in\mathbb{Q}_{p}\left(K^{\prime}\right)$.

[TABLE]

Let us focus on the first term. After having observed that $P(q)=q$, we apply the stability property proved in Theorem 6.10 and, again, triangular inequality

[TABLE]

We choose $q=P\left(v\right)$. Note that the second term of (21) converges as expected by the approximation properties of the $L^{2}$-projection. Plugging (22) into (21):

[TABLE]

*We conclude by summing over $K\in\mathcal{G}_{h}$. *

7 A priori error estimate

The preparatory results of Propositions 6.6 and 6.12 were needed in order to prove the following convergence theorem.

Theorem 7.1**.**

Let $\frac{1}{2}<k\leq p$. There exists $\overline{\beta}>0$ such that, for every $\beta\geq\overline{\beta}$, if $u\in H^{k+1}(\Omega)$ is the solution to (2) and $u_{h}\in\widetilde{V}_{h}$ solution to (7), then

[TABLE]

where $S_{h}$ is the strip of width $ch$, $c\geq 1$, such that $S_{h}\supseteq\bigcup_{K\in\mathcal{M}_{h}^{b}}\left(\widetilde{K}\cup\widetilde{K}^{\prime}\right)$, and $K^{\prime}\in\mathcal{M}_{h}^{g}\cap\mathcal{N}(K)$. Moreover, (24) holds for every $r$ such that:

• •

$0\leq r<p-\frac{1}{2}$* with the stabilization in the parametric domain of Definition 6.1;*

• •

$0\leq r\leq p$* with the stabilization in the physical domain of Definition 6.2.*

Proof 7.2**.**

From Theorems 5.2, 6.4 and 6.10 we know that $\overline{a}_{h}(\cdot,\cdot)$ is coercive w.r.t. $\left\lVert\cdot\right\rVert_{1,h,\Omega}$, i.e. there exists $\alpha>0$ such that for every $u_{h}\in\widetilde{V}_{h}$

[TABLE]

Let $v_{h}\in\widetilde{V}_{h}$. Using the triangular inequality and coercivity, we get:

[TABLE]

Then, recalling that $u_{h}$ solves (7), we get

[TABLE]

Since $u$ solves (2): $\int_{\Omega}fw_{h}=\int_{\Omega}\nabla u\cdot\nabla w_{h}-\int_{\Gamma_{D}}\frac{\partial u}{\partial n}w_{h}$ and {\left.\kern-1.2ptu\vphantom{\big{|}}\right|_{\Gamma_{D}}}=g_{D}, hence:

[TABLE]

Let us now estimate the four terms separately. We will leave $\RomanNumeralCaps 2$ for last since its analysis depends on the choice of the stabilization. Clearly

[TABLE]

where $C>0$ linearly depends on $\beta$. Note that this will not compromise the uniformity of the resulting constant, provided that $\beta$ is chosen as close as possible to $\overline{\beta}$ (see the discussion in Remark 4.3). Using the stability property (8) and taking $K^{\prime}\in\mathcal{N}(K)\cap\mathcal{M}_{h}^{g}$ (if $K$ itself is a good element, then take $K^{\prime}=K$), we get:

[TABLE]

Let us estimate the term $\RomanNumeralCaps 2$. By definition of the norm $\left\lVert\cdot\right\rVert_{1,h,\Omega}$:

[TABLE]

Now, we choose $v_{h}=\Pi_{0}(\widetilde{u})$ and distinguish two cases.

• •

If we use the stabilization in the parametric domain of Definition 6.1, hence apply Proposition 6.6, we get, for any $0\leq r<p-\frac{1}{2}$,

[TABLE]

• •

Employing the stabilization in the physical domain of Definition 6.2, hence apply Proposition 6.12, we obtain

[TABLE]

Therefore, we have that

[TABLE]

where $S_{h}$ is the strip of width $ch$, $c\geq 1$, such that $S_{h}\supseteq\bigcup_{K\in\mathcal{M}_{h}^{b}}\left(\widetilde{K}\cup\widetilde{K}^{\prime}\right)$, and $K^{\prime}\in\mathcal{M}_{h}^{g}\cap\mathcal{N}(K)$ and we can choose any $r$ such that:

• •

$0\leq r<p-\frac{1}{2}$* if we use the stabilization in the parametric domain, hence apply Proposition 6.6;*

• •

$0\leq r\leq p$* if we use the one in the physical domain and use Proposition 6.12.*

As a consequence, we have that

[TABLE]

where in (27), (28) we choose again $v_{h}=\Pi_{0}(\widetilde{u})$.

We now combine the last inequality (29) with (25) and (26) to obtain

[TABLE]

Using approximation results of quasi-interpolants in spline spaces [10], we conclude

[TABLE]

*where $r$ is the same as above. *

Remark 7.3**.**

*As already observed in Remark 6.8, when $u\in H^{k+1}(\Omega)$ with $\frac{1}{2}<k<p-\frac{1}{2}$, both stabilizations give rise to optimal a priori error estimates. When $u\in H^{k+1}(\Omega)$ with $p-\frac{1}{2}\leq k\leq p$ and $k>\frac{1}{2}$, instead, stabilization in Definition 6.1 is sub-optimal. In this case the estimate can be modified and improved using the following result. *

Lemma 7.4**.**

Let $\eps>0$ and $S_{h}$ be defined as in Theorem 7.1. Then, there exists $C>0$ such that

[TABLE]

Proof 7.5**.**

Using the fact that $r<p$, we are able to recover an integer order for the Sobolev norm and so to apply Lemma A.11 with $s=\frac{1}{2}-\eps$:

[TABLE]

*In the last inequality we used the boundedness of the Sobolev-Stein extension operator. *

Proposition 7.6**.**

Let $u\in H^{p+1}(\Omega)$ be the solution to (2) and $u_{h}\in\widetilde{V}_{h}$ solution to (7), obtained using the stabilization in the parametric domain of Definition 6.1. Then, the following error estimate holds:

[TABLE]

Proof 7.7**.**

*It immediately follows combining Theorem 7.1 and Lemma 7.4. *

Remark 7.8**.**

*At the prize of slightly higher regularity request, optimal convergence rate is to be expected also for stabilization in Definition 6.1. *

8 Numerical examples

8.1 Some details about the implementation

For accurate numerical integration, we decompose the trimmed elements into smaller quadrilateral tiles where we compute the integrals. These tiles are reparametrized as Bézier surfaces of the same degree $p$ as the approximation space used to discretize our PDE, see [2] for a detailed explanation. We remark that this reparametrization is also used to compute the boundary integrals.

In order to compute the stabilization terms appearing in (7), first of all for each bad trimmed element $K$ we choose $K^{\prime}$: among all the neighbours of $K$, we choose (the) one with the largest relative overlap $\left\lvert K^{\prime}\cap\Omega\right\rvert/\left\lvert K^{\prime}\right\rvert$. Then we need to locally project functions living in $K^{\prime}$ (or in $Q^{\prime}$) onto the space of polynomials on $K^{\prime}$ (or $Q^{\prime}$) and extend them up to $\Gamma_{K}$. For the stabilization in the parametric domain, by taking as a basis the Bernstein polynomials on $Q^{\prime}$ the projection can be computed by knot insertion, while for the stabilization in the physical domain the $L^{2}$-projection is needed anyhow.

8.2 Validation of stability

Let us repeat the numerical experiment of Section 4 in order to validate the effectiveness of our stabilization technique. Let us solve the eigenvalue problem (5) with the stabilization in Definition 6.1 (since $\mathbf{F}=\mathbf{Id}$, the two proposed stabilizations techniques are equivalent) in the trimmed domain of Figure LABEL:fig:mesh for the same values of $\eps$ used in Section 4. The result is shown in Figure LABEL:fig_eigs_good.

This time we observe that the spectrum remains bounded independently of $\eps$, confirming our method to be stable.

8.3 Validation of the a priori error estimate

In the following we focus on the Poisson problem (7) with the difference that, while we impose Dirichlet boundary conditions weakly on the trimmed parts of the boundary, on the other parts where the mesh is fitted with the boundary we impose them in the strong sense.

Test 1. Let $\Omega=\Omega_{0}\setminus\overline{\Omega}_{1}$ be defined as in Figure 3a, where $\Omega_{0}={\bf F}((0,1)^{2})$ is a quarter of annulus ($\F$ is non linear) constructed with biquadratic NURBS, and $\Omega_{1}$ is the image of a ball in the parametric domain through the isogeometric map, namely $\Omega_{1}=\F(B(0,r))$, with $r=0.76$. We consider as manufactured solution $u_{ex}(x,y)=e^{x}\sin(xy).$ We solve the Poisson problem using the stabilized formulation (7), the stabilization in the parametric domain and the parameters $\beta=1$ and $\theta=0.1$. The results of convergence for different values of $p$, that are displayed Figure 3b, show that we obtain the optimal order of convergence.

Test 2. We now consider the Poisson problem in the L-shaped domain shown in Figure 4a, given by $\Omega=\Omega_{0}\setminus\overline{\Omega}_{1}$, where $\Omega_{0}=(-2,1)\times(-1,2)$ and $\Omega_{1}=(0,1)\times(-1,0)$. The exact solution is chosen as the singular function that, in polar coordinates, reads as $u(r,\varphi)=r^{\frac{2}{3}}\sin\left(\frac{2}{3}\varphi\right)\in H^{\frac{5}{3}-\delta}(\Omega)$, for every $\delta>0$. The function has a singularity at the re-entrant corner in the origin, and the domain is chosen in such a way that the corner is always located in the interior of an element. We employ the formulation (7) together with the stabilization operator in Definition 6.1, noting that since the parametrization is a simple scaling, both stabilizations are equivalent. This time we set the parameters $\theta=1$ and, due to the presence of the singularity, $\beta=(p+1)\cdot 10$. The numerical results of Figure 4b agree with the theory as the method converges with order $\frac{2}{3}$, and the sub-optimal behaviour is due to the low regularity of the reference solution.

Test 3. The goal of this test is to show that, when the regularity of the mapping $\F$ is low between trimmed elements and their neighbors, the stabilization in the physical domain is more effective than the ones based on polynomial extensions in the parametric domain (as it is the case for our stabilization in the parametric domain, but also for the method proposed in [30]). Let us consider again as the domain $\Omega_{0}$ the quarter of annulus, this time parametrized with a different map ${\bf F}$: starting from the standard biquadratic NURBS parametrization, we perform knot insertion adding the knot $\xi=0.75$, with multiplicity $2$, in the direction corresponding to the angular coordinate, that corresponds to the thick black line in Figure 5a. In order to get a geometry of class $C^{0}$, we set the second coordinate of one control point, highlighted in Figure 5b, equal to $0.5$ in homogeneous coordinates.

Note that the new parametrization is only of class $C^{0}$ in correspondence of the knot line given by $\F(\{(x,y):x\in(0,1),\,y=0.75\})$. To ensure that this knot line is located between $K$ and $K^{\prime}$, we define the trimmed domain as $\Omega=\F\left((0,1)\times(0,0.75+\eps)\right)$, with $\eps=10^{-8}$. Here we set $\theta=1$ and, because of the lower regularity of the parametrization, $\beta=(p+1)\cdot 25$. We know from Remarks 6.8 and 6.9 that the convergence rate deriving from the stabilization in Definition 6.1 (and any stabilization based on polynomial extensions in the parametric domain) may suffer of sub-optimality. In particular, from Figure 6a, we see that the error with the stabilization in the parametric domain is converging just as $h^{\frac{1}{2}}$ for any degree $p$, while in Figure 6b we observe that the desired convergence rates are reached when using the stabilization in the physical domain.

8.4 Conditioning

Even if an exhaustive discussion about the conditioning of the stiffness matrix in trimmed geometries is beyond the scope of this work (for a more detailed discussion on the topic see, for instance, [18, 19]), we would like to present some numerical experiments for the sake of completeness. We focus again on the formulation (7) of the Poisson problem. Again, we impose Dirichlet boundary conditions weakly on the trimmed parts of the boundary, and strongly on the fitted parts.

Test 1. Let us come back to the quarter of annulus with a hole and, as above, we employ B-splines of degree $p=3$. In Figure 7a we show that our stabilization coupled with a simple diagonal scaling, which can be interpreted as a left-right Jacobi preconditioner, is able to solve the conditioning issue. In Figure 7b we compare the effectiveness of the diagonal rescaling with and without the stabilization, and we observe that the effect of the stabilization is marginal with respect to the one of the diagonal preconditioner. The stabilization used is the one in the parametric domain with $\beta=1$ and $\theta=0.1$.

Test 2. Let us consider the same configuration as in the test of Figure LABEL:domain_eps, for which we notice again that the two stabilizations are equivalent. Let us take B-splines of degree $p=3$, as mesh size $h=2^{-5}$, and set the penalization parameter $\beta=1$. After a simple diagonal rescaling as preconditioner, we compare the condition number of the stiffness matrix, as a function of $\eps$, obtained for the non-stabilized ($\theta=0$) and the stabilized ($\theta=1$) formulations. Note that as the ratio in Definition 5.1 is the same for all cut elements, it is sufficient to consider only these two values of $\theta$. The results in Figure 8a show the diagonal rescaling is acting as a robust preconditioner with respect to the size of the trimming. Then, we perform uniform dyadic refinement and we plot the condition number as a function of the mesh-size $h$, obtaining the plots in Figures 8b and 8c. The results suggest a better behaviour of the condition number when a stabilized formulation is employed to solve the problem.

Test 3. This test is inspired by [18]. Let us embed $\Omega=\left(0.19,0.78\right)\times\left(0.22,0.78\right)$ in the untrimmed domain $\Omega_{0}=(0,1)^{2}$ with un underlying mesh of size $h=2^{-3}$. We consider B-splines of degree $p=2$. Now, let us rotate $\Omega$ around its barycenter for different angles $\alpha$ (see Figure 9a). For each $\alpha=i\frac{\pi}{200}$, $i=0,\dots,100$ we face a specific trimming configuration where there may appear B-splines whose support intersects in a “pathological way” the domain $\Omega$. Let us denote the “smallest volume fraction” $\eta:=\min_{K\in\mathcal{G}_{h}}\left\lvert\Omega\cap K\right\rvert$. In Figure 9b we plot the condition number of the stiffness matrix against the smallest volume fraction, in order to compare the non stabilized case with the stabilized (with parameter $\theta=0.5$) and diagonally rescaled one. Let us observe that even if the behaviour of the condition number appears to be much better after stabilization and diagonal rescaling, it is still strongly affected from the way the mesh is cut by the trimming boundary. In this regard this is a counter-example to the fact that diagonal rescaling, together with our stabilization, is a robust preconditioner with respect to the trimming operation. This will be object of further investigations by the authors in the future.

Appendix A Auxiliary theoretical results

Lemma A.1**.**

There exists $C>0$ depending on $\Gamma_{D}$, but independent of the mesh-boundary intersection, such that for every $K\in\mathcal{M}_{0,h}$

[TABLE]

Proof A.2**.**

*It follows straightforward from Lemma 3 in [23]. *

Corollary A.3**.**

There exists $C>0$ depending on $\Gamma_{D}$, but independent of the mesh-boundary intersection, such that for every $K\in\mathcal{M}_{0,h}$

[TABLE]

Proof A.4**.**

Let us apply Lemma A.1:

[TABLE]

By a standard inverse inequality, see [3], we get

[TABLE]

*where $C>0$ depends on the shape regularity constant of the un-trimmed mesh $\mathcal{M}_{0,h}$. *

Lemma A.5**.**

Let $Q,Q^{\prime}\in\widehat{\mathcal{M}}_{0,h}$ be neighbor elements in the sense of Definition 5.1. There exists $C>0$ such that

[TABLE]

*where $C$ depends on $k$, on the shape regularity of the mesh and on the distance between $Q$ and $Q^{\prime}$. *

Proof A.6**.**

*The proof follows by a standard scaling argument (see [39] for instance). *

The next one says that the $L^{2}$ norm on the cut portion of an element $Q$ controls the $L^{\infty}$ (and hence any other) norm on the whole element with an equivalence constant depending on the relative measure of the cut portion.

Lemma A.7**.**

Let $\theta\in(0,1]$. There exists $C>0$ such that for every $Q\in\widehat{\mathcal{M}}_{0,h}$ and every $S\subset Q$ measurable such that $\left\lvert S\right\rvert\geq\theta\left\lvert Q\right\rvert$, we have

[TABLE]

*where $C$ depends only on $\theta$, $k$ and the mesh regularity. *

Proof A.8**.**

*See Proposition 1 in [22]. *

Lemma A.9** (Hardy’s inequality, [9]

).**

Let $\Omega\subset\R^{d}$ be a bounded open set of class $C^{1}$. Then there is a constant $C>0$ such that

[TABLE]

*where $\text{d}(x):=\operatorname{dist}(x,\Gamma)$. *

Remark A.10**.**

*Viceversa, it is possible to characterize functions in $H^{1}_{0}(\Omega)$ as functions in $H^{1}(\Omega)$ such that $\frac{u}{d}\in L^{2}(\Omega)$ ([9]). *

Lemma A.11**.**

Let $\Omega_{1}\subset\Omega$ with boundary $\Gamma_{1}$ such that $\Omega_{1}=\{x\in\Omega:\operatorname{dist}(x,\Gamma)\geq Ch\}$, where $C\geq 1$ fixed and $\operatorname{dist}(\Gamma,\Gamma_{1})\leq Ch$. It holds that

[TABLE]

*where the interpolation space $H^{s}_{i}(\Omega)$ or $\left(H^{1}_{0}(\Omega),L^{2}(\Omega)\right)_{s,2}$ is isomorphic to $H^{s}(\Omega)$ for $0\leq s<\frac{1}{2}$, to $H^{\frac{1}{2}}_{00}(\Omega)$ for $s=\frac{1}{2}$ and to $H^{s}_{0}(\Omega)$ for $\frac{1}{2}<s\leq 1$ (see [46]). *

Proof A.12**.**

We prove the following (like in [29]):

[TABLE]

We define $d(x):=\operatorname{dist}(x,\Gamma)$ $\forall\ x\in\Omega\setminus\overline{\Omega}_{1}$. By assumption $d(x)\leq Ch$, hence $1\leq\frac{Ch^{2}}{\left\lvert d(x)\right\rvert^{2}}$.

[TABLE]

where we employed Hardy’s inequality from Lemma A.9. Moreover:

[TABLE]

At this point, let us interpolate estimates (31) and (32), getting

[TABLE]

Acknowledgements

We would like to thank Pablo Antolín who provided us a tool to perform integration on trimmed geometries in GeoPDEs.