Cut Finite Element Methods for Elliptic Problems on Multipatch Parametric Surfaces

TL;DR

This paper introduces a cut finite element method for solving elliptic PDEs on complex, multipatch parametric surfaces, effectively handling trimmed patches and interface continuity.

Contribution

It develops a novel cut finite element approach with Nitsche's method for elliptic problems on multipatch surfaces, including stability and error analysis.

Findings

The method achieves optimal convergence rates in energy and L2 norms.

Numerical examples confirm theoretical stability and accuracy.

The approach effectively manages trimmed patches and interface conditions.

Abstract

We develop a finite element method for the Laplace--Beltrami operator on a surface described by a set of patchwise parametrizations. The patches provide a partition of the surface and each patch is the image by a diffeomorphism of a subdomain of the unit square which is bounded by a number of smooth trim curves. A patchwise tensor product mesh is constructed by using a structured mesh in the reference domain. Since the patches are trimmed we obtain cut elements in the vicinity of the interfaces. We discretize the Laplace--Beltrami operator using a cut finite element method that utilizes Nitsche's method to enforce continuity at the interfaces and a consistent stabilization term to handle the cut elements. Several quantities in the method are conveniently computed in the reference domain where the mappings impose a Riemannian metric. We derive a priori estimates in the energy and $L^2$ norm and also present several numerical examples confirming our theoretical results.

Tables

Table 1. Table 1 : Polynomial degree for nested 1D integrals and resulting number of quadrature points for each boundary segment assuming no optimizations, such as removing zero weight points, are used.

$(f,{\hat{\gamma }}_k)$	${\hat{x}}_2$ integrand	${\hat{x}}_1$ integrand	points/seg.
$(Q_2,P_1)$	$P_2$	$P_5$	$2\times 3=6$
$(Q_4,P_1)$	$P_4$	$P_9$	$3\times 5=15$
$(Q_6,P_1)$	$P_6$	$P_{13}$	$4\times 7=28$
$(Q_2,P_2)$	$P_2$	$P_{11}$	$2\times 6=12$
$(Q_4,P_2)$	$P_4$	$P_{19}$	$3\times 10=30$
$(Q_6,P_2)$	$P_6$	$P_{27}$	$4\times 14=56$

Equations

$$T_x({\Omega }_i)=\mathrm{span}\{{\hat{\partial }}_1{F}_i,{\hat{\partial }}_2{F}_i\}$$

$$a={\hat{a}}_1{\hat{\partial }}_1{F}_i+{\hat{a}}_2{\hat{\partial }}_2{F}_i=\left[\begin{array}{cc}{\hat{\partial }}_1{F}_i& {\hat{\partial }}_2{F}_i\end{array}\right]\left[\begin{array}{c}{\hat{a}}_1\\ {\hat{a}}_2\end{array}\right]=\hat{D}{F}_i\hat{a}$$

$$\hat{D}{F}_i=F_i\otimes \hat{\nabla }=\left[\begin{array}{cc}{\hat{\partial }}_1{F}_i& {\hat{\partial }}_2{F}_i\end{array}\right]$$

$${\hat{g}}_{\hat{x}}(\hat{a},\hat{b})=a\cdot b$$

$${\hat{g}}_{\hat{x}}(\hat{a},\hat{b})=\left[\begin{array}{cc}{\hat{a}}_1& {\hat{a}}_2\end{array}\right]\left[\begin{array}{cc}{\hat{g}}_{11}& {\hat{g}}_{12}\\ {\hat{g}}_{21}& {\hat{g}}_{22}\end{array}\right]\left[\begin{array}{c}{\hat{b}}_1\\ {\hat{b}}_2\end{array}\right]={\hat{a}}^T\hat{G}\hat{b}$$

$$\Vert \hat{v}{\Vert }_{{\hat{g}}_{\hat{x}}}^2={\hat{g}}_{\hat{x}}(\hat{v},\hat{v})$$

$$\Vert a{\Vert }_{{\mathbb{R}}^d}^2=a\cdot a={\hat{g}}_{\hat{x}}(\hat{a},\hat{a})=\Vert \hat{a}{\Vert }_{{\hat{g}}_{\hat{x}}}^2$$

$$\cos \theta =\frac{a\cdot b}{\Vert a{\Vert }_{{\mathbb{R}}^d}\Vert b{\Vert }_{{\mathbb{R}}^d}}$$

$$\hat{a}={\hat{G}}^{-1}\hat{D}{F}_i^{T}a$$

$${\hat{g}}_{\hat{x}}(\hat{a},\hat{b})=a\cdot b=a\cdot (\hat{D}{F}_i\hat{b})=(\hat{D}{F}_i^{T}a)\cdot \hat{b}\forall \hat{b}\in T_{\hat{x}}({\hat{\Omega }}_i)$$

$$c\le {\lambda }_{\min }(\hat{G})\le {\lambda }_{\max }(\hat{G})\le C\forall \hat{x}\in {\hat{\Omega }}_i$$

$$c\Vert \hat{a}{\Vert }_{{\hat{g}}_{\hat{x}}}^2\le \Vert \hat{a}{\Vert }_{{\mathbb{R}}^2}^2\le C\Vert \hat{a}{\Vert }_{{\hat{g}}_{\hat{x}}}^2\hat{x}\in {\hat{\Omega }}_i,\hat{a}\in T_{\hat{x}}({\hat{\Omega }}_i)$$

$${\int }_{\omega }fdx$$

$${\int }_{\gamma }fd\gamma ={\int }_{\hat{\gamma }}\hat{f}\Vert \hat{\tau }{\Vert }_{\hat{g}}d\hat{\gamma }={\int }_0^l\hat{f}\circ \hat{\gamma }(s)\Vert \hat{\tau }\circ \hat{\gamma }(s){\Vert }_{{\hat{g}}_{\hat{\gamma }(s)}}ds$$

$$(v,w{)}_{\omega }$$

$$\hat{\nabla }=\left[\begin{array}{c}{\hat{\partial }}_1\\ {\hat{\partial }}_2\end{array}\right]$$

$$\nabla u=\hat{D}{F}_i\widehat{\nabla u}$$

$$\hat{\nabla }\hat{u}=\hat{\nabla }(u\circ F_i)$$

$$\widehat{\nabla u}$$

$$(\text{div}v,w{)}_{\Omega }=-(v,\nabla w{)}_{\Omega }$$

$$\widehat{\text{div}v}=|\widehat{G}|^{-1/2}\widehat{\nabla}\cdot\Big{(}|\widehat{G}|^{1/2}\widehat{v}\Big{)}$$

$$\Delta v=\text{div}(\nabla v)$$

$$\widehat{\Delta v}=|\hat{G}{|}^{-1/2}\hat{\nabla }\cdot (|\hat{G}{|}^{1/2}{\hat{G}}^{-1}\hat{\nabla }\hat{v})$$

$$-(\Delta v,w{)}_{\omega }=(\nabla v,\nabla w{)}_{\omega }-(\nabla v,nw{)}_{\partial \omega }$$

$$(\hat{v},\hat{w}{)}_{H^k({\hat{\Omega }}_i)}=\sum _{|a|\le k}({\hat{D}}^a\hat{v},{\hat{D}}^a\hat{w}{)}_{L^2({\hat{\Omega }}_i)},\Vert \hat{v}{\Vert }_{H^s({\hat{\Omega }}_i)}^2=(\hat{v},\hat{v}{)}_{H^k({\hat{\Omega }}_i)}$$

$$H^k({\Omega }_i)=H^k({\hat{\Omega }}_i)\circ F_i^{-1}=\{v:\Omega \to \mathbb{R}|v=\hat{v}\circ F_i^{-1},\hat{v}\in H^k({\hat{\Omega }}_i)\}$$

$$(v,w{)}_{H^k({\Omega }_i)}=(\hat{v},\hat{w}{)}_{H^k({\hat{\Omega }}_i)},\Vert v{\Vert }_{H^k({\Omega }_i)}=\Vert \hat{v}{\Vert }_{H^k({\hat{\Omega }}_i)}$$

$$-\Delta u$$

$$\displaystyle\left\llbracket u\right\rrbracket$$

$$n_i\cdot \nabla u_i+n_j\cdot \nabla u_j$$

$$\llbracket w\rrbracket\big{|}_{\Gamma_{ij}}=w_{i}-w_{j}$$

Full text

Subject Classification Codes:

