# A posteriori error control for stationary coupled bulk-surface equations

A posteriori error control for stationary coupled bulk-surface equations Abstract We consider a system of two coupled elliptic equations, one defined on a bulk domain and the other one on the boundary surface. Problems of this kind are relevant for applications in engineering, chemistry and in biology, e.g., biological signal transduction. For the a posteriori error control of the coupled system, a residual error estimator is derived which takes into account the approximation errors due to the finite element discretization in space as well as the polyhedral approximation of the surface. An adaptive refinement algorithm controls the overall error. Numerical experiments illustrate the performance of the a posteriori error estimator and the proposed adaptive algorithm with several benchmark examples. 1. Introduction Coupled reaction diffusion processes in the bulk and on the surface of some domain $${\it{\Omega}}\subset\mathbb{R}^d,\,d\in\{2,3\}$$, have recently attracted interest from an analytical point of view (see, e.g., Elliott & Ranner, 2013; Burman et al., 2015b; Egger et al., 2015) in different application areas such as biology (see, e.g., Novak et al., 2007; Rätz & Röger, 2012; Giese et al., 2015). In all these problems one has to simultaneously account for transport in normal as well as in tangential direction to the boundary surface, whereas most often either only normal or only tangential phenomena are considered. We study a stationary prototype problem in domain $${\it{\Omega}}$$ with piecewise smooth boundary $${\it{\Gamma}}:=\partial{\it{\Omega}}$$, which can be decomposed into a finite set of $$N_{\it{\Gamma}}$$ smooth surface patches $$\{{\it{\Gamma}}^i\}_{i=1}^{N_{\it{\Gamma}}}$$. We seek the solution $$u:{\it{\Omega}}\to\mathbb{R}$$ and $$v:{\it{\Gamma}}\to\mathbb{R}$$ of the stationary coupled diffusion-reaction problem   −Δu+u =f in Ω, (1.1a)  (αu−βv)+∂nu =0 on ∪i=1NΓΓi, (1.1b)  −ΔΓv+v+∂nu =g on ∪i=1NΓΓi, (1.1c)  ∇Γv∣Γi⋅nΓi+∇Γv∣Γj⋅nΓj =0 on ∂Γi∩∂Γj. (1.1d) For each point on a surface patch $${\it{\Gamma}}^i$$, $$n$$ is the outer normal of $${\it{\Omega}}$$, and $$\nabla_{{\it{\Gamma}}^i}$$ and $${\it{\Delta}}_{\it{\Gamma}}$$ denote the tangential gradient and the Laplace–Beltrami operator, respectively. In the compatibility condition (1.1d), $$n_{\it{\Gamma}}^i$$ denotes the outer normal of the patch $${\it{\Gamma}}^i$$, which is often called the co-normal to $${\it{\Gamma}}^i$$. For each point in $$\partial{\it{\Gamma}}^i$$, it lies in the tangent plane to $${\it{\Gamma}}^i$$. The equation system (1.1) describes diffusive transport and reaction of a bulk species $$u$$ and a surface species $$v$$, which are coupled by some binding–unbinding process that transforms $$u$$ to $$v$$ by binding it to $${\it{\Gamma}}$$ and vice versa. For the sake of a convenient presentation, the diffusion coefficients and reaction rates are assumed to be $$1$$, while the binding and unbinding rates are given by the positive constants $$\alpha$$ and $$\beta$$. The existence and uniqueness of the solution to (1.1) have been proved by Elliott & Ranner (2013) for globally $$C^2$$ boundaries. Moreover, optimal a priori error estimates for finite element approximations of arbitrary order were derived for sufficiently smooth boundaries and data. For the time-dependent version of (1.1), convergence to equilibrium and a priori error estimates were shown by Egger et al. (2015). In Burman et al. (2015b), a CutFEM method that uses the same volume grid to discretize both the volume and the surface equation is analysed. We complement the analysis of Elliott & Ranner (2013) by the derivation of an a posteriori error estimator for the error $$e:=(u-{u_h},v-{v_h})$$ of a lowest order continuous FEM approximation $$({u_h},{v_h})$$ of the solution $$(u,v)$$ of (1.1) measured in the energy norm. For this, we split the overall residual into an equilibration, a consistency and a data approximation residual. We transfer the classical residual error estimators for volume domains and for surfaces to the setting of the coupled system. While the a posteriori analysis for problems on polyhedral domains is rather mature (see, e.g., Babuška & Rheinboldt, 1978; Verfürth, 1996; Ainsworth & Oden, 2000), and a priori analysis for finite elements on curved domains is well established (see, e.g., Scott, 1973; Zláamal, 1973; Bernardi, 1989), much less can be found for the a posteriori analysis on curved domains. Most notably, Dörfler & Rumpf (1998) put a focus on a very coarse domain approximation, where large parts of the domain may not be discretized at all and the refinement algorithm has to explore and detect the necessary parts of the domain which are important for the discretization. This is somehow opposite to the setting we are concerned with since, due to the coupling to the surface equation, we expect the complete domain to be required for an accurate solution. For the solution of PDE problems on surfaces, various numerical approaches can be used, cf. the recent review of Dziuk & Elliott (2013). As a first approach, finite element methods on (surface-)meshes that approximate the surface have been devised. For smooth surfaces, they are analysed by Dziuk (1988), and Demlow & Dziuk (2007) derived an a posteriori error estimate which splits the error into contributions from the discretization, geometric errors and data approximations. For an IP-dG scheme, an error estimator of similar structure for a reaction diffusion equation on a smooth surface was analysed by Dedner & Madhavan (2016), in which case the geometric part was shown to be of a higher order. Dassi et al. (2015) transferred ideas of Demlow & Dziuk (2007) into the context of anisotropic meshes and an $$L^1$$ interpolation error estimate is derived together with heuristic control of the geometric error. A posteriori estimates with ZZ-type reconstruction of the solution are derived by Dassi (2014), where again the extension to anisotropic meshes is considered. Higher order methods for pure surface problems were analysed by Demlow (2009). $$L^2$$ and pointwise estimates for $${\it{\Gamma}}\in C^3$$ are given by Camacho & Demlow (2015). For less smooth surfaces, Bonito et al. (2013) derived a convergent refinement algorithm for parametric finite elements, and quasi-optimality of the resulting discretizations is shown under certain assumptions. (Bonito et al., 2016) contains the extension to higher order schemes. Our error estimator extends some of the results by Bonito et al. (2013) for the coupled bulk-surface equations. One main difficulty and contribution of this paper is the analysis of the geometry error of nonconformity caused by the polyhedral approximation of the surface manifold. As an alternative approach, unfitted volume meshes can be used to solve PDE problems on lower-dimensional surfaces. On the one hand, one can extend the equations from the surface to the volume (cf. Bertalmío et al., 2001; Burger, 2009). This approach is in particular well suited for surfaces that are implicitly defined as the level set of a smooth function. Additional computational costs for solving in higher space dimension can be reduced by narrow band techniques as shown by Deckelnick et al. (2010). On the other hand, in TraceFEM or CutFEM the trial functions are defined as traces of functions from spaces related to the volume mesh (see, e.g., Olshanskii et al., 2009; Chernyshenko & Olshanskii, 2015; Burman et al., 2015a). For coupled volume-surface problems with smooth surfaces, TraceFEM and CutFEM are applied and analysed by Gross et al. (2015) and Burman et al. (2015b), respectively. The outline of the paper is as follows: In Section 2, we introduce the description of curved elements and recall some basic results from differential geometry which include integral transformations used later on. In particular, nonconformity estimates are derived which provide bounds for the error caused by the domain approximations. Then, the weak formulation and FEM discretization of the problem (1.1) is presented in Section 3. Section 4 is dedicated to the derivation of the residual a posteriori error estimator for the overall error which is the main result of this paper. It consists of a discretization part and a geometric nonconformity part, both of which can be controlled a posteriori by computable indicators. Numerical experiments in Section 5 illustrate the performance of the derived a posteriori error estimator with several benchmark problems. In order to simplify the notation, we write $$a\lesssim b$$ if $$a\leq C\,b$$ with some mesh-independent positive constant $$C$$. Moreover, the common notation for Lebesgue and Sobolev spaces is employed. 2. Domain approximation A finite element method for (1.1) involves an approximation of the curved domain $${\it{\Omega}}$$ by some interpolating polyhedral domain $${\it{\Omega}}_h$$. If $${\it{\Omega}}$$ is a fixed piecewise smoothly bounded domain, all nonsmooth edges or points on the boundary are known a priori, and no degeneration of the smooth patches $${\it{\Gamma}}^i \subseteq \partial{\it{\Omega}}$$ can happen, in contrast to problems on evolving domains. It is thus reasonable to define $${\it{\Omega}}$$ as the deformation of some polyhedral domain $${\it{\Omega}}_h^0$$ with an associated triangulation $$\mathcal{T}_h^0$$ into simplices. We refer to $$\mathcal{T}_h^0$$ as the macro triangulation and denote the deformation by $${\it{\Phi}}^0$$. We assume that $${\it{\Phi}}^0$$ is globally continuous, and the restriction to each boundary face or boundary edge is a $$C^1$$ function. Moreover, the generality of $${\it{\Phi}}^0$$ is limited by the requirement that there are no self intersections of $${\it{\Gamma}}$$ and that $${\it{\Gamma}}$$ does not intersect the inner facets of $$\mathcal{T}_h^0$$, as detailed below. In the following, we consider different triangulations $$\mathcal{T}_h$$ with their respective polyhedral domain $${\it{\Omega}}_h$$ and functions $${\it{\Phi}}$$ which describe the deformation of $$\partial {\it{\Omega}}_h$$ to match $$\partial {\it{\Omega}}$$. We always assume that $$\mathcal{T}_h$$ has been obtained from $$\mathcal{T}_h^0$$ by several steps of refinement by bisection without indicating the refinement level. Each refinement step is followed by a deformation such that the respective polygonal domain $$\partial {\it{\Omega}}_h$$ interpolates $$\partial {\it{\Omega}}$$. 2.1 Simplicial triangulations Let $${\it{\Omega}}_h\subset\mathbb{R}^{d}$$ be a polyhedral domain with the boundary $${{\it{\Gamma}}_{\!h}} = \partial{\it{\Omega}}_h$$, which interpolates $${\it{\Gamma}}$$. Suppose $$\mathcal{T}_h$$ is a regular partition of $${\it{\Omega}}_h$$ into simplices. For $$d=3$$, we denote the sets of nodes, edges and faces of the triangulation $$\mathcal{T}_h$$ by $$\mathcal{N}_h$$, $$\mathcal{E}_h$$ and $${\mathcal{F}}_h$$, respectively. If $$d=2$$ we let $${\mathcal{F}}_h$$ denote the set of edges and let $$\mathcal{E}_h$$ coincide with $$\mathcal{N}_h$$. For each $$T_h\in\mathcal{T}_h$$, $$F_h\in{\mathcal{F}}_h$$, $$E_h\in\mathcal{E}_h$$, we set $$h_{T} = \operatorname{diam}(T_h)$$, $$h_{F} = \operatorname{diam}(F_h)$$, $$h_{E} = \operatorname{diam}(E_h)$$ and define by this the (global) mesh-size functions $$h_{\mathcal{T}}$$, $$h_{{\mathcal{F}}}$$, $$h_{\mathcal{E}}$$. The jump of some (piecewise) function $$u\in L^2({\it{\Omega}}_h;\mathbb{R}^d)$$ over a face $$F_h\in{\mathcal{F}}_h$$ is defined by $$[u]_{F_h}:= u\cdot n_+ + u\cdot n_-$$ where $$n_+, n_-$$ denote the outer normals on $$F_h=T_+\cap T_-$$ with respect to the elements $$T_+,T_-\in\mathcal{T}_h$$. The subset $${\mathcal{F}}_h^\partial := \{ F_h \in {\mathcal{F}}_h : \ F_h \subset {{\it{\Gamma}}_{\!h}} \}$$ defines a triangulation of $${{\it{\Gamma}}_{\!h}}$$ into $$(d-1)$$ dimensional simplices. Analogous to the volume, we introduce the sets of boundary edges and boundary nodes, denoted by $$\mathcal{E}_h^\partial$$ and $$\mathcal{N}_h^\partial$$. For $$E_h\in\mathcal{E}_h^\partial$$ and $$F_+,F_-\in{\mathcal{F}}_h^\partial$$ such that $$E_h = F_+\cap F_-$$, the jump of some (piecewise) function $$v\in L^2(\partial{\it{\Omega}}_h;\mathbb{R}^d)$$ over the boundary edge $$E_h$$ is defined by $$[[v]]_{E_h} := v\cdot n_{{\it{\Gamma}}_{\!h}}^+ + v\cdot n_{{\it{\Gamma}}_{\!h}}^-$$, where $$n_{{\it{\Gamma}}_{\!h}}^+, n_{{\it{\Gamma}}_{\!h}}^-$$ denote the outer co-normals on $$E_h$$ with respect to the boundary facets $$F_+,F_-\in{\mathcal{F}}_h^\partial$$. Note that $$n_{{\it{\Gamma}}_{\!h}}^+, n_{{\it{\Gamma}}_{\!h}}^-$$ are not necessarily aligned. For the analysis below, we need the following assumptions on the triangulation: (A1) $$\mathcal{T}_h$$ interpolates $${\it{\Gamma}}$$, i.e., the set of boundary nodes satisfies $$\mathcal{N}_h^\partial \subset {\it{\Gamma}}$$. (A2) Each element $$T_h \in \mathcal{T}_h$$ has at most one boundary face $$F_h \subseteq \partial{\it{\Omega}}_h \cap \partial T_h$$. If $$T_h$$ has no boundary face, then there is at most one boundary edge $$E_h \subseteq \partial{\it{\Omega}}_h \cap \partial T_h$$. (A3) All refinements of $$\mathcal{T}_h$$ by bisection and subsequent deformation of $$\partial{\it{\Omega}}_h$$ to interpolate $$\partial {\it{\Omega}}$$ lead to a shape regular family of triangulations. On fixed polyhedral volume domains $${\it{\Omega}}$$, it is known that appropriate bisection algorithms guarantee the shape regularity, and thus assumption (A3) is not necessary in this case. For the case of finite elements on curved surfaces, we refer to the discussion by Demlow & Dziuk (2007) and Bonito et al. (2013). Let $$\hat T$$ denote the $$d$$-dimensional unit simplex. Given an element $$T_h\in\mathcal{T}_h$$, there is an affine reference transformation which we refer to as $$J:\hat{T} \to T_h$$, without any further index related to $$T_h$$. If $$T_h$$ has one boundary face $$F_h\subset{\mathcal{F}}_h^\partial$$ then we assume that $$\hat{F} := \{ \hat x \in \hat T: \hat x_d = 0 \}$$ is mapped by $$J$$ to the boundary segment $$F_h\subset {{\it{\Gamma}}_{\!h}}$$ as depicted in Fig. 1. Analogously, if $$T_h\cap{\it{\Gamma}}_h$$ consists only of one boundary edge $$E_h$$, we assume that $$E_h = J(\hat{E})$$, where $$\hat{E} := \{ \hat x \in \hat T: \hat x_d =\hat x_{d-1} = 0 \}$$. Due to the shape regularity assumption (A3), there is a constant $$L_h>0$$ such that   hLh−1|z|≤‖∇^Jz‖L∞(T^)≤hLh|z| for all z∈Rd. (2.1) Fig. 1. View largeDownload slide Transformation $$\mathcal{X}$$ and $$J$$ from reference element $$\hat T$$ to a curved element $$T$$ and to its polyhedral approximation $$T_h$$, respectively. Fig. 1. View largeDownload slide Transformation $$\mathcal{X}$$ and $$J$$ from reference element $$\hat T$$ to a curved element $$T$$ and to its polyhedral approximation $$T_h$$, respectively. 2.2 Exact triangulation of curved domains Let $${\it{\Omega}}_h$$ and $$\mathcal{T}_h$$ be as in Section 2.1 satisfying (A1)-(A3). The assumed refinement process from $$\mathcal{T}_h^0$$ to $$\mathcal{T}_h$$ also yields an according boundary deformation $${\it{\Phi}}$$ from $${\it{\Phi}}^0$$, which moves $${\it{\Gamma}}_h$$ to $${\it{\Gamma}}$$. (A1) implies that $${\it{\Phi}} = 0$$ in all boundary nodes. We assume that the refinement process is such that (A4) $${\it{\Phi}}: {{\it{\Gamma}}_{\!h}} \to \mathbb{R}^d$$ is globally continuous and the restriction to each boundary face $$F_h \in{\mathcal{F}}_h^\partial$$ or boundary edge $$E_h\in\mathcal{E}_h^\partial$$ is a $$C^1$$ function. At first, we restrict our considerations to a single boundary element $$T_h \in \mathcal{T}_h$$ such that $$T_h\cap{{\it{\Gamma}}_{\!h}} = F_h \in{\mathcal{F}}_h^\partial$$ or $$T_h\cap{{\it{\Gamma}}_{\!h}} = E_h\in\mathcal{E}_h^\partial$$. We denote the restriction of the boundary deformation to $$F_h$$ or $$E_h$$ again by $${\it{\Phi}}$$. With the pull-back to the reference element $$\hat{T}$$, we can define an extension to $$\hat{T}$$ as   Φ^T(x^) ={Φ(J(x^1,…,x^d−1,0)ρ(0,x^d))if Th∩Γh=Fh,Φ(J(x^1,…,x^d−2,0,0)ρ(x^d−1,x^d))if Th∩Γh=Eh,  (2.2) with some monotonous $$C^1$$-function $$\rho$$ which interpolates between $$0$$ and $$1$$ such that $$\rho(0,0)=1$$. The simplest choice is $$\rho(s,t) = 1 -s -t$$. Lemma 2.1 Let $$T_h\in\mathcal{T}_h$$ be a boundary element and $${\it{\Phi}} \in C^1(T_h\cap{\it{\Gamma}}_h, \mathbb{R}^d)$$ be a boundary deformation that vanishes in all boundary nodes. Let $$\hat{{\it{\Phi}}}_T$$ be defined according to (2.2) and set $$1 \leq L_\rho := \| \hat{\nabla} \rho \|_{L^\infty}$$. Suppose there is $$L_T>0$$ such that (a)$$\| \nabla_{\!F_h} {\it{\Phi}} \|_{L^\infty(F_h)} \leq (1+L_\rho)^{-1} L_T$$, if $$T_h\cap{{\it{\Gamma}}_{\!h}} = F_h\in{\mathcal{F}}_h^\partial$$, or (b)$$\| \nabla_{\!E_h} {\it{\Phi}} \|_{L^\infty(E_h)} \leq (1+L_\rho)^{-1} L_T$$, if $$T_h\cap{{\it{\Gamma}}_{\!h}} = E_h\in\mathcal{E}_h^\partial$$. In either case it holds   ‖∇^Φ^T‖L∞(T^) ≤LT‖∇^J‖L∞(T^). (2.3) Proof. We assume $$T_h \cap {{\it{\Gamma}}_{\!h}} = F_h$$. The other case follows from analogous arguments. From the construction of $$\hat{{\it{\Phi}}}_T$$ in (2.2) and abbreviating $$\hat{x}_{\hat{F}} = ( \hat{x}_1, \dots, \hat{x}_{d-1} , 0 ) \in \hat{F}$$, we check   ∇^Φ^T(x^) =( ρ(0,x^d)∇FhΦ(J(x^F^))∇^J,ρ′(0,x^d)Φ(J(x^F^)) ). (2.4) By the triangle inequality and, since $$0\leq\rho\leq 1$$, we infer   ‖∇^Φ^T‖L∞(T^) ≤‖∇FhΦ‖L∞(Fh)‖∇^J‖L∞(T^)+maxx^∈T^|ρ′(0,x^d)|‖Φ‖L∞(Fh). (2.5) Because $${\it{\Phi}}(J(0)) = 0$$, and using that $${\it{\Phi}}(J(\hat{x}_{\hat F})) = \int_{0}^1 \frac{d}{d s} {\it{\Phi}}( J (s \cdot \hat{x}_{\hat{F}})) \,\mathrm{d}s$$, we conclude $$\| {\it{\Phi}} \|_{L^\infty(F_h)} \leq \| \nabla_{\!F_h}{\it{\Phi}} \|_{L^\infty(F_h)} \| \hat\nabla J \|_{L^\infty(\hat T)}$$ to confirm the assertion. □ Next, we define for all $$\hat{x} \in \hat{T}$$ the nonlinear transformation   X(x^) =J(x^)+Φ^T(J(x^)) (2.6) and call $$T := {\mathcal{X}}(\hat{T})$$ the curved element related to $$T_h$$. Note that by the above construction, $$T$$ can have more than one curved face if $$d=3$$. If $$L_T<1$$, we conclude with Lemma 2.1 that $${\mathcal{X}}$$ is a $$C^1$$-diffeomorphism and the Jacobi matrix $$\hat \nabla {\mathcal{X}}$$ is positive definite (cf. Ciarlet & Raviart, 1972). In particular, it holds   (1−LT)‖∇^J‖L∞(T^)≤‖∇^X‖L∞(T^)≤(1+LT)‖∇^J‖L∞(T^). (2.7) Upon setting $$L := L_h \, \max((1+L_T) , (1-L_T)^{-1})$$, we get from (2.1) and (2.7) the bi-Lipschitz property of the nonlinear transformation   hL−1|z|≤‖∇^Xz‖L∞(T^)≤hL|z| for all z∈Rd. (2.8) Finally, we set $$\hat{{\it{\Phi}}}_T=0$$ on all inner elements and assume that (A5) there is a global constant $$L_T<1$$ such that for all boundary elements $$T_h\in\mathcal{T}_h$$ the extensions $$\hat{\it{\Phi}}_T$$ of $${\it{\Phi}}$$ according to (2.2) satisfies $$\| \hat{\nabla} \hat{{\it{\Phi}}}_T \|_{L^\infty(\hat{T})} \leq L_T \, \| \hat{\nabla} J \|_{L^\infty(\hat{T})}$$. Then, for a given triangulation $$\mathcal{T}_h$$ of $${\it{\Omega}}_h$$, we introduce the exact triangulation$$\mathcal{T}$$ of $${\it{\Omega}}$$ consisting of curved elements such that $${\it{\Gamma}}=\partial{\it{\Omega}}$$ is matched exactly. (A4) provides that $$\mathcal{T}$$ is a covering of $${\it{\Omega}}$$ and (A5) is related to the quasi-monotonicity of the geometric indicator of Bonito et al. (2013, 2016) and also implies that the approximation of $${\it{\Omega}}$$ by $${\it{\Omega}}_h$$ has reached the saturated state according to Dörfler & Rumpf (1998). Example 2.2 (a) Let $${\it{\Omega}} = B_1(0) \subset \mathbb{R}^2$$ the unit circle. To approximate $${\it{\Omega}}$$, we choose $${\it{\Omega}}_h^0 = \{(x,y)\in\mathbb{R}^2 : |x|+|y| \leq 1\}$$. A macro triangulation $$\mathcal{T}_h^0$$ of $${\it{\Omega}}_h^0$$ can be constructed from four triangles that all share the midpoint as a common node; see Fig. 2(left). One possible choice for a deformation of $${\it{\Omega}}_h^0$$ to $${\it{\Omega}}$$ is   Φ0(x,y)=(|x|+|y|R(x,y)−1)(xy), (2.9) where $$R(x,y) = \sqrt{x^2+y^2}$$. On the element $$T_h \in \mathcal{T}_h^0$$ indicated in Fig. 2(left), we easily verify that   ∇ΓhΦ0(x,y)=12(y/R(x,y)3−1−x/R(x,y)3+1) on Fh=∂Th∩Γh (2.10) such that on $$F_h$$ we have $$| \nabla_{{\it{\Gamma}}_{\!h}}{\it{\Phi}}^0 |^2 = 1 +\frac{1}{2} \frac{1 -2R(x,y)}{R(x,y)^4}$$. With $$R(x,y) \geq 1/\sqrt{2}$$ it holds $$| \nabla_{{\it{\Gamma}}_{\!h}}{\it{\Phi}}^0 | \leq \sqrt{2}-1 < 1/2$$ and Lemma 2.1 can be applied with $$L_\rho =1$$ and $$L_T<1$$ in order to satisfy (A5). Let $$\mathcal{T}_h^*$$ be a refinement of $$\mathcal{T}_h$$ by bisection and denote by $$\mathcal{I}_h^\ast {\it{\Phi}}^0$$ the piecewise affine interpolation of $${\it{\Phi}}^0$$ on $$\mathcal{T}_h^*$$. Then $$\mathcal{T}_h = \left( \mathcal{I}_h^\ast {\it{\Phi}}^0 \right)(\mathcal{T}_h^*)$$ is an interpolating triangulation of $${\it{\Omega}}$$ defining the approximating domain $${\it{\Omega}}_h$$ as shown in Fig. 2(right). $${\it{\Omega}}$$ is obtained from $${\it{\Omega}}_h$$ by applying the boundary deformation $${\it{\Phi}} = {\it{\Phi}}^0 \circ \left( \mathcal{I}_h^\ast {\it{\Phi}}^0 \right)^{-1}$$. (b) A simple example of a domain $${\it{\Omega}}$$ with only piecewise smooth boundary can be obtained from setting $${\it{\Phi}}^0(x,y)= 0$$ for $$x\leq 0$$. (c) The deformation $${\it{\Phi}}^0$$ according to (2.9) can also be used to map the rectangular cuboid $${\it{\Omega}}_h^0 = \{ (x,y,z) \in \mathbb{R}^3 : |x|+|y|\leq 1 , 0\leq z\leq 1 \}$$ to a cylindrical domain $${\it{\Omega}}$$. Note that on the top and bottom side of $${\it{\Omega}}_h$$, the deformation lies in the tangential plane to $${\it{\Gamma}}$$ and $${{\it{\Gamma}}_{\!h}}$$. For the cylindrical domain $${\it{\Omega}}$$, the deformation $${\it{\Phi}}^0$$ cannot be derived from a normal projection within a narrow band around $${\it{\Gamma}}$$. Fig. 2. View largeDownload slide Left: globally smoothly bounded domain $${\it{\Omega}}$$ and its polygonal approximation $${\it{\Omega}}_h^0$$ with a macro triangulation $$\mathcal{T}_h^0$$. Middle: $$\mathcal{T}_h^*$$ obtained from refinement of $$\mathcal{T}_h^0$$ by three steps of bisection. Right: interpolating polygonal domain $${\it{\Omega}}_h$$ with the triangulation $$\mathcal{T}_h$$. Fig. 2. View largeDownload slide Left: globally smoothly bounded domain $${\it{\Omega}}$$ and its polygonal approximation $${\it{\Omega}}_h^0$$ with a macro triangulation $$\mathcal{T}_h^0$$. Middle: $$\mathcal{T}_h^*$$ obtained from refinement of $$\mathcal{T}_h^0$$ by three steps of bisection. Right: interpolating polygonal domain $${\it{\Omega}}_h$$ with the triangulation $$\mathcal{T}_h$$. 2.3 Approximation of curved elements To quantify the deviation of a curved boundary face or edge from the polygonal interpolation, we define the geometric element indicator  λΓ(Th) :=h−1{‖∇^F^(X−J)‖L∞(F^) if Th∩Γh=Fh,‖∇^E^(X−J)‖L∞(E^) if Th∩Γh=Eh,  (2.11) where $$\hat\nabla_{\!\hat{F}}({\mathcal{X}} -J)$$ denotes the tangential gradient of $$({\mathcal{X}} -J)$$ in $$\hat F$$ given by the first $$d-1$$ partial derivatives. Likewise, $$\hat\nabla_{\!\hat{E}}({\mathcal{X}} -J)$$ is the tangential gradient in $$\hat E$$ subject to the first $$d-2$$ partial derivatives. Although $$\lambda_{\it{\Gamma}}(T_h)$$ is only defined on a boundary face or edge by (2.11), it in fact already characterizes the difference between $$T$$ and $$T_h$$ completely, i.e., from Lemma 2.1 and (2.1) we deduce   ‖∇^(X−J)‖L∞(T^) ≲hλΓ(Th), (2.12)  λΓ(Th) ≲LhLT. (2.13) The partial derivatives of $${\mathcal{X}}$$ or $$J$$ with respect to $$\hat x_1, \dots ,\hat x_d$$ build a complete set of linear independent tangential vectors to $${\it{\Gamma}}$$ or $${\it{\Gamma}}_h$$, respectively. For $$\hat{x} \in \hat{F}$$ we define   T(x^) :=∇^F^X(x^)=[∂^1X(x^),⋯,∂^d−1X(x^)], (2.14a)  Th(x^) :=∇^F^J=[∂^1J,⋯,∂^d−1J]. (2.14b) Then, the first fundamental form of $${\it{\Gamma}}$$ and $${\it{\Gamma}}_h$$ is given by the symmetric positive definite matrix   G=TTTandGh=ThTTh, (2.15) respectively. Given $$\hat{x}\in\hat{T}$$ and $$u\in H^1({\it{\Omega}})$$, we set $$\hat{u}(\hat{x}) := u\circ{\mathcal{X}}(\hat{x}) = u(x)$$. Likewise, for $$\hat{x}_{\hat{F}}\in\hat{F}$$ and $$v\in H^1(\partial{\it{\Gamma}} \cap T)$$, we set $$\hat{v}(\hat{x}_{\hat{F}}) := v\circ{\mathcal{X}}(\hat{x}_{\hat{F}}) = v(x_{F})$$. The gradients on $$\hat T$$ and $$\hat F$$ are related to respective gradients on $${\it{\Omega}}$$ and $${\it{\Gamma}}$$ by   ∇^u^=∇u∇^X and ∇^F^v^=∇ΓvT. (2.16) With the change of variables formulas   ∫T^u^det(∇X)dx^=∫X(T^)udx and ∫F^v^det(G)ds^=∫X(F^)vds (2.17) we get the identities   ∫T^∇^u^⋅∇^X−1∇^X−T∇^u^ det(∇X)dx^ =∫X(T^)∇u⋅∇udx, (2.18a)  ∫F^∇^F^v^⋅G−1∇^F^v^det(G)ds^ =∫X(F^)∇Γv⋅∇Γvds. (2.18b) Next, we want to get bounds for the perturbation of the transformation in terms of the geometric element indicator. Lemma 2.3 Let $$T_h\in\mathcal{T}_h$$ and $$T\in\mathcal{T}$$ be the corresponding curved element. Let $${\mathcal{X}}$$, $$\mathbb{T}$$, $$G$$, and $$J$$, $$\mathbb{T}_h$$, $$G_h$$ be defined as above. Suppose there is one face $$F\in{\mathcal{F}}^\partial$$ with $$T\cap{\it{\Gamma}} = F$$. Then,   ‖G−Gh‖ ≲h2λΓ(Th), (2.19a)  |det(G)−det(Gh)| ≲h2d−2λΓ(Th), (2.19b) where $$\|\cdot\|$$ denotes a compatible matrix norm. Moreover, it holds   ‖det(G)det(Gh)−1‖L∞(F^) ≲λΓ(Th), (2.19c)  ‖Th(detGdetGhG−1−Gh−1)ThT‖L∞(F^) ≲λΓ(Th). (2.19d) Proof. For any $$z\in\mathbb{R}^d$$ it holds $$z^T G z = \| \hat \nabla_{\!\hat F} {\mathcal{X}}\,z \|^2$$, and hence the eigenvalues of $$G$$ are bounded by $$h^{-2}L^{-2}|z|^2 \leq z^T G z \leq h^{-2}L^2|z|^2$$ due to (2.8). Analogously, we conclude that the eigenvalues of $$G_h$$ lie in $$[h L_h^{-1},h L_h]$$. We thus infer $$\| \mathbb{T} \|_{L^2(\hat F)}, \| \mathbb{T}_h \|_{L^2(\hat F)} \lesssim h$$ and $$\det(G), \det(G_h) \lesssim h^{2d-2}$$. By the definition of $$G$$, $$G_h$$ and $$\lambda_{\it{\Gamma}}(T_h)$$, we obtain   ‖G−Gh‖ =‖(T−Th)TT+ThT(T−Th)‖ ≤hλΓ(Th)(‖T‖+‖Th‖), (2.20) which yields the first assertion. We use a Taylor expansion to deduce that there exists $$0 \leq \theta \leq 1$$ such that   det(G)−det(Gh) =θ(Ddet(Gh))(G−Gh) (2.21)   =θdet(Gh)trace⁡(Gh−1(G−Gh)). (2.22) By Hölder’s inequality we conclude that   |det(G)−det(Gh)|≲|det(Gh)|‖Gh−1‖‖G−Gh‖ (2.23) to infer (2.19b). For (2.19c), we note that   det(G)−det(Gh)=det(G)−det(Gh)det(G)+det(Gh) (2.24) and that all powers of $$h$$ cancel after division by $$\det(G_h)$$. Finally, using the abbreviations $$q:=\sqrt{\det(G)}$$ and $$q_h:=\sqrt{\det(G_h)}$$ leads to   Th(qqhG−1−Gh−1)ThT=Th(q−qhqhG−1+G−1(Gh−G)Gh−1)ThT. (2.25) The application of (2.19a) and (2.19c) now concludes the proof. □ Corollary 2.4 Let $$T_h\in\mathcal{T}_h$$ and $$T\in\mathcal{T}$$ be the corresponding curved element. Let $${\mathcal{X}}$$ and $$J$$ be defined as above. We set $$A := \hat\nabla{\mathcal{X}}^T \hat\nabla{\mathcal{X}}$$ and $$A_h := \hat\nabla J^T \hat\nabla J$$. Suppose there is one face $$F\in{\mathcal{F}}^\partial$$ or edge $$E\in\mathcal{E}_h^\partial$$ such that $$T\cap{\it{\Gamma}} = F$$ or $$T\cap{\it{\Gamma}} = E$$, respectively. Then,   ‖det(∇^X)det(∇^J)−1‖L∞(T^) ≲λΓ(Th), (2.26a)  ‖∇^J(|det∇^X||det∇^J|A−1−Ah−1)∇^JT‖L∞(T^) ≲λΓ(Th). (2.26b) Proof. The assertions follow from the same arguments as in Lemma 2.3. □ 2.4 Lifts and geometric nonconformity estimates We first define the lifting of functions defined on the curved domain or surface onto their corresponding polyhedral or polygonal discretization and vice versa. Definition 2.5 For functions $$u_h:{\it{\Omega}}_h\to\mathbb{R}$$ and $$u:{\it{\Omega}}\to\mathbb{R}$$ we define the lift and inverse lift by   uhℓ(X∘J−1(x)):=uh(x) and u−ℓ(x):=uh(X∘J−1(x)) for x∈Ωh. (2.27) Likewise, for $$v_h:{\it{\Gamma}}_h\to\mathbb{R}$$ and $$v:{\it{\Gamma}}\to\mathbb{R}$$, we define   vhℓ(X∘J−1(x)):=vh(x) and v−ℓ(x):=v(X∘J−1(x)) for x∈Γh. (2.28) Next, we estimate the error which arises due to integration on perturbed domains in order to compare functions on the curved domain and on its polyhedral approximation. Lemma 2.6 For $$u_h, \phi^{-\ell}\in H^1({\it{\Omega}}_h)$$ it holds   |∫Ωuhℓϕdx−∫Ωhuhϕ−ℓdx| ≲∑Th∈ThλΓ(Th)‖uh‖L2(Th)‖ϕ−ℓ‖L2(Th), (2.29a)  |∫Ω∇uhℓ⋅∇ϕdx−∫Ωh∇uh⋅∇ϕ−ℓdx| ≲∑Th∈ThλΓ(Th)‖∇uh‖L2(Th)‖∇ϕ−ℓ‖L2(Th). (2.29b) Likewise, for $$v_h, \psi^{-\ell}\in H^1({{\it{\Gamma}}_{\!h}})$$ it holds   |∫Γvhℓψds−∫Γhvhψ−ℓds| ≲∑Th∈ThλΓ(Th)‖vh‖L2(Fh)‖ψ−ℓ‖L2(Fh), (2.29c)  |∫Γ∇Γvhℓ⋅∇Γψds−∫Γh∇Γhvh⋅∇Γhψ−ℓds| ≲∑Th∈ThλΓ(Th)‖∇Γvh‖L2(Fh)‖∇Γψ−ℓ‖L2(Fh). (2.29d) Proof. The identities follow from transformation to the reference domain and a pull-back. We start with the $$L^2$$ cases. By a localization to elements and with (2.17) we derive   ∫Ωuhℓϕdx−∫Ωhuhϕ−ℓdx =∑Thi∈Th∫T^u^h(|det∇^Xi|−|det∇^Ji|)ϕ^dx^ =∑Thi∈Th∫Thiuh(|det∇^Xi||det∇^Ji|−1)ϕ−ℓdx. Application of (2.26a) yields the assertions (2.29a). A similar argument and (2.19c) can be employed with the surface integral to get (2.29c). One can proceed in an analogous way for the $$H^1$$ cases. With (2.17) and (2.18), we obtain   ∫Ω∇uhℓ⋅∇ϕdx −∫Ωh∇uh⋅∇ϕ−ℓdx =∑Thi∈Th∫T^∇^u^h((Ai)−1 |det∇^Xi|−(Ahi)−1 |det∇^Ji|)⋅∇^ϕ^dx^ =∑Thi∈Th∫Thi∇uh∇^Ji(|det∇^X||det∇^J|(Ai)−1−(Ahi)−1)(∇^Ji)T⋅∇ϕ−ℓdx. Application of (2.26b) yields the assertions (2.29b) and, again, similar arguments and (2.19d) can be employed with the surface integrals to get (2.29d). □ For sufficiently fine triangulations, it is known that there is a norm equivalence for the lifting of functions in the volume (cf. Elliott & Ranner, 2013) and on the surface (see Dziuk, 1988; Demlow, 2009). Here, we can get the equivalence from Lemma 2.6. Lemma 2.7 Assume the triangulation $$\mathcal{T}_h$$ is sufficiently fine such that $$L_h L_T <1$$. Then the global geometry estimator satisfies $$\lambda_{\it{\Gamma}} := \max_{T_h\in\mathcal{T}_h}\lambda_{\it{\Gamma}}(T_h) <1$$ and for $$u\in H^1({\it{\Omega}})$$ and $$v\in H^1({\it{\Gamma}})$$, it holds   ‖u‖L2(Ω) ≈‖u−ℓ‖L2(Ωh),‖∇u‖L2(Ω) ≈‖∇u−ℓ‖L2(Ωh), (2.30a)  ‖v‖L2(Γ) ≈‖v−ℓ‖L2(Γh),‖∇Γv‖L2(Γ) ≈‖∇Γhv−ℓ‖L2(Γh), (2.30b) where $$a\approx b$$ means $$a\lesssim b$$ and $$b\lesssim a$$. Proof. Choosing $$\phi = u_h^{\ell}$$ in Lemma 2.6 yields   |‖uhℓ‖L2(Ω)2−‖uh‖L2(Ωh)2|≲λΓ‖uh‖L2(Ωh)2. (2.31) It follows $$\| u_h^\ell \|_{L^2({\it{\Omega}})}^2 \lesssim (1 +\lambda_{{\it{\Gamma}}}) \, \| u_h \|_{L^2({\it{\Omega}}_h)}^2$$ and $$\| u_h \|_{L^2({\it{\Omega}}_h)}^2 \lesssim (1 -\lambda_{{\it{\Gamma}}})^{-1} \, \| u_h^\ell \|_{L^2({\it{\Omega}})}^2$$. In an analogous way one can show (2.30b). □ 3. Problem formulation The finite element method for (1.1) is based on the weak formulation introduced next. Subsequently, we state the corresponding discrete problem for which we apply the notion of an exact triangulation and its simplicial approximation as described in the previous section. 3.1 Weak formulation of the continuous problem Let $${\it{\Omega}}_h^0 \subseteq \mathbb{R}^d$$ be a polyhedral domain with a triangulation $$\mathcal{T}_h^0({\it{\Omega}}_h^0)$$ that satisfies (A1)-(A3). Let $${\it{\Omega}}$$ be the domain resulting from a deformation of $${\it{\Omega}}_h^0$$ by the continuous mapping $${\it{\Phi}}^0$$ that satisfies (A4)-(A5). For each boundary facet $$F_h^i$$ of $$\mathcal{T}_h^{0}$$, the restriction of $${\it{\Phi}}^0$$ to $$F_h^i$$ is a $$C^1$$ function due to (A4). We set $${\it{\Gamma}}^i := {\it{\Phi}}^0(F_h^i)$$. Then $$\partial {\it{\Omega}} = \cup_i {\it{\Gamma}}^i$$ and each boundary patch $${\it{\Gamma}}^i$$ has a piecewise $$C^1$$ boundary $$\partial {\it{\Gamma}}^i$$. We introduce the spaces   U :=H1(Ω), (3.1a)  V :={v∈L2(Γ):∇Γv∣Γi∈L2(Γi), v∣Γi=v∣Γj on ∂Γi∩∂Γj} (3.1b) and denote their dual spaces by $$U^*$$ and $$V^*$$, respectively. Moreover, we define the weak form of the Laplace–Beltrami operator $$-{\it{\Delta}}_{{\it{\Gamma}}} : V \to \mathbb{R}$$ such that for each $$v\in V$$ it holds   ∫Γ−ψΔΓvds:=∫Γ∇Γψ⋅∇Γvds−∑Γi=1NΓ∫∂Γiψ∇Γv⋅nΓidσfor all ψ∈V. (3.2) Note that in general the sum in (3.2) does not vanish, because on $$\bar {\it{\Gamma}}^i \cap \bar {\it{\Gamma}}^j$$ ($$i\neq j$$) the co-normal vectors $$n_{\it{\Gamma}}^i$$ and $$n_{\it{\Gamma}}^j$$ are not parallel unless $${\it{\Gamma}}$$ is a global $$C^1$$ surface. Therefore, the additional condition (1.1d) is necessary to recover from (3.2) the strong form of the Laplace–Beltrami operator for piecewise $$C^2$$ surfaces. The weak form of (1.1) reads: Given $$f\in U^*$$ and $$g\in V^*$$, find $$(u,v) \in U\times V$$ such that, for all $$(\phi,\psi) \in U\times V$$,   ∫Ω∇ϕ⋅∇u+ϕudx+∫Γϕ(αu−βv)ds =∫Ωϕfdx, (3.3a)  ∫Γ∇Γψ⋅∇Γv+vψds−∫Γψ(αu−βv)ds =∫Γψgds. (3.3b) This can also be written in the standard variational form   a((u,v),(ϕ,ψ))=ℓ((ϕ,ψ))for all (ϕ,ψ)∈U×V, (3.4) with the bilinear and linear forms   a((u,v),(ϕ,ψ)) :=α∫Ω∇ϕ⋅∇u+ϕudx+β∫Γ∇Γψ⋅∇Γv+ψvds +∫Γ(αϕ−βψ)(αu−βv)ds, (3.5a)  ℓ((ψ,ϕ)) :=α∫Ωϕfdx+β∫Γψgds. (3.5b) In an analogous way to Elliott & Ranner (2013), we easily verify that $$a(\cdot,\cdot)$$ is a continuous and coercive bilinear form on $$U\times V$$. Hence, the Lax–Milgram theorem provides existence of a unique solution to (3.4). 3.2 Discrete problem Let $$\mathcal{T}_h$$ be a triangulation obtained from repeated refinement of $$\mathcal{T}_h^0$$ by bisection and subsequent deformation such that (A1)-(A5) are satisfied. $$\mathcal{T}_h$$ defines a polyhedral domain $${\it{\Omega}}_h\subset\mathbb{R}^{d}$$ with the boundary $${{\it{\Gamma}}_{\!h}} = \partial{\it{\Omega}}_h$$. We introduce the spaces of continuous finite element functions   Uh :={uh∈C(Ωh¯) :uh|Th∈P1(Th)∀Th∈Th}, (3.6a)  Vh :={vh∈C(Γh) : vh|Fh∈P1(Fh)∀Fh∈Fh∂}, (3.6b) where $$P_1$$ is the space of piecewise polynomials of maximal degree one. The discrete variational problem reads: Given $$f_h \in U_h$$ and $$g_h \in V_h$$ approximating $$f$$ and $$g$$, respectively, find $$(u_h,v_h)\in U_h\times V_h$$ such that   ah((uh,vh),(ϕh,ψh))=ℓh((ϕh,ψh))for all (ϕh,ψh)∈Uh×Vh (3.7) with the discrete bilinear and linear forms   ah((uh,vh),(ϕh,ψh)) :=α∫Ωh∇ϕh⋅∇uh+ϕhuhdx+β∫Γh∇Γhψh⋅∇Γhvh+ψhvhds +∫Γh(αϕh−βψh)(αuh−βvh)ds, (3.8a)  ℓh((ϕh,ψh)) :=α∫Ωhϕhfhdx+β∫Γhψhghds. (3.8b) Existence and uniqueness of solutions to (3.7) can be proven by the Lax–Milgram Theorem (cf. Elliott & Ranner, 2013). Note that the finite element method defined above is nonconforming since $$U_h$$ and $$V_h,$$ in general, are not subspaces of $$U$$ and $$V$$. Moreover, $$a_h$$ and $$\ell_h$$ differ from $$a$$ and $$\ell$$ since the integrals are defined over different domains. 4. Residual error estimation This section is devoted to the derivation of a residual estimator of the error measured in the energy norm for the coupled problem. We introduce the energy norm on $$U \times V$$ by   |||(ϕ,ψ)|||2:=a((ϕ,ψ),(ϕ,ψ))=α‖ϕ‖U2+β‖ψ‖V2+‖αϕ−βψ‖L2(Γ)2. (4.1) Given some finite element function $$(u_h,v_h) \in U_h\times V_h$$, the residual $$R$$ with regard to the variational problem (3.4) is defined by   ⟨R,(ϕ,ψ)⟩:=a((uhℓ,vhℓ),(ϕ,ψ))−ℓ(ϕ,ψ) for all (ϕ,ψ)∈U×V. (4.2) Proposition 4.1 Let $$(u,v)\in U\times V$$ be the solution of (3.4). Given $$(u_h,v_h) \in U_h\times V_h$$, the energy norm of the error satisfies   |||(u−uhℓ,v−vhℓ)|||=|||R|||∗. (4.3) Proof. We set $$e:= (u_h^{\ell} -u, v_h^{\ell}-v)$$ and easily check that for all $$(\phi,\psi) \in U\times V$$ we have $$\langle R,(\phi,\psi) \rangle = a(e,(\phi,\psi))$$. Hölder inequality directly implies $$\langle R,(\phi,\psi)\rangle \leq \left\lvert\!\left\lvert\!\left\lvert{e}\right\rvert\!\right\rvert\!\right\rvert\left\lvert\!\left\lvert\!\left\lvert{(\phi,\psi)}\right\rvert\!\right\rvert\!\right\rvert$$ for all $$(\phi,\psi) \in U\times V$$, and thus $$\left\lvert\!\left\lvert\!\left\lvert{R}\right\rvert\!\right\rvert\!\right\rvert_\ast \leq \left\lvert\!\left\lvert\!\left\lvert{e}\right\rvert\!\right\rvert\!\right\rvert$$. On the other hand, since $$\langle R,e\rangle = \left\lvert\!\left\lvert\!\left\lvert{e}\right\rvert\!\right\rvert\!\right\rvert^2$$, it follows   |||e|||=⟨R,e⟩|||e|||≤supw∈(U×V)∖{0}⟨R,w⟩|||w|||=|||R|||∗. □ We split the residual into $$R = R_h +R_\textrm{nc} +R_\textrm{osc}$$ with   ⟨Rh,(ϕ,ψ)⟩ :=ah((uh,vh),(ϕ−ℓ,ψ−ℓ))−ℓh(ϕ−ℓ,ψ−ℓ), (4.4a)  ⟨Rnc,(ϕ,ψ)⟩ :=a((uhℓ,vhℓ),(ϕ,ψ))−ah((uh,vh),(ϕ−ℓ,ψ−ℓ)) −ℓhℓ(ϕ,ψ) +ℓh(ϕ−ℓ,ψ−ℓ), (4.4b)  ⟨Rosc,(ϕ,ψ)⟩ :=ℓhℓ(ϕ,ψ)−ℓ(ϕ,ψ). (4.4c) $$R_h$$ is the discretization part of the residual, $$R_\textrm{nc}$$ is the nonconformity part which quantifies the quality of the polyhedral domain approximation $${\it{\Omega}}_h$$ with respect to the exact domain $${\it{\Omega}}$$ and $$R_\textrm{osc}$$ is the data approximation part of the residual, where $$\ell_h^\ell(\phi,\psi) := \alpha \int_{{\it{\Omega}}} f_h^{\ell} \phi \,\mathrm{d}x + \beta\int_{{\it{\Gamma}}} g_h^{\ell} \psi \,\mathrm{d}s$$. In the subsequent sections, computable error bounds for $$R_h$$ and $$R_\textrm{nc}$$ are derived and then combined to yield the total a posteriori error estimate. 4.1 Discretization residual The derivation of the upper bound for the discretization residual $$R_h$$ largely follows the classical approach of residual error estimators (see, e.g., Verfürth, 1996; Ainsworth & Oden, 2000). A recent account on the universal application of the principle is given by Carstensen et al. (2012). In the proof, we employ the stable interpolation operator $${\it{\Pi}}_h:U^{-\ell} = H^1({\it{\Omega}}_h)\to U_h$$ of Scott & Zhang (1990). For $$\phi\in H^1({\it{\Omega}}_h)$$ and any $$T_h\in\mathcal{T}_h$$, $$F_h\in{\mathcal{F}}_h$$, one verifies   ‖ϕ−Πhϕ‖L2(Th) ≲hT‖∇ϕ‖L2(ωTh), (4.5a)  ‖ϕ−Πhϕ‖L2(Fh) ≲hF1/2‖∇ϕ‖L2(ωFh), (4.5b) where and $$\omega_{F_h}$$ denotes the volume patch consisting of those two elements which share the common face $$F_h$$, and $$\omega_{T_h}$$ is the volume patch containing $$T_h$$ as well as all neighbouring elements sharing at least on common node with $$T_h$$. In an analogous way, we introduce the stable interpolation operator $${\it{\Pi}}_h^\partial:V^{-\ell} \to V_h$$ on the polyhedral boundary $${\it{\Gamma}}_h$$. It holds   ‖ψ−Πh∂ψ‖L2(Fh) ≲hF‖∇Γhψ‖L2(ωFh∂), (4.6a)  ‖ψ−Πh∂ψ‖L2(Eh) ≲hE1/2‖∇Γhψ‖L2(ωEh∂), (4.6b) with surface patches $$\omega_{F_h}^\partial \subset {{\it{\Gamma}}_{\!h}}$$ and $$\omega_{E_h}^\partial \subset {{\it{\Gamma}}_{\!h}}$$. Lemma 4.2 Let $$(u_h,v_h) \in U_h\times V_h$$ be a solution of the discrete problem (3.7). For all $$(\phi,\psi) \in U\times V$$, the discretization residual $$R_h$$ satisfies the estimate   ⟨Rh,(ϕ,ψ)⟩≲α(ηB+ηC)‖∇ϕ‖L2(Ω)+βηS‖∇Γψ‖L2(Γ), (4.7) where the error indicators on the right-hand side are given by   ηC :=‖hFh∂1/2(αuh−βvh+∂nuh)‖L2(Fh∂), (4.8a)  ηB :=‖hTh(Δuh−uh+fh)‖L2(Th) +‖hFh1/2[∂nuh]Fh‖L2(Fh∖Fh∂), (4.8b)  ηS :=‖hFh∂(ΔΓhvh−(1+β)vh+αuh+gh)‖L2(Fh∂) +‖hEh∂1/2[[∂nvh]]E∂‖L2(Eh∂). (4.8c) Here, the discrete norms are defined as summation over the elements or faces, respectively. By $$[\partial_n u_h]_{{\mathcal{F}}_h}$$ we denote the usual jump of the normal derivative at element faces, and $$[[ \partial_n v_h ]]_{\mathcal{E}^\partial}$$ is the corresponding jump quantity on the surface defined as in (1.1d). Proof. Let $$(\phi,\psi) \in U\times V$$, then $$({\it{\Pi}}_h\phi^{-\ell},{\it{\Pi}}^\partial_h\psi^{-\ell})\in U_h\times V_h$$ is in the kernel of $$R_h$$. Thus, a splitting of $$R_h$$ into element contributions and an integration by parts result in   ⟨Rh,(ϕ,ψ)⟩ =α∑Th∈Th(∫Th(−Δuh+uh−fh)(ϕ−ℓ−Πhϕ−ℓ)dx +∑Fh∈∂Th∫Fh∂nuh(ϕ−ℓ−Πhϕ−ℓ)ds) +α∑Fh∈Fh∂∫Fh(αuh−βvh)(ϕ−ℓ−Πhϕ−ℓ)ds +β∑Fh∈Fh∂(∫Fh(−ΔΓhvh+(1+β)vh−αuh−gh)(ψ−ℓ−Πh∂ψ−ℓ)ds +∑Eh∈∂Fh∫Eh∂nvh(ψ−ℓ−Πh∂ψ−ℓ)dσ). (4.9) In the sum over all $$T_h\in\mathcal{T}_h$$ we collect all contributions over the inner faces $$F_h \in {\mathcal{F}}_h\setminus{\mathcal{F}}_h^\partial$$ and move the boundary face parts to the coupling term. The Cauchy–Schwarz inequality yields   ⟨Rh,(ϕ,ψ)⟩ ≲α∑Th∈Th‖Δuh−uh+fh‖L2(Th)‖ϕ−ℓ−Πhϕ−ℓ‖L2(Th) +α∑Fh∈Fh∖Fh∂‖[∂nuh]Fh‖L2(Fh)‖ϕ−ℓ−Πhϕ−ℓ‖L2(Fh) +α∑Fh∈Fh∂‖αuh−βvh+∂nuh‖L2(Fh)‖ϕ−ℓ−Πhϕ−ℓ‖L2(Fh) +β∑Fh∈Fh∂‖ΔΓhvh−(1+β)vh+αuh+gh‖L2(Fh)‖ψ−ℓ−Πh∂ψ−ℓ‖L2(Fh) +β∑Eh∈Eh∂‖[[∂nvh]]Eh∂‖L2(Eh)‖ψ−ℓ−Πh∂ψ−ℓ‖L2(Eh). (4.10) Then the properties of the interpolation (4.5) and (4.6), and the norm equivalence of the lifted functions according to Lemma 2.7 complete the proof. □ 4.2 Nonconformity residual The geometric element indicator $$\lambda(T)$$ and the comparison results of Section 2 lead to an upper bound for $$R_\textrm{nc}$$ which is due to the polyhedral approximation of the curved domain $${\it{\Omega}}$$. Lemma 4.3 Assume that the triangulation $$\mathcal{T}_h$$ is sufficiently fine such that Lemma 2.7 holds. For any $$(\phi,\psi)\in U\times V$$, the nonconformity residual $$R_{\textrm{nc}}$$ satisfies   ⟨Rnc,(ϕ,ψ)⟩≲ηncB,0‖ϕ‖L2(Ω)+ηncB,1‖∇ϕ‖L2(Ω) +ηncC‖αϕ−βψ‖L2(Γ) +ηncS,0‖ψ‖L2(Γ)+ηncS,1‖∇Γψ‖L2(Γ), where the element indicators are given by   ηncB,0 :=α∑Th∈ThλΓ(Th)‖uh−fh‖L2(Th), (4.11a)  ηncB,1 :=α∑Th∈ThλΓ(Th)‖∇uh‖L2(Th), (4.11b)  ηncC :=∑Fh∈Fh∂λΓ(Fh)‖αuh−βvh‖L2(Fh), (4.11c)  ηncS,0 :=β∑Fh∈Fh∂λΓ(Fh)‖vh−gh‖L2(Fh), (4.11d)  ηncS,1 :=β∑Fh∈Fh∂λΓ(Fh)‖∇Γhvh‖L2(Fh). (4.11e) Proof. Let $$(\phi,\psi)\in U\times V$$. We rearrange terms to split $$\langle R_\textrm{nc}, (\phi,\psi) \rangle = \mathcal{I}_B +\mathcal{I}_S +\mathcal{I}_C$$ with   IB =α(∫Ω∇uhℓ⋅∇ϕ+(uhℓ−fhℓ)ϕdx−∫Ωh∇uh⋅∇ϕ−ℓ+(uh−fh)ϕ−ℓdx),IS =β(∫Γ∇Γvhℓ⋅∇Γψ+(vhℓ−ghℓ)ψds−∫Γh∇Γhvh⋅∇Γhψ−ℓ+(vh−gh)ψ−ℓds),IC =(∫Γ(αuhℓ−βvhℓ)(αϕ−βψ)ds−∫Γh(αuh−βvh)(αϕ−ℓ−βψ−ℓ)ds). To estimate the coupling term $$\mathcal{I}_C$$, we split the integral into its contributions from the surface elements and apply (2.17) and Lemma 2.6. This yields   IC =α∑Fh∈Fh∂∫Fh(det(G)det(Gh)−1)(αuh−βvh)(αϕ−ℓ−βψ−ℓ)ds ≲α∑Fh∈Fh∂λΓ(Th)‖αuh−βvh‖L2(Fh)‖αϕ−ℓ−βψ−ℓ‖L2(Γh). In an analogous way, we estimate the surface part $$\mathcal{I}_S$$, where we also use (2.18) for the gradient terms. It follows   IS =β∑Fh∈Fh∂∫Fh∇vh(det(G)det(Gh)T(G−1−Gh−1)TT)⋅∇ψ−ℓds +β∑Fh∈Fh∂∫Fh(det(G)det(Gh)−1)(vh−gh)ψ−ℓds ≲β∑Fh∈Fh∂λΓ(Th)(‖∇Fhvh‖L2(Fh)‖∇Fhψ−ℓ‖L2(Γh)+‖vh−gh‖L2(Fh)‖ψ−ℓ‖L2(Γh)). Likewise, the bulk term can be bounded by   IB ≲α∑Th∈ThλΓ(Th)(‖∇uh‖L2(Th)‖∇ϕ−ℓ‖L2(Ωh)+‖uh−fh‖L2(Th)‖ϕ−ℓ‖L2(Ωh)). Application of Lemma 2.7 completes the proof. □ 4.3 Overall error estimator As a direct consequence of the preceding Sects. 4.1 and 4.2, we devise an a posteriori estimator for the overall error of a finite element solution of the coupled problem, and prove its reliability and efficiency. In the subsequent error analysis, we employ piecewise constant approximations of the data. We define $$f_h:{\it{\Omega}}_h \to \mathbb{R}$$ by $$f_h\vert_{T_h} = \left\lvert{T}\right\rvert^{-1}\,\int_{T} f \,\mathrm{d}x$$ for each $$T_h\in \mathcal{T}_h$$, where $$T\in\mathcal{T}$$ is the corresponding curved element. In an analogous way, for each $$F_h\in {\mathcal{F}}_h^\partial$$, we set $$g_h\vert_{F_h} := \left\lvert{F}\right\rvert^{-1}\,\int_{F} g \,\mathrm{d}s$$, where $$F$$ is the corresponding curved boundary face. The oscillations of the data $$f$$ and $$g$$ are defined by   osc(f,Th) :=hTh‖f−ℓ−fh‖L2(Th)andosc(f,Th):=∑Th∈Thosc(f,Th), (4.12a)  osc(g,Fh) :=hFh‖g−ℓ−gh‖L2(Fh)andosc(g,Fh∂):=∑Fh∈Fh∂osc(g,Fh). (4.12b) Theorem 4.4 Assume that the triangulation $$\mathcal{T}_h$$ is sufficiently fine such that Lemma 2.7 holds. For the solution $$(u_h,v_h)\in U_h\times V_h$$ of the discrete problem (3.7), it holds   |||(u−uhℓ,v−vhℓ)|||≲ηh+ηnc+osch, (4.13)  ηh≲|||(u−uhℓ,v−vhℓ)|||+ηnc+osch, (4.14) with the discretization error indicator of (4.8) and the nonconformity indicator according to (4.11)   ηh :=ηB+ηC+ηS,ηnc :=hTηncB,0+ηncB,1+hF∂ηncS,0+ηncS,1+ηncC, and data oscillations $$\text{osc}_h := \text{osc}(\,f,{\mathcal{T}_h}) + \text{osc}(g,{{\mathcal{F}}_h})$$. Proof. We easily verify that whenever $$\phi$$ or $$\psi$$ is constant on $$T\in\mathcal{T}$$ or $$F\in{\mathcal{F}}^\partial$$, respectively, the local contribution to $$R_\textrm{osc}$$ vanishes. In an analogous way to $$f_h$$ and $$g_h$$, we define the elementwise mean values $$\phi_h$$ and $$\psi_h$$ of $$\phi$$ and $$\psi$$, respectively. Then for all $$(\phi,\psi) \in U\times V$$, it holds   ⟨Rosc,(ϕ,ψ)⟩ ≲∑T∈T‖fhℓ−f‖L2(T)‖ϕ−ϕh‖L2(T)+∑F∈F∂‖ghℓ−g‖L2(T)‖ψ−ψh‖L2(F) ≲∑Th∈ThhTh‖fh−f−ℓ‖L2(Th)‖∇ϕ−ℓ‖L2(Th) +∑Fh∈Fh∂hFh‖gh−g−ℓ‖L2(Th)‖∇Γψ−ℓ‖L2(Fh) ≲∑Th∈Thosc(f,Th)‖∇ϕ‖L2(Ω)+∑Fh∈Fh∂osc(g,Fh)‖∇Γψ‖L2(Γ). Then, the first part (reliability) is an immediate consequence of Lemmas 4.2 and 4.3. To prove the efficiency, we apply the standard technique due to Verfürth. We have to consider all terms contributing to the indicator $$\eta_h$$, and thus carry out the proof in several steps. (i) Efficiency of$$\boldsymbol{\eta_B}$$ For any $$T_h\in\mathcal{T}_h$$ and for any $$F_h\in{\mathcal{F}}_h$$, we introduce the (volume) bubble functions $$b_{T_h}$$ and $$b_{F_h}$$ with $$\text{supp} b_{T_h} = T_h$$ and $$\text{supp} b_{F_h} = \omega_{F_h}$$, where $$\omega_{F_h}$$ consists of the elements adjacent to $$F_h$$. The bubble functions are normalized such that $$\max_{x\in T_h} b_{T_h}(x)=1$$ and $$\max_{x\in F_h} b_{F_h}(x)=1$$. For any $$T_h\in\mathcal{T}_h$$ and $$F_h\in{\mathcal{F}}_h$$, let   rThB :=Δuh−uh+fh,rFhB :=[∂nuh],ϕTh :=bThrThB,ϕFh :=bFhrFhB. By the equivalence of norms in finite dimensional spaces and the properties of the bubble functions, one can show that   ‖rThB‖L2(Th)2 ≈∫ThrThBϕThdxand‖rFhB‖L2(Fh)2 ≈∫FhrFhBϕFhds, (4.15)  ‖ϕTh‖L2(Th) ≲‖rThB‖L2(Th)and‖ϕFh‖L2(ωFh) ≲hFh1/2‖rFhB‖L2(Fh). (4.16) Details can be found in Verfürth, 1996, 2013; Brenner & Carstensen, 2004. To estimate $${r^B_{T_h}}$$, we use (1.1) on the curved element $$T\in\mathcal{T}$$ corresponding to $$T_h\in\mathcal{T}_h$$. When integrating by parts, the boundary terms vanish due to the bubble functions. We get   ‖rThB‖L2(Th)2 ≲∫Th−∇uh⋅∇ϕTh−uhϕThdx+∫ThfhϕThdx +∫T∇u⋅∇ϕThℓ+uϕThℓdx−∫TfϕThℓdx. For the integrals over $$T_h$$, we insert the respective integrals over $$T$$ and apply Lemma 2.6   ‖rThB‖L2(Th)2 ≲∫T∇(u−uhℓ)⋅∇ϕThℓ+(u−uhℓ)ϕThℓdx+∫T(fhℓ−f)ϕThℓdx +λΓ(T)(‖∇uh‖L2(Th)‖∇ϕTh‖L2(Th)+‖uh−fh‖L2(Th)‖ϕTh‖L2(Th)). We apply the Cauchy–Schwarz inequality and then use Lemma 2.7 for the inverse lifting of $$\left\lVert{ \phi_{T_h}^{\ell} }\right\rVert_{L^2(T)}$$. Then, an inverse inequality and (4.16) yield   ‖rThB‖L2(Th) ≲hTh−1(‖u−uhℓ‖H1(T)+osc(f,T)) +λΓh(T)(hTh−1‖∇uh‖L2(Th)+‖uh−fh‖L2(Th)). On the inner facets $$F\in{\mathcal{F}}_h\setminus{\mathcal{F}}_h^\partial$$ we proceed similarly. After partial integration, we add (1.1), i.e.,   ‖rFhB‖L2(Fh)2 ≲∫ωFhΔuhϕFhdx+∫ωFh∇uh⋅∇ϕFhdx ≲∫ωFh(rThB+uh−fh)ϕFhdx+∫ωFh∇uh⋅∇ϕFhdx +∫ωF(f−u)ϕFhℓdx−∫ωF∇u⋅∇ϕFhℓdx. Lemma 2.6 leads to   ‖rFhB‖L2(Fh)2 ≲∫ωFhrThBϕFhdx +∫ωF(f−fhℓ)ϕFhℓdx +∫ωF∇(uhℓ−u)⋅∇ϕFhℓ+(uhℓ−u)ϕFhℓdx +λΓ(T)(‖∇uh‖L2(ωFh)‖∇ϕFh‖L2(ωFh)+‖uh−fh‖L2(ωFh)‖ϕFh‖L2(ωFh)) ≲‖rThB‖L2(ωFh)‖ϕFh‖L2(ωFh) +‖f−fhℓ‖L2(ωF)‖ϕFhℓ‖L2(ωF) +‖uhℓ−u‖H1(ωF)‖ϕFhℓ‖H1(ωF) +λΓ(T)(‖∇uh‖L2(ωFh)‖∇ϕFh‖L2(ωFh)+‖uh−fh‖L2(ωFh)‖ϕFh‖L2(ωFh)). With the Cauchy–Schwarz inequality, Lemma 2.7, an inverse inequality and (4.16), we deduce   ‖rFhB‖L2(Fh) ≲‖hTh1/2rThB‖L2(ωFh)+‖hTh−1/2(u−ℓ−uh)‖H1(ωF)+‖hTh1/2(f−ℓ−fh)‖L2(ωFh) +∑Th∈ωFhλΓ(T)(hTh−1/2‖∇uh‖L2(Th)+hTh1/2‖uh−fh‖L2(Th)). Since $$\eta_B = \left\lVert{h_{\mathcal{T}_h} {r^B_{T_h}}}\right\rVert_{L^2({\mathcal{T}_h})} + \left\lVert{h_{\mathcal{T}_h}^{1/2} {r^B_{F_h}}}\right\rVert_{L^2({{\mathcal{F}}_h\setminus {\mathcal{F}}_h^\partial})}$$, combining the previous estimates yields   ηB ≲‖u−uhℓ‖H1(Ω)+osc(f,T)+∑Th∈ThλΓ(T)(‖∇uh‖L2(Th)+hTh‖uh−fh‖L2(Th)). ≲|||(u−uhℓ,v−vhℓ)|||+osch+hTηncB,0+ηncB,1 (ii) Efficiency of$$\eta_C$$ For any boundary face $$F_h\in {\mathcal{F}}_h^\partial$$, let   rFhC:=αuh−βvh+∂nuhandϕFh:=bFhrFhC. Then, we estimate as before   ‖rFhC‖L2(Fh)2 ≲∫FhrFhCϕFhds−∫F(αu−βv+∂nu)ϕFhℓds ≲∫Th(rThB+uh−fh)ϕFhdx+∫Th∇uh⋅∇ϕFhdx+∫Fh(αuh−βvh)ϕFhds −∫T(u−f)ϕFhℓdx−∫T∇u⋅∇ϕFhℓdx−∫F(αu−βvh)ϕFhℓds ≲(‖rThB‖L2(Th)+‖f−fhℓ‖L2(T))‖ϕFh‖L2(T) +‖u−uhℓ‖H1(T)‖ϕFh‖H1(T)  +‖α(uhℓ−u)−β(vhℓ−v)‖L2(F)‖ϕFhℓ‖L2(Fh) +λΓ(F)(‖uh−fh‖L2(Th)‖ϕFh‖L2(Th)+‖∇uh‖L2(Th)‖∇ϕFh‖L2(Th) +‖(αuh−βvh)‖L2(Fh)‖ϕFh‖L2(Fh)). A trace inequality and an inverse estimate yield $$\left\lVert{\phi_{F_h}}\right\rVert_{L^2(F_h)} \lesssim h_{F_h}^{-1/2} \left\lVert{\phi_{F_h}}\right\rVert_{L^2(T_h)}$$. We conclude as before   ‖rFhC‖L2(Fh) ≲‖hTh1/2rThB‖L2(Th)+‖hTh1/2(f−ℓ−fh)‖L2(T)+‖hTh−1/2(u−ℓ−uh)‖H1(T) +‖α(u−uhℓ)−β(v−vhℓ)‖L2(F) +λΓ(F)(hFh−1/2‖∇uh‖L2(Th)+hFh1/2‖uh−fh‖L2(Th)+‖αuh−βvh‖L2(Fh)). Due to the properties of the bubble functions, this directly results in the estimate   ηC =‖hFh∂1/2rFhC‖L2(Fh∂) ≲‖u−uhℓ‖H1(Ω)+osc(f,T)+‖hF∂1/2(α(u−uhℓ)−β(v−vhℓ))‖L2(Γ) +∑Fh∈Fh∂λΓ(F)(‖∇uh‖L2(Th)+hFh‖uh−fh‖L2(Th)+hFh1/2‖αuh−βvh‖L2(Fh)) ≲|||(u−uhℓ,v−vhℓ)|||+osch+hTηncB,0+ηncB,1+hF∂1/2ηncC. (iii) Efficiency of$$\boldsymbol{\eta_S}$$ We introduce for $$F_h\in {\mathcal{F}}_h^\partial$$ and $$E_h\in\mathcal{E}_h^\partial$$ the (surface) bubble functions $$b_{F_h}^\partial$$ with $$\text{supp} b_{F_h}^\partial = F_h$$ and $$b_{E_h}^\partial$$ with $$\text{supp} b_{E_h}^\partial = \omega_{E_h}^\partial$$. For any $$F\in{{\mathcal{F}}_h}^\partial$$ and for any $$E\in{\mathcal{E}_h}^\partial$$, let   rFhS :=ΔΓhvh+αuh−(1+β)vh+gh,rEhS :=[[∂nvh]],ψFh :=bFh∂rFhS,ψEh :=bEh∂rEhS. From the properties of the bubble functions, we deduce   ‖rFhS‖L2(Fh)2 ≈∫FhrFhSψFhdsand‖rEhS‖L2(Eh)2 ≈∫EhrEhSψEhdσ,‖ψFh‖L2(Fh) ≲‖rFhS‖L2(Fh)and‖ψEh‖L2(ωEh∂) ≲hEh1/2‖rEhS‖L2(Eh). As before, we integrate (1.1b) and (1.1c) over $$F$$, and integrate by parts to get   ‖rFhS‖L2(Fh)2 ≲∫Fh∇Γhvh⋅∇ΓhψFh+(αuh−(1+β)vh+gh)ψFhds −∫F∇Γv⋅∇ΓψFhℓ+(αu−(1+β)v+g)ψFhℓds. With Lemma 2.6, the properties of the bubble function $$b_{F_h}^\partial$$ and an inverse inequality, we deduce   ‖rFhS‖L2(Fh) ≲hFh−1(‖vhℓ−v‖H1(Fh) +osc(g,F))+‖α(uhℓ−u)−β(vhℓ−v)‖L2(Fh) +λΓ(F)(hF−1‖∇Γhvh‖L2(Fh)+‖vh−gh‖L2(Fh)+‖αuh−βvh‖L2(Fh)). Moreover, with the properties of $$b_{E_h}^\partial$$ and Lemma 2.6,   ‖rEhS‖L2(Eh)2 ≲∫ωEh∂rFhSψEhds +∫ωEh∂∇Γhvh⋅∇ΓhψEh+(vh−gh)ψEhds +∫ωEh∂(βvh−αuh)ψEhds ≲∫ωEh∂rFhSψEhds +∫ωE∂∇Γ(vhℓ−v)⋅∇ΓψEhℓ+(vhℓ−v)ψEhℓds +∫ωE∂(β(vhℓ−v)−α(uhℓ−u))ψEhℓds +∫ωE∂(g−ghℓ)ψEhℓds +∑Fh∈ωEh∂λΓ(Fh)(‖∇Γhvh‖L2(ωEh∂)‖∇ΓhψEh‖L2(ωEh∂) +(‖vh−gh‖L2(ωEh∂)+‖αuh−βvh‖L2(ωEh∂))‖ψEh‖L2(ωEh∂)). It immediately follows   ‖rEhS‖L2(Eh) ≲hF1/2‖rFhS‖L2(ωEh∂)+hF−1/2‖vhℓ−v‖H1(ωE∂)+hF−1/2osc(g,ωEh∂) +hF1/2‖α(uhℓ−u)−β(vhℓ−vh)‖L2(ωE∂) +∑Fh∈ωEh∂λΓ(Fh)(hFh−1/2‖∇Γhvh‖L2(Fh)+hFh1/2‖vh−gh‖L2(Fh) +hFh1/2‖αuh−βvh‖L2(Fh)). With $$\eta_S = \left\lVert{h_{{\mathcal{F}}_h} {r^S_{F_h}}}\right\rVert_{L^2({\mathcal{F}}_h^\partial)} + \left\lVert{h_{{\mathcal{F}}_h}^{1/2}\, {r^S_{E_h}}}\right\rVert_{L^2(\mathcal{E}_h^\partial)}$$ and the finite overlap of patch elements,   ηS ≲‖v−vhℓ‖H1(Γ)+‖hF(α(uhℓ−u)−β(vhℓ−vh))‖L2(Γ)+osc(g,Fh∂) +∑Fh∈Fh∂λΓ(Fh)(‖∇Γhvh‖L2(Fh)+hFh‖βvh−gh‖L2(Fh)+hFh‖αuh−βvh‖L2(Fh)). This verifies the claimed statement and completes the proof. □ 5. Numerical experiments We implemented a finite element method based on the concepts outlined by Alberty et al. (1999). The volume grid is described in two arrays: c4n contains the coordinates of the nodes and thereby defined the global node indices. n4e lists for each element $$T_h$$ the global nodes indices of the corners, and thereby defines the local node and face indices. For the surface representation, we introduce the additional array n4f containing for each surface element $$F\in {\mathcal{F}}^\partial$$ the global node indices of its corners. In the implementation, we do not follow the above convention that within an element the boundary face has to be labelled last. Instead, there is an array bdy_face containing for each $$T_h\in\mathcal{T}_h$$ either the value $$0$$ for interior elements or the local face index $$i\in\{ 1,\ldots,(d+1) \}$$ of the boundary face $$F_h \subset T_h$$. To compute the $$P1$$ stiffness matrix in the volume, on each element $$T_h\in\mathcal{T}_h$$ we need to determine the constant gradient vectors $$\nabla \phi_i$$ for $$i=1,\ldots,(d+1)$$. On each surface element $$F_h\in{\mathcal{F}}_h^\partial$$ we have to evaluate $$\nabla \psi_j$$ for $$j=1,\ldots,d$$. Suppose $$F_h \in {\mathcal{F}}_h^\partial$$ and $$T_h\in\mathcal{T}_h$$ such that $$F_h \subset T_h$$. Let $$\sigma\in\{ 1,\ldots,(d+1) \}$$ denote the local face index of the boundary face $$F_h$$. Then, the outer normal to $${\it{\Gamma}}_h$$ on $$F_h$$ is   nFh=−∇ϕσT/|∇ϕσ|. (5.1) Within $$T_h$$, the corner nodes of $$F_h\subset T_h$$ have the local indices $$i\in\{1,\ldots,d+1\} \setminus \sigma$$, and the gradients $$\nabla \psi_j$$, of the surface finite element shape functions are given by   ∇ψj=∇ϕi(j)T−(∇ϕi(j)T⋅nFh)nFh for j=1,…,d. (5.2) In the adaptive algorithm, elements are marked for refinement using the indicator   ηTh:=αηC,Th+αηB,Th+14βηS,Th, (5.3) where the terms on the right-hand side refer to the error contributions of element $$T_h\in \mathcal{T}_h$$ to the corres-ponding estimators defined in (4.8). The elements with the largest values of $$\eta_{T_h}$$ are marked for refinement until the marked elements contribute to three-fourth of the total estimated error. For the mesh refinement, we apply the bisection algorithm of Bartels & Schreier (2012). Note that marking a boundary element $$T_h\in\mathcal{T}_h$$ due to a large surface or coupling indicators on $$F_h\in{\mathcal{F}}_h^\partial$$ or $$F_h\subset T_h$$ might not lead to an immediate refinement of the surface element because the refinement edge of $$T_h$$ might be in the interior of $${\it{\Omega}}_h$$. For uniformly refined grids, the number of degrees of freedom (DoF) in the volume is proportional to $$h_T^{-3}$$ and the number of DoFs on the surface is proportional to $$h_F^{-2}$$. In the following we will call a quantity proportional to $$h_T$$ or $$h_F$$ whenever it shows a respective dependence on the DoFs in the volume or on the surface, respectively. We set $$\alpha=\beta=1$$ in the following experiments. 5.1 Convergence to known regular solution As a first benchmark, we verify our method with the example from (Elliott & Ranner, 2013). On the unit sphere $${\it{\Omega}} =B_{1}(0) \subset\mathbb{R}^3$$ with boundary $${\it{\Gamma}}=\partial{\it{\Omega}}$$, the data $$f$$ and $$g$$ are prescribed such that the exact solution of the problem (3.3) is   u(x,y,z) =exp⁡(−x(x−1)−y(y−1)), (5.4a)  v(x,y,z) =exp⁡(−x(x−1)−y(y−1))[1+x(1−2x)+y(1−2y)]. (5.4b) The computed solution $$(u_h,v_h)$$ is plotted in Fig. 3. From Fig. 4, we see that both the uniform and the adaptive refinement yield the same asymptotic behaviour with respect to the number of DoFs. On uniformly refined grids, the $$H^1({\it{\Omega}})$$-error of $$u_h$$ decreases proportionally to $$h_T$$. Likewise, the $$H^1({\it{\Gamma}})$$-error of $$v_h$$ is proportional to $$h_F$$. Examining Fig. 5 confirms that the estimator $$\eta_h$$ and its contributions $$\eta_B$$ and $$\eta_S$$ related to the volume and surface correctly follow the reduction of the overall error and the individual $$H^1$$ errors of $$u_h$$ and $$v_h$$, respectively. Moreover, we observe that during the first refinement steps, the estimator $$\eta_S$$ of the surface error is dominant, leading to some extra refinement of $${\it{\Gamma}}_h$$ in the adaptive algorithm, cf. Fig. 3. Due to the high regularity of the prescribed solution, the later refinement steps are almost uniform. The estimate $$\eta_C$$ of the coupling error is at least about one order of magnitude smaller than the other error contributions, and thus does not have any influence on the refinement. Additional evaluation of the $$L^2$$ errors revealed that uniform refinement leads to $$L^2$$ error reduction with optimal order $$O(h_T^2)$$ in the volume and $$O(h_F^2)$$ on the surface. When refining the mesh according to the $$H^1$$ estimator $$\eta_{T_h}$$, we still observe the $$L^2$$ error reduction proportional to the $$h_T^2$$ in the volume, but slightly lower than $$h_F^2$$ on the surface. Fig. 3. View largeDownload slide Experiment of Section 5.1: cross-section of the solution $$u_h$$ and the adaptively refined grid (left) and solution $$v_h$$ on the surface (right). Fig. 3. View largeDownload slide Experiment of Section 5.1: cross-section of the solution $$u_h$$ and the adaptively refined grid (left) and solution $$v_h$$ on the surface (right). Fig. 4. View largeDownload slide Experiment of Section 5.1: bulk, surface and overall error with respect to the exact solution in the $$H^1$$-norm over the number of unknowns for uniform refinement (left) and adaptive refinement (right). Fig. 4. View largeDownload slide Experiment of Section 5.1: bulk, surface and overall error with respect to the exact solution in the $$H^1$$-norm over the number of unknowns for uniform refinement (left) and adaptive refinement (right). Fig. 5. View largeDownload slide Experiment of Section 5.1: discretization error estimator $$\eta_h$$ with contributions $$\eta_B, \eta_S,\eta_C$$ over the number of unknowns with uniform refinement (left) and adaptive refinement (right). Fig. 5. View largeDownload slide Experiment of Section 5.1: discretization error estimator $$\eta_h$$ with contributions $$\eta_B, \eta_S,\eta_C$$ over the number of unknowns with uniform refinement (left) and adaptive refinement (right). 5.2 Nonsmooth boundary and grid adaption due to data In the next experiment with $${\it{\Omega}} = [-1,1]^3$$, we want to examine the effect of a domain $${\it{\Omega}}$$ with nonsmooth boundary $${\it{\Gamma}}$$ and the ability of the adaptive algorithm to resolve singularities of the solution. Then, the surface $${\it{\Gamma}}$$ is only piecewise smooth. Since all faces are flat (i.e., affine), the discretization of the surface $${\it{\Gamma}}_h = {\it{\Gamma}}$$ is exact, and hence $$\lambda_{\it{\Gamma}} = 0$$. We set   xf :=(0.7,0.6,0.5),f(x) :=(|x−xf|2+10−4)−1, (5.5a)  xg :=(0.1,−1.0,−0.3),g(x) :=(|x−xg|3/2+10−5)−1. (5.5b) By this, we introduce ‘weak singularities’ in $${\it{\Omega}}$$ and on $${\it{\Gamma}}$$. The computed solution is plotted in Figure 6 with a cross-section of the bulk solution $$u_h$$ on the left-hand side. Figure 7 shows that adaptive refinement leads to reduction rates of $$\eta_B$$ and $$\eta_S$$ proportional to $$O(h_T)$$ and $$O(h_F)$$, respectively, whereas for uniform refinement the reduction rates are lower and do not reach the same asymptotic behaviour. In particular, we see in Fig. 7, that using uniform refinement the reduction of the surface estimate $$\eta_S$$ is considerably slower than $$O(h_F)$$, while for $$\eta_B$$ the difference to the $$O(h_T)$$ rate is less pronounced. Keeping in mind that the volume grid has more DoFs than the according surface grid, we observe that in each uniform refinement step, the indicator $$\eta_S$$ on the surface is larger than $$\eta_B$$ in the volume. As before, the coupling error $$\eta_C$$ is only marginal with regard to $$\eta$$. Moreover, we see from Fig. 6 that the adaptive algorithm has refined the grid towards the weak singularities of the data. These singularities are well separated from the nonsmooth edges of $${\it{\Gamma}}$$ where no extra refinement of the grid can be observed. Fig. 6. View largeDownload slide Experiment of Section 5.2: Solution for data specified by (5.5): Cross-section of $$u_h$$ with the locally refined grid (left) and $$v_h$$ (right). Fig. 6. View largeDownload slide Experiment of Section 5.2: Solution for data specified by (5.5): Cross-section of $$u_h$$ with the locally refined grid (left) and $$v_h$$ (right). Fig. 7. View largeDownload slide Experiment of Section 5.2: Uniform refinement leads to a reduction rate of the estimated error $$\eta_h$$ slower than $$O(h_T)$$ in the volume and $$O(h_F)$$ on the surface (left), whereas for the adaptively refined grids $$\eta_h$$ meets the optimal rates (right). Fig. 7. View largeDownload slide Experiment of Section 5.2: Uniform refinement leads to a reduction rate of the estimated error $$\eta_h$$ slower than $$O(h_T)$$ in the volume and $$O(h_F)$$ on the surface (left), whereas for the adaptively refined grids $$\eta_h$$ meets the optimal rates (right). 5.3 Domain with corner singularities In the final experiment we examine the effect of corner singularities. Note that the classical example of a re-entrant corner and $$u(r,\phi) = r^{\alpha} \sin(\alpha \phi)$$ cannot be transferred to the coupled problem (1.1). Let the domain $${\it{\Omega}}$$ be given by a ball of radius 1 where one octant is removed (see Fig. 8). We expect $$u$$ and $$v$$ to be of low regularity at the reentrant corner points of the domain $${\it{\Omega}}$$ and the surface $${\it{\Gamma}}$$, respectively. We set $$x_1 = \frac{1}{2} (-1, 1, 1)$$, $$x_2 = \frac{1}{2} ( 1,-1, 1)$$, $$x_3 = \frac{1}{2} ( 1, 1,-1)$$, and choose   f(x) =0, (5.6a)  g(x) =20−20(|x|−exp⁡(−4|x−x1|2)−exp⁡(−4|x−x2|2)+exp⁡(−4|x−x3|2)). (5.6b) Fig. 8. View largeDownload slide Experiment of Section 5.3: The domain $${\it{\Omega}}$$ with the surface patches $${\it{\Gamma}}^i$$, $$i=1,2,3,4$$ (left). Refinement of the grid towards the corner at the north pole $$(0,0,1)$$ due to the adaptive algorithm (right). Fig. 8. View largeDownload slide Experiment of Section 5.3: The domain $${\it{\Omega}}$$ with the surface patches $${\it{\Gamma}}^i$$, $$i=1,2,3,4$$ (left). Refinement of the grid towards the corner at the north pole $$(0,0,1)$$ due to the adaptive algorithm (right). The computed solution $$(u_h,v_h)$$ is plotted in Fig. 9. Fig. 9. View largeDownload slide Experiment of Section 5.3: Solution $$u_h$$ (left) and $$v_h$$ (right). Fig. 9. View largeDownload slide Experiment of Section 5.3: Solution $$u_h$$ (left) and $$v_h$$ (right). Similar to the experiment of Section 5.2, we conclude from Fig. 10 that adaptive refinement leads to reduction rates of $$\eta_B$$ and $$\eta_S$$ proportional to $$O(h_T)$$ and $$O(h_F)$$, respectively, whereas for uniform refinement the reduction rates are lower and do not reach the same asymptotic behaviour. However, when compared to Fig. 7, the deviation of $$\eta_S$$ from the $$O(h_F)$$ line is considerably less pronounced in this experiment. This indicates that for the coupled system the geometrically induced singularities are weaker than those known in the pure volume or the pure surface problem. At least, it is not possible to trigger stronger singularities by using data (5.6). Fig. 10. View largeDownload slide Experiment of Section 5.3: Uniform refinement leads to a reduction rate of the estimated error $$\eta_h$$ and the volume and surface contributions $$\eta_B$$ and $$\eta_S$$ slower than $$O(h_T)$$ and $$O(h_F)$$, respectively (left). For the adaptively refined grids $$\eta_h$$, the optimal rates are attained (right). Fig. 10. View largeDownload slide Experiment of Section 5.3: Uniform refinement leads to a reduction rate of the estimated error $$\eta_h$$ and the volume and surface contributions $$\eta_B$$ and $$\eta_S$$ slower than $$O(h_T)$$ and $$O(h_F)$$, respectively (left). For the adaptively refined grids $$\eta_h$$, the optimal rates are attained (right). References Ainsworth, M. & Oden, J. ( 2000) A Posteriori Error Estimation in Finite Element Analysis . Pure and Applied Mathematics (New York). New York: Wiley-Interscience [John Wiley & Sons], pp. xx+240. Google Scholar CrossRef Search ADS   Alberty, J., Carstensen, C. & Funken, S. ( 1999) Remarks around 50 lines of matlab: short finite element implementation. Numer. Algorithms , 20, 117– 137. Google Scholar CrossRef Search ADS   Babuška, I. & Rheinboldt, W. C. ( 1978) A-posteriori error estimates for the finite element method. Int. J. Numer. Meth. Engng. , 12, 1597– 1615. Google Scholar CrossRef Search ADS   Bartels, S. & Schreier, P. ( 2012) Local coarsening of simplicial finite element meshes generated by bisections. BIT , 52, 559– 569. Google Scholar CrossRef Search ADS   Bernardi, C. ( 1989) Optimal finite-element interpolation on curved domains. SIAM J. Numer. Anal. , 26, 1212– 1240. Google Scholar CrossRef Search ADS   Bertalmío, M., Cheng, L.-T., Osher, S. & Sapiro, G. ( 2001) Variational problems and partial differential equations on implicit surfaces. J. Comput. Phys. , 174, 759– 780. Google Scholar CrossRef Search ADS   Bonito, A., Cascón, J. M., Morin, P. & Nochetto, R. H. ( 2013) AFEM for geometric PDE: the Laplace-Beltrami operator. Analysis and Numerics of Partial Differential Equations  ( Brezzi, F. Franzone, P. C. Gianazza U. & Gilardi G. eds). Milan, Italy: Springer, pp. 257– 306. Google Scholar CrossRef Search ADS   Bonito, A., Cascón, J. M., Morin, P., Mekchay, K. & Nochetto, R. H. ( 2016) High-order AFEM for the Laplace-Beltrami operator: Convergence rates. Found. Comput. Math. , 16, 1473– 1539. Google Scholar CrossRef Search ADS   Brenner, S. & Carstensen, C. ( 2004) Finite element methods. Encyclopedia of computational mechanics  ( Stein, E. De Borst R. & Hughes T. J. R. eds), Chichester: John Wiley & Sons, 73– 117. Burger, M. ( 2009) Finite element approximation of elliptic partial differential equations on implicit surfaces. Comput. Vis. Sci. , 12, 87– 100. Google Scholar CrossRef Search ADS   Burman, E., Claus, S., Hansbo, P., Larson, M. & Massing, A. ( 2015a) CutFEM: Discretizing geometry and partial differential equations. Int. J. Numer. Meth. Engng. , 104, 472– 501. Google Scholar CrossRef Search ADS   Burman, E., Hansbo, P., Larson, M. & Zahedi, S. ( 2015b) Cut finite element methods for coupled bulk–surface problems. Numer. Math. , 133, 203– 231. Google Scholar CrossRef Search ADS   Camacho, F. & Demlow, A. ( 2015) $$L_2$$ and pointwise a posteriori error estimates for FEM for elliptic PDEs on surfaces. IMA J. Num. Anal. , 35, 1199– 1227. Google Scholar CrossRef Search ADS   Carstensen, C., Eigel, M., Hoppe, R. & Löbhard, C. ( 2012) A review of unified a posteriori finite element error control. Numer. Math. Theor. Meth. Appl. , 5, 509– 558. Google Scholar CrossRef Search ADS   Chernyshenko, A. Y. & Olshanskii, M. A. ( 2015) An adaptive octree finite element method for PDEs posed on surfaces. Comput. Methods Appl. Mech. Eng. , 291, 146– 172. Google Scholar CrossRef Search ADS   Ciarlet, P. & Raviart, P.-A. ( 1972) Interpolation theory over curved elements, with applications to finite element methods. Comput. Methods Appl. Mech. Eng. , 1, 217– 249. Google Scholar CrossRef Search ADS   Dassi, F. ( 2014) Finite element techniques for curved boundaries. Ph.D. Thesis , Milan, Italy: Politecnico di Milano. Dassi, F., Perotto, S. & Formaggia, L. ( 2015) A priori anisotropic mesh adaptation on implicitly defined surfaces. SIAM J. Sci. Comp. , 37, A2758– A2782. Google Scholar CrossRef Search ADS   Deckelnick, K., Dziuk, G., Elliott, C. & Heine, C.-J. ( 2010) An h-narrow band finite-element method for elliptic equations on implicit surfaces. IMA J. Num. Anal. , 30, 351– 376. Google Scholar CrossRef Search ADS   Dedner, A. & Madhavan, P. ( 2016) Adaptive discontinuous galerkin methods on surfaces. Numer. Math. , 132, 369– 398. Google Scholar CrossRef Search ADS   Demlow, A. ( 2009) Higher-order finite element methods and pointwise error estimates for elliptic problems on surfaces. SIAM J. Numer. Anal. , 47, 805– 827. Google Scholar CrossRef Search ADS   Demlow, A. & Dziuk, G. ( 2007) An adaptive finite element method for the Laplace-Beltrami operator on implicitly defined surfaces. SIAM J. Numer. Anal. , 45, 421– 442. Google Scholar CrossRef Search ADS   Dörfler, W. & Rumpf, M. ( 1998) An adaptive strategy for elliptic problems including a posteriori controlled boundary approximation. Math. Comp. , 67, 1361– 1382. Google Scholar CrossRef Search ADS   Dziuk, G. ( 1988) Finite elements for the Beltrami operator on arbitrary surfaces. Partial Differential Equations and Calculus of Variations  ( Hildebrand S. & Leis R. eds). Lecture Notes in Mathematics. Heidelberg, Berlin: Springer, pp. 142– 155. Google Scholar CrossRef Search ADS   Dziuk, G. & Elliott, C. ( 2013) Finite element methods for surface PDEs. Acta Numer. , 22, 289– 396. Google Scholar CrossRef Search ADS   Egger, H., Fellner, K., Pietschmann, J.-F. & Tang, B. ( 2015) A finite element method for volume-surface reaction-diffusion systems. Preprint , IGDK-2015-26. Elliott, C. & Ranner, T. ( 2013) Finite element analysis for a coupled bulk-surface partial differential equation. IMA J. Num. Anal. , 33, 377– 402. Google Scholar CrossRef Search ADS   Giese, W., Eigel, M., Westerheide, S., Engwer, C. & Klipp, E. ( 2015) Influence of cell shape, inhomogeneities and diffusion barriers in cell polarization models. Phys. Biol. , 12, 066014. Google Scholar CrossRef Search ADS PubMed  Gross, S., Olshanskii, M. A. & Reusken, A. ( 2015) A trace finite element method for a class of coupled bulk-interface transport problems. ESAIM: M2AN , 49, 1303– 1330. Google Scholar CrossRef Search ADS   Novak, I. L., Gao, F., Choi, Y.-S., Resasco, D., Schaff, J. C. & Slepchenko, B. M. ( 2007) Diffusion on a curved surface coupled to diffusion in the volume: Application to cell biology. JCP , 226, 1271– 1290. Google Scholar CrossRef Search ADS   Olshanskii, M., Reusken, A. & Grande, J. ( 2009) A finite element method for elliptic equations on surfaces. SIAM J. Numer. Anal. , 47, 3339– 3358. Google Scholar CrossRef Search ADS   Rätz, A. & Röger, M. ( 2012) Turing instabilities in a mathematical model for signaling networks. J. Math. Biol. , 65, 1215– 1244. Google Scholar CrossRef Search ADS PubMed  Scott, L. R. ( 1973) Finite element techniques for curved boundaries. Ph.D. Thesis , Massachusetts, Cambridge: Massachusetts Institute of Technology. Scott, L. R. & Zhang, S. ( 1990) Finite element interpolation of nonsmooth functions satisfying boundary conditions. Math. Comp. , 54, 483– 493. Google Scholar CrossRef Search ADS   Verfürth, R. ( 1996) A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques . Chichester, New York, Brisbane, Toronto, Singapore, Stuttgart, Leipzig: Wiley-Teubner. Verfürth, R. ( 2013) A Posteriori Error Estimation Techniques for Finite Element Methods . Oxford, UK: Oxford University Press. Google Scholar CrossRef Search ADS   Zlámal, M. ( 1973) Curved elements in the finite element method. I. SIAM J. Numer. Anal. , 10, 229– 240. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# A posteriori error control for stationary coupled bulk-surface equations

, Volume 38 (1) – Jan 1, 2018
28 pages

/lp/ou_press/a-posteriori-error-control-for-stationary-coupled-bulk-surface-XwgjSyil6a
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drw080
Publisher site
See Article on Publisher Site

### Abstract

Abstract We consider a system of two coupled elliptic equations, one defined on a bulk domain and the other one on the boundary surface. Problems of this kind are relevant for applications in engineering, chemistry and in biology, e.g., biological signal transduction. For the a posteriori error control of the coupled system, a residual error estimator is derived which takes into account the approximation errors due to the finite element discretization in space as well as the polyhedral approximation of the surface. An adaptive refinement algorithm controls the overall error. Numerical experiments illustrate the performance of the a posteriori error estimator and the proposed adaptive algorithm with several benchmark examples. 1. Introduction Coupled reaction diffusion processes in the bulk and on the surface of some domain $${\it{\Omega}}\subset\mathbb{R}^d,\,d\in\{2,3\}$$, have recently attracted interest from an analytical point of view (see, e.g., Elliott & Ranner, 2013; Burman et al., 2015b; Egger et al., 2015) in different application areas such as biology (see, e.g., Novak et al., 2007; Rätz & Röger, 2012; Giese et al., 2015). In all these problems one has to simultaneously account for transport in normal as well as in tangential direction to the boundary surface, whereas most often either only normal or only tangential phenomena are considered. We study a stationary prototype problem in domain $${\it{\Omega}}$$ with piecewise smooth boundary $${\it{\Gamma}}:=\partial{\it{\Omega}}$$, which can be decomposed into a finite set of $$N_{\it{\Gamma}}$$ smooth surface patches $$\{{\it{\Gamma}}^i\}_{i=1}^{N_{\it{\Gamma}}}$$. We seek the solution $$u:{\it{\Omega}}\to\mathbb{R}$$ and $$v:{\it{\Gamma}}\to\mathbb{R}$$ of the stationary coupled diffusion-reaction problem   −Δu+u =f in Ω, (1.1a)  (αu−βv)+∂nu =0 on ∪i=1NΓΓi, (1.1b)  −ΔΓv+v+∂nu =g on ∪i=1NΓΓi, (1.1c)  ∇Γv∣Γi⋅nΓi+∇Γv∣Γj⋅nΓj =0 on ∂Γi∩∂Γj. (1.1d) For each point on a surface patch $${\it{\Gamma}}^i$$, $$n$$ is the outer normal of $${\it{\Omega}}$$, and $$\nabla_{{\it{\Gamma}}^i}$$ and $${\it{\Delta}}_{\it{\Gamma}}$$ denote the tangential gradient and the Laplace–Beltrami operator, respectively. In the compatibility condition (1.1d), $$n_{\it{\Gamma}}^i$$ denotes the outer normal of the patch $${\it{\Gamma}}^i$$, which is often called the co-normal to $${\it{\Gamma}}^i$$. For each point in $$\partial{\it{\Gamma}}^i$$, it lies in the tangent plane to $${\it{\Gamma}}^i$$. The equation system (1.1) describes diffusive transport and reaction of a bulk species $$u$$ and a surface species $$v$$, which are coupled by some binding–unbinding process that transforms $$u$$ to $$v$$ by binding it to $${\it{\Gamma}}$$ and vice versa. For the sake of a convenient presentation, the diffusion coefficients and reaction rates are assumed to be $$1$$, while the binding and unbinding rates are given by the positive constants $$\alpha$$ and $$\beta$$. The existence and uniqueness of the solution to (1.1) have been proved by Elliott & Ranner (2013) for globally $$C^2$$ boundaries. Moreover, optimal a priori error estimates for finite element approximations of arbitrary order were derived for sufficiently smooth boundaries and data. For the time-dependent version of (1.1), convergence to equilibrium and a priori error estimates were shown by Egger et al. (2015). In Burman et al. (2015b), a CutFEM method that uses the same volume grid to discretize both the volume and the surface equation is analysed. We complement the analysis of Elliott & Ranner (2013) by the derivation of an a posteriori error estimator for the error $$e:=(u-{u_h},v-{v_h})$$ of a lowest order continuous FEM approximation $$({u_h},{v_h})$$ of the solution $$(u,v)$$ of (1.1) measured in the energy norm. For this, we split the overall residual into an equilibration, a consistency and a data approximation residual. We transfer the classical residual error estimators for volume domains and for surfaces to the setting of the coupled system. While the a posteriori analysis for problems on polyhedral domains is rather mature (see, e.g., Babuška & Rheinboldt, 1978; Verfürth, 1996; Ainsworth & Oden, 2000), and a priori analysis for finite elements on curved domains is well established (see, e.g., Scott, 1973; Zláamal, 1973; Bernardi, 1989), much less can be found for the a posteriori analysis on curved domains. Most notably, Dörfler & Rumpf (1998) put a focus on a very coarse domain approximation, where large parts of the domain may not be discretized at all and the refinement algorithm has to explore and detect the necessary parts of the domain which are important for the discretization. This is somehow opposite to the setting we are concerned with since, due to the coupling to the surface equation, we expect the complete domain to be required for an accurate solution. For the solution of PDE problems on surfaces, various numerical approaches can be used, cf. the recent review of Dziuk & Elliott (2013). As a first approach, finite element methods on (surface-)meshes that approximate the surface have been devised. For smooth surfaces, they are analysed by Dziuk (1988), and Demlow & Dziuk (2007) derived an a posteriori error estimate which splits the error into contributions from the discretization, geometric errors and data approximations. For an IP-dG scheme, an error estimator of similar structure for a reaction diffusion equation on a smooth surface was analysed by Dedner & Madhavan (2016), in which case the geometric part was shown to be of a higher order. Dassi et al. (2015) transferred ideas of Demlow & Dziuk (2007) into the context of anisotropic meshes and an $$L^1$$ interpolation error estimate is derived together with heuristic control of the geometric error. A posteriori estimates with ZZ-type reconstruction of the solution are derived by Dassi (2014), where again the extension to anisotropic meshes is considered. Higher order methods for pure surface problems were analysed by Demlow (2009). $$L^2$$ and pointwise estimates for $${\it{\Gamma}}\in C^3$$ are given by Camacho & Demlow (2015). For less smooth surfaces, Bonito et al. (2013) derived a convergent refinement algorithm for parametric finite elements, and quasi-optimality of the resulting discretizations is shown under certain assumptions. (Bonito et al., 2016) contains the extension to higher order schemes. Our error estimator extends some of the results by Bonito et al. (2013) for the coupled bulk-surface equations. One main difficulty and contribution of this paper is the analysis of the geometry error of nonconformity caused by the polyhedral approximation of the surface manifold. As an alternative approach, unfitted volume meshes can be used to solve PDE problems on lower-dimensional surfaces. On the one hand, one can extend the equations from the surface to the volume (cf. Bertalmío et al., 2001; Burger, 2009). This approach is in particular well suited for surfaces that are implicitly defined as the level set of a smooth function. Additional computational costs for solving in higher space dimension can be reduced by narrow band techniques as shown by Deckelnick et al. (2010). On the other hand, in TraceFEM or CutFEM the trial functions are defined as traces of functions from spaces related to the volume mesh (see, e.g., Olshanskii et al., 2009; Chernyshenko & Olshanskii, 2015; Burman et al., 2015a). For coupled volume-surface problems with smooth surfaces, TraceFEM and CutFEM are applied and analysed by Gross et al. (2015) and Burman et al. (2015b), respectively. The outline of the paper is as follows: In Section 2, we introduce the description of curved elements and recall some basic results from differential geometry which include integral transformations used later on. In particular, nonconformity estimates are derived which provide bounds for the error caused by the domain approximations. Then, the weak formulation and FEM discretization of the problem (1.1) is presented in Section 3. Section 4 is dedicated to the derivation of the residual a posteriori error estimator for the overall error which is the main result of this paper. It consists of a discretization part and a geometric nonconformity part, both of which can be controlled a posteriori by computable indicators. Numerical experiments in Section 5 illustrate the performance of the derived a posteriori error estimator with several benchmark problems. In order to simplify the notation, we write $$a\lesssim b$$ if $$a\leq C\,b$$ with some mesh-independent positive constant $$C$$. Moreover, the common notation for Lebesgue and Sobolev spaces is employed. 2. Domain approximation A finite element method for (1.1) involves an approximation of the curved domain $${\it{\Omega}}$$ by some interpolating polyhedral domain $${\it{\Omega}}_h$$. If $${\it{\Omega}}$$ is a fixed piecewise smoothly bounded domain, all nonsmooth edges or points on the boundary are known a priori, and no degeneration of the smooth patches $${\it{\Gamma}}^i \subseteq \partial{\it{\Omega}}$$ can happen, in contrast to problems on evolving domains. It is thus reasonable to define $${\it{\Omega}}$$ as the deformation of some polyhedral domain $${\it{\Omega}}_h^0$$ with an associated triangulation $$\mathcal{T}_h^0$$ into simplices. We refer to $$\mathcal{T}_h^0$$ as the macro triangulation and denote the deformation by $${\it{\Phi}}^0$$. We assume that $${\it{\Phi}}^0$$ is globally continuous, and the restriction to each boundary face or boundary edge is a $$C^1$$ function. Moreover, the generality of $${\it{\Phi}}^0$$ is limited by the requirement that there are no self intersections of $${\it{\Gamma}}$$ and that $${\it{\Gamma}}$$ does not intersect the inner facets of $$\mathcal{T}_h^0$$, as detailed below. In the following, we consider different triangulations $$\mathcal{T}_h$$ with their respective polyhedral domain $${\it{\Omega}}_h$$ and functions $${\it{\Phi}}$$ which describe the deformation of $$\partial {\it{\Omega}}_h$$ to match $$\partial {\it{\Omega}}$$. We always assume that $$\mathcal{T}_h$$ has been obtained from $$\mathcal{T}_h^0$$ by several steps of refinement by bisection without indicating the refinement level. Each refinement step is followed by a deformation such that the respective polygonal domain $$\partial {\it{\Omega}}_h$$ interpolates $$\partial {\it{\Omega}}$$. 2.1 Simplicial triangulations Let $${\it{\Omega}}_h\subset\mathbb{R}^{d}$$ be a polyhedral domain with the boundary $${{\it{\Gamma}}_{\!h}} = \partial{\it{\Omega}}_h$$, which interpolates $${\it{\Gamma}}$$. Suppose $$\mathcal{T}_h$$ is a regular partition of $${\it{\Omega}}_h$$ into simplices. For $$d=3$$, we denote the sets of nodes, edges and faces of the triangulation $$\mathcal{T}_h$$ by $$\mathcal{N}_h$$, $$\mathcal{E}_h$$ and $${\mathcal{F}}_h$$, respectively. If $$d=2$$ we let $${\mathcal{F}}_h$$ denote the set of edges and let $$\mathcal{E}_h$$ coincide with $$\mathcal{N}_h$$. For each $$T_h\in\mathcal{T}_h$$, $$F_h\in{\mathcal{F}}_h$$, $$E_h\in\mathcal{E}_h$$, we set $$h_{T} = \operatorname{diam}(T_h)$$, $$h_{F} = \operatorname{diam}(F_h)$$, $$h_{E} = \operatorname{diam}(E_h)$$ and define by this the (global) mesh-size functions $$h_{\mathcal{T}}$$, $$h_{{\mathcal{F}}}$$, $$h_{\mathcal{E}}$$. The jump of some (piecewise) function $$u\in L^2({\it{\Omega}}_h;\mathbb{R}^d)$$ over a face $$F_h\in{\mathcal{F}}_h$$ is defined by $$[u]_{F_h}:= u\cdot n_+ + u\cdot n_-$$ where $$n_+, n_-$$ denote the outer normals on $$F_h=T_+\cap T_-$$ with respect to the elements $$T_+,T_-\in\mathcal{T}_h$$. The subset $${\mathcal{F}}_h^\partial := \{ F_h \in {\mathcal{F}}_h : \ F_h \subset {{\it{\Gamma}}_{\!h}} \}$$ defines a triangulation of $${{\it{\Gamma}}_{\!h}}$$ into $$(d-1)$$ dimensional simplices. Analogous to the volume, we introduce the sets of boundary edges and boundary nodes, denoted by $$\mathcal{E}_h^\partial$$ and $$\mathcal{N}_h^\partial$$. For $$E_h\in\mathcal{E}_h^\partial$$ and $$F_+,F_-\in{\mathcal{F}}_h^\partial$$ such that $$E_h = F_+\cap F_-$$, the jump of some (piecewise) function $$v\in L^2(\partial{\it{\Omega}}_h;\mathbb{R}^d)$$ over the boundary edge $$E_h$$ is defined by $$[[v]]_{E_h} := v\cdot n_{{\it{\Gamma}}_{\!h}}^+ + v\cdot n_{{\it{\Gamma}}_{\!h}}^-$$, where $$n_{{\it{\Gamma}}_{\!h}}^+, n_{{\it{\Gamma}}_{\!h}}^-$$ denote the outer co-normals on $$E_h$$ with respect to the boundary facets $$F_+,F_-\in{\mathcal{F}}_h^\partial$$. Note that $$n_{{\it{\Gamma}}_{\!h}}^+, n_{{\it{\Gamma}}_{\!h}}^-$$ are not necessarily aligned. For the analysis below, we need the following assumptions on the triangulation: (A1) $$\mathcal{T}_h$$ interpolates $${\it{\Gamma}}$$, i.e., the set of boundary nodes satisfies $$\mathcal{N}_h^\partial \subset {\it{\Gamma}}$$. (A2) Each element $$T_h \in \mathcal{T}_h$$ has at most one boundary face $$F_h \subseteq \partial{\it{\Omega}}_h \cap \partial T_h$$. If $$T_h$$ has no boundary face, then there is at most one boundary edge $$E_h \subseteq \partial{\it{\Omega}}_h \cap \partial T_h$$. (A3) All refinements of $$\mathcal{T}_h$$ by bisection and subsequent deformation of $$\partial{\it{\Omega}}_h$$ to interpolate $$\partial {\it{\Omega}}$$ lead to a shape regular family of triangulations. On fixed polyhedral volume domains $${\it{\Omega}}$$, it is known that appropriate bisection algorithms guarantee the shape regularity, and thus assumption (A3) is not necessary in this case. For the case of finite elements on curved surfaces, we refer to the discussion by Demlow & Dziuk (2007) and Bonito et al. (2013). Let $$\hat T$$ denote the $$d$$-dimensional unit simplex. Given an element $$T_h\in\mathcal{T}_h$$, there is an affine reference transformation which we refer to as $$J:\hat{T} \to T_h$$, without any further index related to $$T_h$$. If $$T_h$$ has one boundary face $$F_h\subset{\mathcal{F}}_h^\partial$$ then we assume that $$\hat{F} := \{ \hat x \in \hat T: \hat x_d = 0 \}$$ is mapped by $$J$$ to the boundary segment $$F_h\subset {{\it{\Gamma}}_{\!h}}$$ as depicted in Fig. 1. Analogously, if $$T_h\cap{\it{\Gamma}}_h$$ consists only of one boundary edge $$E_h$$, we assume that $$E_h = J(\hat{E})$$, where $$\hat{E} := \{ \hat x \in \hat T: \hat x_d =\hat x_{d-1} = 0 \}$$. Due to the shape regularity assumption (A3), there is a constant $$L_h>0$$ such that   hLh−1|z|≤‖∇^Jz‖L∞(T^)≤hLh|z| for all z∈Rd. (2.1) Fig. 1. View largeDownload slide Transformation $$\mathcal{X}$$ and $$J$$ from reference element $$\hat T$$ to a curved element $$T$$ and to its polyhedral approximation $$T_h$$, respectively. Fig. 1. View largeDownload slide Transformation $$\mathcal{X}$$ and $$J$$ from reference element $$\hat T$$ to a curved element $$T$$ and to its polyhedral approximation $$T_h$$, respectively. 2.2 Exact triangulation of curved domains Let $${\it{\Omega}}_h$$ and $$\mathcal{T}_h$$ be as in Section 2.1 satisfying (A1)-(A3). The assumed refinement process from $$\mathcal{T}_h^0$$ to $$\mathcal{T}_h$$ also yields an according boundary deformation $${\it{\Phi}}$$ from $${\it{\Phi}}^0$$, which moves $${\it{\Gamma}}_h$$ to $${\it{\Gamma}}$$. (A1) implies that $${\it{\Phi}} = 0$$ in all boundary nodes. We assume that the refinement process is such that (A4) $${\it{\Phi}}: {{\it{\Gamma}}_{\!h}} \to \mathbb{R}^d$$ is globally continuous and the restriction to each boundary face $$F_h \in{\mathcal{F}}_h^\partial$$ or boundary edge $$E_h\in\mathcal{E}_h^\partial$$ is a $$C^1$$ function. At first, we restrict our considerations to a single boundary element $$T_h \in \mathcal{T}_h$$ such that $$T_h\cap{{\it{\Gamma}}_{\!h}} = F_h \in{\mathcal{F}}_h^\partial$$ or $$T_h\cap{{\it{\Gamma}}_{\!h}} = E_h\in\mathcal{E}_h^\partial$$. We denote the restriction of the boundary deformation to $$F_h$$ or $$E_h$$ again by $${\it{\Phi}}$$. With the pull-back to the reference element $$\hat{T}$$, we can define an extension to $$\hat{T}$$ as   Φ^T(x^) ={Φ(J(x^1,…,x^d−1,0)ρ(0,x^d))if Th∩Γh=Fh,Φ(J(x^1,…,x^d−2,0,0)ρ(x^d−1,x^d))if Th∩Γh=Eh,  (2.2) with some monotonous $$C^1$$-function $$\rho$$ which interpolates between $$0$$ and $$1$$ such that $$\rho(0,0)=1$$. The simplest choice is $$\rho(s,t) = 1 -s -t$$. Lemma 2.1 Let $$T_h\in\mathcal{T}_h$$ be a boundary element and $${\it{\Phi}} \in C^1(T_h\cap{\it{\Gamma}}_h, \mathbb{R}^d)$$ be a boundary deformation that vanishes in all boundary nodes. Let $$\hat{{\it{\Phi}}}_T$$ be defined according to (2.2) and set $$1 \leq L_\rho := \| \hat{\nabla} \rho \|_{L^\infty}$$. Suppose there is $$L_T>0$$ such that (a)$$\| \nabla_{\!F_h} {\it{\Phi}} \|_{L^\infty(F_h)} \leq (1+L_\rho)^{-1} L_T$$, if $$T_h\cap{{\it{\Gamma}}_{\!h}} = F_h\in{\mathcal{F}}_h^\partial$$, or (b)$$\| \nabla_{\!E_h} {\it{\Phi}} \|_{L^\infty(E_h)} \leq (1+L_\rho)^{-1} L_T$$, if $$T_h\cap{{\it{\Gamma}}_{\!h}} = E_h\in\mathcal{E}_h^\partial$$. In either case it holds   ‖∇^Φ^T‖L∞(T^) ≤LT‖∇^J‖L∞(T^). (2.3) Proof. We assume $$T_h \cap {{\it{\Gamma}}_{\!h}} = F_h$$. The other case follows from analogous arguments. From the construction of $$\hat{{\it{\Phi}}}_T$$ in (2.2) and abbreviating $$\hat{x}_{\hat{F}} = ( \hat{x}_1, \dots, \hat{x}_{d-1} , 0 ) \in \hat{F}$$, we check   ∇^Φ^T(x^) =( ρ(0,x^d)∇FhΦ(J(x^F^))∇^J,ρ′(0,x^d)Φ(J(x^F^)) ). (2.4) By the triangle inequality and, since $$0\leq\rho\leq 1$$, we infer   ‖∇^Φ^T‖L∞(T^) ≤‖∇FhΦ‖L∞(Fh)‖∇^J‖L∞(T^)+maxx^∈T^|ρ′(0,x^d)|‖Φ‖L∞(Fh). (2.5) Because $${\it{\Phi}}(J(0)) = 0$$, and using that $${\it{\Phi}}(J(\hat{x}_{\hat F})) = \int_{0}^1 \frac{d}{d s} {\it{\Phi}}( J (s \cdot \hat{x}_{\hat{F}})) \,\mathrm{d}s$$, we conclude $$\| {\it{\Phi}} \|_{L^\infty(F_h)} \leq \| \nabla_{\!F_h}{\it{\Phi}} \|_{L^\infty(F_h)} \| \hat\nabla J \|_{L^\infty(\hat T)}$$ to confirm the assertion. □ Next, we define for all $$\hat{x} \in \hat{T}$$ the nonlinear transformation   X(x^) =J(x^)+Φ^T(J(x^)) (2.6) and call $$T := {\mathcal{X}}(\hat{T})$$ the curved element related to $$T_h$$. Note that by the above construction, $$T$$ can have more than one curved face if $$d=3$$. If $$L_T<1$$, we conclude with Lemma 2.1 that $${\mathcal{X}}$$ is a $$C^1$$-diffeomorphism and the Jacobi matrix $$\hat \nabla {\mathcal{X}}$$ is positive definite (cf. Ciarlet & Raviart, 1972). In particular, it holds   (1−LT)‖∇^J‖L∞(T^)≤‖∇^X‖L∞(T^)≤(1+LT)‖∇^J‖L∞(T^). (2.7) Upon setting $$L := L_h \, \max((1+L_T) , (1-L_T)^{-1})$$, we get from (2.1) and (2.7) the bi-Lipschitz property of the nonlinear transformation   hL−1|z|≤‖∇^Xz‖L∞(T^)≤hL|z| for all z∈Rd. (2.8) Finally, we set $$\hat{{\it{\Phi}}}_T=0$$ on all inner elements and assume that (A5) there is a global constant $$L_T<1$$ such that for all boundary elements $$T_h\in\mathcal{T}_h$$ the extensions $$\hat{\it{\Phi}}_T$$ of $${\it{\Phi}}$$ according to (2.2) satisfies $$\| \hat{\nabla} \hat{{\it{\Phi}}}_T \|_{L^\infty(\hat{T})} \leq L_T \, \| \hat{\nabla} J \|_{L^\infty(\hat{T})}$$. Then, for a given triangulation $$\mathcal{T}_h$$ of $${\it{\Omega}}_h$$, we introduce the exact triangulation$$\mathcal{T}$$ of $${\it{\Omega}}$$ consisting of curved elements such that $${\it{\Gamma}}=\partial{\it{\Omega}}$$ is matched exactly. (A4) provides that $$\mathcal{T}$$ is a covering of $${\it{\Omega}}$$ and (A5) is related to the quasi-monotonicity of the geometric indicator of Bonito et al. (2013, 2016) and also implies that the approximation of $${\it{\Omega}}$$ by $${\it{\Omega}}_h$$ has reached the saturated state according to Dörfler & Rumpf (1998). Example 2.2 (a) Let $${\it{\Omega}} = B_1(0) \subset \mathbb{R}^2$$ the unit circle. To approximate $${\it{\Omega}}$$, we choose $${\it{\Omega}}_h^0 = \{(x,y)\in\mathbb{R}^2 : |x|+|y| \leq 1\}$$. A macro triangulation $$\mathcal{T}_h^0$$ of $${\it{\Omega}}_h^0$$ can be constructed from four triangles that all share the midpoint as a common node; see Fig. 2(left). One possible choice for a deformation of $${\it{\Omega}}_h^0$$ to $${\it{\Omega}}$$ is   Φ0(x,y)=(|x|+|y|R(x,y)−1)(xy), (2.9) where $$R(x,y) = \sqrt{x^2+y^2}$$. On the element $$T_h \in \mathcal{T}_h^0$$ indicated in Fig. 2(left), we easily verify that   ∇ΓhΦ0(x,y)=12(y/R(x,y)3−1−x/R(x,y)3+1) on Fh=∂Th∩Γh (2.10) such that on $$F_h$$ we have $$| \nabla_{{\it{\Gamma}}_{\!h}}{\it{\Phi}}^0 |^2 = 1 +\frac{1}{2} \frac{1 -2R(x,y)}{R(x,y)^4}$$. With $$R(x,y) \geq 1/\sqrt{2}$$ it holds $$| \nabla_{{\it{\Gamma}}_{\!h}}{\it{\Phi}}^0 | \leq \sqrt{2}-1 < 1/2$$ and Lemma 2.1 can be applied with $$L_\rho =1$$ and $$L_T<1$$ in order to satisfy (A5). Let $$\mathcal{T}_h^*$$ be a refinement of $$\mathcal{T}_h$$ by bisection and denote by $$\mathcal{I}_h^\ast {\it{\Phi}}^0$$ the piecewise affine interpolation of $${\it{\Phi}}^0$$ on $$\mathcal{T}_h^*$$. Then $$\mathcal{T}_h = \left( \mathcal{I}_h^\ast {\it{\Phi}}^0 \right)(\mathcal{T}_h^*)$$ is an interpolating triangulation of $${\it{\Omega}}$$ defining the approximating domain $${\it{\Omega}}_h$$ as shown in Fig. 2(right). $${\it{\Omega}}$$ is obtained from $${\it{\Omega}}_h$$ by applying the boundary deformation $${\it{\Phi}} = {\it{\Phi}}^0 \circ \left( \mathcal{I}_h^\ast {\it{\Phi}}^0 \right)^{-1}$$. (b) A simple example of a domain $${\it{\Omega}}$$ with only piecewise smooth boundary can be obtained from setting $${\it{\Phi}}^0(x,y)= 0$$ for $$x\leq 0$$. (c) The deformation $${\it{\Phi}}^0$$ according to (2.9) can also be used to map the rectangular cuboid $${\it{\Omega}}_h^0 = \{ (x,y,z) \in \mathbb{R}^3 : |x|+|y|\leq 1 , 0\leq z\leq 1 \}$$ to a cylindrical domain $${\it{\Omega}}$$. Note that on the top and bottom side of $${\it{\Omega}}_h$$, the deformation lies in the tangential plane to $${\it{\Gamma}}$$ and $${{\it{\Gamma}}_{\!h}}$$. For the cylindrical domain $${\it{\Omega}}$$, the deformation $${\it{\Phi}}^0$$ cannot be derived from a normal projection within a narrow band around $${\it{\Gamma}}$$. Fig. 2. View largeDownload slide Left: globally smoothly bounded domain $${\it{\Omega}}$$ and its polygonal approximation $${\it{\Omega}}_h^0$$ with a macro triangulation $$\mathcal{T}_h^0$$. Middle: $$\mathcal{T}_h^*$$ obtained from refinement of $$\mathcal{T}_h^0$$ by three steps of bisection. Right: interpolating polygonal domain $${\it{\Omega}}_h$$ with the triangulation $$\mathcal{T}_h$$. Fig. 2. View largeDownload slide Left: globally smoothly bounded domain $${\it{\Omega}}$$ and its polygonal approximation $${\it{\Omega}}_h^0$$ with a macro triangulation $$\mathcal{T}_h^0$$. Middle: $$\mathcal{T}_h^*$$ obtained from refinement of $$\mathcal{T}_h^0$$ by three steps of bisection. Right: interpolating polygonal domain $${\it{\Omega}}_h$$ with the triangulation $$\mathcal{T}_h$$. 2.3 Approximation of curved elements To quantify the deviation of a curved boundary face or edge from the polygonal interpolation, we define the geometric element indicator  λΓ(Th) :=h−1{‖∇^F^(X−J)‖L∞(F^) if Th∩Γh=Fh,‖∇^E^(X−J)‖L∞(E^) if Th∩Γh=Eh,  (2.11) where $$\hat\nabla_{\!\hat{F}}({\mathcal{X}} -J)$$ denotes the tangential gradient of $$({\mathcal{X}} -J)$$ in $$\hat F$$ given by the first $$d-1$$ partial derivatives. Likewise, $$\hat\nabla_{\!\hat{E}}({\mathcal{X}} -J)$$ is the tangential gradient in $$\hat E$$ subject to the first $$d-2$$ partial derivatives. Although $$\lambda_{\it{\Gamma}}(T_h)$$ is only defined on a boundary face or edge by (2.11), it in fact already characterizes the difference between $$T$$ and $$T_h$$ completely, i.e., from Lemma 2.1 and (2.1) we deduce   ‖∇^(X−J)‖L∞(T^) ≲hλΓ(Th), (2.12)  λΓ(Th) ≲LhLT. (2.13) The partial derivatives of $${\mathcal{X}}$$ or $$J$$ with respect to $$\hat x_1, \dots ,\hat x_d$$ build a complete set of linear independent tangential vectors to $${\it{\Gamma}}$$ or $${\it{\Gamma}}_h$$, respectively. For $$\hat{x} \in \hat{F}$$ we define   T(x^) :=∇^F^X(x^)=[∂^1X(x^),⋯,∂^d−1X(x^)], (2.14a)  Th(x^) :=∇^F^J=[∂^1J,⋯,∂^d−1J]. (2.14b) Then, the first fundamental form of $${\it{\Gamma}}$$ and $${\it{\Gamma}}_h$$ is given by the symmetric positive definite matrix   G=TTTandGh=ThTTh, (2.15) respectively. Given $$\hat{x}\in\hat{T}$$ and $$u\in H^1({\it{\Omega}})$$, we set $$\hat{u}(\hat{x}) := u\circ{\mathcal{X}}(\hat{x}) = u(x)$$. Likewise, for $$\hat{x}_{\hat{F}}\in\hat{F}$$ and $$v\in H^1(\partial{\it{\Gamma}} \cap T)$$, we set $$\hat{v}(\hat{x}_{\hat{F}}) := v\circ{\mathcal{X}}(\hat{x}_{\hat{F}}) = v(x_{F})$$. The gradients on $$\hat T$$ and $$\hat F$$ are related to respective gradients on $${\it{\Omega}}$$ and $${\it{\Gamma}}$$ by   ∇^u^=∇u∇^X and ∇^F^v^=∇ΓvT. (2.16) With the change of variables formulas   ∫T^u^det(∇X)dx^=∫X(T^)udx and ∫F^v^det(G)ds^=∫X(F^)vds (2.17) we get the identities   ∫T^∇^u^⋅∇^X−1∇^X−T∇^u^ det(∇X)dx^ =∫X(T^)∇u⋅∇udx, (2.18a)  ∫F^∇^F^v^⋅G−1∇^F^v^det(G)ds^ =∫X(F^)∇Γv⋅∇Γvds. (2.18b) Next, we want to get bounds for the perturbation of the transformation in terms of the geometric element indicator. Lemma 2.3 Let $$T_h\in\mathcal{T}_h$$ and $$T\in\mathcal{T}$$ be the corresponding curved element. Let $${\mathcal{X}}$$, $$\mathbb{T}$$, $$G$$, and $$J$$, $$\mathbb{T}_h$$, $$G_h$$ be defined as above. Suppose there is one face $$F\in{\mathcal{F}}^\partial$$ with $$T\cap{\it{\Gamma}} = F$$. Then,   ‖G−Gh‖ ≲h2λΓ(Th), (2.19a)  |det(G)−det(Gh)| ≲h2d−2λΓ(Th), (2.19b) where $$\|\cdot\|$$ denotes a compatible matrix norm. Moreover, it holds   ‖det(G)det(Gh)−1‖L∞(F^) ≲λΓ(Th), (2.19c)  ‖Th(detGdetGhG−1−Gh−1)ThT‖L∞(F^) ≲λΓ(Th). (2.19d) Proof. For any $$z\in\mathbb{R}^d$$ it holds $$z^T G z = \| \hat \nabla_{\!\hat F} {\mathcal{X}}\,z \|^2$$, and hence the eigenvalues of $$G$$ are bounded by $$h^{-2}L^{-2}|z|^2 \leq z^T G z \leq h^{-2}L^2|z|^2$$ due to (2.8). Analogously, we conclude that the eigenvalues of $$G_h$$ lie in $$[h L_h^{-1},h L_h]$$. We thus infer $$\| \mathbb{T} \|_{L^2(\hat F)}, \| \mathbb{T}_h \|_{L^2(\hat F)} \lesssim h$$ and $$\det(G), \det(G_h) \lesssim h^{2d-2}$$. By the definition of $$G$$, $$G_h$$ and $$\lambda_{\it{\Gamma}}(T_h)$$, we obtain   ‖G−Gh‖ =‖(T−Th)TT+ThT(T−Th)‖ ≤hλΓ(Th)(‖T‖+‖Th‖), (2.20) which yields the first assertion. We use a Taylor expansion to deduce that there exists $$0 \leq \theta \leq 1$$ such that   det(G)−det(Gh) =θ(Ddet(Gh))(G−Gh) (2.21)   =θdet(Gh)trace⁡(Gh−1(G−Gh)). (2.22) By Hölder’s inequality we conclude that   |det(G)−det(Gh)|≲|det(Gh)|‖Gh−1‖‖G−Gh‖ (2.23) to infer (2.19b). For (2.19c), we note that   det(G)−det(Gh)=det(G)−det(Gh)det(G)+det(Gh) (2.24) and that all powers of $$h$$ cancel after division by $$\det(G_h)$$. Finally, using the abbreviations $$q:=\sqrt{\det(G)}$$ and $$q_h:=\sqrt{\det(G_h)}$$ leads to   Th(qqhG−1−Gh−1)ThT=Th(q−qhqhG−1+G−1(Gh−G)Gh−1)ThT. (2.25) The application of (2.19a) and (2.19c) now concludes the proof. □ Corollary 2.4 Let $$T_h\in\mathcal{T}_h$$ and $$T\in\mathcal{T}$$ be the corresponding curved element. Let $${\mathcal{X}}$$ and $$J$$ be defined as above. We set $$A := \hat\nabla{\mathcal{X}}^T \hat\nabla{\mathcal{X}}$$ and $$A_h := \hat\nabla J^T \hat\nabla J$$. Suppose there is one face $$F\in{\mathcal{F}}^\partial$$ or edge $$E\in\mathcal{E}_h^\partial$$ such that $$T\cap{\it{\Gamma}} = F$$ or $$T\cap{\it{\Gamma}} = E$$, respectively. Then,   ‖det(∇^X)det(∇^J)−1‖L∞(T^) ≲λΓ(Th), (2.26a)  ‖∇^J(|det∇^X||det∇^J|A−1−Ah−1)∇^JT‖L∞(T^) ≲λΓ(Th). (2.26b) Proof. The assertions follow from the same arguments as in Lemma 2.3. □ 2.4 Lifts and geometric nonconformity estimates We first define the lifting of functions defined on the curved domain or surface onto their corresponding polyhedral or polygonal discretization and vice versa. Definition 2.5 For functions $$u_h:{\it{\Omega}}_h\to\mathbb{R}$$ and $$u:{\it{\Omega}}\to\mathbb{R}$$ we define the lift and inverse lift by   uhℓ(X∘J−1(x)):=uh(x) and u−ℓ(x):=uh(X∘J−1(x)) for x∈Ωh. (2.27) Likewise, for $$v_h:{\it{\Gamma}}_h\to\mathbb{R}$$ and $$v:{\it{\Gamma}}\to\mathbb{R}$$, we define   vhℓ(X∘J−1(x)):=vh(x) and v−ℓ(x):=v(X∘J−1(x)) for x∈Γh. (2.28) Next, we estimate the error which arises due to integration on perturbed domains in order to compare functions on the curved domain and on its polyhedral approximation. Lemma 2.6 For $$u_h, \phi^{-\ell}\in H^1({\it{\Omega}}_h)$$ it holds   |∫Ωuhℓϕdx−∫Ωhuhϕ−ℓdx| ≲∑Th∈ThλΓ(Th)‖uh‖L2(Th)‖ϕ−ℓ‖L2(Th), (2.29a)  |∫Ω∇uhℓ⋅∇ϕdx−∫Ωh∇uh⋅∇ϕ−ℓdx| ≲∑Th∈ThλΓ(Th)‖∇uh‖L2(Th)‖∇ϕ−ℓ‖L2(Th). (2.29b) Likewise, for $$v_h, \psi^{-\ell}\in H^1({{\it{\Gamma}}_{\!h}})$$ it holds   |∫Γvhℓψds−∫Γhvhψ−ℓds| ≲∑Th∈ThλΓ(Th)‖vh‖L2(Fh)‖ψ−ℓ‖L2(Fh), (2.29c)  |∫Γ∇Γvhℓ⋅∇Γψds−∫Γh∇Γhvh⋅∇Γhψ−ℓds| ≲∑Th∈ThλΓ(Th)‖∇Γvh‖L2(Fh)‖∇Γψ−ℓ‖L2(Fh). (2.29d) Proof. The identities follow from transformation to the reference domain and a pull-back. We start with the $$L^2$$ cases. By a localization to elements and with (2.17) we derive   ∫Ωuhℓϕdx−∫Ωhuhϕ−ℓdx =∑Thi∈Th∫T^u^h(|det∇^Xi|−|det∇^Ji|)ϕ^dx^ =∑Thi∈Th∫Thiuh(|det∇^Xi||det∇^Ji|−1)ϕ−ℓdx. Application of (2.26a) yields the assertions (2.29a). A similar argument and (2.19c) can be employed with the surface integral to get (2.29c). One can proceed in an analogous way for the $$H^1$$ cases. With (2.17) and (2.18), we obtain   ∫Ω∇uhℓ⋅∇ϕdx −∫Ωh∇uh⋅∇ϕ−ℓdx =∑Thi∈Th∫T^∇^u^h((Ai)−1 |det∇^Xi|−(Ahi)−1 |det∇^Ji|)⋅∇^ϕ^dx^ =∑Thi∈Th∫Thi∇uh∇^Ji(|det∇^X||det∇^J|(Ai)−1−(Ahi)−1)(∇^Ji)T⋅∇ϕ−ℓdx. Application of (2.26b) yields the assertions (2.29b) and, again, similar arguments and (2.19d) can be employed with the surface integrals to get (2.29d). □ For sufficiently fine triangulations, it is known that there is a norm equivalence for the lifting of functions in the volume (cf. Elliott & Ranner, 2013) and on the surface (see Dziuk, 1988; Demlow, 2009). Here, we can get the equivalence from Lemma 2.6. Lemma 2.7 Assume the triangulation $$\mathcal{T}_h$$ is sufficiently fine such that $$L_h L_T <1$$. Then the global geometry estimator satisfies $$\lambda_{\it{\Gamma}} := \max_{T_h\in\mathcal{T}_h}\lambda_{\it{\Gamma}}(T_h) <1$$ and for $$u\in H^1({\it{\Omega}})$$ and $$v\in H^1({\it{\Gamma}})$$, it holds   ‖u‖L2(Ω) ≈‖u−ℓ‖L2(Ωh),‖∇u‖L2(Ω) ≈‖∇u−ℓ‖L2(Ωh), (2.30a)  ‖v‖L2(Γ) ≈‖v−ℓ‖L2(Γh),‖∇Γv‖L2(Γ) ≈‖∇Γhv−ℓ‖L2(Γh), (2.30b) where $$a\approx b$$ means $$a\lesssim b$$ and $$b\lesssim a$$. Proof. Choosing $$\phi = u_h^{\ell}$$ in Lemma 2.6 yields   |‖uhℓ‖L2(Ω)2−‖uh‖L2(Ωh)2|≲λΓ‖uh‖L2(Ωh)2. (2.31) It follows $$\| u_h^\ell \|_{L^2({\it{\Omega}})}^2 \lesssim (1 +\lambda_{{\it{\Gamma}}}) \, \| u_h \|_{L^2({\it{\Omega}}_h)}^2$$ and $$\| u_h \|_{L^2({\it{\Omega}}_h)}^2 \lesssim (1 -\lambda_{{\it{\Gamma}}})^{-1} \, \| u_h^\ell \|_{L^2({\it{\Omega}})}^2$$. In an analogous way one can show (2.30b). □ 3. Problem formulation The finite element method for (1.1) is based on the weak formulation introduced next. Subsequently, we state the corresponding discrete problem for which we apply the notion of an exact triangulation and its simplicial approximation as described in the previous section. 3.1 Weak formulation of the continuous problem Let $${\it{\Omega}}_h^0 \subseteq \mathbb{R}^d$$ be a polyhedral domain with a triangulation $$\mathcal{T}_h^0({\it{\Omega}}_h^0)$$ that satisfies (A1)-(A3). Let $${\it{\Omega}}$$ be the domain resulting from a deformation of $${\it{\Omega}}_h^0$$ by the continuous mapping $${\it{\Phi}}^0$$ that satisfies (A4)-(A5). For each boundary facet $$F_h^i$$ of $$\mathcal{T}_h^{0}$$, the restriction of $${\it{\Phi}}^0$$ to $$F_h^i$$ is a $$C^1$$ function due to (A4). We set $${\it{\Gamma}}^i := {\it{\Phi}}^0(F_h^i)$$. Then $$\partial {\it{\Omega}} = \cup_i {\it{\Gamma}}^i$$ and each boundary patch $${\it{\Gamma}}^i$$ has a piecewise $$C^1$$ boundary $$\partial {\it{\Gamma}}^i$$. We introduce the spaces   U :=H1(Ω), (3.1a)  V :={v∈L2(Γ):∇Γv∣Γi∈L2(Γi), v∣Γi=v∣Γj on ∂Γi∩∂Γj} (3.1b) and denote their dual spaces by $$U^*$$ and $$V^*$$, respectively. Moreover, we define the weak form of the Laplace–Beltrami operator $$-{\it{\Delta}}_{{\it{\Gamma}}} : V \to \mathbb{R}$$ such that for each $$v\in V$$ it holds   ∫Γ−ψΔΓvds:=∫Γ∇Γψ⋅∇Γvds−∑Γi=1NΓ∫∂Γiψ∇Γv⋅nΓidσfor all ψ∈V. (3.2) Note that in general the sum in (3.2) does not vanish, because on $$\bar {\it{\Gamma}}^i \cap \bar {\it{\Gamma}}^j$$ ($$i\neq j$$) the co-normal vectors $$n_{\it{\Gamma}}^i$$ and $$n_{\it{\Gamma}}^j$$ are not parallel unless $${\it{\Gamma}}$$ is a global $$C^1$$ surface. Therefore, the additional condition (1.1d) is necessary to recover from (3.2) the strong form of the Laplace–Beltrami operator for piecewise $$C^2$$ surfaces. The weak form of (1.1) reads: Given $$f\in U^*$$ and $$g\in V^*$$, find $$(u,v) \in U\times V$$ such that, for all $$(\phi,\psi) \in U\times V$$,   ∫Ω∇ϕ⋅∇u+ϕudx+∫Γϕ(αu−βv)ds =∫Ωϕfdx, (3.3a)  ∫Γ∇Γψ⋅∇Γv+vψds−∫Γψ(αu−βv)ds =∫Γψgds. (3.3b) This can also be written in the standard variational form   a((u,v),(ϕ,ψ))=ℓ((ϕ,ψ))for all (ϕ,ψ)∈U×V, (3.4) with the bilinear and linear forms   a((u,v),(ϕ,ψ)) :=α∫Ω∇ϕ⋅∇u+ϕudx+β∫Γ∇Γψ⋅∇Γv+ψvds +∫Γ(αϕ−βψ)(αu−βv)ds, (3.5a)  ℓ((ψ,ϕ)) :=α∫Ωϕfdx+β∫Γψgds. (3.5b) In an analogous way to Elliott & Ranner (2013), we easily verify that $$a(\cdot,\cdot)$$ is a continuous and coercive bilinear form on $$U\times V$$. Hence, the Lax–Milgram theorem provides existence of a unique solution to (3.4). 3.2 Discrete problem Let $$\mathcal{T}_h$$ be a triangulation obtained from repeated refinement of $$\mathcal{T}_h^0$$ by bisection and subsequent deformation such that (A1)-(A5) are satisfied. $$\mathcal{T}_h$$ defines a polyhedral domain $${\it{\Omega}}_h\subset\mathbb{R}^{d}$$ with the boundary $${{\it{\Gamma}}_{\!h}} = \partial{\it{\Omega}}_h$$. We introduce the spaces of continuous finite element functions   Uh :={uh∈C(Ωh¯) :uh|Th∈P1(Th)∀Th∈Th}, (3.6a)  Vh :={vh∈C(Γh) : vh|Fh∈P1(Fh)∀Fh∈Fh∂}, (3.6b) where $$P_1$$ is the space of piecewise polynomials of maximal degree one. The discrete variational problem reads: Given $$f_h \in U_h$$ and $$g_h \in V_h$$ approximating $$f$$ and $$g$$, respectively, find $$(u_h,v_h)\in U_h\times V_h$$ such that   ah((uh,vh),(ϕh,ψh))=ℓh((ϕh,ψh))for all (ϕh,ψh)∈Uh×Vh (3.7) with the discrete bilinear and linear forms   ah((uh,vh),(ϕh,ψh)) :=α∫Ωh∇ϕh⋅∇uh+ϕhuhdx+β∫Γh∇Γhψh⋅∇Γhvh+ψhvhds +∫Γh(αϕh−βψh)(αuh−βvh)ds, (3.8a)  ℓh((ϕh,ψh)) :=α∫Ωhϕhfhdx+β∫Γhψhghds. (3.8b) Existence and uniqueness of solutions to (3.7) can be proven by the Lax–Milgram Theorem (cf. Elliott & Ranner, 2013). Note that the finite element method defined above is nonconforming since $$U_h$$ and $$V_h,$$ in general, are not subspaces of $$U$$ and $$V$$. Moreover, $$a_h$$ and $$\ell_h$$ differ from $$a$$ and $$\ell$$ since the integrals are defined over different domains. 4. Residual error estimation This section is devoted to the derivation of a residual estimator of the error measured in the energy norm for the coupled problem. We introduce the energy norm on $$U \times V$$ by   |||(ϕ,ψ)|||2:=a((ϕ,ψ),(ϕ,ψ))=α‖ϕ‖U2+β‖ψ‖V2+‖αϕ−βψ‖L2(Γ)2. (4.1) Given some finite element function $$(u_h,v_h) \in U_h\times V_h$$, the residual $$R$$ with regard to the variational problem (3.4) is defined by   ⟨R,(ϕ,ψ)⟩:=a((uhℓ,vhℓ),(ϕ,ψ))−ℓ(ϕ,ψ) for all (ϕ,ψ)∈U×V. (4.2) Proposition 4.1 Let $$(u,v)\in U\times V$$ be the solution of (3.4). Given $$(u_h,v_h) \in U_h\times V_h$$, the energy norm of the error satisfies   |||(u−uhℓ,v−vhℓ)|||=|||R|||∗. (4.3) Proof. We set $$e:= (u_h^{\ell} -u, v_h^{\ell}-v)$$ and easily check that for all $$(\phi,\psi) \in U\times V$$ we have $$\langle R,(\phi,\psi) \rangle = a(e,(\phi,\psi))$$. Hölder inequality directly implies $$\langle R,(\phi,\psi)\rangle \leq \left\lvert\!\left\lvert\!\left\lvert{e}\right\rvert\!\right\rvert\!\right\rvert\left\lvert\!\left\lvert\!\left\lvert{(\phi,\psi)}\right\rvert\!\right\rvert\!\right\rvert$$ for all $$(\phi,\psi) \in U\times V$$, and thus $$\left\lvert\!\left\lvert\!\left\lvert{R}\right\rvert\!\right\rvert\!\right\rvert_\ast \leq \left\lvert\!\left\lvert\!\left\lvert{e}\right\rvert\!\right\rvert\!\right\rvert$$. On the other hand, since $$\langle R,e\rangle = \left\lvert\!\left\lvert\!\left\lvert{e}\right\rvert\!\right\rvert\!\right\rvert^2$$, it follows   |||e|||=⟨R,e⟩|||e|||≤supw∈(U×V)∖{0}⟨R,w⟩|||w|||=|||R|||∗. □ We split the residual into $$R = R_h +R_\textrm{nc} +R_\textrm{osc}$$ with   ⟨Rh,(ϕ,ψ)⟩ :=ah((uh,vh),(ϕ−ℓ,ψ−ℓ))−ℓh(ϕ−ℓ,ψ−ℓ), (4.4a)  ⟨Rnc,(ϕ,ψ)⟩ :=a((uhℓ,vhℓ),(ϕ,ψ))−ah((uh,vh),(ϕ−ℓ,ψ−ℓ)) −ℓhℓ(ϕ,ψ) +ℓh(ϕ−ℓ,ψ−ℓ), (4.4b)  ⟨Rosc,(ϕ,ψ)⟩ :=ℓhℓ(ϕ,ψ)−ℓ(ϕ,ψ). (4.4c) $$R_h$$ is the discretization part of the residual, $$R_\textrm{nc}$$ is the nonconformity part which quantifies the quality of the polyhedral domain approximation $${\it{\Omega}}_h$$ with respect to the exact domain $${\it{\Omega}}$$ and $$R_\textrm{osc}$$ is the data approximation part of the residual, where $$\ell_h^\ell(\phi,\psi) := \alpha \int_{{\it{\Omega}}} f_h^{\ell} \phi \,\mathrm{d}x + \beta\int_{{\it{\Gamma}}} g_h^{\ell} \psi \,\mathrm{d}s$$. In the subsequent sections, computable error bounds for $$R_h$$ and $$R_\textrm{nc}$$ are derived and then combined to yield the total a posteriori error estimate. 4.1 Discretization residual The derivation of the upper bound for the discretization residual $$R_h$$ largely follows the classical approach of residual error estimators (see, e.g., Verfürth, 1996; Ainsworth & Oden, 2000). A recent account on the universal application of the principle is given by Carstensen et al. (2012). In the proof, we employ the stable interpolation operator $${\it{\Pi}}_h:U^{-\ell} = H^1({\it{\Omega}}_h)\to U_h$$ of Scott & Zhang (1990). For $$\phi\in H^1({\it{\Omega}}_h)$$ and any $$T_h\in\mathcal{T}_h$$, $$F_h\in{\mathcal{F}}_h$$, one verifies   ‖ϕ−Πhϕ‖L2(Th) ≲hT‖∇ϕ‖L2(ωTh), (4.5a)  ‖ϕ−Πhϕ‖L2(Fh) ≲hF1/2‖∇ϕ‖L2(ωFh), (4.5b) where and $$\omega_{F_h}$$ denotes the volume patch consisting of those two elements which share the common face $$F_h$$, and $$\omega_{T_h}$$ is the volume patch containing $$T_h$$ as well as all neighbouring elements sharing at least on common node with $$T_h$$. In an analogous way, we introduce the stable interpolation operator $${\it{\Pi}}_h^\partial:V^{-\ell} \to V_h$$ on the polyhedral boundary $${\it{\Gamma}}_h$$. It holds   ‖ψ−Πh∂ψ‖L2(Fh) ≲hF‖∇Γhψ‖L2(ωFh∂), (4.6a)  ‖ψ−Πh∂ψ‖L2(Eh) ≲hE1/2‖∇Γhψ‖L2(ωEh∂), (4.6b) with surface patches $$\omega_{F_h}^\partial \subset {{\it{\Gamma}}_{\!h}}$$ and $$\omega_{E_h}^\partial \subset {{\it{\Gamma}}_{\!h}}$$. Lemma 4.2 Let $$(u_h,v_h) \in U_h\times V_h$$ be a solution of the discrete problem (3.7). For all $$(\phi,\psi) \in U\times V$$, the discretization residual $$R_h$$ satisfies the estimate   ⟨Rh,(ϕ,ψ)⟩≲α(ηB+ηC)‖∇ϕ‖L2(Ω)+βηS‖∇Γψ‖L2(Γ), (4.7) where the error indicators on the right-hand side are given by   ηC :=‖hFh∂1/2(αuh−βvh+∂nuh)‖L2(Fh∂), (4.8a)  ηB :=‖hTh(Δuh−uh+fh)‖L2(Th) +‖hFh1/2[∂nuh]Fh‖L2(Fh∖Fh∂), (4.8b)  ηS :=‖hFh∂(ΔΓhvh−(1+β)vh+αuh+gh)‖L2(Fh∂) +‖hEh∂1/2[[∂nvh]]E∂‖L2(Eh∂). (4.8c) Here, the discrete norms are defined as summation over the elements or faces, respectively. By $$[\partial_n u_h]_{{\mathcal{F}}_h}$$ we denote the usual jump of the normal derivative at element faces, and $$[[ \partial_n v_h ]]_{\mathcal{E}^\partial}$$ is the corresponding jump quantity on the surface defined as in (1.1d). Proof. Let $$(\phi,\psi) \in U\times V$$, then $$({\it{\Pi}}_h\phi^{-\ell},{\it{\Pi}}^\partial_h\psi^{-\ell})\in U_h\times V_h$$ is in the kernel of $$R_h$$. Thus, a splitting of $$R_h$$ into element contributions and an integration by parts result in   ⟨Rh,(ϕ,ψ)⟩ =α∑Th∈Th(∫Th(−Δuh+uh−fh)(ϕ−ℓ−Πhϕ−ℓ)dx +∑Fh∈∂Th∫Fh∂nuh(ϕ−ℓ−Πhϕ−ℓ)ds) +α∑Fh∈Fh∂∫Fh(αuh−βvh)(ϕ−ℓ−Πhϕ−ℓ)ds +β∑Fh∈Fh∂(∫Fh(−ΔΓhvh+(1+β)vh−αuh−gh)(ψ−ℓ−Πh∂ψ−ℓ)ds +∑Eh∈∂Fh∫Eh∂nvh(ψ−ℓ−Πh∂ψ−ℓ)dσ). (4.9) In the sum over all $$T_h\in\mathcal{T}_h$$ we collect all contributions over the inner faces $$F_h \in {\mathcal{F}}_h\setminus{\mathcal{F}}_h^\partial$$ and move the boundary face parts to the coupling term. The Cauchy–Schwarz inequality yields   ⟨Rh,(ϕ,ψ)⟩ ≲α∑Th∈Th‖Δuh−uh+fh‖L2(Th)‖ϕ−ℓ−Πhϕ−ℓ‖L2(Th) +α∑Fh∈Fh∖Fh∂‖[∂nuh]Fh‖L2(Fh)‖ϕ−ℓ−Πhϕ−ℓ‖L2(Fh) +α∑Fh∈Fh∂‖αuh−βvh+∂nuh‖L2(Fh)‖ϕ−ℓ−Πhϕ−ℓ‖L2(Fh) +β∑Fh∈Fh∂‖ΔΓhvh−(1+β)vh+αuh+gh‖L2(Fh)‖ψ−ℓ−Πh∂ψ−ℓ‖L2(Fh) +β∑Eh∈Eh∂‖[[∂nvh]]Eh∂‖L2(Eh)‖ψ−ℓ−Πh∂ψ−ℓ‖L2(Eh). (4.10) Then the properties of the interpolation (4.5) and (4.6), and the norm equivalence of the lifted functions according to Lemma 2.7 complete the proof. □ 4.2 Nonconformity residual The geometric element indicator $$\lambda(T)$$ and the comparison results of Section 2 lead to an upper bound for $$R_\textrm{nc}$$ which is due to the polyhedral approximation of the curved domain $${\it{\Omega}}$$. Lemma 4.3 Assume that the triangulation $$\mathcal{T}_h$$ is sufficiently fine such that Lemma 2.7 holds. For any $$(\phi,\psi)\in U\times V$$, the nonconformity residual $$R_{\textrm{nc}}$$ satisfies   ⟨Rnc,(ϕ,ψ)⟩≲ηncB,0‖ϕ‖L2(Ω)+ηncB,1‖∇ϕ‖L2(Ω) +ηncC‖αϕ−βψ‖L2(Γ) +ηncS,0‖ψ‖L2(Γ)+ηncS,1‖∇Γψ‖L2(Γ), where the element indicators are given by   ηncB,0 :=α∑Th∈ThλΓ(Th)‖uh−fh‖L2(Th), (4.11a)  ηncB,1 :=α∑Th∈ThλΓ(Th)‖∇uh‖L2(Th), (4.11b)  ηncC :=∑Fh∈Fh∂λΓ(Fh)‖αuh−βvh‖L2(Fh), (4.11c)  ηncS,0 :=β∑Fh∈Fh∂λΓ(Fh)‖vh−gh‖L2(Fh), (4.11d)  ηncS,1 :=β∑Fh∈Fh∂λΓ(Fh)‖∇Γhvh‖L2(Fh). (4.11e) Proof. Let $$(\phi,\psi)\in U\times V$$. We rearrange terms to split $$\langle R_\textrm{nc}, (\phi,\psi) \rangle = \mathcal{I}_B +\mathcal{I}_S +\mathcal{I}_C$$ with   IB =α(∫Ω∇uhℓ⋅∇ϕ+(uhℓ−fhℓ)ϕdx−∫Ωh∇uh⋅∇ϕ−ℓ+(uh−fh)ϕ−ℓdx),IS =β(∫Γ∇Γvhℓ⋅∇Γψ+(vhℓ−ghℓ)ψds−∫Γh∇Γhvh⋅∇Γhψ−ℓ+(vh−gh)ψ−ℓds),IC =(∫Γ(αuhℓ−βvhℓ)(αϕ−βψ)ds−∫Γh(αuh−βvh)(αϕ−ℓ−βψ−ℓ)ds). To estimate the coupling term $$\mathcal{I}_C$$, we split the integral into its contributions from the surface elements and apply (2.17) and Lemma 2.6. This yields   IC =α∑Fh∈Fh∂∫Fh(det(G)det(Gh)−1)(αuh−βvh)(αϕ−ℓ−βψ−ℓ)ds ≲α∑Fh∈Fh∂λΓ(Th)‖αuh−βvh‖L2(Fh)‖αϕ−ℓ−βψ−ℓ‖L2(Γh). In an analogous way, we estimate the surface part $$\mathcal{I}_S$$, where we also use (2.18) for the gradient terms. It follows   IS =β∑Fh∈Fh∂∫Fh∇vh(det(G)det(Gh)T(G−1−Gh−1)TT)⋅∇ψ−ℓds +β∑Fh∈Fh∂∫Fh(det(G)det(Gh)−1)(vh−gh)ψ−ℓds ≲β∑Fh∈Fh∂λΓ(Th)(‖∇Fhvh‖L2(Fh)‖∇Fhψ−ℓ‖L2(Γh)+‖vh−gh‖L2(Fh)‖ψ−ℓ‖L2(Γh)). Likewise, the bulk term can be bounded by   IB ≲α∑Th∈ThλΓ(Th)(‖∇uh‖L2(Th)‖∇ϕ−ℓ‖L2(Ωh)+‖uh−fh‖L2(Th)‖ϕ−ℓ‖L2(Ωh)). Application of Lemma 2.7 completes the proof. □ 4.3 Overall error estimator As a direct consequence of the preceding Sects. 4.1 and 4.2, we devise an a posteriori estimator for the overall error of a finite element solution of the coupled problem, and prove its reliability and efficiency. In the subsequent error analysis, we employ piecewise constant approximations of the data. We define $$f_h:{\it{\Omega}}_h \to \mathbb{R}$$ by $$f_h\vert_{T_h} = \left\lvert{T}\right\rvert^{-1}\,\int_{T} f \,\mathrm{d}x$$ for each $$T_h\in \mathcal{T}_h$$, where $$T\in\mathcal{T}$$ is the corresponding curved element. In an analogous way, for each $$F_h\in {\mathcal{F}}_h^\partial$$, we set $$g_h\vert_{F_h} := \left\lvert{F}\right\rvert^{-1}\,\int_{F} g \,\mathrm{d}s$$, where $$F$$ is the corresponding curved boundary face. The oscillations of the data $$f$$ and $$g$$ are defined by   osc(f,Th) :=hTh‖f−ℓ−fh‖L2(Th)andosc(f,Th):=∑Th∈Thosc(f,Th), (4.12a)  osc(g,Fh) :=hFh‖g−ℓ−gh‖L2(Fh)andosc(g,Fh∂):=∑Fh∈Fh∂osc(g,Fh). (4.12b) Theorem 4.4 Assume that the triangulation $$\mathcal{T}_h$$ is sufficiently fine such that Lemma 2.7 holds. For the solution $$(u_h,v_h)\in U_h\times V_h$$ of the discrete problem (3.7), it holds   |||(u−uhℓ,v−vhℓ)|||≲ηh+ηnc+osch, (4.13)  ηh≲|||(u−uhℓ,v−vhℓ)|||+ηnc+osch, (4.14) with the discretization error indicator of (4.8) and the nonconformity indicator according to (4.11)   ηh :=ηB+ηC+ηS,ηnc :=hTηncB,0+ηncB,1+hF∂ηncS,0+ηncS,1+ηncC, and data oscillations $$\text{osc}_h := \text{osc}(\,f,{\mathcal{T}_h}) + \text{osc}(g,{{\mathcal{F}}_h})$$. Proof. We easily verify that whenever $$\phi$$ or $$\psi$$ is constant on $$T\in\mathcal{T}$$ or $$F\in{\mathcal{F}}^\partial$$, respectively, the local contribution to $$R_\textrm{osc}$$ vanishes. In an analogous way to $$f_h$$ and $$g_h$$, we define the elementwise mean values $$\phi_h$$ and $$\psi_h$$ of $$\phi$$ and $$\psi$$, respectively. Then for all $$(\phi,\psi) \in U\times V$$, it holds   ⟨Rosc,(ϕ,ψ)⟩ ≲∑T∈T‖fhℓ−f‖L2(T)‖ϕ−ϕh‖L2(T)+∑F∈F∂‖ghℓ−g‖L2(T)‖ψ−ψh‖L2(F) ≲∑Th∈ThhTh‖fh−f−ℓ‖L2(Th)‖∇ϕ−ℓ‖L2(Th) +∑Fh∈Fh∂hFh‖gh−g−ℓ‖L2(Th)‖∇Γψ−ℓ‖L2(Fh) ≲∑Th∈Thosc(f,Th)‖∇ϕ‖L2(Ω)+∑Fh∈Fh∂osc(g,Fh)‖∇Γψ‖L2(Γ). Then, the first part (reliability) is an immediate consequence of Lemmas 4.2 and 4.3. To prove the efficiency, we apply the standard technique due to Verfürth. We have to consider all terms contributing to the indicator $$\eta_h$$, and thus carry out the proof in several steps. (i) Efficiency of$$\boldsymbol{\eta_B}$$ For any $$T_h\in\mathcal{T}_h$$ and for any $$F_h\in{\mathcal{F}}_h$$, we introduce the (volume) bubble functions $$b_{T_h}$$ and $$b_{F_h}$$ with $$\text{supp} b_{T_h} = T_h$$ and $$\text{supp} b_{F_h} = \omega_{F_h}$$, where $$\omega_{F_h}$$ consists of the elements adjacent to $$F_h$$. The bubble functions are normalized such that $$\max_{x\in T_h} b_{T_h}(x)=1$$ and $$\max_{x\in F_h} b_{F_h}(x)=1$$. For any $$T_h\in\mathcal{T}_h$$ and $$F_h\in{\mathcal{F}}_h$$, let   rThB :=Δuh−uh+fh,rFhB :=[∂nuh],ϕTh :=bThrThB,ϕFh :=bFhrFhB. By the equivalence of norms in finite dimensional spaces and the properties of the bubble functions, one can show that   ‖rThB‖L2(Th)2 ≈∫ThrThBϕThdxand‖rFhB‖L2(Fh)2 ≈∫FhrFhBϕFhds, (4.15)  ‖ϕTh‖L2(Th) ≲‖rThB‖L2(Th)and‖ϕFh‖L2(ωFh) ≲hFh1/2‖rFhB‖L2(Fh). (4.16) Details can be found in Verfürth, 1996, 2013; Brenner & Carstensen, 2004. To estimate $${r^B_{T_h}}$$, we use (1.1) on the curved element $$T\in\mathcal{T}$$ corresponding to $$T_h\in\mathcal{T}_h$$. When integrating by parts, the boundary terms vanish due to the bubble functions. We get   ‖rThB‖L2(Th)2 ≲∫Th−∇uh⋅∇ϕTh−uhϕThdx+∫ThfhϕThdx +∫T∇u⋅∇ϕThℓ+uϕThℓdx−∫TfϕThℓdx. For the integrals over $$T_h$$, we insert the respective integrals over $$T$$ and apply Lemma 2.6   ‖rThB‖L2(Th)2 ≲∫T∇(u−uhℓ)⋅∇ϕThℓ+(u−uhℓ)ϕThℓdx+∫T(fhℓ−f)ϕThℓdx +λΓ(T)(‖∇uh‖L2(Th)‖∇ϕTh‖L2(Th)+‖uh−fh‖L2(Th)‖ϕTh‖L2(Th)). We apply the Cauchy–Schwarz inequality and then use Lemma 2.7 for the inverse lifting of $$\left\lVert{ \phi_{T_h}^{\ell} }\right\rVert_{L^2(T)}$$. Then, an inverse inequality and (4.16) yield   ‖rThB‖L2(Th) ≲hTh−1(‖u−uhℓ‖H1(T)+osc(f,T)) +λΓh(T)(hTh−1‖∇uh‖L2(Th)+‖uh−fh‖L2(Th)). On the inner facets $$F\in{\mathcal{F}}_h\setminus{\mathcal{F}}_h^\partial$$ we proceed similarly. After partial integration, we add (1.1), i.e.,   ‖rFhB‖L2(Fh)2 ≲∫ωFhΔuhϕFhdx+∫ωFh∇uh⋅∇ϕFhdx ≲∫ωFh(rThB+uh−fh)ϕFhdx+∫ωFh∇uh⋅∇ϕFhdx +∫ωF(f−u)ϕFhℓdx−∫ωF∇u⋅∇ϕFhℓdx. Lemma 2.6 leads to   ‖rFhB‖L2(Fh)2 ≲∫ωFhrThBϕFhdx +∫ωF(f−fhℓ)ϕFhℓdx +∫ωF∇(uhℓ−u)⋅∇ϕFhℓ+(uhℓ−u)ϕFhℓdx +λΓ(T)(‖∇uh‖L2(ωFh)‖∇ϕFh‖L2(ωFh)+‖uh−fh‖L2(ωFh)‖ϕFh‖L2(ωFh)) ≲‖rThB‖L2(ωFh)‖ϕFh‖L2(ωFh) +‖f−fhℓ‖L2(ωF)‖ϕFhℓ‖L2(ωF) +‖uhℓ−u‖H1(ωF)‖ϕFhℓ‖H1(ωF) +λΓ(T)(‖∇uh‖L2(ωFh)‖∇ϕFh‖L2(ωFh)+‖uh−fh‖L2(ωFh)‖ϕFh‖L2(ωFh)). With the Cauchy–Schwarz inequality, Lemma 2.7, an inverse inequality and (4.16), we deduce   ‖rFhB‖L2(Fh) ≲‖hTh1/2rThB‖L2(ωFh)+‖hTh−1/2(u−ℓ−uh)‖H1(ωF)+‖hTh1/2(f−ℓ−fh)‖L2(ωFh) +∑Th∈ωFhλΓ(T)(hTh−1/2‖∇uh‖L2(Th)+hTh1/2‖uh−fh‖L2(Th)). Since $$\eta_B = \left\lVert{h_{\mathcal{T}_h} {r^B_{T_h}}}\right\rVert_{L^2({\mathcal{T}_h})} + \left\lVert{h_{\mathcal{T}_h}^{1/2} {r^B_{F_h}}}\right\rVert_{L^2({{\mathcal{F}}_h\setminus {\mathcal{F}}_h^\partial})}$$, combining the previous estimates yields   ηB ≲‖u−uhℓ‖H1(Ω)+osc(f,T)+∑Th∈ThλΓ(T)(‖∇uh‖L2(Th)+hTh‖uh−fh‖L2(Th)). ≲|||(u−uhℓ,v−vhℓ)|||+osch+hTηncB,0+ηncB,1 (ii) Efficiency of$$\eta_C$$ For any boundary face $$F_h\in {\mathcal{F}}_h^\partial$$, let   rFhC:=αuh−βvh+∂nuhandϕFh:=bFhrFhC. Then, we estimate as before   ‖rFhC‖L2(Fh)2 ≲∫FhrFhCϕFhds−∫F(αu−βv+∂nu)ϕFhℓds ≲∫Th(rThB+uh−fh)ϕFhdx+∫Th∇uh⋅∇ϕFhdx+∫Fh(αuh−βvh)ϕFhds −∫T(u−f)ϕFhℓdx−∫T∇u⋅∇ϕFhℓdx−∫F(αu−βvh)ϕFhℓds ≲(‖rThB‖L2(Th)+‖f−fhℓ‖L2(T))‖ϕFh‖L2(T) +‖u−uhℓ‖H1(T)‖ϕFh‖H1(T)  +‖α(uhℓ−u)−β(vhℓ−v)‖L2(F)‖ϕFhℓ‖L2(Fh) +λΓ(F)(‖uh−fh‖L2(Th)‖ϕFh‖L2(Th)+‖∇uh‖L2(Th)‖∇ϕFh‖L2(Th) +‖(αuh−βvh)‖L2(Fh)‖ϕFh‖L2(Fh)). A trace inequality and an inverse estimate yield $$\left\lVert{\phi_{F_h}}\right\rVert_{L^2(F_h)} \lesssim h_{F_h}^{-1/2} \left\lVert{\phi_{F_h}}\right\rVert_{L^2(T_h)}$$. We conclude as before   ‖rFhC‖L2(Fh) ≲‖hTh1/2rThB‖L2(Th)+‖hTh1/2(f−ℓ−fh)‖L2(T)+‖hTh−1/2(u−ℓ−uh)‖H1(T) +‖α(u−uhℓ)−β(v−vhℓ)‖L2(F) +λΓ(F)(hFh−1/2‖∇uh‖L2(Th)+hFh1/2‖uh−fh‖L2(Th)+‖αuh−βvh‖L2(Fh)). Due to the properties of the bubble functions, this directly results in the estimate   ηC =‖hFh∂1/2rFhC‖L2(Fh∂) ≲‖u−uhℓ‖H1(Ω)+osc(f,T)+‖hF∂1/2(α(u−uhℓ)−β(v−vhℓ))‖L2(Γ) +∑Fh∈Fh∂λΓ(F)(‖∇uh‖L2(Th)+hFh‖uh−fh‖L2(Th)+hFh1/2‖αuh−βvh‖L2(Fh)) ≲|||(u−uhℓ,v−vhℓ)|||+osch+hTηncB,0+ηncB,1+hF∂1/2ηncC. (iii) Efficiency of$$\boldsymbol{\eta_S}$$ We introduce for $$F_h\in {\mathcal{F}}_h^\partial$$ and $$E_h\in\mathcal{E}_h^\partial$$ the (surface) bubble functions $$b_{F_h}^\partial$$ with $$\text{supp} b_{F_h}^\partial = F_h$$ and $$b_{E_h}^\partial$$ with $$\text{supp} b_{E_h}^\partial = \omega_{E_h}^\partial$$. For any $$F\in{{\mathcal{F}}_h}^\partial$$ and for any $$E\in{\mathcal{E}_h}^\partial$$, let   rFhS :=ΔΓhvh+αuh−(1+β)vh+gh,rEhS :=[[∂nvh]],ψFh :=bFh∂rFhS,ψEh :=bEh∂rEhS. From the properties of the bubble functions, we deduce   ‖rFhS‖L2(Fh)2 ≈∫FhrFhSψFhdsand‖rEhS‖L2(Eh)2 ≈∫EhrEhSψEhdσ,‖ψFh‖L2(Fh) ≲‖rFhS‖L2(Fh)and‖ψEh‖L2(ωEh∂) ≲hEh1/2‖rEhS‖L2(Eh). As before, we integrate (1.1b) and (1.1c) over $$F$$, and integrate by parts to get   ‖rFhS‖L2(Fh)2 ≲∫Fh∇Γhvh⋅∇ΓhψFh+(αuh−(1+β)vh+gh)ψFhds −∫F∇Γv⋅∇ΓψFhℓ+(αu−(1+β)v+g)ψFhℓds. With Lemma 2.6, the properties of the bubble function $$b_{F_h}^\partial$$ and an inverse inequality, we deduce   ‖rFhS‖L2(Fh) ≲hFh−1(‖vhℓ−v‖H1(Fh) +osc(g,F))+‖α(uhℓ−u)−β(vhℓ−v)‖L2(Fh) +λΓ(F)(hF−1‖∇Γhvh‖L2(Fh)+‖vh−gh‖L2(Fh)+‖αuh−βvh‖L2(Fh)). Moreover, with the properties of $$b_{E_h}^\partial$$ and Lemma 2.6,   ‖rEhS‖L2(Eh)2 ≲∫ωEh∂rFhSψEhds +∫ωEh∂∇Γhvh⋅∇ΓhψEh+(vh−gh)ψEhds +∫ωEh∂(βvh−αuh)ψEhds ≲∫ωEh∂rFhSψEhds +∫ωE∂∇Γ(vhℓ−v)⋅∇ΓψEhℓ+(vhℓ−v)ψEhℓds +∫ωE∂(β(vhℓ−v)−α(uhℓ−u))ψEhℓds +∫ωE∂(g−ghℓ)ψEhℓds +∑Fh∈ωEh∂λΓ(Fh)(‖∇Γhvh‖L2(ωEh∂)‖∇ΓhψEh‖L2(ωEh∂) +(‖vh−gh‖L2(ωEh∂)+‖αuh−βvh‖L2(ωEh∂))‖ψEh‖L2(ωEh∂)). It immediately follows   ‖rEhS‖L2(Eh) ≲hF1/2‖rFhS‖L2(ωEh∂)+hF−1/2‖vhℓ−v‖H1(ωE∂)+hF−1/2osc(g,ωEh∂) +hF1/2‖α(uhℓ−u)−β(vhℓ−vh)‖L2(ωE∂) +∑Fh∈ωEh∂λΓ(Fh)(hFh−1/2‖∇Γhvh‖L2(Fh)+hFh1/2‖vh−gh‖L2(Fh) +hFh1/2‖αuh−βvh‖L2(Fh)). With $$\eta_S = \left\lVert{h_{{\mathcal{F}}_h} {r^S_{F_h}}}\right\rVert_{L^2({\mathcal{F}}_h^\partial)} + \left\lVert{h_{{\mathcal{F}}_h}^{1/2}\, {r^S_{E_h}}}\right\rVert_{L^2(\mathcal{E}_h^\partial)}$$ and the finite overlap of patch elements,   ηS ≲‖v−vhℓ‖H1(Γ)+‖hF(α(uhℓ−u)−β(vhℓ−vh))‖L2(Γ)+osc(g,Fh∂) +∑Fh∈Fh∂λΓ(Fh)(‖∇Γhvh‖L2(Fh)+hFh‖βvh−gh‖L2(Fh)+hFh‖αuh−βvh‖L2(Fh)). This verifies the claimed statement and completes the proof. □ 5. Numerical experiments We implemented a finite element method based on the concepts outlined by Alberty et al. (1999). The volume grid is described in two arrays: c4n contains the coordinates of the nodes and thereby defined the global node indices. n4e lists for each element $$T_h$$ the global nodes indices of the corners, and thereby defines the local node and face indices. For the surface representation, we introduce the additional array n4f containing for each surface element $$F\in {\mathcal{F}}^\partial$$ the global node indices of its corners. In the implementation, we do not follow the above convention that within an element the boundary face has to be labelled last. Instead, there is an array bdy_face containing for each $$T_h\in\mathcal{T}_h$$ either the value $$0$$ for interior elements or the local face index $$i\in\{ 1,\ldots,(d+1) \}$$ of the boundary face $$F_h \subset T_h$$. To compute the $$P1$$ stiffness matrix in the volume, on each element $$T_h\in\mathcal{T}_h$$ we need to determine the constant gradient vectors $$\nabla \phi_i$$ for $$i=1,\ldots,(d+1)$$. On each surface element $$F_h\in{\mathcal{F}}_h^\partial$$ we have to evaluate $$\nabla \psi_j$$ for $$j=1,\ldots,d$$. Suppose $$F_h \in {\mathcal{F}}_h^\partial$$ and $$T_h\in\mathcal{T}_h$$ such that $$F_h \subset T_h$$. Let $$\sigma\in\{ 1,\ldots,(d+1) \}$$ denote the local face index of the boundary face $$F_h$$. Then, the outer normal to $${\it{\Gamma}}_h$$ on $$F_h$$ is   nFh=−∇ϕσT/|∇ϕσ|. (5.1) Within $$T_h$$, the corner nodes of $$F_h\subset T_h$$ have the local indices $$i\in\{1,\ldots,d+1\} \setminus \sigma$$, and the gradients $$\nabla \psi_j$$, of the surface finite element shape functions are given by   ∇ψj=∇ϕi(j)T−(∇ϕi(j)T⋅nFh)nFh for j=1,…,d. (5.2) In the adaptive algorithm, elements are marked for refinement using the indicator   ηTh:=αηC,Th+αηB,Th+14βηS,Th, (5.3) where the terms on the right-hand side refer to the error contributions of element $$T_h\in \mathcal{T}_h$$ to the corres-ponding estimators defined in (4.8). The elements with the largest values of $$\eta_{T_h}$$ are marked for refinement until the marked elements contribute to three-fourth of the total estimated error. For the mesh refinement, we apply the bisection algorithm of Bartels & Schreier (2012). Note that marking a boundary element $$T_h\in\mathcal{T}_h$$ due to a large surface or coupling indicators on $$F_h\in{\mathcal{F}}_h^\partial$$ or $$F_h\subset T_h$$ might not lead to an immediate refinement of the surface element because the refinement edge of $$T_h$$ might be in the interior of $${\it{\Omega}}_h$$. For uniformly refined grids, the number of degrees of freedom (DoF) in the volume is proportional to $$h_T^{-3}$$ and the number of DoFs on the surface is proportional to $$h_F^{-2}$$. In the following we will call a quantity proportional to $$h_T$$ or $$h_F$$ whenever it shows a respective dependence on the DoFs in the volume or on the surface, respectively. We set $$\alpha=\beta=1$$ in the following experiments. 5.1 Convergence to known regular solution As a first benchmark, we verify our method with the example from (Elliott & Ranner, 2013). On the unit sphere $${\it{\Omega}} =B_{1}(0) \subset\mathbb{R}^3$$ with boundary $${\it{\Gamma}}=\partial{\it{\Omega}}$$, the data $$f$$ and $$g$$ are prescribed such that the exact solution of the problem (3.3) is   u(x,y,z) =exp⁡(−x(x−1)−y(y−1)), (5.4a)  v(x,y,z) =exp⁡(−x(x−1)−y(y−1))[1+x(1−2x)+y(1−2y)]. (5.4b) The computed solution $$(u_h,v_h)$$ is plotted in Fig. 3. From Fig. 4, we see that both the uniform and the adaptive refinement yield the same asymptotic behaviour with respect to the number of DoFs. On uniformly refined grids, the $$H^1({\it{\Omega}})$$-error of $$u_h$$ decreases proportionally to $$h_T$$. Likewise, the $$H^1({\it{\Gamma}})$$-error of $$v_h$$ is proportional to $$h_F$$. Examining Fig. 5 confirms that the estimator $$\eta_h$$ and its contributions $$\eta_B$$ and $$\eta_S$$ related to the volume and surface correctly follow the reduction of the overall error and the individual $$H^1$$ errors of $$u_h$$ and $$v_h$$, respectively. Moreover, we observe that during the first refinement steps, the estimator $$\eta_S$$ of the surface error is dominant, leading to some extra refinement of $${\it{\Gamma}}_h$$ in the adaptive algorithm, cf. Fig. 3. Due to the high regularity of the prescribed solution, the later refinement steps are almost uniform. The estimate $$\eta_C$$ of the coupling error is at least about one order of magnitude smaller than the other error contributions, and thus does not have any influence on the refinement. Additional evaluation of the $$L^2$$ errors revealed that uniform refinement leads to $$L^2$$ error reduction with optimal order $$O(h_T^2)$$ in the volume and $$O(h_F^2)$$ on the surface. When refining the mesh according to the $$H^1$$ estimator $$\eta_{T_h}$$, we still observe the $$L^2$$ error reduction proportional to the $$h_T^2$$ in the volume, but slightly lower than $$h_F^2$$ on the surface. Fig. 3. View largeDownload slide Experiment of Section 5.1: cross-section of the solution $$u_h$$ and the adaptively refined grid (left) and solution $$v_h$$ on the surface (right). Fig. 3. View largeDownload slide Experiment of Section 5.1: cross-section of the solution $$u_h$$ and the adaptively refined grid (left) and solution $$v_h$$ on the surface (right). Fig. 4. View largeDownload slide Experiment of Section 5.1: bulk, surface and overall error with respect to the exact solution in the $$H^1$$-norm over the number of unknowns for uniform refinement (left) and adaptive refinement (right). Fig. 4. View largeDownload slide Experiment of Section 5.1: bulk, surface and overall error with respect to the exact solution in the $$H^1$$-norm over the number of unknowns for uniform refinement (left) and adaptive refinement (right). Fig. 5. View largeDownload slide Experiment of Section 5.1: discretization error estimator $$\eta_h$$ with contributions $$\eta_B, \eta_S,\eta_C$$ over the number of unknowns with uniform refinement (left) and adaptive refinement (right). Fig. 5. View largeDownload slide Experiment of Section 5.1: discretization error estimator $$\eta_h$$ with contributions $$\eta_B, \eta_S,\eta_C$$ over the number of unknowns with uniform refinement (left) and adaptive refinement (right). 5.2 Nonsmooth boundary and grid adaption due to data In the next experiment with $${\it{\Omega}} = [-1,1]^3$$, we want to examine the effect of a domain $${\it{\Omega}}$$ with nonsmooth boundary $${\it{\Gamma}}$$ and the ability of the adaptive algorithm to resolve singularities of the solution. Then, the surface $${\it{\Gamma}}$$ is only piecewise smooth. Since all faces are flat (i.e., affine), the discretization of the surface $${\it{\Gamma}}_h = {\it{\Gamma}}$$ is exact, and hence $$\lambda_{\it{\Gamma}} = 0$$. We set   xf :=(0.7,0.6,0.5),f(x) :=(|x−xf|2+10−4)−1, (5.5a)  xg :=(0.1,−1.0,−0.3),g(x) :=(|x−xg|3/2+10−5)−1. (5.5b) By this, we introduce ‘weak singularities’ in $${\it{\Omega}}$$ and on $${\it{\Gamma}}$$. The computed solution is plotted in Figure 6 with a cross-section of the bulk solution $$u_h$$ on the left-hand side. Figure 7 shows that adaptive refinement leads to reduction rates of $$\eta_B$$ and $$\eta_S$$ proportional to $$O(h_T)$$ and $$O(h_F)$$, respectively, whereas for uniform refinement the reduction rates are lower and do not reach the same asymptotic behaviour. In particular, we see in Fig. 7, that using uniform refinement the reduction of the surface estimate $$\eta_S$$ is considerably slower than $$O(h_F)$$, while for $$\eta_B$$ the difference to the $$O(h_T)$$ rate is less pronounced. Keeping in mind that the volume grid has more DoFs than the according surface grid, we observe that in each uniform refinement step, the indicator $$\eta_S$$ on the surface is larger than $$\eta_B$$ in the volume. As before, the coupling error $$\eta_C$$ is only marginal with regard to $$\eta$$. Moreover, we see from Fig. 6 that the adaptive algorithm has refined the grid towards the weak singularities of the data. These singularities are well separated from the nonsmooth edges of $${\it{\Gamma}}$$ where no extra refinement of the grid can be observed. Fig. 6. View largeDownload slide Experiment of Section 5.2: Solution for data specified by (5.5): Cross-section of $$u_h$$ with the locally refined grid (left) and $$v_h$$ (right). Fig. 6. View largeDownload slide Experiment of Section 5.2: Solution for data specified by (5.5): Cross-section of $$u_h$$ with the locally refined grid (left) and $$v_h$$ (right). Fig. 7. View largeDownload slide Experiment of Section 5.2: Uniform refinement leads to a reduction rate of the estimated error $$\eta_h$$ slower than $$O(h_T)$$ in the volume and $$O(h_F)$$ on the surface (left), whereas for the adaptively refined grids $$\eta_h$$ meets the optimal rates (right). Fig. 7. View largeDownload slide Experiment of Section 5.2: Uniform refinement leads to a reduction rate of the estimated error $$\eta_h$$ slower than $$O(h_T)$$ in the volume and $$O(h_F)$$ on the surface (left), whereas for the adaptively refined grids $$\eta_h$$ meets the optimal rates (right). 5.3 Domain with corner singularities In the final experiment we examine the effect of corner singularities. Note that the classical example of a re-entrant corner and $$u(r,\phi) = r^{\alpha} \sin(\alpha \phi)$$ cannot be transferred to the coupled problem (1.1). Let the domain $${\it{\Omega}}$$ be given by a ball of radius 1 where one octant is removed (see Fig. 8). We expect $$u$$ and $$v$$ to be of low regularity at the reentrant corner points of the domain $${\it{\Omega}}$$ and the surface $${\it{\Gamma}}$$, respectively. We set $$x_1 = \frac{1}{2} (-1, 1, 1)$$, $$x_2 = \frac{1}{2} ( 1,-1, 1)$$, $$x_3 = \frac{1}{2} ( 1, 1,-1)$$, and choose   f(x) =0, (5.6a)  g(x) =20−20(|x|−exp⁡(−4|x−x1|2)−exp⁡(−4|x−x2|2)+exp⁡(−4|x−x3|2)). (5.6b) Fig. 8. View largeDownload slide Experiment of Section 5.3: The domain $${\it{\Omega}}$$ with the surface patches $${\it{\Gamma}}^i$$, $$i=1,2,3,4$$ (left). Refinement of the grid towards the corner at the north pole $$(0,0,1)$$ due to the adaptive algorithm (right). Fig. 8. View largeDownload slide Experiment of Section 5.3: The domain $${\it{\Omega}}$$ with the surface patches $${\it{\Gamma}}^i$$, $$i=1,2,3,4$$ (left). Refinement of the grid towards the corner at the north pole $$(0,0,1)$$ due to the adaptive algorithm (right). The computed solution $$(u_h,v_h)$$ is plotted in Fig. 9. Fig. 9. View largeDownload slide Experiment of Section 5.3: Solution $$u_h$$ (left) and $$v_h$$ (right). Fig. 9. View largeDownload slide Experiment of Section 5.3: Solution $$u_h$$ (left) and $$v_h$$ (right). Similar to the experiment of Section 5.2, we conclude from Fig. 10 that adaptive refinement leads to reduction rates of $$\eta_B$$ and $$\eta_S$$ proportional to $$O(h_T)$$ and $$O(h_F)$$, respectively, whereas for uniform refinement the reduction rates are lower and do not reach the same asymptotic behaviour. However, when compared to Fig. 7, the deviation of $$\eta_S$$ from the $$O(h_F)$$ line is considerably less pronounced in this experiment. This indicates that for the coupled system the geometrically induced singularities are weaker than those known in the pure volume or the pure surface problem. At least, it is not possible to trigger stronger singularities by using data (5.6). Fig. 10. View largeDownload slide Experiment of Section 5.3: Uniform refinement leads to a reduction rate of the estimated error $$\eta_h$$ and the volume and surface contributions $$\eta_B$$ and $$\eta_S$$ slower than $$O(h_T)$$ and $$O(h_F)$$, respectively (left). For the adaptively refined grids $$\eta_h$$, the optimal rates are attained (right). Fig. 10. View largeDownload slide Experiment of Section 5.3: Uniform refinement leads to a reduction rate of the estimated error $$\eta_h$$ and the volume and surface contributions $$\eta_B$$ and $$\eta_S$$ slower than $$O(h_T)$$ and $$O(h_F)$$, respectively (left). For the adaptively refined grids $$\eta_h$$, the optimal rates are attained (right). References Ainsworth, M. & Oden, J. ( 2000) A Posteriori Error Estimation in Finite Element Analysis . Pure and Applied Mathematics (New York). New York: Wiley-Interscience [John Wiley & Sons], pp. xx+240. Google Scholar CrossRef Search ADS   Alberty, J., Carstensen, C. & Funken, S. ( 1999) Remarks around 50 lines of matlab: short finite element implementation. Numer. Algorithms , 20, 117– 137. Google Scholar CrossRef Search ADS   Babuška, I. & Rheinboldt, W. C. ( 1978) A-posteriori error estimates for the finite element method. Int. J. Numer. Meth. Engng. , 12, 1597– 1615. Google Scholar CrossRef Search ADS   Bartels, S. & Schreier, P. ( 2012) Local coarsening of simplicial finite element meshes generated by bisections. BIT , 52, 559– 569. Google Scholar CrossRef Search ADS   Bernardi, C. ( 1989) Optimal finite-element interpolation on curved domains. SIAM J. Numer. Anal. , 26, 1212– 1240. Google Scholar CrossRef Search ADS   Bertalmío, M., Cheng, L.-T., Osher, S. & Sapiro, G. ( 2001) Variational problems and partial differential equations on implicit surfaces. J. Comput. Phys. , 174, 759– 780. Google Scholar CrossRef Search ADS   Bonito, A., Cascón, J. M., Morin, P. & Nochetto, R. H. ( 2013) AFEM for geometric PDE: the Laplace-Beltrami operator. Analysis and Numerics of Partial Differential Equations  ( Brezzi, F. Franzone, P. C. Gianazza U. & Gilardi G. eds). Milan, Italy: Springer, pp. 257– 306. Google Scholar CrossRef Search ADS   Bonito, A., Cascón, J. M., Morin, P., Mekchay, K. & Nochetto, R. H. ( 2016) High-order AFEM for the Laplace-Beltrami operator: Convergence rates. Found. Comput. Math. , 16, 1473– 1539. Google Scholar CrossRef Search ADS   Brenner, S. & Carstensen, C. ( 2004) Finite element methods. Encyclopedia of computational mechanics  ( Stein, E. De Borst R. & Hughes T. J. R. eds), Chichester: John Wiley & Sons, 73– 117. Burger, M. ( 2009) Finite element approximation of elliptic partial differential equations on implicit surfaces. Comput. Vis. Sci. , 12, 87– 100. Google Scholar CrossRef Search ADS   Burman, E., Claus, S., Hansbo, P., Larson, M. & Massing, A. ( 2015a) CutFEM: Discretizing geometry and partial differential equations. Int. J. Numer. Meth. Engng. , 104, 472– 501. Google Scholar CrossRef Search ADS   Burman, E., Hansbo, P., Larson, M. & Zahedi, S. ( 2015b) Cut finite element methods for coupled bulk–surface problems. Numer. Math. , 133, 203– 231. Google Scholar CrossRef Search ADS   Camacho, F. & Demlow, A. ( 2015) $$L_2$$ and pointwise a posteriori error estimates for FEM for elliptic PDEs on surfaces. IMA J. Num. Anal. , 35, 1199– 1227. Google Scholar CrossRef Search ADS   Carstensen, C., Eigel, M., Hoppe, R. & Löbhard, C. ( 2012) A review of unified a posteriori finite element error control. Numer. Math. Theor. Meth. Appl. , 5, 509– 558. Google Scholar CrossRef Search ADS   Chernyshenko, A. Y. & Olshanskii, M. A. ( 2015) An adaptive octree finite element method for PDEs posed on surfaces. Comput. Methods Appl. Mech. Eng. , 291, 146– 172. Google Scholar CrossRef Search ADS   Ciarlet, P. & Raviart, P.-A. ( 1972) Interpolation theory over curved elements, with applications to finite element methods. Comput. Methods Appl. Mech. Eng. , 1, 217– 249. Google Scholar CrossRef Search ADS   Dassi, F. ( 2014) Finite element techniques for curved boundaries. Ph.D. Thesis , Milan, Italy: Politecnico di Milano. Dassi, F., Perotto, S. & Formaggia, L. ( 2015) A priori anisotropic mesh adaptation on implicitly defined surfaces. SIAM J. Sci. Comp. , 37, A2758– A2782. Google Scholar CrossRef Search ADS   Deckelnick, K., Dziuk, G., Elliott, C. & Heine, C.-J. ( 2010) An h-narrow band finite-element method for elliptic equations on implicit surfaces. IMA J. Num. Anal. , 30, 351– 376. Google Scholar CrossRef Search ADS   Dedner, A. & Madhavan, P. ( 2016) Adaptive discontinuous galerkin methods on surfaces. Numer. Math. , 132, 369– 398. Google Scholar CrossRef Search ADS   Demlow, A. ( 2009) Higher-order finite element methods and pointwise error estimates for elliptic problems on surfaces. SIAM J. Numer. Anal. , 47, 805– 827. Google Scholar CrossRef Search ADS   Demlow, A. & Dziuk, G. ( 2007) An adaptive finite element method for the Laplace-Beltrami operator on implicitly defined surfaces. SIAM J. Numer. Anal. , 45, 421– 442. Google Scholar CrossRef Search ADS   Dörfler, W. & Rumpf, M. ( 1998) An adaptive strategy for elliptic problems including a posteriori controlled boundary approximation. Math. Comp. , 67, 1361– 1382. Google Scholar CrossRef Search ADS   Dziuk, G. ( 1988) Finite elements for the Beltrami operator on arbitrary surfaces. Partial Differential Equations and Calculus of Variations  ( Hildebrand S. & Leis R. eds). Lecture Notes in Mathematics. Heidelberg, Berlin: Springer, pp. 142– 155. Google Scholar CrossRef Search ADS   Dziuk, G. & Elliott, C. ( 2013) Finite element methods for surface PDEs. Acta Numer. , 22, 289– 396. Google Scholar CrossRef Search ADS   Egger, H., Fellner, K., Pietschmann, J.-F. & Tang, B. ( 2015) A finite element method for volume-surface reaction-diffusion systems. Preprint , IGDK-2015-26. Elliott, C. & Ranner, T. ( 2013) Finite element analysis for a coupled bulk-surface partial differential equation. IMA J. Num. Anal. , 33, 377– 402. Google Scholar CrossRef Search ADS   Giese, W., Eigel, M., Westerheide, S., Engwer, C. & Klipp, E. ( 2015) Influence of cell shape, inhomogeneities and diffusion barriers in cell polarization models. Phys. Biol. , 12, 066014. Google Scholar CrossRef Search ADS PubMed  Gross, S., Olshanskii, M. A. & Reusken, A. ( 2015) A trace finite element method for a class of coupled bulk-interface transport problems. ESAIM: M2AN , 49, 1303– 1330. Google Scholar CrossRef Search ADS   Novak, I. L., Gao, F., Choi, Y.-S., Resasco, D., Schaff, J. C. & Slepchenko, B. M. ( 2007) Diffusion on a curved surface coupled to diffusion in the volume: Application to cell biology. JCP , 226, 1271– 1290. Google Scholar CrossRef Search ADS   Olshanskii, M., Reusken, A. & Grande, J. ( 2009) A finite element method for elliptic equations on surfaces. SIAM J. Numer. Anal. , 47, 3339– 3358. Google Scholar CrossRef Search ADS   Rätz, A. & Röger, M. ( 2012) Turing instabilities in a mathematical model for signaling networks. J. Math. Biol. , 65, 1215– 1244. Google Scholar CrossRef Search ADS PubMed  Scott, L. R. ( 1973) Finite element techniques for curved boundaries. Ph.D. Thesis , Massachusetts, Cambridge: Massachusetts Institute of Technology. Scott, L. R. & Zhang, S. ( 1990) Finite element interpolation of nonsmooth functions satisfying boundary conditions. Math. Comp. , 54, 483– 493. Google Scholar CrossRef Search ADS   Verfürth, R. ( 1996) A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques . Chichester, New York, Brisbane, Toronto, Singapore, Stuttgart, Leipzig: Wiley-Teubner. Verfürth, R. ( 2013) A Posteriori Error Estimation Techniques for Finite Element Methods . Oxford, UK: Oxford University Press. Google Scholar CrossRef Search ADS   Zlámal, M. ( 1973) Curved elements in the finite element method. I. SIAM J. Numer. Anal. , 10, 229– 240. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations