Graded Parametric CutFEM and CutIGA for Elliptic Boundary Value Problems in Domains with Corners

TL;DR

This paper introduces a graded parametric CutFEM and CutIGA approach for elliptic boundary value problems with corner singularities, achieving optimal stability and convergence without needing to identify the corner angle.

Contribution

The method uses a mesh grading strategy via a suitable mapping that does not require knowing the corner's opening angle, enhancing robustness and applicability.

Findings

Proves stability and optimal convergence of the method.

Achieves accurate solutions near corners without angle identification.

Employs structured meshes in the reference domain.

Abstract

We develop a parametric cut finite element method for elliptic boundary value problems with corner singularities where we have weighted control of higher order derivatives of the solution to a neighborhood of a point at the boundary. Our approach is based on identification of a suitable mapping that grades the mesh towards the singularity. In particular, this mapping may be chosen without identifying the opening angle at the corner. We employ cut finite elements together with Nitsche boundary conditions and stabilization in the vicinity of the boundary. We prove that the method is stable and convergent of optimal order in the energy norm and $L^2$ norm. This is achieved by mapping to the reference domain where we employ a structured mesh.

Figures

Equations

$$u=\mathcal{O}(r^{\alpha }),0<\alpha <1,\text{as r→0},$$

$$-\Delta u$$

$$u$$

$$u=u_{\mathrm{r}}+u_{\mathrm{s}}$$

$$\Vert u_{\mathrm{r}}{\Vert }_{H^s(\Omega )}\lesssim \Vert f{\Vert }_{H^{s-2}(\Omega )}+\Vert g{\Vert }_{H^{s-1/2}(\partial \Omega )}$$

$$u_{\mathrm{s}}(r,\theta )=r^{\pi /{\omega }_c}\sin (\theta \pi /{\omega }_c)$$

$$|D^{k}v{|}^2=({\partial }_r^{k}v{)}^2+\sum _{\begin{array}{c}m+n\le k,\\ m\ge 0,n\ge 1\end{array}}{\left(r^{m-k}{\partial }_r^m{\partial }_{\theta }^{n}v\right)}^2$$

$$\Vert u_{\mathrm{s}}{\Vert }_{H^{k,\alpha }(\Omega )}^2=\sum _{j=1}^k\Vert r^{{\alpha }_j}{D}^j{u}_{\mathrm{s}}{\Vert }_{\Omega }^2<\infty$$

$$j-1-\frac{\pi }{{\omega }_c}<{\alpha }_j$$

$$|{\partial }_r^j{u}_{\mathrm{s}}|$$

$$|r^{m-j}{\partial }_r^m{\partial }_{\theta }^n{u}_{\mathrm{s}}|$$

$$\Vert r^{{\alpha }_j}{D}^j{u}_{\mathrm{s}}{\Vert }_{\Omega }^2$$

$$={\int }_0^R{r}^{2{\alpha }_j+2(\pi /{\omega }_c)-2j+1}dr$$

$$={\left[\frac{r^{2{\alpha }_j+2(\pi /{\omega }_c)-2(j-1)}}{2{\alpha }_j+2(\pi /{\omega }_c)-2(j-1))}\right]}_0^R$$

$$(j-1)-(\pi /{\omega }_c)<{\alpha }_j$$

$$F_{\gamma }:\left[\begin{array}{c}\hat{r}\\ \hat{\theta }\end{array}\right]\mapsto \left[\begin{array}{c}{\hat{r}}^{\gamma }\\ \hat{\theta }\end{array}\right]=\left[\begin{array}{c}r\\ \theta \end{array}\right]$$

$$\hat{\Omega }=F_{\gamma }^{-1}(\Omega )$$

$$v(r)=\hat{v}(\hat{r}):\hat{\Omega }\to \mathbb{R}$$

$$\Vert \nabla u{\Vert }_{\Omega }\sim \Vert \hat{\nabla }\hat{u}{\Vert }_{\hat{\Omega }}$$

$$\Vert {\hat{D}}^k\hat{u}{\Vert }_{\hat{\Omega }}^2\lesssim \sum _{j=1}^k\Vert r^{{\alpha }_j}{D}^{j}u{\Vert }_{\Omega }^2$$

$${\alpha }_j=(j-k)+{\alpha }_k\text{with}{\alpha }_k=(k-1)\frac{\gamma -1}{\gamma }$$

$$r={\hat{r}}^{\gamma },\hat{r}=r^{1/\gamma }$$

$$\frac{dr}{d\hat{r}}=\gamma {\hat{r}}^{\gamma -1}=\gamma r^{(\gamma -1)/\gamma },\frac{d\hat{r}}{dr}={\gamma }^{-1}{r}^{(1-\gamma )/\gamma }={\gamma }^{-1}{\hat{r}}^{(1-\gamma )}$$

$$rdrd\theta =\gamma {\hat{r}}^{2(\gamma -1)}\hat{r}d\hat{r}d\hat{\theta },\hat{r}d\hat{r}d\hat{\theta }={\gamma }^{-1}{r}^{(1-\gamma )/\gamma }drd\theta$$

$$\frac{du}{dr}=\frac{d\hat{u}}{d\hat{r}}\frac{d\hat{r}}{dr}=\frac{d\hat{u}}{d\hat{r}}{\gamma }^{-1}{\hat{r}}^{(1-\gamma )}$$