65M60, 65M85.

Keywords:

cut finite elements, fictitious domain, Nitsche’s method, a priori error estimates, multipatch surface, Laplace-Beltrami operator

1. 1 Introduction

2. 2 The Surface and the Laplace–Beltrami Operator

3. 2.1 Piecewise Parametric Description of the Surface

4. 2.2 Riemannian Metric on the Reference Patches

5. 2.3 Integration and $L^{2}$ Inner Products

6. 2.4 Differential Operators in Reference Coordinates

7. 2.5 Sobolev Spaces

8. 2.6 The Laplace–Beltrami Interface Problem

9. 3 The Finite Element Method

10. 3.1 Construction of the Mesh

11. 3.2 The Finite Element Spaces

12. 3.3 The Method

13. 3.4 Formulation in Reference Coordinates

14. 3.5 Quadrature

15. 4 A Priori Error Estimates

16. 4.1 Norms

17. 4.2 Inverse Inequalities

18. 4.3 Coercivity and Continuity

19. 4.4 Interpolation

20. 4.5 Error Estimates

21. 4.6 Condition Number Estimate

22. 5 Numerical Results

23. 5.1 Model Problems

24. 5.2 Implementation Aspects

25. 5.3 Convergence

26. 5.4 Stability

27. 5.5 Numerical Examples

28. 6 Summary and Future Work

29. A Differential Operators on Surfaces

1 Introduction

Background.

Differential equations on surfaces appear in many applications including transport phenomena on surfaces and elastic membranes and shells. In engineering applications the surface geometry is often described using a CAD model consisting of a partition of the surface into trimmed patches defined by mappings from a reference domain onto the surface. When performing computations on surfaces it is beneficial to directly utilize the available parametric geometry description, in line with the ideas of isogeometric analysis (IGA) **[9]**.

There are various techniques of enforcing interface conditions between patches, for example Lagrange penalty methods or methods based on weak enforcement. In the case of thin shells another approach is the bending strip method **[10]**. The use of Nitsche’s method **[16]**, or variants thereof, to weakly enforce interface conditions between patches is a well established and flexible technique, see for example **[1, 15, 12, 6]** and the references therein. However, constructing a high quality conforming mesh on the trimmed patches is generally a difficult task. In this work we address the problem of conforming mesh construction by allowing the trim curve on each patch to arbitrarily cut the mesh by utilizing a fictitious domain method called the cut finite element method (CutFEM) **[3, 2]**. In the same spirit **[18, 19, 11]** allow cut elements but employ the finite cell method, which is based on a different stabilization mechanism, where a small artificial stiffness is added on the part of the cut element outside of the patch it belongs to.

Contributions.

We develop a general technique for consistent discretization and a framework for analysis of the Laplace–Beltrami operator, which serves as a model second order partial differential operator on a patchwise parametric surface. The patches provide a partition of the surface and each patch is assumed to be the image by a diffeomorphism of a subdomain of the unit square which is bounded by a number of smooth trim curves. A patchwise tensor product mesh is constructed by using a structured mesh in the reference domain. Since the patches are trimmed we obtain cut elements in the vicinity of the interfaces. We discretize the Laplace–Beltrami operator using a cut finite element method that utilizes Nitsche’s method to enforce continuity at the interfaces and a consistent stabilization term to handle the cut elements.

Several quantities in the method are conveniently computed in the reference domain where the mappings impose a Riemannian metric. In particular, the stabilization term only involves derivatives in the reference coordinates, which is convenient since it involves higher order derivatives. We develop a quadrature formula for integration on the cut elements that is applicable to a piecewise smooth boundary. We show that the method is stable and we derive optimal order a priori estimates in the energy and $L^{2}$ norms and also present several numerical examples confirming our theoretical results.

Summarizing the key characteristics of the technique and framework for analysis developed in this work are:

• •

The method for computations on multipatch parametric surfaces is based on a fictitious domain method (CutFEM) which does not require the construction of conforming meshes for the trimmed patches.

• •

A stabilization term is added which allows us to perform a complete stability and error analysis independent of how the trim curves cut the computational mesh. In particular, we include proofs of the basic estimates related to the stabilization term in the case of higher order parametric polynomial spaces, as well as an estimate of the condition number of the stiffness matrix.

• •

Both the method and the analysis are adapted to higher order elements and a quadrature rule for integration of cut higher order elements is suggested. We consider standard Lagrange elements but the analysis may also be applied to the spline spaces used in isogeometric analysis.

Outline.

In Section 2 we define the patchwise parametric surface, recall some basic facts on differential operators on surfaces, and formulate our model problem, in Section 3 we construct the patchwise mesh, the finite element spaces, formulate the finite element method, and provide some details on the implementation including a method for quadrature on cut elements, in Section 4 we prove stability of the method, construct an interpolation operator, prove a priori error estimates in the energy and $L^{2}$ norm and prove an upper bound for the stiffness matrix condition number. In Section 5 we present numerical results confirming our theoretical results. Finally, in Section 6 we summarize our results and comment on possible future developments.

2 The Surface and the Laplace–Beltrami Operator

2.1 Piecewise Parametric Description of the Surface

We define the surface and the piecewise parametrization as follows:

• •

Let $\Omega$ be a piecewise smooth connected surface immersed in $\mathbb{R}^{d},d\geq 2$, which is not necessarily orientable.

• •

For all points $x\in\overline{\Omega}$ let $\gamma_{x}(s)=\{y\in\Omega\,:\,\mathrm{dist}(x,y)=s\}$ define a path on $\Omega$ at a fixed distance $s$ to $x$. The length of this path is denoted $|\gamma_{x}(s)|$.

• •

If $\Omega$ has a boundary $\partial\Omega$ it is assumed to be described by a set of smooth curves. Furthermore, for all points $x\in\partial\Omega$ we require $\lim_{s\rightarrow 0+}\frac{|\gamma_{x}(s)|}{s}=\pi$ which means that while the boundary may include kinks, in an intrinsic sense it is smooth.

• •

For all points $x\in\Omega$ we require $\lim_{s\rightarrow 0+}\frac{|\gamma_{x}(s)|}{s}=2\pi$ which means that while the surface may include sharp edges, in an intrinsic sense it is smooth. More concretely the surface can have sharp edges, like the surface of a cylinder, but is not allowed to have corners, like the surface of a cube.

• •

Let $\mathcal{O}=\{\Omega_{i},i\in\mathcal{I}_{\Omega}\}$ be a partition of $\Omega$ into a finite number of smooth subdomains $\Omega_{i}$ which we denote patches.

• •

We assume that each patch boundary $\Gamma_{i}=\partial\Omega_{i}$ is described by a uniformly bounded number of smooth curves.

• •

The interfaces between the patches in $\mathcal{O}$ are described by the curves in the set $\mathcal{G}=\{\Gamma_{ij}=\overline{\Omega}_{i}\cap\overline{\Omega}_{j},(i,j)\in\mathcal{I}_{\Gamma},i<j\}$ where $\mathcal{I}_{\Gamma}$ is a set of pairs of domain indices $i,j\in\mathcal{I}_{\Omega}$ for neighboring patches.

• •

The boundary is described by the curves in the sets $\mathcal{B}_{\partial\Omega_{D}}$ and $\mathcal{B}_{\partial\Omega_{N}}$ for the Dirichlet and Neumann parts of the boundary, respectively.

• •

For each patch $\Omega_{i}$ we associate a diffeomorphism $F_{i}^{-1}:\Omega_{i}\rightarrow\widehat{\Omega}_{i}\subset I^{2}=[0,1]^{2}$ to the reference domain. We also assume that $F_{i}|_{\widehat{\Omega}_{i}}$, is the restriction of a diffeomorphism $F_{i}:U_{\delta}(\widehat{\Omega}_{i})\rightarrow F_{i}U_{\delta}(\widehat{\Omega}_{i})$, where $U_{\delta}(\widehat{\Omega}_{i})=\{{\widehat{x}}\in I^{2}:d({\widehat{x}},\widehat{\Omega}_{i})<\delta\}$ and $d(x,y)=\|x-y\|_{\mathbb{R}^{2}}$ is the usual Euclidean distance function.

• •

To be able to evaluate functions $u\in H^{p+1}(\Omega)$ on $F_{i}U_{\delta}(\widehat{\Omega}_{i})$ we for patches $\widehat{\Omega}_{i}$ with a smooth interface further assume $F_{i}U_{\delta}(\widehat{\Omega}_{i})\subset\Omega$ while we for patches with a sharp interface define an $H^{p+1}$-extension $E^{*}u$ of $u$ on $F_{i}U_{\delta}(\widehat{\Omega}_{i})$ which is possible as $\partial\Omega_{i}$ in the sharp interface case is a closed smooth curve.

• •

Given $x\in\Omega_{i}$ we denote the corresponding point $F_{i}^{-1}x\in{\widehat{\Omega}}_{i}$ by ${\widehat{x}}$ and given a subset $\omega\subset\Omega_{i}$ we let $\widehat{\omega}=F_{i}(\omega)$. For a function $v:F_{i}U_{\delta}(\widehat{\Omega}_{i})\rightarrow\mathbb{R}$ we let $\widehat{v}:U_{\delta}(\widehat{\Omega}_{i})\rightarrow\mathbb{R}$ denote the pullback $\widehat{v}=v\circ F_{i}$ and for a function $\widehat{v}:U_{\delta}(\widehat{\Omega}_{i})\rightarrow\mathbb{R}$ we let $v:F_{i}U_{\delta}(\widehat{\Omega}_{i})\rightarrow\mathbb{R}$ denote the push forward such that $v\circ F_{i}=\widehat{v}$. Note that the pullback is indeed defined on the slightly larger domain $U_{\delta}(\widehat{\Omega}_{i})$. This property will be convenient when we construct an interpolation operator.

This surface description and notation are illustrated in Figure 1.

2.2 Riemannian Metric on the Reference Patches

Tangent Spaces.

At each point ${\widehat{x}}\in\widehat{\Omega}_{i}$ we let $T_{{\widehat{x}}}(\widehat{\Omega}_{i})$ denote the tangent space of $\widehat{\Omega}_{i}$ and we let $\{\widehat{e}_{1},\widehat{e}_{2}\}$ be a fixed orthonormal basis for $T_{{\widehat{x}}}(\widehat{\Omega}_{i})$, i.e. the same basis is used independent of ${\widehat{x}}$. At each point $x\in\Omega_{i}$ we locally define the tangent space to $\Omega_{i}$ as

[TABLE]

where $\widehat{\partial}_{k}$ is the partial derivative in the $\widehat{e}_{k}$ direction and $\widehat{\partial}_{k}F_{i}=\widehat{\partial}_{k}F_{i}|_{x}$ . Thus, any tangent vector $a\in T_{x}(\Omega_{i})$ can be written

[TABLE]

where $\widehat{a}=[\widehat{a}_{1},\widehat{a}_{2}]^{T}\in T_{{\widehat{x}}}(\widehat{\Omega}_{i})$ and we introduced the notation

[TABLE]

Riemannian Metric.

We equip $T_{x}(\Omega_{i})$ with the Euclidean inner product ”$\cdot$” in $\mathbb{R}^{d}$, i.e. the inner product of the immersing space, and we define the induced inner product ${\widehat{g}}_{{\widehat{x}}}$ on $T_{{\widehat{x}}}(\widehat{\Omega}_{i})$ as follows

[TABLE]

Introducing the symmetric positive definite matrix $\widehat{G}$ with components ${\widehat{g}}_{kl}=\widehat{\partial}_{k}F_{i}\cdot\widehat{\partial}_{l}F_{i}$ the inner product ${\widehat{g}}_{{\widehat{x}}}(\widehat{a},\widehat{b})$ can be written

[TABLE]

This inner product is a Riemannian metric on $\widehat{\Omega}_{i}$ and $\widehat{G}$ is called the metric tensor. We denote the norm on $T_{{\widehat{x}}}(\widehat{\Omega}_{i})$ induced by the Riemannian metric by

[TABLE]

We note that $\widehat{D}F_{i}:T_{{\widehat{x}}}(\widehat{\Omega}_{i})\ni\widehat{a}\mapsto a\in T_{x}(\Omega_{i})$ is an isometry since norms and angles are preserved

[TABLE]

and if $\theta$ ($\widehat{\theta}$) is the angle between $a$ and $b$ ($\widehat{a}$* and $\widehat{b}$) we have*

[TABLE]

We note that given $a\in T_{x}(\Omega_{i})$ we find the corresponding $\widehat{a}\in T_{\widehat{x}}(\widehat{\Omega}_{i})$ using the relation

[TABLE]

since we have the identity

[TABLE]

and thus we conclude that $\widehat{G}\widehat{a}=\widehat{D}F_{i}^{T}a$.

Note that we have a uniform bound on the eigenvalues of $\widehat{G}$, i.e. there exists constants $0<c<C$ such that

[TABLE]

and as a consequence

[TABLE]

We will use the notation $T(\omega)=\sqcup_{x\in\omega}T_{x}(\omega)$ for the tangent bundle over a subset $\omega\subset\Omega$, which is the collection of all the tangent vector spaces at the points in $\omega$. When there is no possibility of confusion we will use the simplified notation $g=g_{x}$, ${\widehat{g}}={\widehat{g}}_{\widehat{x}}$ and correspondingly for the induced norms.

2.3 Integration and ${\bm{L}}^{\bm{2}}$ Inner Products

Integration.

The integral over $\omega\subset\Omega_{i}$ is defined by

[TABLE]

where $|\widehat{G}|=|\text{det}(\widehat{G})|$.

Let $\widehat{\gamma}$ be a curve in $U_{\delta}(\widehat{\Omega}_{i})$ and let $(0,l)\ni s\mapsto\widehat{\gamma}(s)\in U_{\delta}(\widehat{\Omega}_{i})$ be an arclength parametrization and $d\widehat{\gamma}$ the arclength measure. We define the integral over the curve $\gamma=F_{i}\circ\widehat{\gamma}\subset\Omega_{i}$ as follows

[TABLE]

where $\widehat{\tau}=\frac{d\widehat{\gamma}}{ds}$ is the unit tangent vector to $\widehat{\gamma}$ with respect to the Euclidean $\mathbb{R}^{2}$ inner product.

Remark 2.1

In Appendix A we provide some more details on the definition of the integrals and also discuss the extension to higher dimensions.

Inner Products.

We let $(\cdot,\cdot)_{\omega}$ denote the $L^{2}(\omega)$ inner product

[TABLE]

and for curves $\gamma$ we have analogous definitions where the integrals and measures are replaced by integrals over the curve and the appropriate measures.

2.4 Differential Operators in Reference Coordinates

Here we introduce some differential operators and formulate Green’s formula. In Appendix A we also derive these expressions including Green’s formula using basic calculus.

The Gradient.

Let $\widehat{\nabla}$ be the tangential gradient operator in reference coordinates

[TABLE]

The tangential gradient $\nabla u\in T_{p}(\Omega_{i})$ is represented in terms of reference coordinates

[TABLE]

Using the chain rule we obtain the identities

[TABLE]

and thus we conclude that the the gradient representation in reference coordinates is

[TABLE]

The Divergence.

The divergence operator div on $\Omega$ is defined by the identity

[TABLE]

for $w\in C^{\infty}_{0}(\omega)$, and may be expressed in reference coordinates as follows

[TABLE]

The Laplace–Beltrami Operator.

We define the Laplace–Beltrami operator $\Delta$ on the surface $\Omega$ by

[TABLE]

which in reference coordinates is given by the identity

[TABLE]

Green’s Formula.

Green’s formula on $\omega\subset\Omega_{i}$ takes the form

[TABLE]

where $n\in T_{p}(\Omega_{i})$ is the exterior unit normal to the curve $\partial\omega$.

2.5 Sobolev Spaces

In the reference coordinates we let $H^{k}(\widehat{\Omega}_{i})$ denote the usual Sobolev spaces of order $k$ with inner product and norm

[TABLE]

On the surface we define the corresponding spaces of functions $v$ that are liftings of functions $\widehat{v}\in H^{k}(\widehat{\Omega}_{i})$,

[TABLE]

with inner product and norm

[TABLE]

We employ standard notation $L^{2}(\omega)=H^{0}(\omega)$ and $\|\cdot\|_{L^{2}(\omega)}=\|\cdot\|_{\omega}$.

2.6 The Laplace–Beltrami Interface Problem

