Error analysis of an unfitted HDG method
for a class of non-linear elliptic problems.

Given the domain $\Omega$ where (1.1) is posed, we will define a family of polygonal subdomains and admissible triangulations approximating $\Omega$ where we will ultimately pose our discretization. First, consider a family of simply connected domains $\{\Omega_{\alpha}\}_{\alpha>0}$ such that, for every $\alpha$ the following conditions hold: (1) $\Omega_{\alpha}\subseteq\Omega$, (2) the boundary $\Gamma_{\alpha}:=\partial\Omega_{\alpha}$ is a polygon, and (3) for every $\epsilon>0$ there are infinitely many indices $\alpha$ such that $\lambda(\Omega\setminus\Omega_{\alpha})<\epsilon$. In the preceding expression $\lambda(\cdot)$ denotes the Lebesgue measure. These conditions ensure that the family of subdomains $\{\Omega_{\alpha}\}_{\alpha>0}$ will exhaust $\Omega$.

Having built the family $\{\Omega_{\alpha}\}_{\alpha>0}$ satisfying all the conditions above, the next step is to define a family of admissible simplicial triangulations $\{\mathcal{T}_{h}\}_{h>0}$. To be considered admissible, a triangulation $\mathcal{T}_{h}$ must be such that: (1) $\mathcal{T}_{h}$ is a triangulation for at least one $\Omega_{h}\in\{\Omega_{\alpha}\}_{\alpha>0}$ (we will identify both $\mathcal{T}_{h}$ and the respective domain $\Omega_{h}$ with the same subscript $h$, adding copies of $\Omega_{h}$ to $\{\Omega_{\alpha}\}_{\alpha>0}$ if necessary to account for different triangulations of the same domain) (2) it is shape regular, meaning that there exists $\beta>0$ such that for all elements $T\in\mathcal{T}_{h}$ and all $h>0$, $h_{T}/\rho_{T}\leq\beta$, where $h_{T}$ is the diameter of $T$ and $\rho_{T}$ is the diameter of the largest ball contained in $T$, and (3) for every $T\in\mathcal{T}_{h}$ such that $T\cap\Gamma_{h}\neq\varnothing$, the maximum distance between $\boldsymbol{x}\in T\cap\Gamma_{h}$ and $\boldsymbol{y}\in\Gamma$ is of the same order of magnitude as the element diameter $h_{T}$. More precisely, if $d_{loc}:=\max\{d(\boldsymbol{x},\boldsymbol{y}):\boldsymbol{x}\in T\cap\Gamma_{h}\text{ and }\boldsymbol{y}\in\Gamma\}$ then $d_{loc}=\mathcal{O}(h_{T})$. This last requirement, which will be referred to as the local proximity condition and is depicted schematically in Figure 1, is of key importance for the transfer process that will be defined later on.

For every element $T$ in a particular triangulation, we will denote by $\boldsymbol{n}_{T}$ the outward unit normal vector to $T$, or simply $\boldsymbol{n}$ instead of $\boldsymbol{n}_{T}$ whenever the context prevents any confusion. As it is conventional, we will denote the mesh parameter as $\displaystyle h:=\max_{T\in\mathcal{T}_{h}}h_{T}$, which will be assumed to be smaller than one for the sake of simplicity. We will denote by $e$ any face of a simplex and its length by $h_{e}$. Moreover, we will talk about an interior face $e$ if there are two elements $T^{+}$ and $T^{-}$ in the triangulation $\mathcal{T}_{h}$ such that $e=\partial T^{+}\cap\partial T^{-}$. The set of all interior faces will be denoted by $\mathcal{E}_{h}^{\circ}$. In a similar manner, we will talk about a boundary face $e$ if there is an element $T\in\mathcal{T}_{h}$ such that $e=\partial T\cap\Gamma_{h}$; the set of boundary faces will be denoted by $\mathcal{E}_{h}^{\partial}$. Note that with these definitions, the entirety of the faces of the triangulation denoted by $\mathcal{E}_{h}$ (often referred to as the skeleton of the mesh) can then be decomposed as $\mathcal{E}_{h}=\mathcal{E}_{h}^{\circ}\cup\mathcal{E}_{h}^{\partial}$.