$$\Vert \nabla u{\Vert }_{\Omega }^2$$

$$\displaystyle=\int_{\widehat{\Omega}}\Big{(}|\hat{\partial}_{\hat{r}}\hat{u}|^{2}(\gamma^{-1}\hat{r}^{(1-\gamma)})^{2}+\hat{r}^{-2\gamma}|\hat{\partial}_{\hat{\theta}}\hat{u}|^{2}\Big{)}\gamma\hat{r}^{2(\gamma-1)}\hat{r}\,d\hat{r}d\hat{\theta}$$

$$\displaystyle=\int_{\widehat{\Omega}}\Big{(}\gamma^{-1}|\hat{\partial}_{\hat{r}}\hat{u}|^{2}+\gamma\hat{r}^{-2}|\hat{\partial}_{\hat{\theta}}\hat{u}|^{2}\Big{)}\hat{r}\,d\hat{r}d\hat{\theta}$$

$$\displaystyle\sim\int_{\widehat{\Omega}}\Big{(}|\hat{\partial}_{\hat{r}}\hat{u}|^{2}+\hat{r}^{-2}|\hat{\partial}_{\hat{\theta}}\hat{u}|^{2}\Big{)}\hat{r}\,d\hat{r}d\hat{\theta}$$

$$=\Vert \hat{\nabla }\hat{u}{\Vert }_{\hat{\Omega }}^2$$

$$|{\hat{D}}^k\hat{v}{|}^2=\underset{I}{\underbrace{({\hat{\partial }}_{\hat{r}}^k\hat{v}{)}^2}}+\stackrel{II}{\overbrace{\sum _{\begin{array}{c}m\ge 0,n\ge 1\\ k\ge m+n\end{array}}{\left({\hat{r}}^{m-k}{\hat{\partial }}_{\hat{r}}^m{\hat{\partial }}_{\hat{\theta }}^n\hat{v}\right)}^2}}$$

$${\hat{\partial }}_{\hat{r}}^k\hat{v}\sim \sum _{j=1}^k\sum _{\{q_{\ell }\}\in {\mathcal{Q}}_j^k}{\partial }_r^{j}v\cdot \prod _{\ell =1}^j{\hat{\partial }}_{\hat{r}}^{q_{\ell }}r$$

$${\hat{\partial }}_{\hat{r}}^{q_{\ell }}r=\left(\prod _{i=1}^{q_{\ell }}(\gamma -(i-1))\right){\hat{r}}^{\gamma -q_{\ell }}\sim r^{(\gamma -q_{\ell })/\gamma }$$

$$\sum _{\{q_{\ell }\}\in {\mathcal{Q}}_j^k}{\partial }_r^{j}v\cdot \prod _{\ell =1}^j{\hat{\partial }}_{\hat{r}}^{q_{\ell }}r$$

Full text

1 Introduction

A classical issue in the use and development of finite element methods (FEMs) are problems where the exact solution may be singular at certain points, due to the presence of nonconvex corners or jumps in data. In such cases standard FEMs perform poorly as the low regularity of the solution causes loss of convergence. By utilizing known information about the singularities, various methods that regain optimal order convergence have been devised. In this contribution we consider problems with nonconvex corners in the context of an unfitted finite element method – the cut finite element method (CutFEM) – and propose a CutFEM for which we prove stability and optimal order error estimates.

CutFEM and CutIGA.

CutFEM is a framework for finite element methods where the physical domain may cut the computational mesh arbitrarily while retaining the optimal approximation and stability properties of standard FEMs, see [6, 5]. This framework is applicable also to isogeometric analysis (IGA), which combines powerful spline approximation spaces on structured computational meshes with high precision geometry descriptions via CAD, see, e.g., [12, 7, 19]. In general, IGA generates a cut computational mesh as the CAD consists of patchwise parametric mappings and trim curves. It is therefore natural to utilize the mathematically rigorous CutFEM framework also in IGA. We denote this combination of techniques cut isogeometric analysis (CutIGA), which we have previously explored in [8, 9].

Previous Work.

Over the years considerable efforts have been made to construct special finite element methods that deal with singularities arising from domains with (nonconvex) corners or jumps in data, see the very extensive survey [17] and the references therein. A great deal about these singularities is actually known, especially in the planar domain case, and this information is also typically needed to construct a method that efficiently deals with singularities. The survey [17] classifies methods into the following three categories:

1. I.

Methods involving local refinement. These are based on the assumption that the location of the singular point is known, and that the behavior of the solution when approaching the singular point is

[TABLE]

where $r$ is the distance to the singular point. Most methods belong to this category, as does the method in the present work, and it includes standard approaches such as $h/p$-FEM, and various parametric techniques. If the location of the singular point is unknown it is possible to instead use an adaptive procedure where the approximation space is iteratively tuned to the problem. Such procedures are typically based on a posteriori error indicators and the tuning can consist of local mesh refinement ($h$-refinement), locally increasing the polynomial order ($p$-refinement), and relocation of the mesh ($r$-refinement), see, e.g., [1, 2, 14]. 2. II.

Methods supplementing the approximation space with singular functions. Here the assumption is that the leading singular functions in the expansion of the solution when $r\to 0$ are known. These can then be supplemented to the usual approximation space. Typically, one or two singular functions is sufficient to resolve the singular effects. In this category we also find methods such as XFEM/GFEM, see [3, 10]. 3. III.

Combined methods. This final category is based on the assumption that the complete expansion of the solution, i.e., both the singular and analytical parts, is locally known in a singular subdomain $r\leq r_{0}$.

In the context of IGA there have been some recent contributions dealing with corner singularities. A method belonging to the first category was proposed in [13] where the parametric mappings describing the geometry were modified to grade the approximation space appropriately towards the singularity. In another contribution [20] a method related to the first and second category was presented in which the approximation space was enriched by certain basis functions constructed using push forward operations with a mapping similar to the one in [13].

The present work belongs to the first category and is based on a mapping in the form of simple radial scaling, which has the effect of a graded mesh. As the CutFEM uses a geometric description of the domain that is independent of the computational mesh, the mapping only affects the computational mesh, not the domain. The simple expression for the mapping allows us to conveniently reformulate the method in a way that avoids numerical instabilities due to large derivatives of the mapping close to the singular point. Also, using a mapping which is smooth everywhere (except at the singular point) means that the regularity of the discrete approximation space will not be affected, making the method suitable also for IGA approximation spaces of arbitrarily high regularity.

New Contributions.

We develop a CutFEM for elliptic problems for domains with corners, including treatment of singularities arising at nonconvex corners. Our approach is based on identification of a suitable mapping from a reference domain, where a quasi-uniform mesh is used, that grades the mesh towards the singular corner. The mapping is bijective but has zero derivative in the corner. We employ weak enforcement of Dirichlet conditions, allow cut elements at the boundary and stabilize the method using Ghost penalty, see [4, 18]. We prove that the method is stable by mapping to the reference domain and utilizing the additional stability provided by the Ghost penalty. Furthermore, we can prove optimal order a priori error estimates in the energy and $L^{2}$ norm using a bound on derivatives of order $k$ in the reference domain in terms of a weighted norm in the physical domain which is bounded for the singular solution. The analysis holds for polygonal domains and piecewise smooth boundaries. Furthermore, the method is not sensitive to the choice of grading parameter as long as the grading is strong enough. As a consequence, the grading that works for a corner with opening angle close to $2\pi$ works for all angles. This observation also indicates that the approach may be extended to three dimensional situations, which we plan on investigating in future work.

Outline.

The remainder of this work is organized as follows: In Section 2 we formulate the model problem and describe its regularity properties in terms of weighted norms and then we introduce the grading mapping; in Section 3 we present the parametric cut finite element method, expressed both in the physics domain and in the reference domain; in Section 4 we prove our error estimates; in Section 5 we present numerical examples using tensor product elements of various order and regularity; and finally, in Section 6 we give some concluding remarks.

2 The Model Problem and Regularity Properties

2.1 Model Problem

Let $\Omega$ be a domain in $\mathbb{R}^{2}$ with boundary $\partial\Omega$ and consider the elliptic problem: find $u:\Omega\rightarrow\mathbb{R}$ such that

[TABLE]

We assume that the boundary $\partial\Omega$ consists of a finite number of smooth curve segments $\{\Gamma_{i}\}$ that meet in corners $\mathcal{C}=\{c_{j}\}$, some of which are nonconvex. Such a domain is illustrated in Figure 1(a). For brevity let us from now on consider the situation that we have one corner $c$ which is nonconvex, as illustrated in Figure 1(b). Since the singular behavior is local, the extension to several singular corners is straightforward.

2.2 The Regularity of the Solution

The solution $u$ to the above problem with a single nonconvex corner can be decomposed into a regular and singular part

[TABLE]

where the regular part fulfills the elliptic shift property

[TABLE]

In the vicinity of the corner $c$ the singular part of the solution takes the form

[TABLE]

where the parameter $\omega_{c}\in(\pi,2\pi)$ is the opening angle, see Figure 1(b). Here we use polar coordinates $\{r,\theta\}$ centered at $c$ and with the angular component $\theta$ measured counterclockwise from the tangent line originating in $c$ when passing the boundary in the counterclockwise direction. Since $1>\pi/\omega_{c}>1/2$ the solution is always in $H^{1}(\Omega)$ but not in $H^{2}(\Omega)$. More precisely the regularity can be described using weighted norms, see e.g. [16].

Remark 2.1 (More General Problems).

For brevity we focus on the Dirichlet problem (2.1)–(2.2). In literature more general boundary conditions are studied, for example mixed boundary conditions with both non-homogeneous Neumann and Dirichlet parts that gives singular parts of the solution different from (2.5). Another source than nonconvex corners for inducing singular solutions pertains to jumps in data. However, in most cases the method presented below will still be applicable as the essential assumption is (1.1) rather than the complete expression for the singular part of the solution.

Definition 2.1 (Total Derivative Magnitude).

The magnitude of the total derivative of order $k$ in polar coordinates is given by

[TABLE]

That this definition is reasonable is readily seen by expanding $(\partial_{r}+r^{-1}\partial_{\theta})^{k}v$, evaluating the derivatives of all coefficients of the form $r^{-\ell}$ using the product rule, and noting that the total derivative magnitude contains all terms in the expansion modulo integer coefficients.

Lemma 2.1 (Weighted Sobolev Space).

The singular part of the solution satisfies, for $k=1,2,\dots$,

[TABLE]

where the powers $\{\alpha_{j}\}$ must satisfy the condition

[TABLE]

• **Proof.**Observe that the derivatives appearing in $D^{j}u_{\mathrm{s}}$ have the bounds

[TABLE]

where $m,n$ are given in Definition 2.1 above. Let $R$ be the radius of a ball, which contains $\Omega$ and is centered at $c$. Using (2.9)–(2.10) we compute the following bound on the weighted norm

[TABLE]

and we note that this bound is finite for

[TABLE]

which is precisely the condition (2.8). ∎

2.3 Weighted Sobolov Spaces and Parametric Mapping

We will now construct a mapping from a reference domain to $\Omega$ such that the pullback of the singular part of the solution $u_{\mathrm{s}}$ has full regularity in terms of the reference coordinates.

The Parametric Mapping.

For convenience we here work with polar coordinates, both for the physical coordinates and for the reference coordinates. Consider the bijective mapping $F_{\gamma}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ defined by

[TABLE]

for $\gamma>0$, and define the domain in reference coordinates

[TABLE]

For $v:\Omega\to\mathbb{R}$ we define the pullback

[TABLE]

where we for notational brevity in the calculations below only include the radial coordinate as the angular coordinate is unchanged in the mapping (2.15).

Theorem 2.1 (Equivalence of Norms).

(1) For $\gamma>0$, it holds

[TABLE]

*where the hidden constants depend only on $\gamma$.

(2) For $k=1,2,\dots$, and $\gamma>0$, it holds*

[TABLE]

where $\widehat{D}^{k}\hat{u}$ is the total derivative of order $k$ in the reference domain and

[TABLE]

• Proof.The First Order Case (2.18). Consider the mapping

[TABLE]

with derivatives

[TABLE]

and measures

[TABLE]

Using that $u(r)=\hat{u}(\hat{r})$ we by the chain rule have

[TABLE]

Taking the norm of the gradient and changing coordinates we find the equivalence

[TABLE]

with constants only dependent on $\gamma$ and $\gamma^{-1}$.

The General Case (2.19).

By Definition 2.1 the total derivative magnitude is

[TABLE]

Term $\bm{I}$.

Using the chain rule we get the following expression with derivatives in terms of physical coordinates instead of reference coordinates

[TABLE]

where $\mathcal{Q}_{j}^{k}$ consists of all sets $\{q_{\ell}\}_{\ell=1}^{j}$ of positive integers such that $\sum_{\ell=1}^{j}q_{\ell}=k$. Taking $q_{\ell}$ derivatives of $r$ with respect to $\hat{r}$ we get

[TABLE]

with a constant depending only on $\gamma$. This gives the equivalence

[TABLE]

where we in the second equivalence used that $\sum_{\ell=1}^{j}q_{\ell}=k$. In summary, we have the equivalence

[TABLE]

For the measure we have

[TABLE]

Integrating Term $I$ over the reference domain, using the equivalence (2.35), and the Cauchy–Schwarz inequality we get

[TABLE]

where we identify

[TABLE]

Term $\bm{II}$.

Again using the chain rule we get the following expression with derivatives in terms of physical coordinates instead of reference coordinates

[TABLE]

where we recall that by definition $\prod_{\ell=1}^{0}\hat{\partial}_{\hat{r}}^{q_{\ell}}r=1$. By (2.32) we have the equivalence

[TABLE]

where we in the second equivalence used that $\sum_{\ell=1}^{j}q_{\ell}=m$. In summary, we have the equivalence

[TABLE]

Integrating Term $II$ over the reference domain gives

[TABLE]

with $\alpha_{k}$ as in (2.39). Here we in (2.47) used the equivalence (2.43)–(2.45) and the Cauchy-Schwarz inequality to move the sum over $j$ outside the square; in (2.48) we used the fact that all index pairs $\{j,n\}$ generated by the inner sum over $j$ are also generated by the outer sum as index pairs $\{m,n\}$; and in (2.49) we used the transformation of the measure (2.36).

Conclusion.

From (2.37)–(2.38) and (2.46)–(2.50) we have

[TABLE]

with $\alpha_{j}$ given by (2.39). ∎

Remark 2.2 (The Scaling Parameter $\bm{\gamma}$).

Combining the condition (2.8) for $u_{\mathrm{s}}$ to be in the weighted Sobolev space with the expression for $\alpha_{j}$ (2.20) gives the inequality

[TABLE]

Choosing the scaling parameter $\gamma$ according to this inequality means that the pullback of the solution $\hat{u}$, including the singular part, is in $H^{k}(\widehat{\Omega})$.

3 The Parametric Cut Finite Element Method

We will now construct a parametric cut finite element method based on the map $F_{\gamma}$. Since the derivative of the mapping is zero in the origin for $\gamma>1$ we define the mesh parameters in the Nitsche forms using the well defined mesh parameter in the reference domain where we use a quasi-uniform mesh.

3.1 The Mesh and Finite Element Space

• •

Let $F:\widehat{\Omega}\rightarrow\Omega$ with $\widehat{\Omega}\subset[-1,1]^{2}\subset\mathbb{R}^{2}$ where $\widehat{\Omega}$ is the reference domain such that the center $(0,0)\in\partial\widehat{\Omega}$ corresponds to the singular point. For simplicity we here assume $F=F_{\gamma}$ but in general the mapping $F$ is a composite mapping $F=F_{*}\circ F_{\gamma}$ where the additional bijective mapping $F_{*}$ is readily incorporated in the method below using standard techniques.

• •

Let $\widehat{\mathcal{K}}_{h,0}$ be a quasi-uniform mesh on $[-1,1]^{2}$ consisting of shape regular elements with mesh parameter $h$. Typically we employ uniform quadrilaterals. Let $\widehat{\mathcal{K}}_{h}=\{K\in\widehat{\mathcal{K}}_{h,0}:K\cap\Omega\neq\emptyset\}$ be the active mesh in the reference domain. Let $\mathcal{K}_{h}=F(\widehat{\mathcal{K}}_{h})$ be the induced active mesh in the physical domain.

• •

Let $\widehat{\mathcal{F}}_{h}$ be the set of interior faces in $\widehat{\mathcal{K}}_{h}$ such that each face in $\widehat{\mathcal{F}}_{h}$ is adjacent to at least one element $K\in\widehat{\mathcal{K}}_{h}$ that is cut by the boundary $\partial\widehat{\Omega}$, i.e., $\overline{K}\cap\partial\widehat{\Omega}\neq\emptyset$.

• •

Let $\widehat{V}_{h,0}$ be a finite element space consisting of continuous piecewise polynomials of order $p$ defined on $\widehat{\mathcal{K}}_{h,0}$ and let $\widehat{V}_{h}=\widehat{V}_{h,0}|_{\widehat{\mathcal{K}}_{h}}$ be the active finite element space in the reference domain. The active finite element space in the physical domain is given by $V_{h}=\widehat{V}_{h}\circ F^{-1}$.

• •

It is noted in the analysis below, see Remark 4.1, that by choosing $\gamma=2p$ in the mapping we obtain optimal order convergence for our method regardless of the opening angle $\omega_{c}\in(\pi,2\pi)$.

3.2 The Method

The Method in Physical Cartesian Coordinates.

Find $u_{h}\in V_{h}$ such that

[TABLE]

where

[TABLE]

Here $[\cdot]$ denotes the jump between neighboring elements; $\pi_{h}$ is a projection operator onto $V_{h}$, which we define in the analysis below; and $h_{\Omega}\sim h\hat{r}^{\gamma-1}=hr^{\frac{\gamma-1}{\gamma}}$ is a mesh function defined in such a way that

[TABLE]

where the reference forms $\widehat{a}_{h}$ and $\widehat{l}_{h}$ are given below. Note that since $\pi_{h}$ is a projection, i.e., $\pi_{h}v=v$ for $v\in V_{h}$, the method (3.1) is unaffected by its presence. The stability form $s_{h}$ is naturally defined in the reference domain since we, in the analysis below, prove coercivity in the reference domain. It is furthermore easier to evaluate higher-order derivatives in the reference domain.

Galerkin Orthogonality.

Due to the inclusion of the projection operator $\pi_{h}$ the method is not consistent yielding the perturbed Galerkin orthogonality

[TABLE]

Remark 3.1 (Inconsistency).

This inconsistency is artificial in the sense that for sufficiently regular $u$ the interpolation operator $\pi_{h}$ is not needed and $s_{h}(u,v)=0$. However, we include $\pi_{h}$ in $s_{h}$ so that the form $A_{h}$ is well defined also for exact solutions of lower regularity, which is the case of the solution to the dual problem in the proof of the $L^{2}$ estimate below.

The Method in Reference Cartesian Coordinates.

Transforming back to Euclidian coordinates $(\hat{x},\hat{y})=(\hat{r}\cos\hat{\theta},\hat{r}\sin\hat{\theta})$ in the reference domain we obtain the method

[TABLE]

where $\widehat{A}_{h}(\hat{v},\hat{w})=\widehat{a}_{h}(\hat{v},\hat{w})+\hat{s}_{h}(\hat{v},\hat{w})$. The reference forms are defined by (3.4) and

[TABLE]

where $\hat{n}$ is the outward pointing unit normal to $\partial\widehat{\Omega}$ in Cartesian reference coordinates and we employed the notation

[TABLE]

See Appendix A for details of the derivation. We note that $B$ is symmetric and positive definite

[TABLE]

with $c_{\gamma}=\min(\gamma,\gamma^{-1})$, and uniformly bounded

[TABLE]

with $C_{\gamma}=\max(\gamma,\gamma^{-1})$. Furthermore, since $S_{\hat{\theta}}$ is a rotation matrix the identity $S_{\hat{\theta}}S_{\hat{\theta}}^{T}=I$ holds and also $B^{m}$ is uniformly bounded

[TABLE]

Remark 3.2 (Numerical Stability).

By implementing the method in reference coordinates according to (3.8) we conveniently avoid numerical issues induced by the mapping $F_{\gamma}$ having unbounded derivatives when approaching the singular point. Even though the forms in the physical respectively reference domains are mathematically identical, evaluation of the forms in the physical domain involves products between terms tending to infinity and to zero when approaching the singular point, causing numerical accuracy and stability issues. Thanks to simplification, this issue does not exist in the reference domain forms (3.9) and (3.10).

3.3 The Method with Multiple Nonconvex Corners

We consider situations where the domain $\Omega$ is partitioned into a finite set of patches $\{\Omega_{i}\}_{i=0}^{n}$ in such a way that each singular corner belongs to one patch as illustrated in Figure 2. The scaling can then be applied locally to each patch $\Omega_{i}$ containing a singular corner. The coupling between the patches is done weakly using Nitsche’s method which does not require the meshes to match across the interface, see the parametric multipatch method in [15]. This approach enables us to use the grading locally in the vicinity of each corner without difficulty.