Let $n_{i}\in T_{x}(\Omega_{i})$ be the outward pointing normal to the patch boundary $\partial\Omega_{i}$ as illustrated over a patch interface $\Gamma_{ij}$ in Figure 2. We formulate our model problem for a surface without boundary: Given $f$ such that $(f,1)_{\Omega}=0$, find $u:\Omega\rightarrow\mathbb{R}$, with $(u,1)_{\Omega}=0$, such that

[TABLE]

where the jump operator is defined by

[TABLE]

Remark 2.2

We pose our model problem on a surface without boundary to simplify the presentation in the analysis. However, we comment in Section 3.3 on how the method is easily adapted to boundary conditions and in Section 5 we also present a numerical example with boundary conditions.

Remark 2.3

We use the interface formulation since we have parametric mappings defined on the partition $\mathcal{O}$ of $\Omega$ in contrast to the standard manifold description which is based on a partition of unity and compatibility conditions between the local parametrizations.

Let $V$ be the Hilbert space

[TABLE]

Then we have following weak formulation of (2.28),

[TABLE]

for all $v\in V$. It then follows from the Lax–Milgram lemma that (2.28) has a unique solution in $V$ for $f\in V^{\prime}$, the dual of $V$, such that $\int_{\Omega}f\,dx=0$. Furthermore, we also have the elliptic regularity result

[TABLE]

3 The Finite Element Method

3.1 Construction of the Mesh

Let $\widehat{\mathcal{K}}_{h,i,0}$ be a uniform structured tensor product mesh on the unit square $I^{2}=[0,1]^{2}$ consisting of elements ${\widehat{K}}$ with mesh size $h$. For each patch $\Omega_{i}$ we define the active background mesh in the reference domain as

[TABLE]

and the corresponding mesh on the surface is obtained by

[TABLE]

Let $\widehat{\mathcal{F}}_{h,i}$ be the set of interior faces belonging to elements in $\widehat{\mathcal{K}}_{h,i}$ that intersects the boundary $\partial\widehat{\Omega}_{i}$. The mesh and the set of edges $\widehat{\mathcal{F}}_{h,i}$ are illustrated in Figure 3. Finally, the collection of meshes

[TABLE]

provides a mesh on the surface with cut elements in the vicinity of the interfaces.

3.2 The Finite Element Spaces

Let $\widehat{V}_{h,0}$ be a finite element space on $I^{2}=[0,1]^{2}$ of continuous piecewise tensor product polynomials of order $p$ defined on the mesh $\widehat{\mathcal{K}}_{h,i,0}$. We define the spaces

[TABLE]

and we note that the functions in $V_{h}$ are discontinuous across the interfaces $\Gamma_{ij}\in\mathcal{G}$.

Remark 3.1

Since each patch boundary $\widehat{\Gamma}_{i}=\partial\widehat{\Omega}_{i}$ can arbitrarily cut the mesh $\widehat{\mathcal{K}}_{h,i,0}$ we typically get cut elements on both sides of an interface and on the boundary.

3.3 The Method

Using the formulation in **[8]** we introduce the Nitsche bilinear form $a_{h}$ given by

[TABLE]

where $\beta$ is a positive parameter, and the linear functional $l$ is defined

[TABLE]

The finite element method takes the form: find $u_{h}\in V_{h}/\mathbb{R}$ such that

[TABLE]

where the bilinear form $A_{h}$ is defined

[TABLE]

The stabilization form $j_{h}$ is defined

[TABLE]

where $j_{h,i}$ is short for $j_{h,{\widehat{\mathcal{F}}_{h,i}}}$ which is defined

[TABLE]

where $\widehat{\mathcal{F}}$ is a set of faces, $\{\gamma_{k}\}_{k=1}^{p}$ are positive parameters, and $\widehat{D}^{k}_{n}v$ on a face $\widehat{F}\in\widehat{\mathcal{F}}$ is the $k$:th derivative in the face normal direction to $\widehat{F}$ with respect to the Euclidean $\mathbb{R}^{2}$ inner product. We recall the definition of the jump and define the normal flux average over interfaces

[TABLE]

In the context of unfitted finite elements the stabilization form (3.12) was first analysed for linear elements in **[3]** and extended to higher order elements in **[14]**.

Note that, for $u|_{\Omega_{i}}\in H^{p+1}(\Omega_{i})$, the finite element method (3.9) is consistent and thus the error $u-u_{h}$ satisfies the Galerkin orthogonality

[TABLE]

Remark 3.2** (Penalty Parameter)**

The parameter $\beta$ is chosen large enough as in standard Nitsche type methods and suitable choices of the parameters $\{\gamma_{k}\}_{k=1}^{p}$ are provided in Section 5 below.

Remark 3.3** (Adaptation to Boundary Conditions)**

While formulated above for a surface without boundary the method (3.9) is easily adapted to boundary conditions. For non-homogeneous Dirichlet and Neumann boundary conditions

[TABLE]

where $\partial\Omega_{D}\cup\partial\Omega_{D}=\partial\Omega$ and $\partial\Omega_{D}\cap\partial\Omega_{D}=\emptyset$ we introduce the modified Nitsche forms

[TABLE]

and the resulting method reads: find $u_{h}\in V_{h}$ such that

[TABLE]

Note that for a problem with Dirichlet boundary we no longer need $(u_{h},1)_{\Omega}=(f,1)_{\Omega}=0$.

Remark 3.4** (Adaptation to Convection–Diffusion)**

The method is also easily extended to cover convection–diffusion operators such as

[TABLE]

where $\epsilon$ is a constant and $b|_{\Omega_{i}}=b_{i}:\Omega_{i}\rightarrow T_{p}(\Omega_{i})$ is a tangential vector field. In this case we get some additional terms and the bilinear form reads

[TABLE]

where $\gamma$ is an upwind parameter and $\{\cdot\}$ for the non-flux terms is the usual average $\{w\}=\frac{w_{i}+w_{j}}{2}$. Note that we require the vector field $b$ to be consistent across interfaces in the sense that $n_{i}\cdot b_{i}+n_{j}\cdot b_{j}=0$ on $\Gamma_{ij}$.

3.4 Formulation in Reference Coordinates

In order to assemble the load vector and stiffness matrix we have the following expressions in reference coordinates

[TABLE]

For the interface terms on $\Gamma_{ij}$ we note that $p_{ij}=F_{j}^{-1}F_{i}:\widehat{\Gamma}_{ij}^{i}\rightarrow\widehat{\Gamma}_{ij}^{j}$ is a bijection. Here we introduced the notation $\widehat{\Gamma}_{ij}^{i}=F_{i}^{-1}(\Gamma_{ij})$. We may use $p_{ij}$ to pull back values from $\widehat{\Gamma}_{ij}^{j}$ on $\widehat{\Gamma}_{ij}^{i}$, i.e. fetch values from the reference representation of $\Gamma_{ij}$ on $\widehat{\Omega}_{j}$ based on coordinates in the reference representation of $\Gamma_{ij}$ on $\widehat{\Omega}_{i}$. We obtain the identity

[TABLE]

and for the other interface term we get

[TABLE]

where

[TABLE]

and $\widehat{\nu}_{i}$ is the unit normal to $\widehat{\tau}_{i}$ with respect to the Euclidean $\mathbb{R}^{2}$ inner product in $T(\widehat{\Omega}_{i})$. To verify (3.29) we used the identity (2.19) for the gradient $\widehat{\nabla v}_{i}=\widehat{G}^{-1}\widehat{\nabla}\widehat{v}_{i}$, and the fact that the reference coordinates $\widehat{n}_{i}$ of the normal $n_{i}$ may be expressed in terms of $\widehat{\nu}_{i}$,

[TABLE]

*which follows from the fact that $\frac{\widehat{G}^{-1}\widehat{\nu}}{\|\widehat{G}^{-1}\widehat{\nu}\|_{{\widehat{g}}}}$ is a unit vector with respect to the metric inner product ${\widehat{g}}$ which is also *${\widehat{g}}-$orthogonal to the tangent vector $\widehat{\tau}$,

[TABLE]

3.5 Quadrature

Quadrature on Cut Elements.

To compute the terms implied from the above forms we generate a quadrature scheme for evaluation of integrals in the reference patches on the form

[TABLE]

where the domain of integration is the intersection between a reference patch $\widehat{\Omega}_{i}$ and a finite element ${\widehat{K}}$, and the integrand $f$ stems from tensor product polynomials of arbitrary order. We denote this intersection $\widehat{\omega}=\widehat{\Omega}_{i}\cap{\widehat{K}}$ and assume its boundary can be described by the union of $N$ non-overlapping curves $\widehat{\gamma}_{k}:[0,1]\rightarrow\mathbb{R}^{2}$, i.e.