Throughout the paper, we will be working with functions that in general will not be continuous across mesh elements. For a scalar-valued function we will use the symbol $[\![{w}]\!]:=w^{+}-w^{-}$ to refer to its jump across any given interior face. At the boundary faces, the jump will be defined as $[\![{w}]\!]:=w-\varphi_{h}$, where $\varphi_{h}$ is the approximation of the boundary data at $\Gamma_{h}$ that will be defined later. In the case of vector-valued functions $\boldsymbol{v}$, we will be interested in the discontinuity of its normal component across interior faces, which will be denoted by $[\![{\boldsymbol{v}}]\!]:=\boldsymbol{v}^{+}\cdot\boldsymbol{n}^{+}+\boldsymbol{v}^{-}\cdot\boldsymbol{n}^{-}$.

A remark on the local proximity condition and mesh refinement: The local proximity condition limits the minimum size that the elements near the boundary of an admissible triangulation can attain. Therefore, mesh refinement in this context must be understood as moving through a sequence of computational domains in the set $\{\Omega_{h}\}_{h>0}$ and their corresponding admissible triangulations in $\{\mathcal{T}_{h}\}_{h>0}$ as the parameter $h\to 0$. As it will be shown later, the error estimates will not depend on the particular domain $\Omega_{h}$ or triangulation $\mathcal{T}_{h}$ as long as the three requirements on the mesh stated above are satisfied. Possible ways of building sequences of admissible triangulations and computational domians have been detailed in [7] for uniform meshes and in [18] for adaptively refined triangulations.

Having defined the family of polygonal subdomains and admissible triangulations on which the discretization will be performend, we will now proceed to detail the process through which the boundary information will be transferred from the boundary into the computational domain. In order to do that we will have to tessellate the region enclosed between the two boundaries $\Gamma$ and $\Gamma_{h}$ as follows.

Given a triangulation $\mathcal{T}_{h}$ of the computational domain $\Omega_{h}$ and a boundary face $e\in\mathcal{E}_{h}^{\partial}$, we will denote by $T^{e}$ the unique element of $\mathcal{T}_{h}$ such that $e\cap\overline{T^{e}}=e$. To every point $\boldsymbol{x}\in e$, we will associate—in a smooth fashion—a point $\overline{\boldsymbol{x}}\in\Gamma$ and set $l(\boldsymbol{x})=|\boldsymbol{x}-\overline{\boldsymbol{x}}|$. We will define the extension patch $T^{e}_{ext}$ as

where $\boldsymbol{t}=\boldsymbol{t}(\boldsymbol{x})$ is the unit vector anchored at $\boldsymbol{x}$ and pointing in the direction of $\overline{\boldsymbol{x}}$. With this notation, the line segment connecting $\boldsymbol{x}$ to $\overline{\boldsymbol{x}}$ can be parameterized by

The point $\overline{\boldsymbol{x}}\in\Gamma$ and therefore the vector $\boldsymbol{t}(\boldsymbol{x})$ can be specified in several ways. Here, we will consider that the point has been determined in such a way that

This assumption is made with the sole purpose of making the analysis simpler, it can in fact be relaxed to the existence of a constant $a_{0}$ such that $0<a_{0}\leq\boldsymbol{t}(\boldsymbol{x})\cdot\boldsymbol{n}$ for every $\boldsymbol{x}$ belonging to a boundary edge. The numerical method described here is remarkably robust with respect to the method used to choose $\boldsymbol{t}(\boldsymbol{x})$. Previously, the direction had been determined using the algorithm proposed by [7], which assigns $\overline{\boldsymbol{x}}$ in such a way that the three following conditions are satisfied: (1) $\overline{\boldsymbol{x}}$ is unique, (2) any two different line segments $\sigma_{\boldsymbol{t}}$ do not intersect each other inside $T^{e}_{ext}$, and (3) the segments $\sigma_{\boldsymbol{t}}$ do not intersect the interior of $\Omega_{h}$. As it is proven in the aforementioned reference, these three conditions guarantee that the union of $T^{e}_{ext}$ completely covers $\Omega_{h}^{c}:=\Omega\setminus\overline{\Omega_{h}}$. An alternate method was used and tested numerically in [17, 18], where $\boldsymbol{t}(\boldsymbol{x})$ was determined using a weighted average of the normal vectors from neighboring boundary edges.

