A Nitsche Method for Elliptic Problems on Composite Surfaces

TL;DR

This paper introduces a Nitsche finite element method for solving elliptic PDEs on complex composite surfaces composed of multiple intersecting surfaces, providing stability and error estimates.

Contribution

It develops a novel Nitsche-based finite element approach for elliptic problems on composite surfaces with interfaces, including stability analysis and numerical implementations.

Findings

Method achieves stability under certain geometric assumptions

Error estimates are derived for exact geometry representation

Numerical examples demonstrate method effectiveness

Abstract

We develop a finite element method for elliptic partial differential equations on so called composite surfaces that are built up out of a finite number of surfaces with boundaries that fit together nicely in the sense that the intersection between any two surfaces in the composite surface is either empty, a point, or a curve segment, called an interface curve. Note that several surfaces can intersect along the same interface curve. On the composite surface we consider a broken finite element space which consists of a continuous finite element space at each subsurface without continuity requirements across the interface curves. We derive a Nitsche type formulation in this general setting and by assuming only that a certain inverse inequality and an approximation property hold we can derive stability and error estimates in the case when the geometry is exactly represented. We discuss several different realizations, including so called cut meshes, of the method. Finally, we present numerical examples.

Figures

Equations

$$-1<C\le {\nu }^{+}(x)\cdot {\nu }^{-}(x)$$

$$\Omega =\bigcup _{i\in {\mathcal{I}}_{\Omega }}{\Omega }_i$$

$${\Omega }_i\cap {\Omega }_j=\partial {\Omega }_i\cap \partial {\Omega }_j={\Gamma }_{ij}$$

$$\bigcup _{k\in {\mathcal{I}}_{{\Gamma }_I}}{\Gamma }_k=\bigcup _{i,j\in {\mathcal{I}}_{\Omega }}{\Gamma }_{ij}$$

$$0<c\le \mu (x)\le C$$

$$\sigma (u)=\mu {\nabla }_{\Omega }u$$

$$-{\nabla }_{\Omega }\cdot \sigma (u)$$

$$\sum _{k\in {\mathcal{I}}_{\Omega }(j)}{\nu }_k\cdot {\sigma }_k(u_k)$$

$$u_k$$

$$\nu \cdot \sigma (u)$$

$$u$$

$$V_{g_D}(\Omega )=\left\{v:\Omega \to \mathbb{R}:\begin{array}{rlrl}& v_i\in H^1({\Omega }_i),& & \forall i\in {\mathcal{I}}_{\Omega }\\ & v_k=v_l\text{ on Γj},& & k,l\in {\mathcal{I}}_{\Omega }(j),\forall j\in {\mathcal{I}}_{{\Gamma }_I}\\ & v=g_D\text{ on ΓD}& & \end{array}\right\}$$

$$(v,w{)}_{H^1(\Omega )}=\sum _{i\in {\mathcal{I}}_{\Omega }}({\nabla }_{{\Omega }_i}{v}_i,{\nabla }_{{\Omega }_i}{v}_i{)}_{{\Omega }_i}+(v,w{)}_{{\Omega }_i},\Vert v{\Vert }_{H^1(\Omega )}^2=\sum _{i\in {\mathcal{I}}_{\Omega }}\Vert v_i{\Vert }_{H^1({\Omega }_i)}^2$$

$$(f,v{)}_{\Omega }$$

$$=\sum _{i\in {\mathcal{I}}_{\Omega }}(-{\nabla }_{\Omega }\cdot \sigma (u),v{)}_{{\Omega }_i}$$

$$=\sum _{i\in {\mathcal{I}}_{\Omega }}(\sigma (u),{\nabla }_{\Omega }v{)}_{{\Omega }_i}-(\nu \cdot \sigma (u),v{)}_{\partial {\Omega }_i}$$

$$=\sum _{i\in {\mathcal{I}}_{\Omega }}(\sigma (u),{\nabla }_{\Omega }v{)}_{{\Omega }_i}-(\nu \cdot \sigma (u),v{)}_{{\Gamma }_N}-\sum _{j\in {\mathcal{I}}_{{\Gamma }_I}}\sum _{k\in {\mathcal{I}}_{\Omega }(j)}({\nu }_k\cdot {\sigma }_k(u_k),v_k{)}_{{\Gamma }_j}$$

$$=\sum _{i\in {\mathcal{I}}_{\Omega }}(\sigma (u),{\nabla }_{\Omega }v{)}_{{\Omega }_i}-(g_N,v{)}_{{\Gamma }_N}$$

