UNFITTED FINITE ELEMENT METHODS USING ... - Semantic Scholar
SIAM J. NUMER. ANAL. Vol. 52, No. 4, pp. 2137–2162
c 2014 Society for Industrial and Applied Mathematics
UNFITTED FINITE ELEMENT METHODS USING BULK MESHES FOR SURFACE PARTIAL DIFFERENTIAL EQUATIONS∗ KLAUS DECKELNICK† , CHARLES M. ELLIOTT‡ , AND THOMAS RANNER§ Abstract. In this paper, we define new unfitted finite element methods for numerically approximating the solution of surface partial differential equations using bulk finite elements. The key idea is that the n-dimensional hypersurface, Γ ⊂ Rn+1 , is embedded in a polyhedral domain in Rn+1 consisting of a union, Th , of (n + 1)-simplices. The unifying feature of the methodological approach is that the finite element approximating space is based on continuous piecewise linear finite element functions on the bulk triangulation Th which is independent of Γ. Our first method is a sharp interface method (SIF ) which uses the bulk finite element space in an approximating weak formulation obtained from integration on a polygonal approximation, Γh , of Γ. The full gradient is used rather than the projected tangential gradient and it is this which distinguishes SIF from the method of [M. A. Olshanskii, A. Reusken, and J. Grande, SIAM J. Numer. Anal., 47 (2009), pp. 3339–3358]. The second method is a narrow band method (NBM ) in which the region of integration is a narrow band of width O(h). NBM is similar to the method of [K. Deckelnick et al., IMA J. Numer. Anal., 30 (2010), pp. 351–376] but again the full gradient is used in the discrete weak formulation. The a priori error analysis in this paper shows that the methods are of optimal order in the surface L2 and H 1 norms and have the advantage that the normal derivative of the discrete solution is small and converges to zero. Our third method combines bulk finite elements, discrete sharp interfaces, and narrow bands in order to give an unfitted finite element method for parabolic equations on evolving surfaces. We show that our method is conservative so that it preserves mass in the case of an advection-diffusion conservation law. Numerical results are given which illustrate the rates of convergence. Key words. unfitted finite elements, cut cells, error analysis, narrow band, sharp interface, elliptic and parabolic surface equations AMS subject classifications. 35R01, 65N30, 65N15, 65M60 DOI. 10.1137/130948641
1. Introduction. In this article we propose and analyze numerical methods based on bulk finite element meshes for the following model elliptic equation on a stationary surface. Model elliptic equation on stationary surface. Let Γ be a smooth hypersurface in Rn+1 and f ∈ L2 (Γ). We seek solutions u : Γ → R of (1.1)
−ΔΓ u + u = f
on Γ.
The methods can be extended in natural ways to deal with variable coefficients and nonlinearities. The approach may be extended to the following advection-diffusion equation on a moving surface. ∗ Received
by the editors December 10, 2013; accepted for publication (in revised form) May 30, 2014; published electronically August 19, 2014. http://www.siam.org/journals/sinum/52-4/94864.html † Institut f¨ ur Analysis und Numerik, Otto-von-Guericke-Universit¨ at Magdeburg, 39106 Magdeburg, Germany ([email protected]). ‡ Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK (C.M.Elliott@warwick. ac.uk). § Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK. Current address: School of Computing, University of Leeds, LS2 9JT, UK ([email protected]). The work of this author was supported by an EPSRC Ph.D. studentship (grant EP/P504333/1 and EP/P50516X/1) and the Warwick Impact Fund. 2137
2138
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
Model parabolic equation on evolving surface. Let {Γ(t)} be a family of smooth hypersurfaces in Rn+1 for t ∈ [0, T ]. Denoting by ∂ • u the material derivative of u and v the velocity of Γ(t) (see section 5 for notation), we seek solutions u : t Γ(t) × {t} of the advection-diffusion equation (1.2a) Γ(t) × {t}, ∂ • u + u∇Γ · v − ΔΓ u = f on t∈(0,T )
(1.2b)
u(·, 0) = u0
on Γ(0).
Surface partial differential equations or partial differential equations (PDEs) on manifolds arise in a wide variety of applications in materials science, fluid dynamics, and biology, [5, 32, 24, 29, 25, 2, 28, 26]. Computational approaches include surface finite elements on triangulated surfaces [17, 19, 18, 15, 14, 23, 22], bulk finite element or finite difference meshes for the approximation of implicit surface formulations [9, 30, 10, 21, 11], bulk finite element or finite difference meshes on narrow bands [12, 38], and bulk finite element meshes and sharp interface weak forms [34, 33, 16]. An important feature of the methods cited above is the avoidance of charts both in the problem formulation and the numerical methods. For example, the surface finite element method is based simply on triangulated surfaces and requires the geometry solely through the knowledge of the vertices of the triangulation, whereas methods based on implicit surfaces require only the level set function Φ which encodes all the necessary geometry. Another feature of some of these methods is the use of unfitted bulk meshes. Here we use the terminology unfitted finite element methods (sometimes called cut cell methods) when the underlying meshes that form the computational domain are not fitted to the domain in which the PDE holds. The motivation for using finite element spaces on meshes not fitting to the domain came from the desire to solve free or moving boundary problems. Such methods were introduced in [3, 4] for elliptic equations in curved domains; see also [31, 8, 27]. In this setting we are concerned with bulk meshes independent of the surface. The new methods. The new unfitted finite element methods for surface elliptic equations proposed in this paper are variants of the bulk finite element approaches using a sharp interface or a narrow band. The new scheme for advection diffusion on an evolving surface is a hybrid of these. In the following we sketch the main ideas of these methods describing the details in sections 3–5. Sharp interface method (SIF ). Given an interpolation Γh of Γ, we use a bulk finite element space VhI of the form VhI = {φh ∈ C 0 (UhI ) | φh|T ∈ P1 (T ) for each T ∈ ThI }, where ThI is a set of elements which intersect Γh and UhI = T ∈T I T ; see section 3. h The discrete scheme approximating the model elliptic equation (1.1) is find uh ∈ VhI such that ∇uh · ∇φh + uh φh dσh = f e φh dσh for all φh ∈ VhI , (1.3) Γh
Γh
where f e is an extension of f . The method is related to the following method of Olshanskii, Reusken, and Grande, introduced in [34]: find uh ∈ VhΓ such that (1.4) ∇Γh uh · ∇Γh φh + uh φh dσh = f e φh dσh for all φh ∈ VhΓ . Γh
Γh
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
2139
There are two significant differences. Note the use of the full gradient in (1.3) as opposed to the tangential gradient. This gives control over the normal derivative of the finite element solution which is lacking in (1.4). Another difference relates to the use of the finite element space VhΓ , which essentially consists of the traces on Γh of elements in VhI . While VhI has a natural basis, this does not seem to be the case for VhΓ . The “standard basis” of finite element hat functions is only a spanning set for VhΓ . Narrow band method (NBM ). We use the bulk finite element space VhB on the triangulation ThB , VhB = {φh ∈ C 0 (UhB ) | φh|T ∈ P1 (T ) for each T ∈ ThB }. Here ThB consists of those triangles intersecting a narrow band domain D h defined by the ±h level sets of an interpolated level set function Ih Φ and UhB = T ∈T B T . h The discrete scheme approximating the model elliptic equation (1.1) is find uh ∈ VhB such that (1.5)
∇uh · ∇φh + uh φh |∇Ih Φ| dx =
Dh
f e φh |∇Ih Φ| dx Dh
for all φh ∈ VhB .
This is similar to the method in [12] except that NBM uses the full instead of projected gradients thus avoiding the resulting degeneracy. As a result we are able to prove an optimal L2 -error bound which was not obtained for the method in [12]. It is also the case that the normal derivative of the discrete solution converges to zero. Hybrid unfitted evolving surface method. The discrete problem approxim+1 m ∈ Vhm+1 such that mating (1.2) is given um h ∈ Vh , m = 0, . . . , N − 1, find uh
(1.6)
Γm+1 h
um+1 φh dσh − h
τm + 2h
m+1 Dh
∇um+1 h
Γm h
e,m+1 um ) dσh h φh (· + τm v
· ∇φh ∇Ih Φm+1 dx = τm
Γm+1 h
f e,m+1 φh dσh
for all φh ∈ Vhm+1 . Here v e,m denotes an extension of the surface velocity at time level m. We use time step labeled analogues of the notation for NBM; see section 5 for the details. Here, u0h is appropriate initial data. Because of this combination of narrow band and sharp interface discretization, under some mild constraints on the discretization parameters (see section 5) our numerical scheme preserves the important property that solutions of (1.2) conserve mass in the case that f ≡ 0. Outline. The paper is organized as follows: In section 2 we introduce our notation and collect some auxiliary results. In sections 3 and 4 we present and analyze unfitted methods for the model elliptic equation (1.1). In section 5 we describe how a combination of these two approaches can be used to calculate solutions of the advection-diffusion equation on evolving hypersurfaces, (1.2). Details of the implementation and several numerical examples illustrating the orders of convergence are presented in section 6.
2140
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
2. Preliminaries. 2.1. Surface calculus. Let Γ be a connected compact smooth hypersurface embedded in Rn+1 (n = 1, 2). We assume that there exists a smooth function Φ : U → R such that Γ = {x ∈ U | Φ(x) = 0}
∇Φ(x) = 0, x ∈ U,
and
where U is an open neighborhood of Γ. For a function z : Γ → R we define its tangential gradient by (2.1) ∇Γ z(p) := ∇ z (p) − ∇ z (p) · ν(p) ν(p), p ∈ Γ, where z : U → R is an arbitrary smooth extension of z to U and ν(x) =
∇Φ(x) |∇Φ(x)|
is a unit vector to the level sets of Φ. It can be shown that ∇Γ z(p) is independent of the particular choice of z. We denote by Di z, 1 ≤ i ≤ n + 1, the components of ∇Γ z. Furthermore, we let ΔΓ z = ∇Γ · ∇Γ z =
n+1
Di Di z
i=1
be the Laplace–Beltrami operator of z. In what follows it will be convenient to use special coordinates which are adapted to Φ. Consider for p ∈ Γ the system of ODEs γp (s) =
(2.2)
∇Φ(γp (s)) , |∇Φ(γp (s))|2
γp (0) = p.
It can be shown that there exists δ > 0 so that the solution γp of (2.2) exists uniquely on (−δ, δ) uniformly in p ∈ Γ, so that we can define the mapping F : Γ × (−δ, δ) → Rn+1 by (2.3)
F (p, s) := γp (s),
p ∈ Γ, |s| < δ.
d Since ds Φ(γp (s)) = 1 and γp (0) = p ∈ Γ, we infer that Φ(γp (s)) = s, |s| < δ, and hence that x = F (p, s) implies that |Φ(x)| < δ. As a result, we deduce that F is a diffeomorphism of Γ × (−δ, δ) onto Uδ := {x ∈ U | |Φ(x)| < δ} with inverse
F −1 (x) = (p(x), Φ(x)),
(2.4)
x ∈ Uδ ,
where p : Uδ → Rn+1 satisfies p(x) ∈ Γ, x ∈ Uδ . For later purposes it is convenient to expand p and its derivatives in terms of Φ. Let us fix x ∈ Uδ and define the function η(τ ) := F (p(x), (1 − τ )Φ(x)), τ ∈ [0, 1]. Since
∂F ∂s
(p, s) = γp (s) we have
∇Φ(γp(x) ((1 − τ )Φ(x))) η (τ ) = −Φ(x)γp(x) ((1 − τ )Φ(x)) = −Φ(x) . ∇Φ(γp(x) ((1 − τ )Φ(x))) 2
2141
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
Observing that γp(x) (Φ(x)) = F (p(x), Φ(x)) = x and using similar arguments to calculate η (τ ) we find that ηk (0) = −Φ(x)
Φxk (x)
2, |∇Φ(x)|
n+1 2Φxk (x)Φxr (x) Φxl (x)Φxl xr (x) 2 ηk (0) = Φ(x) , δkr − 2 4 |∇Φ(x)| |∇Φ(x)| l,r=1
k = 1, . . . , n + 1. Since η(1) = p(x), η(0) = x we deduce with the help of Taylor’s theorem that (2.5)
n+1 Φxk (x) 1 2Φxk (x)Φxr (x) Φxl (x)Φxl xr (x) 2 pk (x) = xk − Φ(x) + Φ(x) δkr − |∇Φ(x)|2 2 |∇Φ(x)|2 |∇Φ(x)|4 l,r=1
+ Φ(x)3 rk (x),
k = 1, . . . , n + 1,
where rk are smooth functions. In a similar way we may write ∇Φ(x) = ∇Φ(p(x)) + Φ(x)G(x),
(2.6)
1 where G(x) = 0 D2 Φ(F (p(x), τ Φ(x))) ∂F ∂s (p(x), τ Φ(x)) dτ . Let us next use the function p in order to define a particular extension of z : Γ → R: x ∈ Uδ .
z e (x) := z(p(x)),
(2.7)
Since p(F (p(x), s)) = p(x) we deduce that s → z e (F (p(x), s)) is independent of s and thus ∇z e (x) · ν(x) = 0,
(2.8)
x ∈ Uδ .
In order to express the derivatives of z e in terms of the tangential derivatives of z we first deduce from (2.5) that pk,xi (x) =δik −
Φxk (x)Φxi (x) |∇Φ(x)|
+ Φ(x)Φxi (x)
2
n+1
−
Φ(x)Φxk xi (x)
|∇Φ(x)|
δkr −
l,r=1
2
+ 2Φ(x)Φxk (x)
2Φxk (x)Φxr (x) |∇Φ(x)|
2
n+1
Φxl (x)Φxl xi (x)
l=1
|∇Φ(x)|
Φxl (x)Φxl xr (x) |∇Φ(x)|4
4
+ Φ(x)2 αik (x).
Combining this relation with (2.6) we deduce that (2.9) pk,xi (x) = δik − νi (p(x))νk (p(x)) + aik (x)Φ(x), pk,xi xj (x) = − (2.10)
+
Φxi (x)Φxk xj (x) Φxj (x)Φxk xi (x) − |∇Φ(x)|2 |∇Φ(x)|2 n+1 Φxi (x)Φxj (x) Φxl (x)Φxk xl (x) + βkij (x)νk (p(x)) + γkij (x)Φ(x), |∇Φ(x)|2 |∇Φ(x)|2 l=1
2142
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
where aik , βkij , γkij are smooth functions. Differentiating (2.7) and using (2.9), (2.10) n+1 as well as the fact that k=1 Dk z(p(x))νk (p(x)) = 0 we obtain (2.11) ∇z e (x) = I + Φ(x)A(x) ∇Γ z(p(x)), (2.12) 1 ∇ · |∇Φ(x)|∇z e (x) |∇Φ(x)|
⎛
= (ΔΓ z)(p(x)) + Φ(x) ⎝
n+1
blk (x)D l Dk z(p(x)) +
k,l=1
n+1
⎞ ck (x)D k z(p(x))⎠ ,
k=1
where A = (aik ), blk , and ck are again smooth. 2.2. Bulk finite element space and inequalities. In what follows we assume that U is polyhedral. Let (Th )0 0. As a consequence,
Φ − Ih Φ L∞ (U) + h ∇(Φ − Ih Φ) L∞ (U) ≤ Ch2 ,
so that we may assume that there exist constants c0 , c1 such that c0 ≤ |∇Ih Φ(x)| ≤ c1 ,
(2.16)
x ∈ U, 0 < h ≤ h0 .
Let us next define (2.17)
Γh := {x ∈ U | Ih Φ(x) = 0}
and
Dh := {x ∈ U | |Ih Φ(x)| < h}
as approximations of the given hypersurface Γ and the neighborhood Dh := {x ∈ U | |Φ(x)| < h} ; see Figure 1, for example. Note that Γh is a polygon whose facets are line segments if n = 1 and a polyhedral surface whose facets consist of triangles or quadrilaterals if n = 2. The corresponding decomposition of Γh is in general not shape regular and can have arbitrary small elements. Furthermore, we introduce Fh : U → Rn+1 by Fh (x) := F (p(x), Ih Φ(x)),
x ∈ U,
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
2143
Fig. 1. A cartoon of the domains of the sharp interface (left) and the narrow band (right) method. The surface Γ, resp., the set D h is displayed in red, the approximations Γh , resp., Dh in blue, and the domains UhI , UhB in gray.
where F was defined in (2.3). From the properties of F we infer that (2.18)
p(Fh (x)) = p(x)
and
Φ(Fh (x)) = Ih Φ(x) if Fh (x) ∈ Uδ , Fh (x) = p(x) if x ∈ Γh .
(2.19)
Lemma 2.1. There exists 0 < h1 ≤ h0 such that for 0 < h ≤ h1 the mapping Fh : Dh → Dh := {x ∈ U | |Φ(x)| < h} is bi-Lipschitz with Fh (Γh ) = Γ. Furthermore,
Fh − Id L∞ (U) + h DFh − I L∞ (U) ≤ ch2 , |det DFh | − |∇Ih Φ| ≤ ch2 . |∇Φ| ∞
(2.20) (2.21)
L
(U)
Proof. Since F (p(x), Φ(x)) = x we deduce with the help of (2.15) |Fh (x) − x| = |F (p(x), Ih Φ(x)) − F (p(x), Φ(x))| ≤ c|Ih Φ(x) − Φ(x)| ≤ ch2 . Differentiating the relation Fi (p(x), Φ(x)) = xi , i = 1, . . . , n + 1, we obtain n+1
Dk Fi (p(x), Φ(x))pk,xj (x) +
k=1
∂Fi (p(x), Φ(x))Φxj (x) = δij , ∂s
i, j = 1, . . . , n + 1,
and hence Fhi,xj (x) =
n+1
Dk Fi (p(x), Ih Φ(x))pk,xj (x) +
k=1
∂Fi (p(x), Ih Φ(x))(Ih Φ)xj (x) ∂s
∂Fi (p(x), Φ(x)) Ih Φ − Φ xj (x) ∂s n+1 + Dk Fi (p(x), Ih Φ(x)) − Dk Fi (p(x), Φ(x)) pk,xj (x)
= δij + (2.22)
k=1
∂Fi ∂Fi (p(x), Ih Φ(x)) − (p(x), Φ(x)) (Ih Φ)xj (x) ∂s ∂s Φxi (x) Ih Φ − Φ xj (x) + rij (x), = δij + |∇Φ(x)|2 +
2144
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
where |rij (x)| ≤ ch2 in view of (2.15). This implies (2.20). In particular we deduce that Fh is bi-Lipschitz provided that h is sufficiently small, whereas the properties Fh (Dh ) = Dh and Fh (Γh ) = Γ follow from (2.18). Finally we deduce from (2.22) that ∇Φ ∇Φ · ∇Ih Φ · ∇(Ih Φ − Φ) + ch = + ch |∇Φ|2 |∇Φ|2 2 ∇Φ |∇Ih Φ| |∇Ih Φ| 1 ∇Ih Φ |∇Ih Φ| − − + ch = + dh , = |∇Φ| 2 |∇Ih Φ| |∇Φ| |∇Φ| |∇Φ|
|det DFh | = 1 +
where |ch | , |dh | ≤ ch2 proving (2.21). We introduce μh : Γh → R via dσ(p(x)) = μh (x) dσh (x). It is well known (see Proposition 2.1 in [15], (3.37) in [34]) that |1 − μh | ≤ ch2
(2.23)
on Γh .
Using the properties of Fh together with the coarea formula and (2.9),(2.10), (2.11), (2.23) one can prove the following result on the equivalence of certain norms. Lemma 2.2. There exist constants c1 , c2 > 0 which are independent of h, such that for all z ∈ H 1 (Γ) c1 z e L2 (Γh ) ≤ z L2 (Γ) ≤ c2 z e L2 (Γh ) , 1 1 c1 √ z e L2 (Dh ) ≤ z L2 (Γ) ≤ c2 √ z e L2 (Dh ) , h h e c1 ∇z L2 (Γh ) ≤ ∇Γ z L2 (Γ) ≤ c2 ∇z e L2 (Γh ) , 1 1 c1 √ ∇z e L2 (Dh ) ≤ ∇Γ z L2 (Γ) ≤ c2 √ ∇z e L2 (Dh ) . h h If in addition z ∈ H 2 (Γ) then 1 c1 √ D2 z e L2 (D ) ≤ z H 2 (Γ) . h h 2.3. Variational form of elliptic equation and Strang’s second lemma. It is well known [1] that for every f ∈ L2 (Γ) there exists a unique solution u ∈ H 2 (Γ) of (1.1) which satisfies
u H 2 (Γ) ≤ c f L2 (Γ) .
(2.24)
Let us write (1.1) in weak form: (2.25)
a(u, ϕ) = l(ϕ)
where
a(w, ϕ) = Γ
for all ϕ ∈ H 1 (Γ),
∇Γ w · ∇Γ ϕ + wϕ dσ,
l(ϕ) =
f ϕ dσ. Γ
Next, suppose that Vh is a finite–dimensional space and V e := {v e | v ∈ H 1 (Γ)}. Assume that ah : (Vh + V e ) × (Vh + V e ) → R is a symmetric, positive semidefinite bilinear form which is, in addition, positive definite on Vh × Vh . Furthermore, let lh : Vh → R be linear. Then the approximate problem (2.26)
ah (uh , vh ) = lh (vh )
for all vh ∈ Vh
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
has a unique solution uh ∈ Vh . Introducing
v h := ah (v, v),
2145
v ∈ Vh + V e ,
we have by Strang’s second lemma (2.27)
ue − uh h ≤ 2 inf ue − vh h + sup vh ∈Vh
φh ∈Vh
|ah (ue , φh ) − lh (φh )| .
φh h
3. Sharp interface method (SIF ). 3.1. Setting up the method. Let us begin by observing that if T ∈ Th satisfies H n (T ∩ Γh ) > 0, then the following two cases can occur: (1) Γh ∩ int(T ) = ∅, in which case H n (∂T ∩ Γh ) = 0; (2) T ∩ Γh = ∂T ∩ Γh in which case T ∩ Γh is the face between two elements. We may now define a unique subset ThI ⊂ Th by taking all elements satisfying case 1 and in case 2 taking just one of the two elements T . The numerical method does not depend on which element is chosen. We may therefore conclude that there exists N ⊂ Γh with H n (N ) = 0 and a subset ThI ⊂ Th such that every x ∈ Γh \ N belongs to exactly one T ∈ ThI . We then define UhI = T. T ∈ThI
Clearly UhI ⊆ Uδ if h is small enough. We define the finite element space VhI by VhI = {φh ∈ C 0 (UhI ) | φh|T ∈ P1 (T ) for each T ∈ ThI }. Note that ∇φh is defined on Γh \ N in view of the definition of ThI . In particular the unit normal νh to Γh is given by νh =
(3.1)
∇Ih Φ |∇Ih Φ|
on Γh \ N,
and we use (3.1) in order to extend νh to UhI . Let us next turn to the approximation ¯δ ) so error for the space VhI . Note that for a function z ∈ H 2 (Γ) we have z e ∈ C 0 (U e that Ih z is well–defined. Lemma 3.1. Let z ∈ H 2 (Γ). Then (3.2)
z e − Ih z e L2 (Γh ) + h ∇(z e − Ih z e ) L2 (Γh ) ≤ ch2 z H 2 (Γ) .
Proof. We first observe that Theorem 3.7 in [34] yields (3.3)
z e − Ih z e L2 (Γh ) + h ∇Γh (z e − Ih z e ) L2 (Γh ) ≤ ch2 z H 2 (Γ) .
Hence, it remains to bound ∇(z e −Ih z e )·νh L2 (Γh ) . To do so, we start by considering an element T ∈ ThI . Then we see that 2 |∇(z e − Ih z e ) · νh | dσh T ∩Γh 2 2 e ≤2 |∇z · νh | dσh + 2 |∇(Ih z e ) · νh | dσh T ∩Γh T ∩Γh 2 2 e −1 ≤2 |∇z · (νh − ν)| dσh + ch(T ) |∇(Ih z e ) · νh | dx =: I1 + I2 , T ∩Γh
T
2146
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
in view of (2.8) and the fact that H n (T ∩ Γh ) ≤ ch(T )−1 H n+1 (T ). Note that by (3.1) and (2.15)
ν − νh L∞ (T )
(3.4)
∇Φ ∇Ih Φ − = ≤ ch(T ) |∇Φ| |∇Ih Φ| L∞ (T )
so that I1 ≤ ch2
T ∩Γh
2
|∇z e | dσh .
Furthermore, recalling (2.14) and using again (3.4) I2 ≤ ch(T )−1
T
2 2 2 |∇z e · (νh − ν)| + |∇(z e − Ih z e )| dx ≤ ch z e H 2 (T ) .
We use the bounds for I1 , I2 and sum over all elements T ∈ ThI , then apply Lemma 2.2 to see |∇(z e − Ih z e ) · νh |2 dσh ≤ ch2 ∇z e 2L2 (Γh ) + ch z e 2H 2 (Dc h ) ≤ ch2 z 2H 2 (Γ) , 1
Γh
since T ⊂ Dc1 h for all T ∈ ThI in view of (2.16). 3.2. The method. Let us write (1.3) in the form find uh ∈ VhI such that (3.5)
ah (uh , φh ) = lh (φh )
for all φh ∈ VhI ,
where ah (wh , φh ) =
Γh
∇wh · ∇φh + wh φh dσh ,
lh (φh ) =
Γh
f e φh dσh .
In order to verify that the symmetric bilinear form ah is positive definite on VhI × VhI we note that ah (φh , φh ) = 0 implies that Γh ∩T
2 |∇φh | + φ2h dσh = 0
for all T ∈ ThI .
Since H n (T ∩ Γh ) > 0 for T ∈ ThI we infer that ∇φh = 0 and hence φh is constant on these elements. Using again that H n (T ∩ Γh ) > 0 we deduce that φh = 0 on each T ∈ ThI so that φh ≡ 0 in VhI . Hence (3.5) has a unique solution uh ∈ VhI and (3.6)
1
uh h = ∇uh 2L2 (Γh ) + uh 2L2 (Γh ) 2 ≤ c f e L2 (Γh ) ≤ c f L2 (Γ) .
3.3. Error analysis. The following error bounds hold. Theorem 3.2. Let u be the solution of (1.1) and uh the solution of the finite element scheme (3.5). Then (3.7)
ue − uh L2 (Γh ) + h ∇(ue − uh ) L2 (Γh ) ≤ ch2 f L2 (Γ) .
2147
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
Proof. In view of the definition of · h , (2.27), and Lemma 3.1 we have for eh := ue − uh 1 2 2
eh L2 (Γh ) + ∇eh L2 (Γh ) 2 1 2 2 ≤ 2 ue − Ih ue L2 (Γh ) + ∇(ue − Ih ue ) L2 (Γh ) 2
(3.8)
+ sup φh ∈VhI
|ah (ue , φh ) − lh (φh )|
φh h
≤ ch u H 2 (Γ) + sup
φh ∈VhI
|ah (ue , φh ) − lh (φh )| .
φh h
In order to estimate the second term we choose φh ∈ VhI , then for ϕh := φh ◦ Fh−1 , ah (ue , φh ) − lh (φh ) = ah (ue , φh ) − a(u, ϕh ) + l(ϕh ) − lh (φh ) ≡ I + II. Using the transformation rule and (2.11) we obtain ∇Γ u · ∇Γ ϕh + uϕh dσ = (∇Γ u) ◦ p · (∇Γ ϕh ) ◦ p + (u ◦ p) (ϕh ◦ p) μh dσh Γ Γ h (I + ΦA)−1 ∇ue · (I + ΦA)−1 ∇ϕeh + ue ϕeh μh dσh . (3.9) = Γh
Since ϕeh (x) = ϕh (p(x)) = φh (Fh−1 (p(x))) we derive ∇ϕeh (x) = [Dp(x)]T [DFh−1 (p(x))]T ∇φh (Fh−1 (p(x))) = [Dp(x)]T [DFh (Fh−1 (p(x)))]−T ∇φh (Fh−1 (p(x))). We infer from (2.22) that (3.10)
(DFh )−T = I −
1 ∇ηh ⊗ ν + Bh |∇Φ|
with |Bh | ≤ ch2 ,
where ηh = Ih Φ − Φ. It follows from (2.19) that Fh−1 (p(x)) = x, x ∈ Γh , which together with (2.9) implies
1 e ∇ϕh = (I − ν ⊗ ν) I − ∇ηh ⊗ ν ∇φh + qh on Γh , |qh | ≤ ch2 |∇φh |. |∇Φ| Taking into account that ∇ue · ν = 0 we therefore have ∇ue · ∇ϕeh = ∇ue · ∇φh −
1 (∇ue · ∇ηh )(∇φh · ν) + ∇ue · qh |∇Φ|
on Γh .
Inserting this relation into (3.9) and recalling the definition of ah we find that e ∇u · ∇φh − μh (I + ΦA)−T (I + ΦA)−1 ∇ue · ∇ϕeh + |(μh − 1)ue ϕeh | dσh |I| ≤ Γh 2 e e ≤ ch u h φh h + ch u h |(∇φh · ν)| dσh , Γh
2148
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
where we used (2.15), (2.23), and the fact that ϕeh = φh on Γh . Similarly, |II| = f ϕh dσ − f e φh dσh ≤ |1 − μh | |f e | |φh | dσh Γ
Γh
Γh
≤ ch2 f e L2 (Γh ) φh h ≤ ch2 f L2 (Γ) φh h . Combining these estimates with (2.24) we have (3.11)
2
|ah (u , φh ) − lh (φh )| ≤ ch f L2 (Γ) φh h + ch f L2(Γ) e
Γh
|∇φh · ν| dσh
for all φh ∈ VhI , which inserted into (3.8) yields
eh L2 (Γh ) + ∇eh L2 (Γh ) ≤ ch f L2 (Γ) .
(3.12)
In order to improve the L2 -error bound we employ the usual Aubin–Nitsche argument. Denote by w ∈ H 2 (Γ) the solution of the dual problem a(ϕ, w) = eh ϕ dσ for all ϕ ∈ H 1 (Γ) with eh = eh ◦ Fh−1 , Γ
which satisfies
w H 2 (Γ) ≤ c eh L2 (Γ) .
(3.13) We have in view of (1.3)
(3.14)
2 eh , w) = a( eh , w) − ah (eh , we )
eh L2 (Γ) = a( + ah (eh , we − Ih we ) + ah (ue , Ih we ) − lh (Ih we ) ≡ I + II + III.
Similarly as above we deduce with the help of (3.12) and Lemma 2.2 |I| ≤ ch eh h we h ≤ ch2 f L2 (Γ) w H 1 (Γ) . Next, Lemma 3.1 and (3.12) imply |II| ≤ eh h we − Ih we h ≤ ch2 f L2 (Γ) w H 2 (Γ) . Finally, (3.11), the fact that ∇we · ν = 0, and Lemma 3.1 yield |∇(Ih we − we ) · ν| dσh |III| ≤ ch2 f L2 (Γ) Ih we h + ch f L2 (Γ) Γh
2
≤ ch f L2 (Γ) w H 2 (Γ) . Inserting the above estimates into (3.14) and recalling (3.13) we obtain
eh L2 (Γ) ≤ ch2 f L2 (Γ) , which together with Lemma 2.2 completes the proof since eeh = eh on Γh .
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
2149
4. Narrow band method (NBM). 4.1. Setting up the method. Recalling the definition of Dh (2.17), we consider ThB = {T ∈ Th | H n+1 (T ∩ Dh ) > 0} and UhB = T. T ∈ThB
We define the finite element space VhB on the triangulation ThB by VhB = {φh ∈ C 0 (UhB ) | φh|T ∈ P1 (T ) for each T ∈ ThB }. Let us first examine the approximation error for the space VhB . Lemma 4.1. We have for each function z ∈ H 2 (Γ) that Ih z e ∈ VhB satisfies (4.1)
√ 1 √ z e − Ih z e L2 (Dh ) + h ∇(z e − Ih z e ) L2 (Dh ) ≤ ch2 z H 2 (Γ) . h
Proof. We infer from (2.14) and Lemma 2.2 that 1 e
z − Ih z e 2L2 (Dh ) + h ∇(z e − Ih z e ) 2L2 (Dh ) h
1 ≤
z e − Ih z e 2L2 (T ) + h ∇(z e − Ih z e ) 2L2 (T ) h T ∩Dh =∅ ≤ ch3
z e 2H 2 (T ) ≤ ch3 z e 2H 2 (D(1+c )h ) ≤ ch4 z 2H 2 (Γ) , 1
T ∩Dh =∅
since T ⊂ D(1+c1 )h for all T ∩ Dh = ∅ in view of (2.16). 4.2. The method. Let us write (1.5) in the form find uh ∈ VhB such that ah (uh , φh ) = lh (φh )
(4.2) where
for all φh ∈ VhB ,
1 ∇wh · ∇φh + wh φh |∇Ih Φ| dx, ah (wh , φh ) = 2h Dh 1 lh (φh ) = f e φh |∇Ih Φ| dx. 2h Dh
Note that the factors h1 in each of the above terms are there to aid the notation for the error analysis. In a similar way as for SIF one can verify that ah is positive definite on VhB × VhB . Hence, the finite element scheme (4.2) has a unique solution uh ∈ VhB which satisfies (4.3)
uh h =
1 2h
Dh
12 2 |∇uh | + u2h |∇Ih Φ| dx ≤ c f L2 (Γ) .
4.3. Error analysis. Before we prove our main error bound we formulate a technical lemma which will be helpful in the error analysis. Lemma 4.2. Suppose that u ∈ H 2 (Γ) is a solution of (1.1). Then, 1 1 |∇Ih Φ| −1 e e ah (u , φ) = dx + S, φ f φ ◦ Fh |∇Φ| dx + (∇ue · ∇ηh )(∇φ · ν) 2h Dh 2h Dh |∇Φ|
2150
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
for all φ ∈ H 1 (Dh ), where ηh = Ih Φ − Φ and |S, φ| ≤ Ch2 u H 2 (Γ) φ h . Proof. To begin, we derive from (2.12) and (1.1) that −
(4.4) where
1 ∇ · |∇Φ| ∇ue + ue = f e + R |∇Φ| ⎛
R(x) = −Φ(x) ⎝
(4.5)
n+1
blk (x)D l Dk u(p(x)) +
k,l=1
n+1
in Uδ , ⎞ ck (x)D k u(p(x))⎠ .
k=1
We multiply (4.4) by φ ◦ Fh−1 |∇Φ|, φ ∈ H 1 (Dh ), and integrate over Dh . Since ∂u h ∂ν = 0 on ∂D we obtain after integration by parts ∇ue · ∇(φ ◦ Fh−1 )|∇Φ| dx + ue φ ◦ Fh−1 |∇Φ| dx h h D D (4.6) −1 e = f φ ◦ Fh |∇Φ| dx + R φ ◦ Fh−1 |∇Φ| dx. e
Dh
Dh
Observing that ∇(φ ◦ Fh−1 ) = [(DFh )−T ◦ Fh−1 ]∇φ ◦ Fh−1 , the transformation rule and Lemma 2.1 imply that −1 e ∇u ·∇(φ◦Fh )|∇Φ| dx = ∇ue ◦Fh ·(DFh )−T ∇φ |∇Φ◦Fh | |det DFh | dx. I := Dh
Dh
Recalling (2.7) and (2.18) we have z e (x) = z(p(x)) = z(p(Fh (x)) = z e (Fh (x)),
(4.7)
from which we deduce by differentiation ∇z e ◦ Fh = (DFh )−T ∇z e ,
(4.8) so that
I=
(DFh )−T ∇ue · (DFh )−T ∇φ |∇Φ ◦ Fh | |det DFh | dx.
Dh
Recalling (3.10), we find with the help of ∇ue · ν = 0 that (DFh )−T ∇ue = ∇ue + Bh ∇ue ,
(DFh )−T ∇φ = ∇φ −
1 (∇φ · ν)∇ηh + Bh ∇φ, |∇Φ|
where |Bh | ≤ ch2 . Furthermore, Lemma 2.1 implies that |∇Φ ◦ Fh | |det DFh | = |∇Ih Φ| + γh , so that in conclusion e ∇u · ∇φ |∇Ih Φ| dx − I= Dh
Dh
where |γh | ≤ ch2 ,
(∇ue · ∇ηh )(∇φ · ν)
|∇Ih Φ| dx + Rh1 , φ, |∇Φ|
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
where
1 Rh , φ ≤ ch2 ∇ue
L2 (Dh )
2151
∇φ L2 (Dh ) ≤ ch3 u H 2 (Γ) φ h ,
in view of Lemma 2.2 and the definition of · h . Similarly, (4.7) and (2.21) yield −1 e u φ ◦ Fh |∇Φ| dx = ue φ |∇Ih Φ| dx + Rh2 , φ Dh
Dh
with Rh2 , φ ≤ ch3 u H 2 (Γ) φ h . Inserting the above identities into (4.6) and dividing by 2h we derive (4.9) 1 1 ah (ue , φ) = f e φ ◦ Fh−1 |∇Φ| dx + R φ ◦ Fh−1 |∇Φ| dx 2h Dh 2h Dh 1 1 |∇Ih Φ| 1 dx − Rh1 , φ − R2 , φ. + (∇ue · ∇ηh )(∇φ · ν) 2h Dh |∇Φ| 2h 2h h In order to rewrite the integral over Dh we note that F (·, s) maps Γ onto Γs = {Φ = s} and that dσs = (1 + O(s)) dσp , where dσs , dσp are the surface elements of Γs , Γ respectively. The coarea formula then yields for integrable g : Dh → R h h 1 dσs ds = g(x) dx = g(x) g(F (p, s))μ(p, s) dσp ds |∇Φ(x)| −h Γ Dh −h Γs 1 ≤ C |s| , |s| < h, p ∈ Γ. μ(p, s) − (4.10) where |∇Φ(F (p, s))| Hence,
Dh
R φ ◦ Fh−1 |∇Φ| dx
h
= −h Γ h
= −h
Γ
R ◦ F φ ◦ Fh−1 ◦ F |∇Φ ◦ F | μ dσp ds R ◦ F φ ◦ Fh−1 ◦ F dσp ds +
h
−h
Γ
r φ ◦ Fh−1 ◦ F dσp ds
≡ T1 + T2 ,
where r(p, s) = R(F (p, s)) μ(p, s) |∇Φ(F (p, s))| − 1 . In order to treat T1 we deduce from (4.5) and the fact that Φ(F (p, s)) = s that ⎛ ⎞ n+1 n+1 R(F (p, s)) = −s ⎝ blk (F (p, s))D l Dk u(p) + ck (F (p, s))D k u(p)⎠ . k,l=1
Since −
h −h h
−h
k=1
s ds = 0, the first term in T1 can be written as
n+1 Γ k,l=1
sDl Dk u(p) blk (F (p, s)) φ ◦ Fh−1 (F (p, s)) − blk (p) φ ◦ Fh−1 (p) dσp ds.
Treating the second term in T1 in the same way and observing that p = F (p, 0) we deduce with the help of the fundamental theorem of calculus that
12 2 2 5 −1 −1 ∇φ ◦ Fh dx + φ ◦ Fh ≤ ch3 u H 2 (Γ) φ h . |T1 | ≤ ch 2 u H 2 (Γ) Dh
2152
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
Next, we infer from (4.5) and (4.10) that | r (p, s)| ≤ cs2 |∇Γ u(p)| + DΓ2 u(p) , so that
5
|T2 | ≤ Ch 2 u H 2 (Γ)
Dh
12 φ ◦ F −1 2 dx ≤ Ch3 u H 2 (Γ) φ h . h
The result now follows from (4.9) together with the bounds on Rh1 and Rh2 . Theorem 4.3. Let u be the solution of (1.1) and uh the solution of the finite element scheme (4.2). Then
u − uh L2 (Γh ) + h e
(4.11)
1 2h
12 |∇(u − uh )| |∇Ih Φ| dx ≤ ch2 f L2 (Γ) . 2
e
Dh
Proof. Let us write eh := ue − uh . We infer from (2.27) and Lemma 4.1 that
eh h ≤ ch u H 2 (Γ) + sup
(4.12)
φh ∈VhB
|ah (ue , φh ) − lh (φh )| .
φh h
The second term on the right-hand side can be estimated with the help of Lemma 4.2. The transformation rule together with (4.7) yields 1 1 f e φh ◦ Fh−1 |∇Φ| dx = f e ◦ Fh φh |∇Φ ◦ Fh | |det DFh | dx, 2h Dh 2h Dh so that we deduce from Lemma 4.2 1 |∇Ih Φ| ah (ue , φh ) − lh (φh ) = (∇ue · ∇ηh ) (∇φh · ν) 2h Dh |∇Φ| 1 + f e φh |∇Φ ◦ Fh | |detDFh | − |∇Ih Φ| dx + S, φh . 2h Dh Using (2.21), (2.15), (2.24), and Lemmas 2.1 and 2.2 we infer that for φh ∈ VhB |ah (u , φh ) − lh (φh )| ≤ c ∇u L2 (Dh ) e
e
12 |∇φh · ν| dx 2
Dh
+ ch f e L2 (Dh ) φh L2 (Dh ) + ch2 u H 2 (Γ) φh h
12 1 2 ≤ ch f L2 (Γ) |∇φh · ν| dx + ch2 f L2 (Γ) φh h , 2h Dh
(4.13)
so that (4.12) implies the following intermediate result: (4.14)
eh h =
1 2h
Dh
12 2 |∇eh | + e2h |∇Ih Φ| dx ≤ ch f L2 (Γ) .
In order to improve the L2 -error bound we define eh := eh ◦ Fh−1 as well as h (p) := 1 E 2h
h
−h
eh (F (p, s)) ds,
p ∈ Γ,
2153
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
with F as above. We denote by w ∈ H 2 (Γ) the unique solution of h −ΔΓ w + w = E which satisfies
w H 2 (Γ) ≤ c E h
(4.15)
on Γ,
L2 (Γ)
.
Similarly to (4.4) the extension we solves −
1 e + R ∇ · |∇Φ| ∇we + we = E h |∇Φ|
in Uδ ,
is obtained from (4.5) by replacing u by w. Using the transformation rule where R together with (4.10) we obtain h h 2 1 1 e ◦ F eh ◦ F |∇Φ ◦ F | μ dσp ds Eh eh ◦ F dσp ds = E Eh 2 = 2h −h Γ 2h −h Γ h L (Γ) h 1 h E + eh ◦ F (1 − |∇Φ ◦ F | μ) dσp ds 2h −h Γ 1 e eh ◦ F −1 |∇Φ| dx E = h 2h Dh h h 1 h + E eh ◦ F (1 − |∇Φ ◦ F | μ) dσp ds. 2h −h Γ The first term can be rewritten with the help of Lemma 4.2 (applied to w instead of u) to give 2 |∇Ih Φ| eh − 1 dx (∇we · ∇ηh )(∇eh · ν) Eh 2 = ah (we , eh ) − S, 2h Dh |∇Φ| L (Γ) h 4 1 h eh ◦ F (1 − |∇Φ ◦ F | μ) dσp ds ≡ + E Ik . 2h −h Γ k=1
In view of Lemma 4.1, (4.13), the fact that ∇we · ν = 0, and (4.14) we have eh | |I1 | + |I2 | ≤ |ah (we − Ih we , eh )| + |ah (ue , Ih we ) − lh (Ih we )| + |S, ≤ ch w H 2 (Γ) eh h + ch2 f L2 (Γ) Ih we h
12 1 e e 2 + ch f L2 (Γ) |∇(Ih w − w ) · ν| dx + ch2 w H 2 (Γ) eh h 2h Dh ≤ ch2 f L2 (Γ) w H 2 (Γ) . Furthermore, (2.15), (4.10), and (4.14) imply
e 2 |I3 |+|I4 | ≤ ch w h + Eh
eh h ≤ ch f L2 (Γ) w H 2 (Γ) + E h L2 (Γ)
so that we obtain together with (4.15) Eh 2 ≤ ch2 f L2 (Γ) . L (Γ)
L2 (Γ)
,
2154
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
Next, since F (p, 0) = p we may write for p ∈ Γ h (p) − eh (p) = 1 E 2h
h
−h
s
∇ eh (F (p, τ )) ·
0
∂F (p, τ ) dτ ds, ∂s
and hence we obtain with the help of (4.14) Eh − eh
L2 (Γ)
√ ≤c h
12 |∇ eh | dx ≤ ch ∇eh h ≤ ch2 f L2 (Γ) . 2
Dh
In conclusion we deduce that
eh L2 (Γh ) ≤ c eh L2 (Γ) ≤ E h h−e
L2 (Γ)
+ E h
L2 (Γ)
≤ ch2 f L2 (Γ)
and the theorem is proved. 5. A hybrid method for equations on evolving surfaces. 5.1. The setting. The aim of this section is to combine ideas employed in sections 3 and 4 for the stationary problem in order to develop a finite element method for an advection–diffusion equation on an evolving hypersurface in which there is an underlying conservation law with a diffusive flux, [18]. More precisely, let (Γ(t))t∈[0,T ] be a family of compact, connected smooth hypersurfaces embedded in Rn+1 for n = 1, 2. We suppose that Γ(t) = {x ∈ N (t) | Φ(x, t) = 0},
where ∇Φ(x, t) = 0, x ∈ N (t),
and N (t) is an open neighborhood of Γ(t). We assume that N (t) is chosen so small that we can construct the function p(·, t) as in section 2.1. Given a velocity field v(·, t) : Γ(t) → Rn+1 , not necessarily in the normal direction, we then consider the following initial value problem (5.1a) Γ(t) × {t}, ∂ • u + u∇Γ · v − ΔΓ u = f on t∈(0,T )
(5.1b)
u(·, 0) = u0
on Γ(0).
Here, ∂ • η denotes the material derivative of a function η : which is given by
t∈(0,T )
Γ(t) × {t} → R
∂ • η = ∂t η + v · ∇η if η is extended into a neighborhood of t∈(0,T ) Γ(t) × {t}. 5.2. The method. In order to discretize the above problem we choose a partition 0 = t0 < t1 < · · · < tN = T of [0, T ] with τm := tm+1 − tm , m = 0, . . . , N − 1, and τ := maxm=0,...,N −1 τm . Also, let Th be an unfitted regular triangulation with mesh size h of a region containing N (t), t ∈ [0, T ]. For m = 0, 1, . . . , N we set Γm h = {x ∈ N (tm ) | Ih Φ(x, tm ) = 0}
and
Dhm = {x ∈ N (tm ) | |Ih Φ(x, tm )| < h},
as well as Thm := {T ∈ Th | H n+1 (T ∩ Dhm ) > 0}
and
Uhm :=
T ∈Thm
T.
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
2155
Here we assume that 0 < h ≤ h0 , where h0 is chosen so small that there exists c0 , c1 > 0 such that N (t) × {t}. c0 ≤ |∇Ih Φ(x, t)| ≤ c1 , (x, t) ∈ t∈(0,T )
Finally, we introduce Vhm = {φh ∈ C 0 (Uhm ) | φh |T ∈ P1 (T ) for each T ∈ Thm }. In what follows we shall frequently use the abbreviation z m (x) := z(x, tm ). In order to motivate our method we fix m ∈ {0, 1, . . . , N − 1} and let Ψ be the solution of Ψt (x, t) + DΨ(x, t)v e (x, t) = 0,
Ψ(x, tm+1 ) = x,
where v e (x, t) := v(p(x, t), t). For a sufficiently smooth function ϕ : N (tm+1 ) → R we define η(x, t) := ϕ(Ψ(x, t)). Clearly, η(·, tm+1 ) = ϕ and a short calculation shows that ∂ • η = 0. Assuming that u is a solution of (5.1a) we obtain with the help of the Leibniz formula and integration by parts • d ∂ (uη) + uη∇Γ · v dσ uη dσ|t=tm+1 = dt Γ(t) Γ(t ) m+1 ϕ ∂ • u + u∇Γ · v dσ = Γ(t ) m+1 ϕΔΓ u + ϕf dσ = Γ(tm+1 )
=−
Γ(tm+1 )
∇Γ u · ∇Γ ϕ dσ +
f ϕ dσ. Γ(tm+1 )
Since Ψ(·, tm+1 ) ≡ id, a Taylor expansion shows that Ψ(x, tm ) = Ψ(x, tm+1 − τm ) ≈ x − τm Ψt (x, tm+1 ) = x + τm DΨ(x, tm+1 )v e,m+1 (x) = x + τm v e,m+1 (x). Thus we may approximate the left-hand side of the above relation by d 1 uη dσ|t=tm+1 ≈ um+1 ϕ dσ − um ϕ(· + τm v e,m+1 ) dσ . dt Γ(t) τm Γ(tm+1 ) Γ(tm ) The above calculations motivate the following scheme, in which we use the narrow m band approach in order to discretize the elliptic part. Given um h ∈ Vh , m = 1, . . . , N − m+1 m+1 1, find uh ∈ Vh such that e,m+1 um+1 φ dσ − um ) dσh h h h φh (· + τm v h (5.2)
Γm+1 h
τm + 2h
Γm h
m+1 Dh
∇um+1 h
· ∇φh |∇Ih Φ
m+1
| dx = τm
Γm+1 h
f e,m+1 φh dσh
for all φh ∈ Vhm+1 . Here, u0h = Ih u0 . Existence and uniqueness of um+1 follows in a h similar way as for NBM in the elliptic case.
2156
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
5.3. Mass conservation. An important property of solutions of (5.1a) is conservation of mass in the case that Γ(t) f (·, t) dσ = 0. The following lemma shows that our numerical scheme preserves this property under some mild constraints on the discretization parameters. In fact this discrete conservation law is a remarkable property of the scheme and relies on the use of both a sharp interface and a narrow band approach. m Lemma 5.1. Let um h ∈ Vh , m = 1, . . . , N , be the solutions of (5.2) in the case √ f e,m+1 dσh = 0, m = 0, . . . , N −1. Then provided that 0 < h ≤ h1 and τ ≤ γ h Γm+1 h
(5.3)
Γm h
um h dσh =
Γ0h
u0h dσh .
Proof. Let us first observe that m+1 {x + τm v e,m+1 (x) | x ∈ Γm , h } ⊂ Uh
(5.4)
m = 0, . . . , N − 1,
provided that h, τ are sufficiently small. To see this, let x ∈ Γm h and choose an element T ∈ Th such that x ∈ T . Then, Φm+1 (x + τm v e,m+1 (x)) = Φm (x) + τm ∇Φm (x) · v e,m+1 (x) + τm Φt (x, tm ) + Rm (x), 2 where |Rm (x)| ≤ cτm . Observing that Φt + ∇Φ · v = 0 on
t∈(0,T )
Γ(t) × {t} we write
∇Φm (x) · v e,m+1 (x) + Φt (x, tm ) = ∇Φm (x) · v m+1 (pm+1 (x)) + Φt (x, tm ) = ∇Φm (x) − ∇Φm+1 (pm+1 (x)) · v m+1 (pm+1 (x)) + Φt (x, tm ) − Φt (pm+1 (x), tm+1 ), so that m+1 2 Φ (x + τm v e,m+1 (x)) ≤ |Φm (x)| + cτm x − pm+1 (x) + cτm 2 2 ≤ |Φm (x)| + cτm Φm+1 (x) + cτm ≤ c(h2 + τm ), in view of (2.5) and since |Φm (x)| = |Φm (x) − Ih Φm (x)| ≤ ch2 . As a result, (Ih Φm+1 )(x + τm v e,m+1 (x)) ≤ Ih Φm+1 − Φm+1
L∞
2 2 + c(h2 + τm ) ≤ c(h2 + τm )
c 2014 Society for Industrial and Applied Mathematics
UNFITTED FINITE ELEMENT METHODS USING BULK MESHES FOR SURFACE PARTIAL DIFFERENTIAL EQUATIONS∗ KLAUS DECKELNICK† , CHARLES M. ELLIOTT‡ , AND THOMAS RANNER§ Abstract. In this paper, we define new unfitted finite element methods for numerically approximating the solution of surface partial differential equations using bulk finite elements. The key idea is that the n-dimensional hypersurface, Γ ⊂ Rn+1 , is embedded in a polyhedral domain in Rn+1 consisting of a union, Th , of (n + 1)-simplices. The unifying feature of the methodological approach is that the finite element approximating space is based on continuous piecewise linear finite element functions on the bulk triangulation Th which is independent of Γ. Our first method is a sharp interface method (SIF ) which uses the bulk finite element space in an approximating weak formulation obtained from integration on a polygonal approximation, Γh , of Γ. The full gradient is used rather than the projected tangential gradient and it is this which distinguishes SIF from the method of [M. A. Olshanskii, A. Reusken, and J. Grande, SIAM J. Numer. Anal., 47 (2009), pp. 3339–3358]. The second method is a narrow band method (NBM ) in which the region of integration is a narrow band of width O(h). NBM is similar to the method of [K. Deckelnick et al., IMA J. Numer. Anal., 30 (2010), pp. 351–376] but again the full gradient is used in the discrete weak formulation. The a priori error analysis in this paper shows that the methods are of optimal order in the surface L2 and H 1 norms and have the advantage that the normal derivative of the discrete solution is small and converges to zero. Our third method combines bulk finite elements, discrete sharp interfaces, and narrow bands in order to give an unfitted finite element method for parabolic equations on evolving surfaces. We show that our method is conservative so that it preserves mass in the case of an advection-diffusion conservation law. Numerical results are given which illustrate the rates of convergence. Key words. unfitted finite elements, cut cells, error analysis, narrow band, sharp interface, elliptic and parabolic surface equations AMS subject classifications. 35R01, 65N30, 65N15, 65M60 DOI. 10.1137/130948641
1. Introduction. In this article we propose and analyze numerical methods based on bulk finite element meshes for the following model elliptic equation on a stationary surface. Model elliptic equation on stationary surface. Let Γ be a smooth hypersurface in Rn+1 and f ∈ L2 (Γ). We seek solutions u : Γ → R of (1.1)
−ΔΓ u + u = f
on Γ.
The methods can be extended in natural ways to deal with variable coefficients and nonlinearities. The approach may be extended to the following advection-diffusion equation on a moving surface. ∗ Received
by the editors December 10, 2013; accepted for publication (in revised form) May 30, 2014; published electronically August 19, 2014. http://www.siam.org/journals/sinum/52-4/94864.html † Institut f¨ ur Analysis und Numerik, Otto-von-Guericke-Universit¨ at Magdeburg, 39106 Magdeburg, Germany ([email protected]). ‡ Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK (C.M.Elliott@warwick. ac.uk). § Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK. Current address: School of Computing, University of Leeds, LS2 9JT, UK ([email protected]). The work of this author was supported by an EPSRC Ph.D. studentship (grant EP/P504333/1 and EP/P50516X/1) and the Warwick Impact Fund. 2137
2138
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
Model parabolic equation on evolving surface. Let {Γ(t)} be a family of smooth hypersurfaces in Rn+1 for t ∈ [0, T ]. Denoting by ∂ • u the material derivative of u and v the velocity of Γ(t) (see section 5 for notation), we seek solutions u : t Γ(t) × {t} of the advection-diffusion equation (1.2a) Γ(t) × {t}, ∂ • u + u∇Γ · v − ΔΓ u = f on t∈(0,T )
(1.2b)
u(·, 0) = u0
on Γ(0).
Surface partial differential equations or partial differential equations (PDEs) on manifolds arise in a wide variety of applications in materials science, fluid dynamics, and biology, [5, 32, 24, 29, 25, 2, 28, 26]. Computational approaches include surface finite elements on triangulated surfaces [17, 19, 18, 15, 14, 23, 22], bulk finite element or finite difference meshes for the approximation of implicit surface formulations [9, 30, 10, 21, 11], bulk finite element or finite difference meshes on narrow bands [12, 38], and bulk finite element meshes and sharp interface weak forms [34, 33, 16]. An important feature of the methods cited above is the avoidance of charts both in the problem formulation and the numerical methods. For example, the surface finite element method is based simply on triangulated surfaces and requires the geometry solely through the knowledge of the vertices of the triangulation, whereas methods based on implicit surfaces require only the level set function Φ which encodes all the necessary geometry. Another feature of some of these methods is the use of unfitted bulk meshes. Here we use the terminology unfitted finite element methods (sometimes called cut cell methods) when the underlying meshes that form the computational domain are not fitted to the domain in which the PDE holds. The motivation for using finite element spaces on meshes not fitting to the domain came from the desire to solve free or moving boundary problems. Such methods were introduced in [3, 4] for elliptic equations in curved domains; see also [31, 8, 27]. In this setting we are concerned with bulk meshes independent of the surface. The new methods. The new unfitted finite element methods for surface elliptic equations proposed in this paper are variants of the bulk finite element approaches using a sharp interface or a narrow band. The new scheme for advection diffusion on an evolving surface is a hybrid of these. In the following we sketch the main ideas of these methods describing the details in sections 3–5. Sharp interface method (SIF ). Given an interpolation Γh of Γ, we use a bulk finite element space VhI of the form VhI = {φh ∈ C 0 (UhI ) | φh|T ∈ P1 (T ) for each T ∈ ThI }, where ThI is a set of elements which intersect Γh and UhI = T ∈T I T ; see section 3. h The discrete scheme approximating the model elliptic equation (1.1) is find uh ∈ VhI such that ∇uh · ∇φh + uh φh dσh = f e φh dσh for all φh ∈ VhI , (1.3) Γh
Γh
where f e is an extension of f . The method is related to the following method of Olshanskii, Reusken, and Grande, introduced in [34]: find uh ∈ VhΓ such that (1.4) ∇Γh uh · ∇Γh φh + uh φh dσh = f e φh dσh for all φh ∈ VhΓ . Γh
Γh
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
2139
There are two significant differences. Note the use of the full gradient in (1.3) as opposed to the tangential gradient. This gives control over the normal derivative of the finite element solution which is lacking in (1.4). Another difference relates to the use of the finite element space VhΓ , which essentially consists of the traces on Γh of elements in VhI . While VhI has a natural basis, this does not seem to be the case for VhΓ . The “standard basis” of finite element hat functions is only a spanning set for VhΓ . Narrow band method (NBM ). We use the bulk finite element space VhB on the triangulation ThB , VhB = {φh ∈ C 0 (UhB ) | φh|T ∈ P1 (T ) for each T ∈ ThB }. Here ThB consists of those triangles intersecting a narrow band domain D h defined by the ±h level sets of an interpolated level set function Ih Φ and UhB = T ∈T B T . h The discrete scheme approximating the model elliptic equation (1.1) is find uh ∈ VhB such that (1.5)
∇uh · ∇φh + uh φh |∇Ih Φ| dx =
Dh
f e φh |∇Ih Φ| dx Dh
for all φh ∈ VhB .
This is similar to the method in [12] except that NBM uses the full instead of projected gradients thus avoiding the resulting degeneracy. As a result we are able to prove an optimal L2 -error bound which was not obtained for the method in [12]. It is also the case that the normal derivative of the discrete solution converges to zero. Hybrid unfitted evolving surface method. The discrete problem approxim+1 m ∈ Vhm+1 such that mating (1.2) is given um h ∈ Vh , m = 0, . . . , N − 1, find uh
(1.6)
Γm+1 h
um+1 φh dσh − h
τm + 2h
m+1 Dh
∇um+1 h
Γm h
e,m+1 um ) dσh h φh (· + τm v
· ∇φh ∇Ih Φm+1 dx = τm
Γm+1 h
f e,m+1 φh dσh
for all φh ∈ Vhm+1 . Here v e,m denotes an extension of the surface velocity at time level m. We use time step labeled analogues of the notation for NBM; see section 5 for the details. Here, u0h is appropriate initial data. Because of this combination of narrow band and sharp interface discretization, under some mild constraints on the discretization parameters (see section 5) our numerical scheme preserves the important property that solutions of (1.2) conserve mass in the case that f ≡ 0. Outline. The paper is organized as follows: In section 2 we introduce our notation and collect some auxiliary results. In sections 3 and 4 we present and analyze unfitted methods for the model elliptic equation (1.1). In section 5 we describe how a combination of these two approaches can be used to calculate solutions of the advection-diffusion equation on evolving hypersurfaces, (1.2). Details of the implementation and several numerical examples illustrating the orders of convergence are presented in section 6.
2140
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
2. Preliminaries. 2.1. Surface calculus. Let Γ be a connected compact smooth hypersurface embedded in Rn+1 (n = 1, 2). We assume that there exists a smooth function Φ : U → R such that Γ = {x ∈ U | Φ(x) = 0}
∇Φ(x) = 0, x ∈ U,
and
where U is an open neighborhood of Γ. For a function z : Γ → R we define its tangential gradient by (2.1) ∇Γ z(p) := ∇ z (p) − ∇ z (p) · ν(p) ν(p), p ∈ Γ, where z : U → R is an arbitrary smooth extension of z to U and ν(x) =
∇Φ(x) |∇Φ(x)|
is a unit vector to the level sets of Φ. It can be shown that ∇Γ z(p) is independent of the particular choice of z. We denote by Di z, 1 ≤ i ≤ n + 1, the components of ∇Γ z. Furthermore, we let ΔΓ z = ∇Γ · ∇Γ z =
n+1
Di Di z
i=1
be the Laplace–Beltrami operator of z. In what follows it will be convenient to use special coordinates which are adapted to Φ. Consider for p ∈ Γ the system of ODEs γp (s) =
(2.2)
∇Φ(γp (s)) , |∇Φ(γp (s))|2
γp (0) = p.
It can be shown that there exists δ > 0 so that the solution γp of (2.2) exists uniquely on (−δ, δ) uniformly in p ∈ Γ, so that we can define the mapping F : Γ × (−δ, δ) → Rn+1 by (2.3)
F (p, s) := γp (s),
p ∈ Γ, |s| < δ.
d Since ds Φ(γp (s)) = 1 and γp (0) = p ∈ Γ, we infer that Φ(γp (s)) = s, |s| < δ, and hence that x = F (p, s) implies that |Φ(x)| < δ. As a result, we deduce that F is a diffeomorphism of Γ × (−δ, δ) onto Uδ := {x ∈ U | |Φ(x)| < δ} with inverse
F −1 (x) = (p(x), Φ(x)),
(2.4)
x ∈ Uδ ,
where p : Uδ → Rn+1 satisfies p(x) ∈ Γ, x ∈ Uδ . For later purposes it is convenient to expand p and its derivatives in terms of Φ. Let us fix x ∈ Uδ and define the function η(τ ) := F (p(x), (1 − τ )Φ(x)), τ ∈ [0, 1]. Since
∂F ∂s
(p, s) = γp (s) we have
∇Φ(γp(x) ((1 − τ )Φ(x))) η (τ ) = −Φ(x)γp(x) ((1 − τ )Φ(x)) = −Φ(x) . ∇Φ(γp(x) ((1 − τ )Φ(x))) 2
2141
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
Observing that γp(x) (Φ(x)) = F (p(x), Φ(x)) = x and using similar arguments to calculate η (τ ) we find that ηk (0) = −Φ(x)
Φxk (x)
2, |∇Φ(x)|
n+1 2Φxk (x)Φxr (x) Φxl (x)Φxl xr (x) 2 ηk (0) = Φ(x) , δkr − 2 4 |∇Φ(x)| |∇Φ(x)| l,r=1
k = 1, . . . , n + 1. Since η(1) = p(x), η(0) = x we deduce with the help of Taylor’s theorem that (2.5)
n+1 Φxk (x) 1 2Φxk (x)Φxr (x) Φxl (x)Φxl xr (x) 2 pk (x) = xk − Φ(x) + Φ(x) δkr − |∇Φ(x)|2 2 |∇Φ(x)|2 |∇Φ(x)|4 l,r=1
+ Φ(x)3 rk (x),
k = 1, . . . , n + 1,
where rk are smooth functions. In a similar way we may write ∇Φ(x) = ∇Φ(p(x)) + Φ(x)G(x),
(2.6)
1 where G(x) = 0 D2 Φ(F (p(x), τ Φ(x))) ∂F ∂s (p(x), τ Φ(x)) dτ . Let us next use the function p in order to define a particular extension of z : Γ → R: x ∈ Uδ .
z e (x) := z(p(x)),
(2.7)
Since p(F (p(x), s)) = p(x) we deduce that s → z e (F (p(x), s)) is independent of s and thus ∇z e (x) · ν(x) = 0,
(2.8)
x ∈ Uδ .
In order to express the derivatives of z e in terms of the tangential derivatives of z we first deduce from (2.5) that pk,xi (x) =δik −
Φxk (x)Φxi (x) |∇Φ(x)|
+ Φ(x)Φxi (x)
2
n+1
−
Φ(x)Φxk xi (x)
|∇Φ(x)|
δkr −
l,r=1
2
+ 2Φ(x)Φxk (x)
2Φxk (x)Φxr (x) |∇Φ(x)|
2
n+1
Φxl (x)Φxl xi (x)
l=1
|∇Φ(x)|
Φxl (x)Φxl xr (x) |∇Φ(x)|4
4
+ Φ(x)2 αik (x).
Combining this relation with (2.6) we deduce that (2.9) pk,xi (x) = δik − νi (p(x))νk (p(x)) + aik (x)Φ(x), pk,xi xj (x) = − (2.10)
+
Φxi (x)Φxk xj (x) Φxj (x)Φxk xi (x) − |∇Φ(x)|2 |∇Φ(x)|2 n+1 Φxi (x)Φxj (x) Φxl (x)Φxk xl (x) + βkij (x)νk (p(x)) + γkij (x)Φ(x), |∇Φ(x)|2 |∇Φ(x)|2 l=1
2142
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
where aik , βkij , γkij are smooth functions. Differentiating (2.7) and using (2.9), (2.10) n+1 as well as the fact that k=1 Dk z(p(x))νk (p(x)) = 0 we obtain (2.11) ∇z e (x) = I + Φ(x)A(x) ∇Γ z(p(x)), (2.12) 1 ∇ · |∇Φ(x)|∇z e (x) |∇Φ(x)|
⎛
= (ΔΓ z)(p(x)) + Φ(x) ⎝
n+1
blk (x)D l Dk z(p(x)) +
k,l=1
n+1
⎞ ck (x)D k z(p(x))⎠ ,
k=1
where A = (aik ), blk , and ck are again smooth. 2.2. Bulk finite element space and inequalities. In what follows we assume that U is polyhedral. Let (Th )0 0. As a consequence,
Φ − Ih Φ L∞ (U) + h ∇(Φ − Ih Φ) L∞ (U) ≤ Ch2 ,
so that we may assume that there exist constants c0 , c1 such that c0 ≤ |∇Ih Φ(x)| ≤ c1 ,
(2.16)
x ∈ U, 0 < h ≤ h0 .
Let us next define (2.17)
Γh := {x ∈ U | Ih Φ(x) = 0}
and
Dh := {x ∈ U | |Ih Φ(x)| < h}
as approximations of the given hypersurface Γ and the neighborhood Dh := {x ∈ U | |Φ(x)| < h} ; see Figure 1, for example. Note that Γh is a polygon whose facets are line segments if n = 1 and a polyhedral surface whose facets consist of triangles or quadrilaterals if n = 2. The corresponding decomposition of Γh is in general not shape regular and can have arbitrary small elements. Furthermore, we introduce Fh : U → Rn+1 by Fh (x) := F (p(x), Ih Φ(x)),
x ∈ U,
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
2143
Fig. 1. A cartoon of the domains of the sharp interface (left) and the narrow band (right) method. The surface Γ, resp., the set D h is displayed in red, the approximations Γh , resp., Dh in blue, and the domains UhI , UhB in gray.
where F was defined in (2.3). From the properties of F we infer that (2.18)
p(Fh (x)) = p(x)
and
Φ(Fh (x)) = Ih Φ(x) if Fh (x) ∈ Uδ , Fh (x) = p(x) if x ∈ Γh .
(2.19)
Lemma 2.1. There exists 0 < h1 ≤ h0 such that for 0 < h ≤ h1 the mapping Fh : Dh → Dh := {x ∈ U | |Φ(x)| < h} is bi-Lipschitz with Fh (Γh ) = Γ. Furthermore,
Fh − Id L∞ (U) + h DFh − I L∞ (U) ≤ ch2 , |det DFh | − |∇Ih Φ| ≤ ch2 . |∇Φ| ∞
(2.20) (2.21)
L
(U)
Proof. Since F (p(x), Φ(x)) = x we deduce with the help of (2.15) |Fh (x) − x| = |F (p(x), Ih Φ(x)) − F (p(x), Φ(x))| ≤ c|Ih Φ(x) − Φ(x)| ≤ ch2 . Differentiating the relation Fi (p(x), Φ(x)) = xi , i = 1, . . . , n + 1, we obtain n+1
Dk Fi (p(x), Φ(x))pk,xj (x) +
k=1
∂Fi (p(x), Φ(x))Φxj (x) = δij , ∂s
i, j = 1, . . . , n + 1,
and hence Fhi,xj (x) =
n+1
Dk Fi (p(x), Ih Φ(x))pk,xj (x) +
k=1
∂Fi (p(x), Ih Φ(x))(Ih Φ)xj (x) ∂s
∂Fi (p(x), Φ(x)) Ih Φ − Φ xj (x) ∂s n+1 + Dk Fi (p(x), Ih Φ(x)) − Dk Fi (p(x), Φ(x)) pk,xj (x)
= δij + (2.22)
k=1
∂Fi ∂Fi (p(x), Ih Φ(x)) − (p(x), Φ(x)) (Ih Φ)xj (x) ∂s ∂s Φxi (x) Ih Φ − Φ xj (x) + rij (x), = δij + |∇Φ(x)|2 +
2144
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
where |rij (x)| ≤ ch2 in view of (2.15). This implies (2.20). In particular we deduce that Fh is bi-Lipschitz provided that h is sufficiently small, whereas the properties Fh (Dh ) = Dh and Fh (Γh ) = Γ follow from (2.18). Finally we deduce from (2.22) that ∇Φ ∇Φ · ∇Ih Φ · ∇(Ih Φ − Φ) + ch = + ch |∇Φ|2 |∇Φ|2 2 ∇Φ |∇Ih Φ| |∇Ih Φ| 1 ∇Ih Φ |∇Ih Φ| − − + ch = + dh , = |∇Φ| 2 |∇Ih Φ| |∇Φ| |∇Φ| |∇Φ|
|det DFh | = 1 +
where |ch | , |dh | ≤ ch2 proving (2.21). We introduce μh : Γh → R via dσ(p(x)) = μh (x) dσh (x). It is well known (see Proposition 2.1 in [15], (3.37) in [34]) that |1 − μh | ≤ ch2
(2.23)
on Γh .
Using the properties of Fh together with the coarea formula and (2.9),(2.10), (2.11), (2.23) one can prove the following result on the equivalence of certain norms. Lemma 2.2. There exist constants c1 , c2 > 0 which are independent of h, such that for all z ∈ H 1 (Γ) c1 z e L2 (Γh ) ≤ z L2 (Γ) ≤ c2 z e L2 (Γh ) , 1 1 c1 √ z e L2 (Dh ) ≤ z L2 (Γ) ≤ c2 √ z e L2 (Dh ) , h h e c1 ∇z L2 (Γh ) ≤ ∇Γ z L2 (Γ) ≤ c2 ∇z e L2 (Γh ) , 1 1 c1 √ ∇z e L2 (Dh ) ≤ ∇Γ z L2 (Γ) ≤ c2 √ ∇z e L2 (Dh ) . h h If in addition z ∈ H 2 (Γ) then 1 c1 √ D2 z e L2 (D ) ≤ z H 2 (Γ) . h h 2.3. Variational form of elliptic equation and Strang’s second lemma. It is well known [1] that for every f ∈ L2 (Γ) there exists a unique solution u ∈ H 2 (Γ) of (1.1) which satisfies
u H 2 (Γ) ≤ c f L2 (Γ) .
(2.24)
Let us write (1.1) in weak form: (2.25)
a(u, ϕ) = l(ϕ)
where
a(w, ϕ) = Γ
for all ϕ ∈ H 1 (Γ),
∇Γ w · ∇Γ ϕ + wϕ dσ,
l(ϕ) =
f ϕ dσ. Γ
Next, suppose that Vh is a finite–dimensional space and V e := {v e | v ∈ H 1 (Γ)}. Assume that ah : (Vh + V e ) × (Vh + V e ) → R is a symmetric, positive semidefinite bilinear form which is, in addition, positive definite on Vh × Vh . Furthermore, let lh : Vh → R be linear. Then the approximate problem (2.26)
ah (uh , vh ) = lh (vh )
for all vh ∈ Vh
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
has a unique solution uh ∈ Vh . Introducing
v h := ah (v, v),
2145
v ∈ Vh + V e ,
we have by Strang’s second lemma (2.27)
ue − uh h ≤ 2 inf ue − vh h + sup vh ∈Vh
φh ∈Vh
|ah (ue , φh ) − lh (φh )| .
φh h
3. Sharp interface method (SIF ). 3.1. Setting up the method. Let us begin by observing that if T ∈ Th satisfies H n (T ∩ Γh ) > 0, then the following two cases can occur: (1) Γh ∩ int(T ) = ∅, in which case H n (∂T ∩ Γh ) = 0; (2) T ∩ Γh = ∂T ∩ Γh in which case T ∩ Γh is the face between two elements. We may now define a unique subset ThI ⊂ Th by taking all elements satisfying case 1 and in case 2 taking just one of the two elements T . The numerical method does not depend on which element is chosen. We may therefore conclude that there exists N ⊂ Γh with H n (N ) = 0 and a subset ThI ⊂ Th such that every x ∈ Γh \ N belongs to exactly one T ∈ ThI . We then define UhI = T. T ∈ThI
Clearly UhI ⊆ Uδ if h is small enough. We define the finite element space VhI by VhI = {φh ∈ C 0 (UhI ) | φh|T ∈ P1 (T ) for each T ∈ ThI }. Note that ∇φh is defined on Γh \ N in view of the definition of ThI . In particular the unit normal νh to Γh is given by νh =
(3.1)
∇Ih Φ |∇Ih Φ|
on Γh \ N,
and we use (3.1) in order to extend νh to UhI . Let us next turn to the approximation ¯δ ) so error for the space VhI . Note that for a function z ∈ H 2 (Γ) we have z e ∈ C 0 (U e that Ih z is well–defined. Lemma 3.1. Let z ∈ H 2 (Γ). Then (3.2)
z e − Ih z e L2 (Γh ) + h ∇(z e − Ih z e ) L2 (Γh ) ≤ ch2 z H 2 (Γ) .
Proof. We first observe that Theorem 3.7 in [34] yields (3.3)
z e − Ih z e L2 (Γh ) + h ∇Γh (z e − Ih z e ) L2 (Γh ) ≤ ch2 z H 2 (Γ) .
Hence, it remains to bound ∇(z e −Ih z e )·νh L2 (Γh ) . To do so, we start by considering an element T ∈ ThI . Then we see that 2 |∇(z e − Ih z e ) · νh | dσh T ∩Γh 2 2 e ≤2 |∇z · νh | dσh + 2 |∇(Ih z e ) · νh | dσh T ∩Γh T ∩Γh 2 2 e −1 ≤2 |∇z · (νh − ν)| dσh + ch(T ) |∇(Ih z e ) · νh | dx =: I1 + I2 , T ∩Γh
T
2146
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
in view of (2.8) and the fact that H n (T ∩ Γh ) ≤ ch(T )−1 H n+1 (T ). Note that by (3.1) and (2.15)
ν − νh L∞ (T )
(3.4)
∇Φ ∇Ih Φ − = ≤ ch(T ) |∇Φ| |∇Ih Φ| L∞ (T )
so that I1 ≤ ch2
T ∩Γh
2
|∇z e | dσh .
Furthermore, recalling (2.14) and using again (3.4) I2 ≤ ch(T )−1
T
2 2 2 |∇z e · (νh − ν)| + |∇(z e − Ih z e )| dx ≤ ch z e H 2 (T ) .
We use the bounds for I1 , I2 and sum over all elements T ∈ ThI , then apply Lemma 2.2 to see |∇(z e − Ih z e ) · νh |2 dσh ≤ ch2 ∇z e 2L2 (Γh ) + ch z e 2H 2 (Dc h ) ≤ ch2 z 2H 2 (Γ) , 1
Γh
since T ⊂ Dc1 h for all T ∈ ThI in view of (2.16). 3.2. The method. Let us write (1.3) in the form find uh ∈ VhI such that (3.5)
ah (uh , φh ) = lh (φh )
for all φh ∈ VhI ,
where ah (wh , φh ) =
Γh
∇wh · ∇φh + wh φh dσh ,
lh (φh ) =
Γh
f e φh dσh .
In order to verify that the symmetric bilinear form ah is positive definite on VhI × VhI we note that ah (φh , φh ) = 0 implies that Γh ∩T
2 |∇φh | + φ2h dσh = 0
for all T ∈ ThI .
Since H n (T ∩ Γh ) > 0 for T ∈ ThI we infer that ∇φh = 0 and hence φh is constant on these elements. Using again that H n (T ∩ Γh ) > 0 we deduce that φh = 0 on each T ∈ ThI so that φh ≡ 0 in VhI . Hence (3.5) has a unique solution uh ∈ VhI and (3.6)
1
uh h = ∇uh 2L2 (Γh ) + uh 2L2 (Γh ) 2 ≤ c f e L2 (Γh ) ≤ c f L2 (Γ) .
3.3. Error analysis. The following error bounds hold. Theorem 3.2. Let u be the solution of (1.1) and uh the solution of the finite element scheme (3.5). Then (3.7)
ue − uh L2 (Γh ) + h ∇(ue − uh ) L2 (Γh ) ≤ ch2 f L2 (Γ) .
2147
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
Proof. In view of the definition of · h , (2.27), and Lemma 3.1 we have for eh := ue − uh 1 2 2
eh L2 (Γh ) + ∇eh L2 (Γh ) 2 1 2 2 ≤ 2 ue − Ih ue L2 (Γh ) + ∇(ue − Ih ue ) L2 (Γh ) 2
(3.8)
+ sup φh ∈VhI
|ah (ue , φh ) − lh (φh )|
φh h
≤ ch u H 2 (Γ) + sup
φh ∈VhI
|ah (ue , φh ) − lh (φh )| .
φh h
In order to estimate the second term we choose φh ∈ VhI , then for ϕh := φh ◦ Fh−1 , ah (ue , φh ) − lh (φh ) = ah (ue , φh ) − a(u, ϕh ) + l(ϕh ) − lh (φh ) ≡ I + II. Using the transformation rule and (2.11) we obtain ∇Γ u · ∇Γ ϕh + uϕh dσ = (∇Γ u) ◦ p · (∇Γ ϕh ) ◦ p + (u ◦ p) (ϕh ◦ p) μh dσh Γ Γ h (I + ΦA)−1 ∇ue · (I + ΦA)−1 ∇ϕeh + ue ϕeh μh dσh . (3.9) = Γh
Since ϕeh (x) = ϕh (p(x)) = φh (Fh−1 (p(x))) we derive ∇ϕeh (x) = [Dp(x)]T [DFh−1 (p(x))]T ∇φh (Fh−1 (p(x))) = [Dp(x)]T [DFh (Fh−1 (p(x)))]−T ∇φh (Fh−1 (p(x))). We infer from (2.22) that (3.10)
(DFh )−T = I −
1 ∇ηh ⊗ ν + Bh |∇Φ|
with |Bh | ≤ ch2 ,
where ηh = Ih Φ − Φ. It follows from (2.19) that Fh−1 (p(x)) = x, x ∈ Γh , which together with (2.9) implies
1 e ∇ϕh = (I − ν ⊗ ν) I − ∇ηh ⊗ ν ∇φh + qh on Γh , |qh | ≤ ch2 |∇φh |. |∇Φ| Taking into account that ∇ue · ν = 0 we therefore have ∇ue · ∇ϕeh = ∇ue · ∇φh −
1 (∇ue · ∇ηh )(∇φh · ν) + ∇ue · qh |∇Φ|
on Γh .
Inserting this relation into (3.9) and recalling the definition of ah we find that e ∇u · ∇φh − μh (I + ΦA)−T (I + ΦA)−1 ∇ue · ∇ϕeh + |(μh − 1)ue ϕeh | dσh |I| ≤ Γh 2 e e ≤ ch u h φh h + ch u h |(∇φh · ν)| dσh , Γh
2148
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
where we used (2.15), (2.23), and the fact that ϕeh = φh on Γh . Similarly, |II| = f ϕh dσ − f e φh dσh ≤ |1 − μh | |f e | |φh | dσh Γ
Γh
Γh
≤ ch2 f e L2 (Γh ) φh h ≤ ch2 f L2 (Γ) φh h . Combining these estimates with (2.24) we have (3.11)
2
|ah (u , φh ) − lh (φh )| ≤ ch f L2 (Γ) φh h + ch f L2(Γ) e
Γh
|∇φh · ν| dσh
for all φh ∈ VhI , which inserted into (3.8) yields
eh L2 (Γh ) + ∇eh L2 (Γh ) ≤ ch f L2 (Γ) .
(3.12)
In order to improve the L2 -error bound we employ the usual Aubin–Nitsche argument. Denote by w ∈ H 2 (Γ) the solution of the dual problem a(ϕ, w) = eh ϕ dσ for all ϕ ∈ H 1 (Γ) with eh = eh ◦ Fh−1 , Γ
which satisfies
w H 2 (Γ) ≤ c eh L2 (Γ) .
(3.13) We have in view of (1.3)
(3.14)
2 eh , w) = a( eh , w) − ah (eh , we )
eh L2 (Γ) = a( + ah (eh , we − Ih we ) + ah (ue , Ih we ) − lh (Ih we ) ≡ I + II + III.
Similarly as above we deduce with the help of (3.12) and Lemma 2.2 |I| ≤ ch eh h we h ≤ ch2 f L2 (Γ) w H 1 (Γ) . Next, Lemma 3.1 and (3.12) imply |II| ≤ eh h we − Ih we h ≤ ch2 f L2 (Γ) w H 2 (Γ) . Finally, (3.11), the fact that ∇we · ν = 0, and Lemma 3.1 yield |∇(Ih we − we ) · ν| dσh |III| ≤ ch2 f L2 (Γ) Ih we h + ch f L2 (Γ) Γh
2
≤ ch f L2 (Γ) w H 2 (Γ) . Inserting the above estimates into (3.14) and recalling (3.13) we obtain
eh L2 (Γ) ≤ ch2 f L2 (Γ) , which together with Lemma 2.2 completes the proof since eeh = eh on Γh .
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
2149
4. Narrow band method (NBM). 4.1. Setting up the method. Recalling the definition of Dh (2.17), we consider ThB = {T ∈ Th | H n+1 (T ∩ Dh ) > 0} and UhB = T. T ∈ThB
We define the finite element space VhB on the triangulation ThB by VhB = {φh ∈ C 0 (UhB ) | φh|T ∈ P1 (T ) for each T ∈ ThB }. Let us first examine the approximation error for the space VhB . Lemma 4.1. We have for each function z ∈ H 2 (Γ) that Ih z e ∈ VhB satisfies (4.1)
√ 1 √ z e − Ih z e L2 (Dh ) + h ∇(z e − Ih z e ) L2 (Dh ) ≤ ch2 z H 2 (Γ) . h
Proof. We infer from (2.14) and Lemma 2.2 that 1 e
z − Ih z e 2L2 (Dh ) + h ∇(z e − Ih z e ) 2L2 (Dh ) h
1 ≤
z e − Ih z e 2L2 (T ) + h ∇(z e − Ih z e ) 2L2 (T ) h T ∩Dh =∅ ≤ ch3
z e 2H 2 (T ) ≤ ch3 z e 2H 2 (D(1+c )h ) ≤ ch4 z 2H 2 (Γ) , 1
T ∩Dh =∅
since T ⊂ D(1+c1 )h for all T ∩ Dh = ∅ in view of (2.16). 4.2. The method. Let us write (1.5) in the form find uh ∈ VhB such that ah (uh , φh ) = lh (φh )
(4.2) where
for all φh ∈ VhB ,
1 ∇wh · ∇φh + wh φh |∇Ih Φ| dx, ah (wh , φh ) = 2h Dh 1 lh (φh ) = f e φh |∇Ih Φ| dx. 2h Dh
Note that the factors h1 in each of the above terms are there to aid the notation for the error analysis. In a similar way as for SIF one can verify that ah is positive definite on VhB × VhB . Hence, the finite element scheme (4.2) has a unique solution uh ∈ VhB which satisfies (4.3)
uh h =
1 2h
Dh
12 2 |∇uh | + u2h |∇Ih Φ| dx ≤ c f L2 (Γ) .
4.3. Error analysis. Before we prove our main error bound we formulate a technical lemma which will be helpful in the error analysis. Lemma 4.2. Suppose that u ∈ H 2 (Γ) is a solution of (1.1). Then, 1 1 |∇Ih Φ| −1 e e ah (u , φ) = dx + S, φ f φ ◦ Fh |∇Φ| dx + (∇ue · ∇ηh )(∇φ · ν) 2h Dh 2h Dh |∇Φ|
2150
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
for all φ ∈ H 1 (Dh ), where ηh = Ih Φ − Φ and |S, φ| ≤ Ch2 u H 2 (Γ) φ h . Proof. To begin, we derive from (2.12) and (1.1) that −
(4.4) where
1 ∇ · |∇Φ| ∇ue + ue = f e + R |∇Φ| ⎛
R(x) = −Φ(x) ⎝
(4.5)
n+1
blk (x)D l Dk u(p(x)) +
k,l=1
n+1
in Uδ , ⎞ ck (x)D k u(p(x))⎠ .
k=1
We multiply (4.4) by φ ◦ Fh−1 |∇Φ|, φ ∈ H 1 (Dh ), and integrate over Dh . Since ∂u h ∂ν = 0 on ∂D we obtain after integration by parts ∇ue · ∇(φ ◦ Fh−1 )|∇Φ| dx + ue φ ◦ Fh−1 |∇Φ| dx h h D D (4.6) −1 e = f φ ◦ Fh |∇Φ| dx + R φ ◦ Fh−1 |∇Φ| dx. e
Dh
Dh
Observing that ∇(φ ◦ Fh−1 ) = [(DFh )−T ◦ Fh−1 ]∇φ ◦ Fh−1 , the transformation rule and Lemma 2.1 imply that −1 e ∇u ·∇(φ◦Fh )|∇Φ| dx = ∇ue ◦Fh ·(DFh )−T ∇φ |∇Φ◦Fh | |det DFh | dx. I := Dh
Dh
Recalling (2.7) and (2.18) we have z e (x) = z(p(x)) = z(p(Fh (x)) = z e (Fh (x)),
(4.7)
from which we deduce by differentiation ∇z e ◦ Fh = (DFh )−T ∇z e ,
(4.8) so that
I=
(DFh )−T ∇ue · (DFh )−T ∇φ |∇Φ ◦ Fh | |det DFh | dx.
Dh
Recalling (3.10), we find with the help of ∇ue · ν = 0 that (DFh )−T ∇ue = ∇ue + Bh ∇ue ,
(DFh )−T ∇φ = ∇φ −
1 (∇φ · ν)∇ηh + Bh ∇φ, |∇Φ|
where |Bh | ≤ ch2 . Furthermore, Lemma 2.1 implies that |∇Φ ◦ Fh | |det DFh | = |∇Ih Φ| + γh , so that in conclusion e ∇u · ∇φ |∇Ih Φ| dx − I= Dh
Dh
where |γh | ≤ ch2 ,
(∇ue · ∇ηh )(∇φ · ν)
|∇Ih Φ| dx + Rh1 , φ, |∇Φ|
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
where
1 Rh , φ ≤ ch2 ∇ue
L2 (Dh )
2151
∇φ L2 (Dh ) ≤ ch3 u H 2 (Γ) φ h ,
in view of Lemma 2.2 and the definition of · h . Similarly, (4.7) and (2.21) yield −1 e u φ ◦ Fh |∇Φ| dx = ue φ |∇Ih Φ| dx + Rh2 , φ Dh
Dh
with Rh2 , φ ≤ ch3 u H 2 (Γ) φ h . Inserting the above identities into (4.6) and dividing by 2h we derive (4.9) 1 1 ah (ue , φ) = f e φ ◦ Fh−1 |∇Φ| dx + R φ ◦ Fh−1 |∇Φ| dx 2h Dh 2h Dh 1 1 |∇Ih Φ| 1 dx − Rh1 , φ − R2 , φ. + (∇ue · ∇ηh )(∇φ · ν) 2h Dh |∇Φ| 2h 2h h In order to rewrite the integral over Dh we note that F (·, s) maps Γ onto Γs = {Φ = s} and that dσs = (1 + O(s)) dσp , where dσs , dσp are the surface elements of Γs , Γ respectively. The coarea formula then yields for integrable g : Dh → R h h 1 dσs ds = g(x) dx = g(x) g(F (p, s))μ(p, s) dσp ds |∇Φ(x)| −h Γ Dh −h Γs 1 ≤ C |s| , |s| < h, p ∈ Γ. μ(p, s) − (4.10) where |∇Φ(F (p, s))| Hence,
Dh
R φ ◦ Fh−1 |∇Φ| dx
h
= −h Γ h
= −h
Γ
R ◦ F φ ◦ Fh−1 ◦ F |∇Φ ◦ F | μ dσp ds R ◦ F φ ◦ Fh−1 ◦ F dσp ds +
h
−h
Γ
r φ ◦ Fh−1 ◦ F dσp ds
≡ T1 + T2 ,
where r(p, s) = R(F (p, s)) μ(p, s) |∇Φ(F (p, s))| − 1 . In order to treat T1 we deduce from (4.5) and the fact that Φ(F (p, s)) = s that ⎛ ⎞ n+1 n+1 R(F (p, s)) = −s ⎝ blk (F (p, s))D l Dk u(p) + ck (F (p, s))D k u(p)⎠ . k,l=1
Since −
h −h h
−h
k=1
s ds = 0, the first term in T1 can be written as
n+1 Γ k,l=1
sDl Dk u(p) blk (F (p, s)) φ ◦ Fh−1 (F (p, s)) − blk (p) φ ◦ Fh−1 (p) dσp ds.
Treating the second term in T1 in the same way and observing that p = F (p, 0) we deduce with the help of the fundamental theorem of calculus that
12 2 2 5 −1 −1 ∇φ ◦ Fh dx + φ ◦ Fh ≤ ch3 u H 2 (Γ) φ h . |T1 | ≤ ch 2 u H 2 (Γ) Dh
2152
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
Next, we infer from (4.5) and (4.10) that | r (p, s)| ≤ cs2 |∇Γ u(p)| + DΓ2 u(p) , so that
5
|T2 | ≤ Ch 2 u H 2 (Γ)
Dh
12 φ ◦ F −1 2 dx ≤ Ch3 u H 2 (Γ) φ h . h
The result now follows from (4.9) together with the bounds on Rh1 and Rh2 . Theorem 4.3. Let u be the solution of (1.1) and uh the solution of the finite element scheme (4.2). Then
u − uh L2 (Γh ) + h e
(4.11)
1 2h
12 |∇(u − uh )| |∇Ih Φ| dx ≤ ch2 f L2 (Γ) . 2
e
Dh
Proof. Let us write eh := ue − uh . We infer from (2.27) and Lemma 4.1 that
eh h ≤ ch u H 2 (Γ) + sup
(4.12)
φh ∈VhB
|ah (ue , φh ) − lh (φh )| .
φh h
The second term on the right-hand side can be estimated with the help of Lemma 4.2. The transformation rule together with (4.7) yields 1 1 f e φh ◦ Fh−1 |∇Φ| dx = f e ◦ Fh φh |∇Φ ◦ Fh | |det DFh | dx, 2h Dh 2h Dh so that we deduce from Lemma 4.2 1 |∇Ih Φ| ah (ue , φh ) − lh (φh ) = (∇ue · ∇ηh ) (∇φh · ν) 2h Dh |∇Φ| 1 + f e φh |∇Φ ◦ Fh | |detDFh | − |∇Ih Φ| dx + S, φh . 2h Dh Using (2.21), (2.15), (2.24), and Lemmas 2.1 and 2.2 we infer that for φh ∈ VhB |ah (u , φh ) − lh (φh )| ≤ c ∇u L2 (Dh ) e
e
12 |∇φh · ν| dx 2
Dh
+ ch f e L2 (Dh ) φh L2 (Dh ) + ch2 u H 2 (Γ) φh h
12 1 2 ≤ ch f L2 (Γ) |∇φh · ν| dx + ch2 f L2 (Γ) φh h , 2h Dh
(4.13)
so that (4.12) implies the following intermediate result: (4.14)
eh h =
1 2h
Dh
12 2 |∇eh | + e2h |∇Ih Φ| dx ≤ ch f L2 (Γ) .
In order to improve the L2 -error bound we define eh := eh ◦ Fh−1 as well as h (p) := 1 E 2h
h
−h
eh (F (p, s)) ds,
p ∈ Γ,
2153
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
with F as above. We denote by w ∈ H 2 (Γ) the unique solution of h −ΔΓ w + w = E which satisfies
w H 2 (Γ) ≤ c E h
(4.15)
on Γ,
L2 (Γ)
.
Similarly to (4.4) the extension we solves −
1 e + R ∇ · |∇Φ| ∇we + we = E h |∇Φ|
in Uδ ,
is obtained from (4.5) by replacing u by w. Using the transformation rule where R together with (4.10) we obtain h h 2 1 1 e ◦ F eh ◦ F |∇Φ ◦ F | μ dσp ds Eh eh ◦ F dσp ds = E Eh 2 = 2h −h Γ 2h −h Γ h L (Γ) h 1 h E + eh ◦ F (1 − |∇Φ ◦ F | μ) dσp ds 2h −h Γ 1 e eh ◦ F −1 |∇Φ| dx E = h 2h Dh h h 1 h + E eh ◦ F (1 − |∇Φ ◦ F | μ) dσp ds. 2h −h Γ The first term can be rewritten with the help of Lemma 4.2 (applied to w instead of u) to give 2 |∇Ih Φ| eh − 1 dx (∇we · ∇ηh )(∇eh · ν) Eh 2 = ah (we , eh ) − S, 2h Dh |∇Φ| L (Γ) h 4 1 h eh ◦ F (1 − |∇Φ ◦ F | μ) dσp ds ≡ + E Ik . 2h −h Γ k=1
In view of Lemma 4.1, (4.13), the fact that ∇we · ν = 0, and (4.14) we have eh | |I1 | + |I2 | ≤ |ah (we − Ih we , eh )| + |ah (ue , Ih we ) − lh (Ih we )| + |S, ≤ ch w H 2 (Γ) eh h + ch2 f L2 (Γ) Ih we h
12 1 e e 2 + ch f L2 (Γ) |∇(Ih w − w ) · ν| dx + ch2 w H 2 (Γ) eh h 2h Dh ≤ ch2 f L2 (Γ) w H 2 (Γ) . Furthermore, (2.15), (4.10), and (4.14) imply
e 2 |I3 |+|I4 | ≤ ch w h + Eh
eh h ≤ ch f L2 (Γ) w H 2 (Γ) + E h L2 (Γ)
so that we obtain together with (4.15) Eh 2 ≤ ch2 f L2 (Γ) . L (Γ)
L2 (Γ)
,
2154
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
Next, since F (p, 0) = p we may write for p ∈ Γ h (p) − eh (p) = 1 E 2h
h
−h
s
∇ eh (F (p, τ )) ·
0
∂F (p, τ ) dτ ds, ∂s
and hence we obtain with the help of (4.14) Eh − eh
L2 (Γ)
√ ≤c h
12 |∇ eh | dx ≤ ch ∇eh h ≤ ch2 f L2 (Γ) . 2
Dh
In conclusion we deduce that
eh L2 (Γh ) ≤ c eh L2 (Γ) ≤ E h h−e
L2 (Γ)
+ E h
L2 (Γ)
≤ ch2 f L2 (Γ)
and the theorem is proved. 5. A hybrid method for equations on evolving surfaces. 5.1. The setting. The aim of this section is to combine ideas employed in sections 3 and 4 for the stationary problem in order to develop a finite element method for an advection–diffusion equation on an evolving hypersurface in which there is an underlying conservation law with a diffusive flux, [18]. More precisely, let (Γ(t))t∈[0,T ] be a family of compact, connected smooth hypersurfaces embedded in Rn+1 for n = 1, 2. We suppose that Γ(t) = {x ∈ N (t) | Φ(x, t) = 0},
where ∇Φ(x, t) = 0, x ∈ N (t),
and N (t) is an open neighborhood of Γ(t). We assume that N (t) is chosen so small that we can construct the function p(·, t) as in section 2.1. Given a velocity field v(·, t) : Γ(t) → Rn+1 , not necessarily in the normal direction, we then consider the following initial value problem (5.1a) Γ(t) × {t}, ∂ • u + u∇Γ · v − ΔΓ u = f on t∈(0,T )
(5.1b)
u(·, 0) = u0
on Γ(0).
Here, ∂ • η denotes the material derivative of a function η : which is given by
t∈(0,T )
Γ(t) × {t} → R
∂ • η = ∂t η + v · ∇η if η is extended into a neighborhood of t∈(0,T ) Γ(t) × {t}. 5.2. The method. In order to discretize the above problem we choose a partition 0 = t0 < t1 < · · · < tN = T of [0, T ] with τm := tm+1 − tm , m = 0, . . . , N − 1, and τ := maxm=0,...,N −1 τm . Also, let Th be an unfitted regular triangulation with mesh size h of a region containing N (t), t ∈ [0, T ]. For m = 0, 1, . . . , N we set Γm h = {x ∈ N (tm ) | Ih Φ(x, tm ) = 0}
and
Dhm = {x ∈ N (tm ) | |Ih Φ(x, tm )| < h},
as well as Thm := {T ∈ Th | H n+1 (T ∩ Dhm ) > 0}
and
Uhm :=
T ∈Thm
T.
UNFITTED FINITE ELEMENT METHODS FOR SURFACE PDES
2155
Here we assume that 0 < h ≤ h0 , where h0 is chosen so small that there exists c0 , c1 > 0 such that N (t) × {t}. c0 ≤ |∇Ih Φ(x, t)| ≤ c1 , (x, t) ∈ t∈(0,T )
Finally, we introduce Vhm = {φh ∈ C 0 (Uhm ) | φh |T ∈ P1 (T ) for each T ∈ Thm }. In what follows we shall frequently use the abbreviation z m (x) := z(x, tm ). In order to motivate our method we fix m ∈ {0, 1, . . . , N − 1} and let Ψ be the solution of Ψt (x, t) + DΨ(x, t)v e (x, t) = 0,
Ψ(x, tm+1 ) = x,
where v e (x, t) := v(p(x, t), t). For a sufficiently smooth function ϕ : N (tm+1 ) → R we define η(x, t) := ϕ(Ψ(x, t)). Clearly, η(·, tm+1 ) = ϕ and a short calculation shows that ∂ • η = 0. Assuming that u is a solution of (5.1a) we obtain with the help of the Leibniz formula and integration by parts • d ∂ (uη) + uη∇Γ · v dσ uη dσ|t=tm+1 = dt Γ(t) Γ(t ) m+1 ϕ ∂ • u + u∇Γ · v dσ = Γ(t ) m+1 ϕΔΓ u + ϕf dσ = Γ(tm+1 )
=−
Γ(tm+1 )
∇Γ u · ∇Γ ϕ dσ +
f ϕ dσ. Γ(tm+1 )
Since Ψ(·, tm+1 ) ≡ id, a Taylor expansion shows that Ψ(x, tm ) = Ψ(x, tm+1 − τm ) ≈ x − τm Ψt (x, tm+1 ) = x + τm DΨ(x, tm+1 )v e,m+1 (x) = x + τm v e,m+1 (x). Thus we may approximate the left-hand side of the above relation by d 1 uη dσ|t=tm+1 ≈ um+1 ϕ dσ − um ϕ(· + τm v e,m+1 ) dσ . dt Γ(t) τm Γ(tm+1 ) Γ(tm ) The above calculations motivate the following scheme, in which we use the narrow m band approach in order to discretize the elliptic part. Given um h ∈ Vh , m = 1, . . . , N − m+1 m+1 1, find uh ∈ Vh such that e,m+1 um+1 φ dσ − um ) dσh h h h φh (· + τm v h (5.2)
Γm+1 h
τm + 2h
Γm h
m+1 Dh
∇um+1 h
· ∇φh |∇Ih Φ
m+1
| dx = τm
Γm+1 h
f e,m+1 φh dσh
for all φh ∈ Vhm+1 . Here, u0h = Ih u0 . Existence and uniqueness of um+1 follows in a h similar way as for NBM in the elliptic case.
2156
K. DECKELNICK, C. M. ELLIOTT, AND T. RANNER
5.3. Mass conservation. An important property of solutions of (5.1a) is conservation of mass in the case that Γ(t) f (·, t) dσ = 0. The following lemma shows that our numerical scheme preserves this property under some mild constraints on the discretization parameters. In fact this discrete conservation law is a remarkable property of the scheme and relies on the use of both a sharp interface and a narrow band approach. m Lemma 5.1. Let um h ∈ Vh , m = 1, . . . , N , be the solutions of (5.2) in the case √ f e,m+1 dσh = 0, m = 0, . . . , N −1. Then provided that 0 < h ≤ h1 and τ ≤ γ h Γm+1 h
(5.3)
Γm h
um h dσh =
Γ0h
u0h dσh .
Proof. Let us first observe that m+1 {x + τm v e,m+1 (x) | x ∈ Γm , h } ⊂ Uh
(5.4)
m = 0, . . . , N − 1,
provided that h, τ are sufficiently small. To see this, let x ∈ Γm h and choose an element T ∈ Th such that x ∈ T . Then, Φm+1 (x + τm v e,m+1 (x)) = Φm (x) + τm ∇Φm (x) · v e,m+1 (x) + τm Φt (x, tm ) + Rm (x), 2 where |Rm (x)| ≤ cτm . Observing that Φt + ∇Φ · v = 0 on
t∈(0,T )
Γ(t) × {t} we write
∇Φm (x) · v e,m+1 (x) + Φt (x, tm ) = ∇Φm (x) · v m+1 (pm+1 (x)) + Φt (x, tm ) = ∇Φm (x) − ∇Φm+1 (pm+1 (x)) · v m+1 (pm+1 (x)) + Φt (x, tm ) − Φt (pm+1 (x), tm+1 ), so that m+1 2 Φ (x + τm v e,m+1 (x)) ≤ |Φm (x)| + cτm x − pm+1 (x) + cτm 2 2 ≤ |Φm (x)| + cτm Φm+1 (x) + cτm ≤ c(h2 + τm ), in view of (2.5) and since |Φm (x)| = |Φm (x) − Ih Φm (x)| ≤ ch2 . As a result, (Ih Φm+1 )(x + τm v e,m+1 (x)) ≤ Ih Φm+1 − Φm+1
L∞
2 2 + c(h2 + τm ) ≤ c(h2 + τm )