Multipatch Meshes and Finite Element Spaces.

For each patch $\Omega_{i}$ we construct a mapping $F^{(i)}$, a mesh $\mathcal{K}_{h}^{(i)}$ and a finite element space $V_{h}^{(i)}$ as outlined in Section 3.1. On the complete multipatch domain $\Omega$ we define the finite element space

[TABLE]

where we note that $V_{h}$ is discontinuous across patch interfaces. This will instead be enforced weakly in the method.

The Multipatch Method in Physical Cartesian Coordinates.

Defining the domain $\Omega=\bigcup_{i=0}^{n}\Omega_{i}$ and the patch interface $\Gamma=\bigcup_{i<j}\overline{\Omega}_{i}\cap\overline{\Omega}_{j}$ we adapt formulation (3.1) to the multi-patch setting by redefining $a_{h}$ and $s_{h}$ as

[TABLE]

In (3.16) we have added Nitsche interface terms that couple the patches and (3.17) now includes Ghost penalty stabilization on the computational mesh of each patch. In the interface terms $\langle\cdot\rangle$ denotes the average between adjoining patches while $[\cdot]$ denotes the jump. Expressing the interface terms in reference Cartesian coordinates is straightforward and we refer to [15] for further details.

4 A Priori Error Estimates

Main Approach.

• •

We start by defining an energy norm associated with the form $\widehat{A}_{h}$ in the reference domain and we establish coercivity and continuity using standard arguments since the mesh is quasi-uniform and $B$ is positive definite and bounded, see (3.12) and (3.13).

• •

Next we define a corresponding energy norm in the physical domain using the mapping which leads to the corresponding coercivity and continuity results for the physical problem.

• •

We define an interpolation operator in the reference domain and recall a standard interpolation error estimate. This then leads to an interpolation error estimate in the physical domain using (2.19).

• •

Using the above results optimal order a priori results follows using standard techniques.

4.1 Basic Results

Extension Operator.

Let $U(\Omega)\subset\mathbb{R}^{2}$ be a neighborhood of $\Omega$ such that $\mathcal{K}_{h}\subset U(\Omega)$ for all $h\in(0,h_{0}]$. Due to the fact that $\omega_{c}<2\pi$ there is an extension operator $\widehat{E}:H^{s}(\widehat{\Omega})\rightarrow H^{s}(\widehat{U(\Omega)})$ that satisfies the stability estimate

[TABLE]

Interpolant.

Let $\widehat{\pi}_{h}:H^{1}(\widehat{U(\Omega)})\rightarrow\widehat{V}_{h}$ be a Scott–Zhang type interpolation operator and recall that we have the error estimate

[TABLE]

In the physical domain we define $\pi_{h}:H^{1}(U(\Omega))\rightarrow V_{h}$ by

[TABLE]

Since the Scott–Zhang interpolation operator is a projection, this is also the case for $\pi_{h}$.

Remark 4.1 (Choice of Scaling Parameter $\bm{\gamma}$).

To obtain optimal order interpolation estimates we need control over $s=p+1$ derivatives in the reference domain. A condition on $\gamma$ for obtaining full regularity of the pullback of the solution in the reference domain was given in Remark 2.2 and choosing $k=p+1$ therein we arrive at the condition

[TABLE]

As $\omega_{c}\in(\pi,2\pi)$ a choice that holds for any opening angle is $\gamma=2p$.

Stabilization Form.

By construction the stabilization form (3.4) satisfies the following abstract properties:

• •

The form is consistent, i.e.,

[TABLE]

• •

The form satisfies the interpolation estimate

[TABLE]

• •

The following inverse inequality holds

[TABLE]

Energy Norm.

Define the following energy norms associated with the Nitsche method

[TABLE]

and the push forward energy norm

[TABLE]

Interpolation in Energy Norm.

For $\hat{v}\in H^{s}(\widehat{\Omega})$ we on every $K\in\widehat{\mathcal{K}}_{h}$ have $(\hat{v}-\widehat{\pi}_{h}\hat{v})|_{K}\in H^{s}(K)$. Assuming $s\geq 2$ and using the element-wise trace inequality

[TABLE]

which holds independent of the position of $\partial\widehat{\Omega}$ in $K$, see [11], and the interpolation estimate (4.2) we obtain

[TABLE]

Using the definition of the push forward energy norm (4.9), the interpolation estimate in the reference domain (4.11), and the estimate (2.19), we obtain the following interpolation estimate in the physical domain

[TABLE]

More precisely we proceeded as follows

[TABLE]

Continuity and Coercivity.

It holds

[TABLE]

and for $\beta>0$ large enough

[TABLE]

These results follows directly using standard arguments for cut finite elements since the mesh is quasi-uniform in the reference domain and $B$ is positive definite and bounded, see (3.12) and (3.13).

• Proof.(4.14). Using the triangle and the Cauchy–Schwarz inequalities, and the bound (3.14) on $B^{2}$ we readily get

[TABLE]

with a hidden constant $c=C_{\gamma}^{2}+\beta$.

(4.15).

Using the Cauchy–Schwarz inequality on the consistency term, and applying a Young’s inequality $2ab\leq\delta a^{2}+\delta^{-1}b^{2}$ with $\delta>0$, we get the bound