[TABLE]

where $\widehat{\gamma}_{k}(s)=\begin{bmatrix}\widehat{\gamma}_{k,1}(s)&\widehat{\gamma}_{k,2}(s)\end{bmatrix}^{T}$ is parametrized such that $\|\widehat{\gamma}_{k}^{\prime}(s)\|_{\mathbb{R}^{2}}\neq 0$ and that $\widehat{\gamma}_{k}(s)$ in the positive direction traverses $\partial\widehat{\omega}$ counter-clockwise. Thus, the exterior unit normal to $\widehat{\omega}$, with respect to the $\mathbb{R}^{2}$ inner product, may be expressed

[TABLE]

We define a vector field

[TABLE]

where $a$ is an arbitrary constant which we choose as $a=\min_{\widehat{x}_{2}\in{\widehat{K}}}\widehat{x}_{2}$ and note that we can express the integrand in (3.32) as $f=\widehat{\nabla}\cdot\widehat{\phi}$ by the fundamental theorem of calculus. We rewrite (3.32) as two nested one dimensional integrals, one in each reference coordinate direction, by applying the divergence theorem in $\mathbb{R}^{2}$ and the following calculations

[TABLE]

where we in (3.41) use (3.34) and in the last equality make a change of integration in the inner 1D integral.

Assuming the integrand $f$ is a tensor product polynomial of degree $p_{f}$, i.e. $f\in Q_{p_{f}}$, and that the boundary representation $\widehat{\gamma}_{k}(s)$ in each dimension is a polynomial of degree $p_{\gamma}$, i.e. $\widehat{\gamma}_{k}\in P_{p_{\gamma}}$, we can deduce the following resulting polynomial degrees of the integrands in the inner and outer 1D integral

[TABLE]

and thus we for each reference dimension can choose the number of Gauss points such that the 1D integrals in (3.41) are evaluated exactly. Let $\{s^{(i)},w^{(i)}\}_{i=1}^{n}$ and $\{\tilde{s}^{(j)},\tilde{w}^{(j)}\}_{j=1}^{\tilde{n}}$ be the set of Gauss quadrature points and weights which exactly integrates polynomials of degree $\leq p_{f}$ and degree $\leq 2p_{f}p_{\gamma}+2p_{\gamma}-1$ on $[0,1]$, respectively. Thus, the resulting quadrature points and weights are

[TABLE]

Note that the sum of the quadrature weights $\{\widehat{w}^{(ijk)}\}$ gives the area of $\widehat{\omega}$. An illustration of the resulting quadrature points for an example intersection is given in Figure 4 and in Table 1 we list the polynomial degree and the number of integration points in the nested 1D integrals depending on the tensor product polynomial of the initial integrand $f$ and the boundary polynomial degree of the boundary representation.

Remark 3.5** (Domain Complexity)**

Note that the formulation of the quadrature rule assumes nothing about the complexity of the integration domain other than that its boundary should be well approximated by piecewise polynomial parametrizations. Thus, complex boundaries or holes pose no problem with this quadrature rule and the resolution of the boundary approximation is independent of the size of the finite elements. On the other hand, allowing arbitrarily complex boundaries within an element means we cannot assume a readily available bulk description, for example a mesh, of the intersection between the element and the domain.

Remark 3.6** (Current Implementation)**

As the tensor product polynomials of our finite element basis functions will be perturbed by the Riemannian metric we compensate for this in the quadrature rule by choosing a higher order rule than indicated by the basis functions alone. Also, in cases where the trimmed patches in the reference domain are not exactly represented by piecewise $P_{1}$ curves we in our current implementation choose a $P_{1}$ representation with a resolution high enough for this error to be negligible. This use of $P_{1}$ representation is however not a limitation of the quadrature rule as seen in the above derivation and an alternative would be to use higher order approximations of the patch boundaries instead.

Remark 3.7** (Negative Quadrature Weights)**

As seen in Figure 4 the quadrature method includes both positive and negative weights which stems from adding and subtracting various parts of the integration domain. This is an undesirable property when considering reduced quadrature as inexact cancellation possibly could lead to loss of coercivity. A possible modification which improves the method in this regard is to replace the constant lower bound $a$ in the integral in (3.35) by a polynomial of degree $\leq p_{\gamma}$ where the polynomial coefficients are chosen such that the number of negative quadrature weights are minimized. If no restriction is placed on the polynomial coefficients, this could lead to some quadrature points being placed slightly outside the element.

However, in the present work we do further not investigate the aspect of reduced integration. We view this as a reference quadrature rule capable of integrating higher order tensor product polynomials and as noted in Remark 3.6 we rather use an increased integration order. In a complicated real world setting it is therefore advisable to chose an alternative quadrature scheme where positive quadrature weights can be guaranteed.

Quadrature on Interfaces.

To compute the interface terms (3.26) and (3.28) we construct a partition of $\widehat{\Gamma}^{i}_{ij}$ which contains both all the intersection points between the curve $\widehat{\Gamma}^{i}_{ij}$ and the mesh $\widehat{\mathcal{K}}_{h,i}$ as well as all the intersection points between the curve $\widehat{\Gamma}^{j}_{ij}$ and the mesh $\widehat{\mathcal{K}}_{h,j}$ mapped back to $\widehat{\Gamma}^{i}_{ij}$ using the mapping $p_{ij}^{-1}=F_{i}^{-1}F_{j}$. Each interval in the partition of $\widehat{\Gamma}^{i}_{ij}$ will thus be associated only with a single element in $\widehat{\mathcal{K}}_{h,i}$ and a single element in $\widehat{\mathcal{K}}_{h,j}$ and we apply a 1D Gauss quadrature rule on each interval. See Figure 5 for an illustration of the partition of the interface.

4 A Priori Error Estimates

Let $a\lesssim b$ denote $a\leq Cb$ with a constant $C$ independent of the mesh parameter $h$.

4.1 Norms

Given a set of faces $\widehat{\mathcal{F}}$ in a mesh let

[TABLE]

We define the following energy norm

[TABLE]

where

[TABLE]

4.2 Inverse Inequalities

On elements which are partially outside the patch domain, as illustrated in Figure 6, we use the following inverse inequality to control a discrete function or its gradient on an element in terms of the gradient on a neighboring element and a suitable face term.

We will below make repeated use of the set of elements cut by the patch boundary $\Gamma_{i}$ and thus we define the set

[TABLE]

and analogously we define $\widehat{\mathcal{K}}_{h,i}(\widehat{\Gamma}_{i})$ in the reference domain.

Lemma 4.1

Let the two elements ${\widehat{K}}_{1},{\widehat{K}}_{2}\in\widehat{\mathcal{K}}_{h,i}$ be neighbors of face $\widehat{F}$ with face normal $n_{\widehat{F}}$. For all $\widehat{v}\in{\widehat{V}}_{h,i}|_{{\widehat{K}}_{1}\cup{\widehat{K}}_{2}}$ the following estimates then hold

[TABLE]

Proof.* We begin by proving estimate (4.5) and we then make use of calculations in this proof when proving the second estimate (4.6).*

Proof of Estimate (4.5).

Let $\widehat{v}_{l}=\widehat{v}|_{{{\widehat{K}}}_{l}}$, $l=1,2.$ Since ${\widehat{v}}_{2}$ is a polynomial we may evaluate $\widehat{v}_{2}$ also on ${\widehat{K}}_{1}$. Using the triangle inequality followed by the inverse estimate $\|\widehat{v}_{2}\|_{{\widehat{K}}_{1}}\lesssim\|\widehat{v}_{2}\|_{{\widehat{K}}_{2}}$, we obtain

[TABLE]

To estimate the second term we note that using Taylor’s formula on $(\widehat{v}_{1}-\widehat{v}_{2})$ at $x_{\widehat{F}}\in{\widehat{F}}$ in the face normal direction $n_{\widehat{F}}$ (with respect to the Euclidean $\mathbb{R}^{2}$ inner product) gives

[TABLE]