$$\sum _{k\in {\mathcal{I}}_{\Omega }(j)}({\nu }_k\cdot {\sigma }_k(u_k),v_k{)}_{{\Gamma }_j}={\left(\sum _{k\in {\mathcal{I}}_{\Omega }(j)}{\nu }_k\cdot {\sigma }_k(u_k),v\right)}_{{\Gamma }_j}=0$$

$$a(u,v)=l(v),\forall v\in V_0$$

$$a(u,v)=\sum _{i\in {\mathcal{I}}_{\Omega }}(\sigma (u),{\nabla }_{\Omega }v{)}_{{\Omega }_i},l(v)=\sum _{i\in {\mathcal{I}}_{\Omega }}(f,v{)}_{{\Omega }_i}+(g_N,v{)}_{{\Gamma }_N}$$

$$\sum _{i\in \mathcal{I}({\Omega }_i)}\Vert u{\Vert }_{H^{{\eta }_i}({\Omega }_i)}^2\lesssim \sum _{i\in \mathcal{I}({\Omega }_i)}\Vert f{\Vert }_{H^{{\eta }_i-2}({\Omega }_i)\cap L^2({\Omega }_i)}^2$$

$${\overline{\Omega }}_i\subseteq {\cup }_{K\in {\mathcal{K}}_{h,i}}\overline{K}$$

$$\text{int}({\Omega }_i)\cap \text{int}(K)\ne \varnothing \forall K\in {\mathcal{K}}_{h,i}$$

$$V_h=\underset{i\in {\mathcal{I}}_{\Omega }}{⨁}{V}_{h,i}$$

$$\Vert v-{\pi }_{h,i}v{\Vert }_{L^2({\mathcal{K}}_{h,i})}+h\Vert {\nabla }_{\Omega }(v-{\pi }_{h,i}v){\Vert }_{L^2({\mathcal{K}}_{h,i})}\lesssim h^s\Vert v{\Vert }_{H^s({\Omega }_i)}1\le s\le p+1$$

$$(f,v{)}_{\Omega }$$

$$=\sum _{i\in {\mathcal{I}}_{\Omega }}(\sigma (u_i),{\nabla }_{\Omega }{v}_i{)}_{{\Omega }_i}-({\nu }_i\cdot \sigma (u_i),v_i{)}_{\partial {\Omega }_i}$$

$$=\sum _{i\in {\mathcal{I}}_{\Omega }}(\sigma (u_i),{\nabla }_{\Omega }{v}_i{)}_{{\Omega }_i}-(\nu \cdot \sigma (u),v{)}_{{\Gamma }_N}-\sum _{j\in {\mathcal{I}}_{{\Gamma }_I}}\underset{★}{\underbrace{\sum _{k\in {\mathcal{I}}_{\Omega }(j)}({\nu }_k\cdot \sigma (u_k),v_k{)}_{{\Gamma }_j}}}$$

$$\sum _{k\in {\mathcal{I}}_{\Omega }(j)}{\nu }_k\cdot \sigma (u_k)=0\text{on Γj}$$

$$\sum _{k\in {\mathcal{I}}_{\Omega }(j)}({\nu }_k\cdot \sigma (u_k),⟨v⟩{)}_{{\Gamma }_j}=0$$

$$⟨v⟩=\sum _{k\in {\mathcal{I}}_{{\Gamma }_I}(j)}{\alpha }_k{v}_k$$

$$★$$

$$\sum _{k\in {\mathcal{I}}_{{\Gamma }_I}(j)}(\beta {\mu }_k{h}^{-1}(u_k-⟨u⟩),v_k-⟨v⟩{)}_{{\Gamma }_j}$$

Full text

1. 1 Introduction

2. 2 Model Problem

3. 2.1 The Composite Surface

4. 2.2 Elliptic Model Problem

5. 2.3 Weak Form

6. 3 The Finite Element Method

7. 3.1 Constructions of Meshes on Composite Surfaces

8. 3.2 The Finite Element Spaces

9. 3.3 The Method

10. 3.4 Relation to the Standard Average-Jump Formulation

11. 3.5 The Stabilization Term

12. 4 Properties of the Finite Element Method

13. 4.1 Coercivity and Continuity

14. 4.2 Interpolation Error Estimate

15. 4.3 Error Estimates

16. 5 Numerical Examples

17. 5.1 Visual Examples

18. 5.2 Convergence

1 Introduction

Background.

Many physical phenomena takes place on geometries that consist of an arrangement of surfaces, for instance transport of surfactants, heat transfer, and flows in cracks. See Figure 1 for examples of surface arrangements. In manufacturing the use of surface arrangements to minimize the amount of needed material, for example in the form of honeycomb sandwich structures, is well established while recent developments of additive manufacturing enable production of even more complex surface structures. The arrangement of surfaces in applications often contain sharp edges, corners, and lines where several surfaces meet. Thus there is significant interest in the development of finite element methods for solving partial differential equations on such general geometries.

New Contributions.

In this contribution we develop a Nitsche method for a diffusion problem on a such an arrangement of surfaces. The key feature is that the formulation can handle interfaces where several surfaces meet at intersecting interfaces, including triple points and sharp edges. The method avoids defining a conormal to each interface and instead the well defined conormal associated with each subsurface is used together with the natural conservation law: the sum of all conormal fluxes is zero at the interface. This conservation law is sometimes refereed to as the Kirchhoff condition. We show that the method is equivalent to the standard Nitsche interface method in flat geometries. The same idea naturally extend to discontinuous Galerkin methods on surfaces, where instead of defining a conormal to each edge which would be needed in a standard discontinuous Galerkin method, see for instance the discussion in [6], the well defined discrete element conormal can be used.

We consider different ways of constructing a mesh on a composite geometry including matching meshes, non matching meshes, and cut meshes. For cut meshes we add a stabilization term that provides control in the vicinity of the interfaces. More precisely we consider a stabilization term, see [2] and [17], that satisfies certain abstract conditions which enable us to prove discrete stability and a priori error estimates. For clarity we restrict our analysis to the case when the geometry is exactly represented. This can be realized using parametric mappings, see [17]. We also give a concrete construction of such a stabilization term based on penalization of jumps in derivatives across faces belonging to elements that intersect the interface. We conclude the paper with some illustrating numerical examples.

Earlier Work.

Since the pioneering work of Dziuk [9] where a continuous Galerkin method for the Laplace-Beltrami operator on a triangulated surface was first proposed there has been several extensions including adaptive and higher order methods, Demlow and Dziuk [8] and Demlow [7], and higher order problems, Clarenz et al. [3] and Larsson and Larson [20]. A standard discontinuous Galerkin method for the Laplace–Beltrami operator on a smooth closed surface was analyzed in [6]. For further extensions including time dependent problems we also refer to the review articles [5, 10]. Models of membranes were considered in [13] and [15]. All of these contributions deals with smooth surfaces and the discontinuous Nitsche formulation proposed in this paper which allows more complex surfaces appears to be new. In [14] we develop a method for plate structures on composite surfaces consisting of plane surfaces with the restriction that only two plates meet at an interface. Here each plate is modeled using a membrane model and a fourth order Kirchhoff model for the bending. Various methods for connecting parametric patches pairwise which also allow for cut meshes have been proposed, see for example [19, 21, 18] and the extension to surfaces in [17].

Outline.

The outline of the remainder of the paper is as follows. In Section 2 we give a short introduction to tangential calculus on surfaces and then we formulate the model problem. In Section 3 we derive the method, study how it relates to the standard discontinuous Galerkin method in the the case of an interface in flat geometry, and introduce the stabilization term. In Section 4 we prove a priori error estimates in the energy and $L^{2}$ norm. Finally, in Section 5 we present some numerical examples illustrating the method on three different test cases.

2 Model Problem

2.1 The Composite Surface

Surface with Boundary.

Let $\Sigma$ be smooth surface embedded in $\mathbb{R}^{3}$ with orientation (or normal) $n$ and boundary $\partial\Sigma$, which has the following properties

• •

There is a smooth closed surface $\widetilde{\Sigma}$ embedded in $\mathbb{R}^{d}$ such that $\Sigma\subset\widetilde{\Sigma}$.

• •

The boundary $\partial\Sigma$ consists of a finite set of smooth curve segments and corner points $N(\partial\Sigma)$. The exterior unit conormal to $\partial\Sigma$ is denoted by $\nu_{\partial\Sigma}$.

• •

At each corner $x\in N(\partial\Sigma)$ on the boundary there is a constant such that

[TABLE]

where $\nu^{\pm}(x)$ denotes the left and right conormal at the corner $x$.

Composite Surface.

We introduce the following notation:

• •

Let $\mathcal{O}=\{\Omega_{i}:i\in\mathcal{I}_{\Omega}\}$ be a set of smooth surfaces with boundaries that satisfy the assumptions above.

• •

A composite surface $\Omega$ is a finite union of surfaces with boundaries

[TABLE]

• •

The intersection between any two surfaces $\Omega_{i}$ and $\Omega_{j}$, $i,j\in\mathcal{I}_{\Omega}$, is either empty or occur at the boundary of the surfaces,

[TABLE]

and $\Gamma_{ij}$ is either a smooth curve segment (an interface curve) or a point (an interface node). Furthermore, given two surfaces $\Omega_{i}$ and $\Omega_{j}$, there is a sequence of surfaces which starts with $\Omega_{i}$ and end at $\Omega_{j}$ such that two consecutive surfaces share a nonempty interface. This means that no surface share only a point with the union of the other surfaces.

• •

We may also assume that $\Gamma_{ij}$ is one of the curve segments in $\partial\Omega_{i}$ and $\partial\Omega_{j}$, if not we may simply modify the boundary description.

• •

Let $\mathcal{G}_{I}=\{\Gamma_{k}:k\in\mathcal{I}_{\Gamma_{I}}\}$ be the set of all non-overlapping interface curves, $\Gamma_{k}\cap\Gamma_{l}=\emptyset$ for $k\neq l$, such that

[TABLE]

Thus the union of the non-overlapping interface curves in $\mathcal{G}_{I}$ includes all interface curves (2.3).

• •

Given an index $j\in\mathcal{I}_{\Gamma_{I}}$ we let $\mathcal{I}_{\Omega}(j)\subseteq\mathcal{I}_{\Omega}$ be the set of indexes corresponding to surfaces that share interface $\Gamma_{j}$.

2.2 Elliptic Model Problem

Notation.

We introduce the following notation:

• •

The tangential gradient $\nabla_{\Omega}$ is defined by $\nabla_{\Omega}v=P_{\Omega}\nabla v$, where $P_{\Omega}=I-n\otimes n$ is the projection onto the tangent plane of $\Omega$.

• •

$\mu:\Omega\rightarrow\mathbb{R}$ is a given function such that $\mu|_{\Omega_{i}}=\mu_{i}$ with $\mu_{i}\in C(\Omega_{i})$, and there are constants such that for all $x\in\Omega$ it holds

[TABLE]

• •

The flux is defined by

[TABLE]

Model Problem.

Consider the problem: find $u:\Omega\rightarrow{\mathbb{R}}$ such that

[TABLE]

and that the following interface and boundary conditions are satisfied:

• •

For each $j\in\mathcal{I}_{\Gamma_{I}}$ the following interface conditions hold

[TABLE]

The first condition is a so called Kirchhoff condition and corresponds to conservation over interfaces while the second condition corresponds to continuity at the interface. Note that (2.8) thus encompasses the case where the interface is a sharp edge such as illustrated in Figure 2.

• •

The following boundary conditions hold

[TABLE]

where the boundary $\partial\Omega$ is decomposed into a Neumann boundary $\Gamma_{N}$ and a Dirichlet boundary $\Gamma_{D}$ such that $\partial\Omega=\Gamma_{N}\cup\Gamma_{D}$, $\Gamma_{N}\cap\Gamma_{D}=\emptyset$ and $\Gamma_{D}\neq\emptyset$.

2.3 Weak Form

Function Spaces.

Let

[TABLE]

be equipped with the $H^{1}$ inner product and associated norm

[TABLE]

Then $V_{0}$ is clearly a Hilbert space.

Derivation of Weak Form.

Starting from (2.7) and using Green’s formula on each of the surfaces $\Omega_{i}$ in $\Omega$ we obtain for $v\in V_{0}$,

[TABLE]

where we used the identity $v=0$ on $\Gamma_{D}$ and $\nu\cdot\sigma(u)=g_{N}$ on $\Gamma_{N}$ as well as the continuity $v_{k}|_{\Gamma_{j}}=v_{l}|_{\Gamma_{j}}=v$ for any $k,l\in\mathcal{I}_{\Omega}(j)$ and the interface condition (2.8) to conclude that

[TABLE]

Thus we arrive at the following weak problem: find $u\in V_{g_{D}}$ such that

[TABLE]

where

[TABLE]

Existence and Uniqueness.

The form $a(\cdot,\cdot)$ is coercive in $V_{0}$ and thus in case $g_{D}=0$ existence and uniqueness of the solution $u\in V_{0}$ follows directly from the Lax-Milgram lemma. If $g_{D}\neq 0$ we let $u_{g_{D}}\in V_{g_{D}}$ be an extension of $g_{D}$, set $u=u_{0}+u_{g_{D}}$, where $u_{0}\in V_{0}$ is determined by $a(u_{0},v)=l(v)-a(u_{g_{D}},v)$ for all $v\in V_{0}$. Here we can again apply the Lax-Milgram lemma.

Regularity of the Solution.

As the assumptions in Section 2.1 allow for very intricate surface configurations it is challenging to give any precise prediction on the regularity of the solution $u$. However, motivated by the regularity properties of elliptic problems in planar domains with nonsmooth boundary, see for example [11], we assume that

[TABLE]

where $\|f\|^{2}_{H^{\eta_{i}-2}(\Omega_{i})\cap L^{2}(\Omega_{i})}=\max(\|f\|_{H^{\eta_{i}-2}(\Omega_{i})},\|f\|_{L^{2}(\Omega_{i})})$.

3 The Finite Element Method

In this section we derive our finite element method on the composite surface $\Omega$. For clarity we consider the situation when $\Omega$ is exactly represented, which may be realized using exact parametric mappings. In the case of domains described by CAD models, the underlying assumption in isogeometric analysis, see [16, 4], the surface is equipped with a parametric finite element space consisting of continuous functions without any continuity requirement across the interface curves. In this setting we can focus on essentials and derive the method together with the basic properties including a stability result and an priori error estimate.

3.1 Constructions of Meshes on Composite Surfaces

There are several different natural ways to construct a mesh on a composite surface:

1. (a)

Each surface is meshed with elements and matching meshes are used across the interface curves. 2. (b)

Each surface is meshed with elements and non matching meshes are used across the interface curves. 3. (c)

A number of surfaces are individually meshed and arranged in such a way that they intersect. In this situation so called cut elements naturally occur close to the interface. 4. (d)

Each of the surfaces is meshed using a cut finite element technique.

Examples of these mesh constructions are illustrated in Figure 3. In the case of cut elements we add a stabilization term which enables us to prove stability and optimal order error estimates.

The Mesh.

To accommodate these different situations we define the mesh as follows. For each $\Omega_{i}$ we assume that there is a family of quasi uniform meshes $\mathcal{K}_{h,i}$ with mesh parameter $h\in(0,h_{0}]$ such that

[TABLE]

and

[TABLE]

The mesh may match $\Omega$ perfectly, i.e, we could have equality in (3.1).

3.2 The Finite Element Spaces

Finite Element Space.

Let $P_{p}(K)$ denote the space of either full or tensor product polynomials of degree less or equal to $p$ on $K$. For each mesh $\mathcal{K}_{h,i}$ there is a family of finite element spaces $V_{h,i}$, $h\in(0,h_{0}]$, such that $V_{h,i}|_{K}\in P_{p}(K)$ for all $K\in\mathcal{K}_{h,i}$. On the composite surface $\Omega$ we define the broken finite element space

[TABLE]

Approximation Property.

For each $i\in\mathcal{I}_{\Omega}$ there is an interpolation operator $\pi_{h,i}:L^{2}(\Omega_{i})\rightarrow V_{h,i}$ such that

[TABLE]

and on the composite surface we have the interpolation operator $\pi_{h}:L^{2}(\Omega)\rightarrow V_{h}$ such that $\pi_{h}|_{V_{h,i}}=\pi_{h,i}$ for $i\in\mathcal{I}_{\Omega}$. See [17] for the construction of such an interpolation operator in the case of cut meshes.

3.3 The Method

Derivation.

For $v\in V_{h}$ we obtain using Green’s formula

[TABLE]

Conservation (2.8) states that for each $j\in\mathcal{I}_{\Gamma_{I}}$,

[TABLE]

which implies

[TABLE]

where $\langle v\rangle$ is the convex combination

[TABLE]

with $0\leq\alpha_{k}\leq 1$ and $\sum_{k\in\in\mathcal{I}_{\Gamma_{I}}(j)}\alpha_{k}=1$. We may thus subtract $\langle v\rangle$ in the interface term $\bigstar$ as follows

[TABLE]

Next for each $j\in\mathcal{I}_{\Gamma_{I}}$ we add the consistent stabilization term

[TABLE]

where $\beta$ is a positive parameter.

We may also symmetrize the method by adding a consistent term and the Dirichlet boundary is taken care of using a standard Nitsche formulation.

The Method.

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

[TABLE]

where

[TABLE]

and

[TABLE]

3.4 Relation to the Standard Average-Jump Formulation

Consider an interface curve $\Gamma_{j}$ which has only two neighboring surfaces and assume that $\mu$ is a positive constant. We may assume that $\mathcal{I}_{\Gamma_{I}}(j)=\{1,2\}$. Then we have

[TABLE]

and the concistency term takes the form

[TABLE]

where $[w]=w_{1}-w_{2}$ is the jump in $w$ across $\Gamma_{j}$ and

[TABLE]

is the average of the normal flux. Note, that in the case when $\Omega_{1}$ and $\Omega_{2}$ are tangent at $\Gamma_{j}$ we have $\nu_{2}=-\nu_{1}$ and we find that