[TABLE]

where we in (4.21) use that $B^{2}$ is bounded via (3.14) and in the last step utilize the inverse inequality (4.7) with hidden constant $c_{I}$, which is possible since $\hat{v}\in\widehat{V}_{h}$. By the positive definiteness of $B$ (3.12) and the bound (4.19)–(4.22) we now obtain

[TABLE]

Choosing $\delta<c_{\gamma}c_{I}^{-1}C_{\gamma}^{-2}$ and $\beta>c_{I}\delta^{-1}$ yield positive constants in (4.25). The proof is completed by transforming part of the $\left(\|\widehat{\nabla}\hat{v}\|_{\widehat{\Omega}}^{2}+\hat{s}_{h}(\hat{v},\hat{v})\right)$ term into a $\|h^{1/2}\widehat{\nabla}\hat{v}\|_{\partial\widehat{\Omega}}^{2}$ term via the inverse inequality (4.7). ∎

Analogously, in physical coordinates it holds

[TABLE]

and for $\beta$ large enough

[TABLE]

These results follows using the definition of the energy norms (4.9) and the corresponding results (4.14) and (4.15) in the reference domain. More precisely

[TABLE]

and

[TABLE]

4.2 Error Estimates

Theorem 4.1.

For $\gamma>p\omega_{c}/\pi$ it holds

[TABLE]

with powers $\{\alpha_{j}\}_{j=1}^{p+1}$ in the weighted Sobolev norm given by (2.39) using $k=p+1$.

• Proof.(4.30). Using the coercivity (4.27), Galerkin orthogonality (3.7), continuity (4.26), and the interpolation estimate (4.12), the estimate follows directly

[TABLE]

(4.31).

For the $L^{2}$ estimate we let $\phi\in H^{1}(\Omega)$ be the solution to the dual problem

[TABLE]

for which we have the weighted regularity estimate

[TABLE]

Setting $\psi=u-u_{h}$ and using the dual problem (4.36), Galerkin orthogonality (3.7), continuity (4.26) and the Cauchy–Schwarz inequality, the interpolation estimates (4.12) and (4.6), and the regularity estimate (4.37), we obtain

[TABLE]

which completes the proof. ∎

5 Numerical Experiments

Implementation.

The method was implemented using the formulation in reference Cartesian coordinates (3.8) and other than that we follow the implementation outlined in [15] including the quadrature. A limitation due to our current implementation of certain geometrical operations, for example computing the intersection between an element and the domain, is that the geometry must be represented as a polygon. However, in our examples this polygon approximation is of very high resolution compared to the mesh size and should have a negligible effect on the results.

Approximation Space.

As approximation space $V_{h}$ we use the space spanned by piecewise quadratic B-splines of maximum regularity, i.e., $V_{h}\subset C^{1}(\Omega)\subset H^{2}(\Omega)$. These are constructed as tensor products of 1D quadratic B-spines on a structured grid, which in general is cut by the domain. For reference purposes we also briefly consider conforming approximation spaces in the form of standard Lagrange triangular elements of various orders.

Parameter Values.

For the mapping $F_{\gamma}$ we choose the radial scaling parameter $\gamma=2p$, as this choice holds for all opening angles $\omega_{c}\in(\pi,2\pi)$, see Remark 4.1. In the method we use the Nitsche penalty parameter $\beta=100$ and the stabilization form parameter $\tau=0.1$.

5.1 Convergence and Stability

Model Problem.

To investigate the convergence and stability of the method we construct a model problem on a domain in the shape of a unit circle sector with various opening angles $\omega_{c}\in(\pi,2\pi)$, see Figure 3. On this domain we manufacture a problem with known singular solution by taking (2.5) as an ansatz for $u$ and deriving the corresponding right hand side $f$ and Dirichlet boundary data $g$.

Convergence.

We here provide numerical studies of convergence using the model problem as extra validation of our theoretical a priori error estimates in Theorem 4.1. First we, as a reference, apply the method to conforming finite element spaces in the form of standard Lagrange triangles of order $p=1,2,3$. These results are presented in Figure 4 where we note that the mapping indeed has the desired effect of recovering the optimal order convergence as implied by our estimates.

Next we apply the method to cut finite element spaces in the form of $C^{1}$ splines, i.e., B-splines of order $p=2$, and the results are presented in Figure 5. In this experiment we note an at least initial loss of optimal convergence rate when the opening angle $\omega_{c}$ approaches $2\pi$. Our explanation for this is that since the finite element space is constructed independently from the geometry the finite element space does not take the opening slit into account. Thus, every basis function that have support on both sides of the slit actually provides an undesirable coupling over the slit. To remedy this we in the next paragraph outline and numerically assess a fix.

Fix for Removing Unwanted Coupling.

As a simple fix for removing the undesired coupling over the slit we propose the following approach: Basis functions whose support consists of two disjoint parts, one above and one below the slit, are associated with two separate degrees of freedom, one connected to elements above the slit and one connected to elements below the slit. This is illustrated in Figure 6(a). Note that this fix will not effect the regularity of the numerical solution. However, in the vicinity of the corner $c$ there will still exist some unwanted coupling as illustrated in Figure 6(b).