for $(x-x_{F})=sn_{\widehat{F}}$ with $s\in\mathbb{R}$. Using an orthonormal coordinate system $\{\widehat{e}_{1},\widehat{e}_{2}\}$ with $\widehat{e}_{1}=n_{\widehat{F}}$, i.e. a coordinate system which is aligned with the face normal and the face itself, we have $x_{\widehat{F}}=[x_{\widehat{F},1},x_{2}]^{T}$ and $n_{\widehat{F}}\cdot(x-x_{\widehat{F}})=s=x_{1}-x_{\widehat{F},1}$. This gives the following expressions for the partial derivatives

[TABLE]

Using the Cauchy–Schwarz inequality for sums we obtain

[TABLE]

and in the same way we have

[TABLE]

where we used an inverse inequality in the last step to remove $\widehat{\partial}_{2}$. In summary, we have

[TABLE]

which concludes the proof of estimate (4.5).

Proof of Estimate (4.6).

Starting in the same way as the proof of (4.5) but without the gradients we have

[TABLE]

As $(\widehat{v}_{1}-\widehat{v}_{2})|_{\widehat{F}}=0$ the Poincaré inequality holds yielding the following estimate

[TABLE]

*and we handle the remaining term as in the proof of (4.5).

---

**

Assumption 4.1** (Patch Geometry)**

For a given element $K$ let $\mathcal{N}_{0}(K)=K$ and for $l=1,2,\dots,$ let $\mathcal{N}_{l}(K)$ be the union of all elements that share a face or a node with an element in $\mathcal{N}_{l-1}(K)$, in other words $\mathcal{N}_{l}(K)$ is the set of elements that are neighbors of distance less or equal to $l$. Assume that there is a positive integer $l$, a maximum mesh parameter $0<h_{0}$, and a positive constant $c$, such that for all $h\in(0,h_{0}]$ and all $K\in\mathcal{K}_{h}(\Gamma)$ there is an element $K^{\prime}\in\mathcal{N}_{l}(K)$ such that $|K^{\prime}\cap\widehat{\Omega}_{i}|\geq c|K|$.

Remark 4.1

This assumption limits the complexity of the reference subdomains $\widehat{\Omega}_{i}$. Note that the assumptions holds for $0<h\leq h_{0}$ with $h_{0}$ small enough when the boundary $\partial\widehat{\Omega}_{i}$ satisfies a cone condition. The assumption does not hold for instance in the vicinity of a cusp. In future work we will return to situations with more general patches including very thin patches and patches with cusps since such patches may occur in CAD models used in practical engineering design.

Lemma 4.2

Given Assumption 4.1 it for ${\widehat{K}}\in\widehat{\mathcal{K}}_{h,i}(\widehat{\Gamma}_{i})$, holds

[TABLE]

where $\mathcal{F}(\mathcal{N}_{l}({\widehat{K}}))$ is the set of interior faces in the neighborhood $\mathcal{N}_{l}({\widehat{K}})$.

Proof.* This estimate follows directly from repeated use of Lemma 4.1 together with the Assumption on Patch Geometry that there is ${\widehat{K}}^{\prime}\in\mathcal{N}_{l}({\widehat{K}})$ such that $|\mathcal{N}_{l}({\widehat{K}})\cap\widehat{\Omega}_{i}|\gtrsim h^{2}$.

---

**

Lemma 4.3

There is a constant such that for all $v\in V_{h,i}$ it holds

[TABLE]

Proof.* We proceed as follows*

[TABLE]

where in (4.22) we used the Cauchy–Schwarz inequality; in (4.23) we divided the integral into element contributions; in (4.24) we mapped to reference coordinates and used the bound

[TABLE]

In (4.25) we used the following inverse trace inequality

[TABLE]

which holds for $\widehat{v}\in\widehat{V}_{h,i}|_{\widehat{K}}$ and we verify below; in (4.26) we used estimate (4.20); in (4.27) we used the bound

[TABLE]

for each of the elements in $\mathcal{N}_{l}({\widehat{K}})$; and finally in (4.28) we used (4.20).

Verification of (4.33).

We have the bounds

[TABLE]

*where we used the fact that the length $\widehat{\Gamma}_{i}\cap{\widehat{K}}$ of the curve segment $|\widehat{\Gamma}_{i}\cap{\widehat{K}}|\lesssim h$ for $h\in(0,h_{0}]$, with $h_{0}$ small enough, which holds since $\Gamma_{i}$ consists of a finite set of smooth curve segments, and at last we used an inverse bound to estimate the $L^{\infty}$ norm in terms of the $L^{2}$ norm.

---

**

4.3 Coercivity and Continuity

Lemma 4.4** (Coercivity and Continuity)**

*The form $A_{h}$ satisfies:

(i) For large enough $\beta$ it holds for all $v\in V_{h}$,*

[TABLE]

(ii) For all $v,w\in V_{h}+H^{3/2+\epsilon}(\Omega)$ it holds

[TABLE]

Proof.* The first statement (i) follows directly from the inverse inequality (4.21) together with standard arguments, see for instance [13]. The second statement (ii) follows directly from the Cauchy–Schwarz inequality.

---

**

4.4 Interpolation

Let $\widehat{\pi}_{h,i}:L^{2}(\widehat{\mathcal{K}}_{h,i})\rightarrow\widehat{V}_{h,i}$ be the Scott–Zhang interpolation operator, see **[20]**. We recall the standard interpolation error estimate

[TABLE]

where $\mathcal{N}({\widehat{K}})\subset\widehat{\mathcal{K}}_{h,i}$ is the set of elements in $\widehat{\mathcal{K}}_{h,i}$ that are neighbors to ${\widehat{K}}$. Next we note that for $0<h\leq h_{0}$, with $h_{0}$ small enough, we have

[TABLE]

where $U_{\delta}(\widehat{\Omega}_{i})=\{x\in I^{2}\,|\,d(x,\widehat{\Omega}_{i})<\delta\}$, see Section 2.1. We define the global interpolation operator $\pi_{h}:L^{2}(\Omega)\rightarrow V_{h}$ as follows

[TABLE]

where we used the fact that $\widehat{v}_{i}$ is defined on $U_{\delta}(\widehat{\Omega}_{i})$ and therefore the right hand side of (4.42) is well defined due to (4.41).

Lemma 4.5** (Interpolation Error Estimate)**

The interpolation operator $\pi_{h}$ defined by (4.42) satisfies

[TABLE]

Proof.* Let $\rho=u-\pi_{h}u$. For each $\widehat{F}\in\widehat{\mathcal{F}}_{h,i}$ there are two neighboring elements $\widehat{K}_{1}$ and $\widehat{K}_{2}$ in $\widehat{\mathcal{K}}_{h,i}$. Using the triangle inequality followed by the trace inequality*

[TABLE]

and the interpolation estimate (4.40) we obtain

[TABLE]

Thus we conclude that

[TABLE]

Next we recall the following trace inequality

[TABLE]

which holds independent of the position of $\widehat{\Gamma}_{i}$ in ${\widehat{K}}$, see **[7]**. Introducing the notation

[TABLE]

we may estimate the remaining terms in the energy norm (4.2) using the triangle inequality in (4.52), the trace inequality (4.50) in (4.54) and the interpolation estimate (4.40) in (4.55) as follows

[TABLE]

*which concludes the proof.

---

**

4.5 Error Estimates

Theorem 4.1** (Energy Error Estimate)**

Let $u\in H^{p+1}(\Omega)$ be the solution to (2.28) and $u_{h}$ the solution to (3.9), then

[TABLE]

Proof.* Adding and subtracting the interpolant we have*

[TABLE]

where we used the interpolation error estimate (4.43). For the second term we have, using the notation $e_{h}=\pi_{h}u-u_{h}$ and Galerkin orthogonality (3.14),

[TABLE]

Thus we conclude that

[TABLE]

*where we finally used the interpolation estimate (4.43) again. Together (4.58) and (4.62) concludes the proof.

---

**

Theorem 4.2** (${\bm{L}}^{\bm{2}}$ Error Estimate)**

Let $u$ be the solution to (2.28) and $u_{h}$ the solution to (3.9), then

[TABLE]

Proof.* Let $e=u-u_{h}$ be the error, and $\phi$ the solution of the dual problem*

[TABLE]

Recall that $\int_{\Omega}u\,dx=\int_{\Omega}u_{h}\,dx=0\,dx$ and thus $\int_{\Omega}e=0\,dx$ and we conclude that the dual problem has a unique solution in $V$, see (2.30), that satisfies

[TABLE]

Multiplying the dual problem by $e$, integrating by parts on each subdomain $\Omega_{i}$, and using the interface conditions on $\phi$ we obtain