[TABLE]

which is the usual average in discontinuous Galerkin methods. Note, in particular, that our formulation thus extends the standard discontinuous Galerkin formulations to surfaces with sharp edges. Finally, we easily find that the penalty term takes the form

[TABLE]

where we used the assumption that $\mu$ is constant and the identities $v_{1}-\langle v\rangle=v_{1}-(\alpha_{1}v_{1}+\alpha_{2}v_{2})=(1-\alpha_{1})v_{1}-\alpha_{2}v_{2}=\alpha_{2}[v]$ and $v_{2}-\langle v\rangle=\alpha_{1}[v]$. We conclude that the penalty term is of the same form as in the standard discontinuous Galerkin interior penalty term.

3.5 The Stabilization Term

Abstract Properties.

In the case of cut elements at an interface we add a stabilization term of the form

[TABLE]

where $s_{h,i}$ is a positive semidefinite bilinear form on $V_{h,i}$, and define the stabilized form

[TABLE]

and the following seminorms on each subsurface $\Omega_{i}$

[TABLE]

We assume that the stabilization term has the following properties:

• •

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

[TABLE]

• •

There is a constant such that for all $i\in\mathcal{I}_{\Omega}$ and $v\in H^{s+1}(\Omega_{i})$ with $0\leq s\leq p$ it holds

[TABLE]

Here and below we use the abbreviated notation $a\lesssim b$ for the inequality $a\leq cb$ where $c$ is a constant independent of the mesh size parameter $h$.

Remark 3.1

In order to prove a bound on the condition number of the stiffness matrix we assume that

[TABLE]

The condition number bound can then be proved using the techniques in [17].

Remark 3.2

Note that we do not assume that the stabilization form is consistent, i.e. for the exact solution $u$, $s_{h}(u,v)=0$ for $v\in V_{h}$, even though this may be the case for a sufficiently regular solution.

Normal Derivative Jump Penalty.

The stabilization form which we will use in this paper takes the form

[TABLE]

where $\mathcal{F}_{h,i}$ is the set of interior faces in the mesh $\mathcal{K}_{h,i}$, which belong to at least one element $K$ such that $\text{int}(K)\cap\partial\Omega_{i}\neq\emptyset$. We note that $s_{h}$ is consistent for sufficiently regular functions, for instance for $v\in C^{p}(\Omega)$ we have $s_{h}(v,w)=0$ for $w\in V_{h,i}$.

We first note that $s_{h}$ is a bilinear positive semidefinite form by construction. To verify (3.29) we recall, see [17] Lemma 4.1, that for two neighboring elements $K_{1}$ and $K_{2}$ in $\mathcal{K}_{h,i}$ which share the face $F$ we have the estimate

[TABLE]

and we may conclude that, for $h\in(0,h_{0}]$ with $h_{0}$ small enough,

[TABLE]

where $\text{int}(\mathcal{K}_{h,i})=\{K\in\mathcal{K}_{h,i}:\text{int}(K)\cap\partial\Omega_{i}=\emptyset\}$ is the set of elements that do not intersect the boundary. Finally, we have the inverse estimate

[TABLE]

where $\mathcal{K}_{h,i}(\partial\Omega_{i})$ is the set of elements that intersect the boundary in a face or in the interior. See [12] for the trace inequality $\|v\|^{2}_{\partial\Omega_{i}\cap K}\lesssim h^{-1}\|v\|^{2}_{K}+h\|\nabla v\|^{2}_{K}$, $v\in H^{1}(K)$, from which (3.35) follows. Combining estimates (3.34) and (3.35) we find that

[TABLE]

To verify (3.30) we consider again the pair of two elements $K_{1}$ and $K_{2}$ sharing the face $F$ and we let $w\in P_{p}(K_{1}\cup K_{2})$. We note that $s_{h,i,F}(w,w)=0$ and thus

[TABLE]

where we used an inverse inequality and the Bramble–Hilbert Lemma to estimate the second term. Summing over all $F\in\mathcal{F}_{h,i}$ and using the approximation property (3.4) give (3.30).

Least Squares Gradient Variation Penalty.

Define

[TABLE]

where $P_{F}:H^{1}(K_{1}\cup K_{2})\rightarrow P_{p}(K_{1}\cup K_{2})$ is the $H^{1}$ projection. This stabilization term is not consistent but it is not difficult to verify that it satisfies the conditions (3.30) and (3.30) using similar arguments as above. See also [1] where related stabilization terms were studied.

4 Properties of the Finite Element Method

4.1 Coercivity and Continuity