For an extended patch $T^{e}_{ext}$ and a mesh element $T^{e}\in\mathcal{T}_{h}$ sharing a boundary face $e\in\mathcal{E}^{\partial}_{h}$, we denote by $h_{e}^{\perp}$ (resp. $H_{e}^{\perp}$) the largest distance between a point inside $T_{e}$ (resp. $T_{ext}^{e}$) and the plane determined by the face $e$. The ratio between these two distances will be denoted by $r_{e}:=H_{e}^{\perp}/h_{e}^{\perp}$ and the maximum such ratio taken over all the boundary edges will be denoted by $\displaystyle R_{\mathcal{T}_{h}}:=\max_{e\in\mathcal{E}_{h}^{\partial}}r_{e}$. Note that the local proximity condition will ensure the existence of a constant R, that we will call the proximity constant, independent of the particular triangulation $\mathcal{T}_{h}$, such that

We will also define the class of non-trivial vector-valued polynomials of degree at most $k$ defined in and across both patches as

We can then introduce, for all those elements with a non-empty intersection with the computational boundary, the element-wise constants

where, in abuse of notation, $\boldsymbol{n}_{e}$ is a constant vector field defined in $T^{e}_{ext}\cup T^{e}\cup\,e$ that coincides with the unit exterior normal vector associated to the face $e$ and pointing in the direction of the extension patch. Above, the norms $\|\cdot\|_{T^{e}_{ext}}$ and $\|\cdot\|_{T^{e}}$ are the standard $L^{2}$ norms supported on the extension patch $T^{e}_{ext}$ and its neighboring element $T^{e}$ respectively. In [5], these constants were bounded in terms of the polynomial degree of the approximation, $k$, and the regularity constant of the mesh, $\beta$, as

where $C_{1}$ and $C_{2}$ depend only on the mesh regularity. As we will see below, these constants will determine the magnitude of the proximity constant and therefore the maximum admissible gap between the computational and physical boundaries.

We will also need to make two technical assumptions relating the proximity constant to the and the diffusion coefficient $\kappa$ and the degree of the polynomial approximation. For each $e\in\mathcal{E}_{h}^{\partial}$ we will require the following to hold:

where $\underline{\kappa}$ and $\overline{\kappa}$ are the lower and upper bounds of $\kappa$, resp., specified in (3.1), whereas $\overline{\tau}$ is the maximum of the stabilization parameter $\tau$ of the HDG scheme. The first of these two conditions, (2.4a), states the well known fact that for small values of the diffusivity, small scale behavior can be expected near the physical boundary, and therefore fine extension patches are required. However, it also provides the additional insight that the distance between the boundaries can be increased at the cost of accepting smaller values of the stabilization factor $\tau$—and hence larger discontinuities in the discrete solution. In a similar vein, (2.4b) relates the range of values of the diffusion coefficient with the maximum separation between the computational and physical boundaries, thus making sure that the external patches are fine enough to resolve possible boundary behavior induced by large variations in diffusivity over the domain. Moreover, it sets a hard upper limit to the mechanism that allows for a larger separation by decreasing $\tau$.

By combining (2.3) with (2.4) it is not hard to show that, for $k>0$, $H^{\perp}_{e}$ must be bounded as

This expression provides insight into the way in which the physics of the problem—through the range of values for $\kappa$—interacts with the discretization—through the parameters $H^{\perp}_{e}$, $h^{\perp}_{e}$, $k$, $\beta$, and $\tau$—and determines the maximum separation between the physical boundary and that of an admissible triangulation. Of particular note is the role played by the polynomial degree of the approximation: for larger values of $k$ the distance between the mesh and the boundary must decrease. The reason for this will become apparent soon, as we will resort to extrapolation to approximate some quantities over the extension patches.

Having established the requirements for an admissible triangulation, and defined the extension patches $T^{e}_{ext}$ in such a way that for each of them there corresponds a single $T^{e}\in\mathcal{T}_{h}$, we can define a way to extend polynomial functions from the computational domain into $\Omega_{h}^{c}$. This extension process will enable us to transfer the boundary condition from $\Gamma$ into the computational boundary $\Gamma_{h}$. Let $p:T^{e}\to\mathbb{R}$ be a polynomial function and $T^{e}_{ext}$ an extension patch associated to $T^{e}$. We will define the extension $\boldsymbol{E}_{h}(p)$ of $p$ to $T^{e}_{ext}$ by extrapolation as follows:

Where, to keep notation simple, a polynomial function $p$ should be understood as its extrapolation $E_{h}(p)$ whenever an evaluation outside of $\Omega_{h}$ is required, which should be clear from the context. For vector-valued polynomial functions, the extension is defined similarly component by component.

To avoid computing on a domain with curved boundary, we wish to pose the boundary vaue problem (1.1) polygonal subdomain $\Omega_{h}\subset\Omega$. This simpler geometry can then be discretized by a uniform triangulation of size satisfying the conditions outlined in sections 2.1 and 2.2.

Recasting (1.1) in mixed form and restricting the resulting equivalent first order system to $\Omega_{h}$ leads to

where the specific relation between $\kappa$, $u$ and $\nabla\,u$ has not been made explicit, and the—a priori unknown—function $\varphi$ encodes the restriction of $u$ to the computational boundary $\Gamma_{h}$. We can recover $\varphi$ following the method proposed by [4] (in one dimension) and extended to higher dimensions by [7]. The idea consists of transferring the Dirichlet data $g$ from $\Gamma$ to $\Gamma_{h}$ along segments called transfer paths by computing a line integral of the flux $\boldsymbol{q}$.

To be precise, given $\boldsymbol{x}\in\Gamma_{h}$ and $\overline{\boldsymbol{x}}\in\Gamma$, equation (2.5a) can be integrated along the segment connecting them. Lets denote by $\boldsymbol{t}(\boldsymbol{x})$ the unit vector anchored at $\boldsymbol{x}$ pointing towards $\overline{\boldsymbol{x}}$, and by $l(\boldsymbol{x})$ the length of the segment connecting them. We then have the following representation for $\varphi$:

Note that $\varphi$ depends on the values of either $u$ or $\nabla\,u$ (through $\kappa^{-1}$), and $\boldsymbol{q}$ over the extended domain $\Omega_{h}^{c}$. As such, we should write $\varphi=\varphi(u,\nabla\,u,\boldsymbol{q},\boldsymbol{x})$ however, to keep notation simple, we will abstain from this and will write simply $\varphi$. In a similar fashion, $g(\overline{\boldsymbol{x}})$ is in fact a function of $\boldsymbol{x}$, since the point $\overline{\boldsymbol{x}}$ varies smoothly with $\boldsymbol{x}$. To avoid the use of cumbersome notation we will write either $g(\overline{\boldsymbol{x}})$ or simply $\overline{g}:=g(\overline{\boldsymbol{x}}(\boldsymbol{x}))$.

To denote spaces of functions we will make use of the standard notation and terminology from Sobolev space theory. Let $\mathcal{O}$ be a domain in $\mathbb{R}^{d}$, and $\Sigma$ be either a Lipschitz curve (if $d=2$) or surface (if $d=3$); for scalar-valued functions and non zero real numbers $s$, we will use the spaces $H^{s}(\mathcal{O})$ and $H^{s}(\Sigma)$ with their usual definition, whereas for the case $s=0$ we will write simply $L^{2}(\mathcal{O})$ and $L^{2}(\Sigma)$. The spaces of vector-valued functions will be denoted in bold face, therefore $\boldsymbol{H}^{s}(\mathcal{O}):=[H^{s}(\mathcal{O})]^{d}$ and $\boldsymbol{H}^{s}(\Sigma):=[H^{s}(\Sigma)]^{d}$.

The $L^{2}$ inner products for both scalar and vector-valued functions on volumes and surfaces will be denoted by $(\cdot,\cdot)_{\mathcal{O}}$ and $\left\langle\cdot,\cdot\right\rangle_{\Sigma}$ respectively. The associated norms will be denoted by $\|\cdot\|_{s,\mathcal{O}}$ and $\|\cdot\|_{s,\Sigma}$ and simply $\|\cdot\|_{\mathcal{O}}$ for the case $s=0$. As is common, will we write $|\cdot|_{s,\mathcal{O}}$ for the $\boldsymbol{H}^{s}$ and $H^{s}$-semi norms.

Given a triangulation $\mathcal{T}_{h}$ we will define the following mesh-dependent inner products over elements and edges

These inner products induce mesh-dependent norms that will be denoted, respectively, by

In the forthcoming analysis, the expression $a\lesssim b$ should be understood as meaning $a\leq Cb$ where $C$ is a positive constant independent of $h$.