We numerically assess the effect of this fix in Figure 7 where we in (a)–(b) note that we with again recover the optimal order convergence implied by our estimates also for opening angles $\omega_{c}$ close to $2\pi$ when using a cut $C^{1}$ spline approximation space. In Figures 7(c)–d we note that the effect of the fix is more pronounced the larger the opening angle. We should however remark that this fix is only necessary in more extreme cases of opening angles or when using a quite large mesh size.

Stability with Respect to Mesh Position.

In the previous examples we have positioned the mesh such that the nonconvex corner $c$ is located in the center of an element when using $C^{1}$ spline approximation spaces. To study how the mesh position in relation to the nonconvex corner $c$ effects the performance of the method we in Figure 8 consider 400 random positions of the mesh in a $C^{1}$ spline approximation space with fixed mesh size $h=0.1$. We note that the method is actually quite insensitive with respect to the mesh position with a standard deviation of the error relative its mean of $0.02$ in $H^{1}(\Omega)$ semi-norm and of $0.04$ in $L^{2}(\Omega)$ norm.

5.2 More Examples

Multiple Nonconvex Corners.

In Figure 9 we present a numerical example on a domain featuring two nonconvex corners. This is solved using the multipatch approach outlined in Section 3.3. The domain is partitioned into three patches where each patch is equipped with its own approximation space and the appropriately chosen map. We note that the finite element solution and gradient magnitude seem to flow nicely over the internal interfaces.

Curved Surface.

Via the parametric map the method naturally handles problems on a curved surface $\Omega_{\mathbb{R}^{3}}\subset\mathbb{R}^{3}$ by the composite map $F=F_{*}\circ F_{\gamma}$ where $F_{*}:\Omega\to\Omega_{\mathbb{R}^{3}}$. As an illustration we present an example solution on a curved surface in Figure 10.

6 Conclusions

We have developed a new parametric higher-order cut finite element method for elliptic boundary value problems with corner singularities. It has the following notable features:

• •

The method is based on a radial map that suitably grades the mesh towards the singularity. This map can be chosen without knowing the exact opening angle $\omega_{c}<2\pi$.

• •

Numerical instabilities due to unbounded derivatives of the map near the corner are avoided by formulating the method in a reference domain.

• •

The method is proven to be stable and to be optimal order convergent in energy and $L^{2}(\Omega)$ norms.

• •

Multiple nonconvex corners are handled by using a previously developed multipatch framework such that each corner can be dealt with individually.

• •

Unwanted coupling over voids induced by the combination of extreme opening angles $\omega_{c}$ and cut approximation spaces is remedied by a proposed fix that restores the initial optimal order convergence to a large extent.

Acknowledgements.

This research was supported in part by the Swedish Foundation for Strategic Research Grant No. AM13-0029, the Swedish Research Council Grants Nos. 2013-4708, 2017-03911 and the Swedish Research Programme Essence.

Appendix A Riemannian Calculus Approach

We first change coordinates from Euclidean to polar and then we change from polar to weighted polar.

A.1 Polar Coordinates

Consider the mapping $F_{P}:B_{P}\rightarrow B\subset\mathbb{R}^{2}$ where $\mathbb{R}^{2}$ is equipped with the Euclidean inner product and

[TABLE]

with the partial derivatives

[TABLE]

The metric tensor $G_{P}$ is defined $(G_{P})_{ij}=\partial_{i}F_{P}\cdot\partial_{j}F_{P}$ which after simplification takes the form

[TABLE]

with inverse

[TABLE]

and determinant

[TABLE]

Gradient.

Let $v:B\rightarrow\mathbb{R}$ and let $v_{P}:B_{P}\rightarrow\mathbb{R}$ be the pullback

[TABLE]

The gradient satisfies

[TABLE]

which gives

[TABLE]

and

[TABLE]

Bilinear Form.

The bilinear form associated with the Laplacian transforms as follows

[TABLE]

where $d\mu=dxdy$ is the standard measure on $\mathbb{R}^{2}$ and $d\mu_{P}$ is the pullback measure.

Consistency Term.

Correspondingly the consistency term transforms as follows

[TABLE]

where $ds$ is the standard line measure on $\mathbb{R}^{2}$ and $(ds)_{P}$ is the pullback measure.

A.2 Radial Mapping in Polar Coordinates

Consider the mapping from $B_{\widehat{P}}\rightarrow B_{P}$ where $B_{P}$ is a disc centered at the origin such that

[TABLE]

with the partial derivatives

[TABLE]

The induced metric has elements $(G_{\widehat{P}})_{ij}=\partial_{i}F_{\widehat{P}}\cdot G_{P}\cdot\partial_{j}F_{\widehat{P}}$, where we use the polar metric inner product,

[TABLE]

with inverse

[TABLE]

and determinant

[TABLE]

Bilinear Form.

The bilinear form associated with the Laplacian transforms as follows

[TABLE]

Thus we obtain a form with a $\gamma$ dependent scaling of the radial and angular derivatives and in particular we note that no radial weight is present. Transforming back to Cartesian coordinates in the reference domain we have the identity

[TABLE]

where

[TABLE]

and we used the identities

[TABLE]

Consistency Term.

Correspondingly the consistency term transforms as follows

[TABLE]

where $ds_{\widehat{P}}$ is the standard line measure on $\mathbb{R}^{2}$. Here we used that

[TABLE]

and also utilized the simplification

[TABLE]

Transforming back to Cartesian coordinates in the reference domain we have the identity

[TABLE]

where we used the identities

[TABLE]