We define the norms

[TABLE]

and

[TABLE]

where

[TABLE]

Lemma 4.1

For large enough penalty parameter $\beta$, there is a constant such that for all $v\in V_{h}$ it holds

[TABLE]

There is a constant such that for $v,w\in V_{h}+V$ it holds

[TABLE]

and

[TABLE]

Proof. For large enough penalty parameter $\beta$ the coercivity follows from the fact that the inverse inequality (3.29) holds together with standard arguments.

4.2 Interpolation Error Estimate

We have the estimate

[TABLE]

This inequality follows from the trace inequality

[TABLE]

for each element $K\in\mathcal{K}_{h,i}$, $i\in\mathcal{I}_{\Omega}(j)$ with a nonempty intersection with $\Gamma_{j}$ and the interpolation estimate (3.4).

4.3 Error Estimates

In the proof of the $L^{2}$ error estimate we will use a duality argument and we now specify the regularity properties needed for our analysis. Let $\phi\in V_{0}$ be the solution to the dual problem

[TABLE]

where $\psi\in L^{2}(\Omega)$. We assume that there are constants $\eta_{i}^{*}\in(3/2,2],i\in\mathcal{I}_{\Omega}$, and a hidden constant such that for all $\psi\in L^{2}(\Omega)$, the regularity estimate

[TABLE]

holds.

Theorem 4.1** (Error Estimates)**

Assume that the exact solution $u$ to (2.20) satisfies the regularity estimate (2.22), then there is a constant such that

[TABLE]

with $\widetilde{\eta}_{i}=\min(\eta_{i},p+1)$

If in addition the solution to the dual problem (4.9) satisfies the regularity estimate (4.10), then there is a constant such that

[TABLE]

where $\eta^{*}=\min_{i\in\mathcal{I}_{\Omega}}\eta^{*}_{i}$

Remark 4.1

Note that it follows from (4.11) and (3.30) that

[TABLE]

Proof. (4.11). Splitting the error into an interpolation error and a discrete error we have

[TABLE]

where we used the energy norm interpolation error estimate (4.7). Next using coercivity (4.4) we have

[TABLE]

Here the numerator may rewritten, using the definition of the method (3.13), $A_{h}(u_{h},v)=l_{h}(v)$, $v\in V_{h}$, and the fact that the unstabilized method is consistent $a_{h}(u,v)=l_{h}(v)$, $v\in V_{h}$, as follows

[TABLE]

Thus we find that we have the bound

[TABLE]

where we used the energy norm interpolation estimate (4.7) and the approximation property (3.30) of $s_{h}$. Combining (4.16) and (4.23), gives

[TABLE]

which together with (4.15) concludes the proof.

(4.12).

Let $\phi$ be the solution to the dual problem (4.9) with $\psi=e$, where $e=u-u_{h}$ is the error. We note that by consistency we have

[TABLE]

Setting $v=e$ we obtain

[TABLE]

Here we used the energy norm bound (4.11) and the energy interpolation estimate (4.7) followed by the elliptic regularity bound (4.10) to conclude that the following estimate holds

[TABLE]

5 Numerical Examples

Geometries.

We demonstrate the method using the three composite surfaces illustrated in Figure 4, where each geometry has different features:

• •

The cube with holes in Figure 4(a) consists of 6 separate surfaces pairwise connected along 12 interface lines. The resulting compound surface feature both sharp edges and corners.

• •

The house of cards in Figure 4(b) consists of 18 separate surfaces and 10 interface lines. Here each interface connects 2-6 surfaces. Also this geometry has sharp edges.

• •

The intersecting cylinders composite surface in Figure 4(c) is constructed by intersecting two cylinder surfaces and cutting the surfaces along their intersection. This construction produces 6 surfaces connected along the interface curves. In this geometry both the surfaces and the interfaces are curved and 4 surfaces meet at each interface.

General Construction.

All examples share the following set-up:

• •

The geometry of each surface is exactly described by a mapping $F\colon\mathbb{R}^{2}\supset\widehat{\Omega}\to\Omega$, where $\widehat{\Omega}_{i}$ is the reference domain and $F(\widehat{\Omega}_{i})=\Omega_{i}$.

• •

The elements used in all examples are parametric quadratic tensor product Lagrange elements ($p=2$) and we allow cut elements at the interfaces. The parametrization is based on the exact map $F$,

[TABLE]

where $\widehat{V}_{h}$ is the finite element space in the reference domain. The stabilization term $s_{h}$ is evaluated in the reference domain, see [17] for further details.

• •

Our Nitsche penalty parameter is chosen to be $\beta=100$ and our stabilization parameter to be $\gamma_{k}=10^{-2}$.