For the discrete formulations that will be introduced in the next sections, we will make use of the following finite dimensional spaces of piece-wise polynomial functions

where, $\mathds{P}_{k}(T)$ denotes the space of polynomials of degree at most $k$ defined in $T\in\mathcal{T}_{h}$. Similarly, $\mathds{P}_{k}(e)$ denotes the space of polynomials of degree at most $k$ defined over a face $e\in\mathcal{E}_{h}$.

We will first consider the case when the coefficient $\kappa$ depends on the solution in the form

For the analysis, we will require the existence of positive constants $\underline{\kappa}$ and $\overline{\kappa}$ such that for all $u\in L^{2}(\Omega)$

Moreover, $\kappa$ will be assumed to be Lipschitz-continuous on $L^{2}(\Omega)$, i.e, there exists $\tilde{L}>0$ such that

The two conditions above, together, imply the existence of constants $\widehat{L}$ and $L$ such that

Note that all these assumptions imply the Lipschitz continuity of $\kappa,\kappa^{1/2}$, and $\kappa^{-1}$ on the subdomain $\Omega_{h}\subset\Omega$ with corresponding Lipschitz constants equal to or smaller than those stated above.

Before introducing the discrete formulation we will recall here the mixed form (2.5), but now we make explicit the dependence $\kappa=\kappa(u)$

The boundary data $\varphi$ on the computational boundary $\Gamma_{h}$ is transferred according to (2.6).

Taking an admissible triangulation $\mathcal{T}_{h}$ of the computational domain $\Omega_{h}$, the HDG discretization of (3.5) reads: Find $(\boldsymbol{q}_{h},u_{h},\widehat{u}_{h})\in\boldsymbol{V}_{h}\times W_{h}\times M_{h}$, such that

In the definition above, we have used the extrapolation $E_{h}\boldsymbol{q}_{h}$ due to the fact that the approximation $\boldsymbol{q}_{h}$ is available only inside of the computational domain $\Omega_{h}$, but the transfer paths along which the integral is computed are defined over the complementary extended region $\Omega_{h}^{c}$.

In this section we employ a Banach fixed-point argument to ensure the well-posedness of the discrete problem (3.6). To that end we will define an operator $\mathcal{J}:W_{h}\to W_{h}$ mapping $\zeta$ to the second component of the triplet $(\boldsymbol{q},u,\widehat{u})\in\boldsymbol{V}_{h}\times W_{h}\times M_{h}$ satisfying, for all $(\boldsymbol{v},w,\mu)\in\boldsymbol{V}_{h}\times W_{h}\times M_{h}$, the HDG system (3.6) where the source has been evaluated at $\zeta$, namely

Above, the term $\varphi(\zeta)$ corresponds to the boundary condition transferred to the computational domain by means of (3.6e). The mapping $\mathcal{J}$ is well defined, as the linearized system (3.7) is uniquely solvable as proven in [5].

The main result of this section—that the mapping $\mathcal{J}$ defined above is a contraction—relies on the validity of a particular inequality—estimate (3.8) below—but is otherwise a simple argument. Since the proof of (3.8) requires a sequence of technical arguments, in the interest of clarity (we will first prove the main theorem assuming that the aforementioned inequality is valid. After having established the well posedness of the discrete problem, the reminder of the section will be devoted to verifying the validity of (3.8). This will be finally established in Lemma 4, after a series of auxiliary results.

Combined with a standard fixed-point argument, the theorem above guarantees the existence and uniqueness of the solution to the discrete HDG system. We will now direct our efforts to showing that inequality (3.8) holds. To that end we will make use of the following auxiliary function and its properties listed in the lemma below—the proof of which can be found on [5, Lemma 5.2].

For the subsequent analysis, we will also make use of the following norm. Given$(\boldsymbol{v},w,\mu)\in\boldsymbol{V}_{h}\times W_{h}\times M_{h}$ and $\xi\in M_{h}$ we define

The proof of (3.8) requires considering the solution $\phi$ to an auxiliary problem (c.f. (B.1) ) and using a duality argument that connects a stability estimate for $\phi$ with our variables of interest. This will be done in Lemma 4. Lemmas 2 and 3 below establish estimates relating the norm (3.11) of the solution $(\boldsymbol{q},u,\widehat{u})$ to problem data that will be used in the final step of Lemma 4.