[TABLE]

where we in (4.68) used the Galerkin orthogonality (3.14) to subtract $A_{h}(e,\pi_{h}\phi)=0$, and we in (4.69) applied the Cauchy–Schwarz inequality and also introduced the norms $\|v\|_{a_{h}}^{2}=a_{h}(v,v)$ and $\|v\|_{j_{h}}^{2}=j_{h}(v,v)$ induced by their respective forms. In the final inequality we applied the energy norm estimate Theorem 4.1.

Estimate of $\bm{I}$.

Using the triangle inequality on the jumps and averages, applying the trace inequality and standard interpolation estimates we obtain

[TABLE]

where we in the last inequality use the elliptic regularity of the dual solution (4.65).

Estimate of $\bm{I}\bm{I}$.

Consider two neighboring elements ${\widehat{K}}_{1},{\widehat{K}}_{2}\in\widehat{\mathcal{K}}_{h,i}$ sharing face $\widehat{F}$. Introducing a patchwise interpolant $\widehat{\pi}_{F}:L^{2}({\widehat{K}}_{1}\cup{\widehat{K}}_{2})\rightarrow Q^{p}({\widehat{K}}_{1}\cup{\widehat{K}}_{2})$, where $Q^{p}$ is the space of tensor product polynomials of degree $\leq p$, we note that we can subtract $\widehat{\pi}_{F}\widehat{\phi}$ inside the stabilization terms as it will give no contribution due to the jump over faces. We proceed as follows

[TABLE]

where in (4.73) we used an inverse trace inequality, in (4.74) we used an inverse inequality, in (4.74) we added and subtracted $\widehat{\phi}$ and used the triangle inequality, and finally in (4.76) we used interpolation estimates. We thus have the estimate

[TABLE]

*which concludes the proof.

---

**

4.6 Condition Number Estimate

To prove an upper bound on the stiffness matrix condition number we follow the approach in **[3, 5]**. Let $\{\varphi_{i}\}_{i=1}^{N}$ be the standard piecewise tensor product polynomial Lagrange basis functions associated with the nodes in $\mathcal{K}_{h}$ and let $\mathcal{A}$ and $\mathcal{M}$ be the stiffness and mass matrices with elements $\mathcal{A}_{ij}=A_{h}(\varphi_{j},\varphi_{i})$ and $\mathcal{M}_{ij}=(\varphi_{j},\varphi_{i})_{\Omega}$, respectively. The condition number for the stiffness matrix is defined by

[TABLE]

where $|\cdot|_{\mathbb{R}^{N}}$ on matrices denotes the operator norm

[TABLE]

and $|\cdot|_{\mathbb{R}^{N}}$ on vectors denotes the Eucledian norm.

Theorem 4.3** (Upper Bound on Condition Number)**

The condition number of the stiffness matrix $\mathcal{A}$ satisfies the estimate

[TABLE]

for all $h\in(0,h_{0}]$ with $h_{0}$ sufficiently small.

Proof.* If $v=\sum_{i=1}^{N}\mathbf{V}_{i}\varphi_{i}$ and $\{\varphi_{i}\}_{i=1}^{N}$ is the usual nodal basis on $\mathcal{K}_{h}$ the following well known estimate holds*

[TABLE]

We will make use of the inverse inequality

[TABLE]

and the discrete Poincaré inequality

[TABLE]

The inverse inequality (4.82) is proven by first applying the triangle inequality to all jump and average terms and then using Lemma 4.3 on the consistency terms, the inverse inequality

[TABLE]

on the jump penalty term and the inverse inequality

[TABLE]

on the stability terms. The discrete Poincaré inequality (4.83) is proven by first separating the cut elements and applying Lemma 4.2 which gives

[TABLE]

where we in the second last inequality apply the standard Poincaré inequality on the first term.

We now turn to estimating $|\mathcal{A}|_{\mathbb{R}^{N}}$ and $|\mathcal{A}^{-1}|_{\mathbb{R}^{N}}$ separately. The product of these estimates will give a bound on the condition number $\kappa(\mathcal{A})$ by its definition (4.78).

Estimate of $|\mathcal{A}|_{\mathbb{R}^{N}}$.

Let $\mathbf{V}\in\widetilde{\mathbb{R}^{N}}$ where $\widetilde{\mathbb{R}^{N}}=\{\mathbf{V}\in\mathbb{R}^{N}\,:\,\mathbf{V}^{T}\mathcal{M}\mathbf{V}=0\}$. In other words $\widetilde{\mathbb{R}^{N}}$ is the space of coefficient vectors corresponding for discrete functions in $V_{h}/\mathbb{R}$. By the definition of the method (in matrix form) and using continuity (4.39) we have

[TABLE]

where we in the last inequality used the inverse estimate (4.82) and (4.81). It follows that

[TABLE]

Estimate of $|\mathcal{A}^{-1}|_{\mathbb{R}^{N}}$.

Let $\mathbf{V}\in\widetilde{\mathbb{R}^{N}}$. Using (4.81), the Poincaré inequality (4.83), coercivity (4.38) and the Cauchy–Schwarz inequality we obtain

[TABLE]

Since $\mathbf{V}$ is arbitrary we by choosing $\mathbf{U}=\mathcal{A}\mathbf{V}\in\widetilde{\mathbb{R}^{N}}$, i.e. $\widetilde{\mathbb{R}^{N}}\ni\mathbf{V}=\mathcal{A}^{-1}\mathbf{U}$, get

[TABLE]

*which in combination with (4.92) concludes the proof.

---

**

5 Numerical Results

In this section we present our numerical experiments to verify convergence rates and the stability of the cut finite element method on patchwise parametrized surfaces. We also provide various numerical examples.

5.1 Model Problems

For our convergence and stability results we choose the same Laplace–Beltrami model problems as in **[17]**; a problem on the unit sphere and a problem on a torus surface. The solutions and load functions to these problems satisfy $(u,1)_{\Omega}=(f,1)_{\Omega}=0$.

Surface and Analytical Solution.

The surfaces and analytical solutions for our two model problems are illustrated in Figure 7 and described below.

• •

Sphere:* The surface $\Omega$ is the unit sphere centered in origo and we use a manufactured problem with analytical solution $u=3x^{2}y-y^{3}$.*

• •

Torus:* The surface $\Omega$ is a torus with inner radius $r=0.6$ and outer radius $R=1$. This surface can be expressed in Cartesian coordinates as the points*

[TABLE]

for $0\leq\theta\leq 2\pi$ and $0\leq\phi<2\pi$ where $\{\theta,\phi\}$ are toroidal coordinates of the surface. We use a manufactured problem with analytical solution $u=\sin(3\phi)\cos(3\theta+\phi)$.

Patchwise Surface Description.

In the presented method the surface $\Omega$ is described by a set of mappings $\{F_{i}\}_{i\in\mathcal{I}_{\Omega}}$ and trimmed patches in reference coordinates $\{\widehat{\Omega}_{i}\}_{i\in\mathcal{I}_{\Omega}}$ such that $\Omega=\bigcup_{i\in\mathcal{I}_{\Omega}}F_{i}(\widehat{\Omega}_{i})$ and $\bigcap_{i\in\mathcal{I}_{\Omega}}F_{i}(\widehat{\Omega}_{i})=\emptyset$. To construct such a description for the two model problems we first create a closed surface approximation of $\Omega$ consisting of a number of polygons $\{T_{i}\}_{i\in\mathcal{I}_{\Omega}}$. For each polygon $T_{i}$ we by a simple affine mapping can create an inverse mapping down to a reference patch $\widehat{\Omega}_{i}$ in $[0,1]^{2}$. To map onto the surface $\Omega$ we from $T_{i}$ use a closest point mapping and by combining the inverse mapping and the closest point mapping we define $F_{i}:\widehat{\Omega}_{i}\rightarrow\Omega_{i}$. The actual patchwise descriptions used for the model problem are illustrated in Figure 8.

5.2 Implementation Aspects

We use tensor product Lagrange finite elements of order $p$ on quadrilaterals in our implementation. In the results below we for the Nitsche interface terms used the parameter $\beta=100$ and for the CutFEM stability terms used parameters $\gamma_{k}=10^{-2}$, $k=1,\dots,p$. The latter choice is numerically investigated in Section 5.4 below. To impose the average constraint $(u_{h},1)_{\Omega}=0$ we use a Lagrange multiplier approach, see for example **[13]**.