• •

We solve the model problem $-\nabla_{\Omega}\cdot\sigma(u)=f$, where $\sigma(u)=\mu\nabla_{\Omega}u$ and $\mu=1$, together with Dirichlet and Neumann boundary conditions.

5.1 Visual Examples

Cube with holes.

The cube illustrated in Figure 4(a) is an example of a composite surface which includes both sharp edges and corners. We consider our model problem with load $f=0$ and non-homogeneous Dirichlet boundary conditions. The boundary in this case consist of the two holes and we let $u=1$ along the top hole and $u=-1$ along the side hole. We present a finite element solution and the magnitude of its gradient in Figure 5 using unfitted meshes on each surface. Note that as implied by the interface conditions (2.8)–(2.9) both the solution and the flux seem to flow continuously over the interfaces. Curiously, the gradient near the cube corners do not seem to include any singularity, as could possibly be expected, but is rather zero at the corner. We highlight this in Figure 6.

House of Cards.

The house of cards composite surface shown in Figure 4(b) includes interfaces joining more than two surfaces. Again we consider our model problem with $f=0$ and we set the non-homogeneous Dirichlet boundary condition $u=1$ on the right side of the top card closest in view and $u=-1$ on the left side of the bottom standing card closest in view. On the remainder of the boundary we use homogeneous Neumann conditions. This leads to singularities in the points where we switch boundary conditions. Note that the singularities spread to the surfaces that couple at the interface. However, as the Kirchhoff condition (2.8) states that the sum of the conormal fluxes should be zero we no longer have pairwise continuity of fluxes over interfaces joining more than two surfaces. We can see this in Figure 8 where we also note the peculiar effect that the singularities due to the switch of boundary condition at the interfaces has a greater effect on interfaces joining more surfaces. While somewhat counter intuitive we isolate this effect in the study presented in Figure 9.

5.2 Convergence

Intersecting Cylinders.

Here we consider the composite surface constructed by intersecting two cylinders with radii $R_{1}<R_{2}$ and cutting the surfaces along the intersection as shown in Figure 4(c). The cylinder axes are parallel to the $xy$-plane, their relative rotation in the $xy$-plane is $\pi/2+\theta$, $0\leq\theta\leq\pi/2$, and their relative offset in $z$-direction is $z_{0}$. The points on each cylinder thus satisfy

[TABLE]

respectively. The cylinder intersection satisfy both these equations which after simplification gives that the intersection is described by

[TABLE]

To describe the example geometry displaced in Figure 4(c) we use parameters $R_{1}=1.5$, $R_{2}=2$, $\theta=2\pi/3$ and $z_{0}=0.15$.

We manufacture a problem with known solution by choosing a solution $u$ which is analytical on each subsurface $\Omega_{i}$ and generating the appropriate right hand side $f=-\Delta_{\Omega}u$ and matching non-homogeneous Dirichlet boundary conditions. For any function $u$ which have a smooth extension to a volumetric neighbourhood to $\Omega_{i}$ the Laplace–Beltrami operator $\Delta_{\Omega}$ can be expressed in terms of the usual $\mathbb{R}^{3}$ Laplacian $\Delta$ and the first and second normal derivatives as

[TABLE]

where $H=\nabla\cdot n$ is the mean curvature of $\Omega_{i}$. In the present example we choose the solution as the restriction of

[TABLE]

to $\Omega$ and by construction (5.5) also is an extension of $u$ to $\mathbb{R}^{3}$. This solution and the magnitude of its gradient is presented in Figure 10 where we separated the subsurfaces to show the solution also on the subsurfaces hidden from view. Note that this solution will also satisfy the Kirchhoff condition (2.8) as $u$ is smooth and that for each subsurface in an interface there is a subsurface with exactly opposing conormal.

Using this manufactured problem we study the convergence of the method (3.13) in $L^{2}(\Omega)$ norm with $p=1$ and $p=2$ finite elements. An example finite element solution using $p=2$ elements is presented in Figure 11 and the convergence results are presented in Figure 12 where we note that we achieve optimal order convergence of $\mathcal{O}(h^{p+1})$. Looking at the $L^{2}$ estimate (4.12) this is expected as the manufactured solution $u$ is smooth by definition and each subsurface in the intersecting cylinders composite surface has a smooth boundary which should give the necessary regularity for the dual solution $\phi$.

Using the same problem set-up we also investigate how the condition number of the stiffness matrix scales with $h$ for $p=1$ elements. As shown in Figure 13 we get the expected scaling of $\mathcal{O}(h^{-2})$.