5.3 Convergence

To confirm our theoretical results in Theorem 4.1 and Theorem 4.2 we present convergence results for the energy norm error and $L^{2}$ norm error in Figure 10 and Figure 11, respectively. Example numerical solutions for the two model problems are displayed in Figure 9. In these studies the geometry representation, i.e. the reference patches $\widehat{\Omega}_{i}$ and mappings $F_{i}$, is kept fixed while the background grid is refined.

5.4 Stability

Patch Position in the Background Mesh.

Depending on how a reference patch $\widehat{\Omega}_{i}$ is positioned in $[0,1]^{2}$ the intersection with the background mesh may produce situations with arbitrary small cut elements. To demonstrate the stability of the method with regard to different cut situations we produce statistical data by randomly rotating each reference patch $\widehat{\Omega}_{i}$ in the background mesh to give random cut situations and repeating the simulation N times. The standard deviation of the energy and $L^{2}$ errors in these simulations are presented in Figure 12.

Condition Number.

The discrete problem can become arbitrarily ill conditioned if the stabilization term $j_{h,i}$ is not included in the form $A_{h}$. To capture this instability we estimate the condition number of the stiffness matrix for numerous patch positions in the reference domain, producing different cut situations. This is illustrated in Figure 13(a) where we estimate the condition number for both the stabilized and unstabilized system in random cut situations. Note that the condition number for the stabilized stiffness matrix scales as $\mathcal{O}(h^{-2})$ in agreement with the bound proven in Theorem 4.3.

Choice of Stability Parameter $\bm{\gamma}$.

For a fixed mesh size $h$ we investigate how the size of the stabilization parameters $\gamma_{k}$, $k=1,\dots,p$, affect numerical stability, i.e. the condition number $\kappa$, and the size of the error in the solution. Assuming all stabilization parameters take on the same value, i.e. $\gamma_{k}=\gamma$, we present a numerical study of this in Figure 13(b) respectively in Figure 14. For small values of $\gamma$ we increasing condition numbers resulting in numerical instabilities, see Figure 13(b), and for large values of $\gamma$ we note that the stabilization term $j_{h,i}$ will start to impact the solution leading to larger errors, see Figure 14. A good middle ground seems to be $\gamma=10^{-2}$.

5.5 Numerical Examples

Surface with Boundary.

With the simple adaption of the method to boundary conditions described in Section 3.3 we give results of a problem with boundary in Figure 15(a), where we have both non-homogeneous Dirichlet conditions and Neumann conditions.

Klein Bottle.

The Klein bottle is a closed non-orientable surface for which it exists no embedding in $\mathbb{R}^{3}$. Let $\Omega$ in Cartesian coordinates be described by the parametrization

[TABLE]

for $0\leq\theta<\pi$ and $0\leq\phi<2\pi$. We manufacture a problem with the analytical solution $u=3\cos^{2}\theta\sin\phi-\sin^{3}\phi$ and the resulting finite element solution is presented in Figure 15(b).

Surface with Sharp Interfaces.

Let $\Omega$ be the closed surface to the half solid torus defined via (5.1) and $\theta\in[0,2\pi],\,\phi\in[\frac{\pi}{2},\frac{3\pi}{2}]$ and $r\leq 0.6$. The resulting geometry consists of half a torus and two circular discs. We manufacture a problem by choosing the same analytical solution as for the torus model problem on the torus part and on the discs we make the ansatz of a single cubic Hermite polynomial in the radial direction with zero solution and radial derivative in the disc center. The analytical solution on the discs are then derived from the interface conditions. In Figure 16 the solution and gradient magnitude of the finite element solution are displayed, and both flow nicely over the interfaces.

6 Summary and Future Work

We have presented and analysed a higher order cut finite element method for elliptic problems on multipatch surfaces. The method has the following fundamental features:

• •

Patches are described by mappings from a reference domain and trim curves.

• •

On each patch a mesh is constructed using structured grids in the reference domain.

• •

The discrete solution is coupled between the patchwise meshes by enforcing interface conditions using Nitsche’s method.

• •

On each patch we handle elements cut by trim curves by adding certain stabilization terms.

• •

The stability and error analysis is independent of how the trim curves cut the mesh.

Real Applications.

While we in this work consider the Laplace–Beltrami operator as a model problem, there are many real problems posed on surfaces to which the same framework for dealing with multipatch surfaces effectively could be applied. For example, there is a great interest in structural mechanics for modeling membranes, plates and shells, and modeling of thin films and lubrication also occur on surfaces.

Extended Analysis.

In the analysis we assume that, at the interface, the trim curves on both patches map exactly onto the same interface curve. However, this is typically not the case when working with geometries extracted from CAD due to the discrete representation of the trim curves. Therefore a useful extension of the analysis would be to consider gaps in the geometry. Another useful extension would be higher order PDE which are common for problems on surfaces and in this setting we can easily construct a tensor product basis with higher order continuity properties, i.e. where the restriction of the finite element space to each patch is a subspace to the proper Hilbert space.

Isogeometric Analysis.

As the presented multipatch method is based on an exact description of the geometry by parametric mappings and features higher order elements it fits perfectly into the framework of isogeometric analysis **[9, 4]**. The CutFEM approach also allows for convenient construction of structured meshes equipped with tensor product spline basis functions.

Appendix A Differential Operators on Surfaces

We provide details for the local forms of the divergence and Laplace–Beltrami operator as well as a derivation of Green’s formula on a surface using only calculus in the reference coordinates and some basic linear algebra.

Divergence.

Starting from the definition

[TABLE]

for $v,w\in C^{\infty}_{0}(\omega)$, for $\omega\subset\Omega_{i}$, we have the identities

[TABLE]

where we integrated using standard Green’s formula in local coordinates. Thus we conclude that

[TABLE]

The Laplace–Beltrami Operator.

Using the definition of the Laplace–Beltrami operator (2.22) we conclude that in local coordinates

[TABLE]

where we used (2.19).

Green’s Formula.

We have the identities

[TABLE]

where we will now verify the identity $\bigstar=\int_{\partial\omega}n\cdot\nabla v\,w\,d\gamma$. We note that $\widehat{\nu}$ is the unit exterior normal to $\partial\omega$, with respect to the usual $\mathbb{R}^{2}$ inner product in $T(\widehat{\Omega}_{i})$, and that

[TABLE]

see (3.30). Thus the term $\bigstar$ on the right hand side in (A.14) may be written in the form

[TABLE]

and we will next study the measure $\|\widehat{G}^{-1}\widehat{\nu}\|_{g}|\widehat{G}|^{1/2}d\widehat{\gamma}$ in more detail.

The Two-Dimensional Case.

Using the identity

[TABLE]

where

[TABLE]

we have $\widehat{\tau}=S\widehat{\nu}$, where $\widehat{\tau}$ is the unit with respect to the Euclidean inner product tangent vector to the curve $\partial\widehat{\omega}$, and

[TABLE]

and therefore we recover the standard curve measure. Thus we conclude that

[TABLE]

which together with (A.14) concludes the derivation of Green’s formula on the surface $\omega$.

The General Case.

A more general approach which also holds in higher dimension is to use the identity

[TABLE]

where $A$ is a square $n\times n$ matrix and $a$ and $b$ are $n$ vectors. We then obtain

[TABLE]

Now we may chose an orthonormal basis in $\mathbb{R}^{n}$ consisting of $\widehat{\nu}$ and $n-1$ tangent vectors $\{t_{i}\}_{i=2}^{n}$. Let $P=I-\widehat{\nu}\otimes\widehat{\nu}$ be the projection onto the tangent plane and then we have the identity

[TABLE]

where $|P\widehat{G}P|$ is the $n-1$ determinant of the tangent part $P\widehat{G}P$ of $\widehat{G}$. This follows directly from the fact that in normal-tangent coordinates $\widehat{G}+\widehat{\nu}\otimes\widehat{\nu}$ takes the form

[TABLE]

with $\widetilde{G}$ the matrix representation of $\widehat{G}$ in a tangent-normal coordinate system. We may thus conclude that

[TABLE]

which is the appropriate measure on $\partial\omega$. Note also that in the two dimensional case $P\widehat{G}P$ has rank one and the determinant equals the absolute value of the scalar $\widehat{\tau}\cdot\widehat{G}\cdot\widehat{\tau}$ and thus

[TABLE]

which is consistent with the definition of the measure based on arclength measure.