Least-Squares versus Partial Least-Squares Finite Element Methods: Robust A Priori and A Posteriori Error Estimates of Augmented Mixed Finite Element Methods

We consider the following elliptic equation with possible discontinuous coefficients (a generalized Darcy’s problem):

(1.1)

$$\left\{\begin{array}[]{lllll}\nabla\cdot\mbox{\boldmath$\sigma$}&=&g,&\mbox{in }\Omega\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{the constitutive equation},\\[2.84526pt] A\nabla u+\mbox{\boldmath$\sigma$}&=&A{\bf f},&\mbox{in }\Omega\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{the equilibrium equation},\\[2.84526pt] u&=&0&\mbox{on }\Gamma_{D}\\ \mbox{\boldmath$\sigma$}\cdot{\bf n}&=&0&\mbox{on }\Gamma_{N},\end{array}\right.$$

The domain $\Omega$ is a bounded, open, connected subset of $\mathbb{R}^{d}(d=2\mbox{ or }3)$ with a Lipschitz continuous boundary $\partial\Omega$. We partition $\partial\Omega$ into two open subsets $\Gamma_{D}$ and $\Gamma_{N}$, such that $\partial\Omega=\overline{\Gamma_{D}}\cup\overline{\Gamma_{N}}$ and $\Gamma_{D}\cap\Gamma_{N}=\emptyset$. For simplicity, we assume that $\Gamma_{D}$ is not empty (i.e., $\mbox{meas}(\Gamma_{D})\neq 0$ ) and is connected. We assume that the diffusion coefficient matrix $A\in L^{\infty}(\Omega)^{d\times d}$ is a given $d\times d$ tensor-valued function; the matrix $A$ is uniformly symmetric positive definite: there exist positive constants $0<\Lambda_{0}\leq\Lambda_{1}$ such that

(1.2)

$$\Lambda_{0}{\bf y}^{T}{\bf y}\leq{\bf y}^{T}A{\bf y}\leq\Lambda_{1}{\bf y}^{T}{\bf y}$$

for all ${\bf y}\in\mathbb{R}^{d}$ and almost all $x\in\Omega$. The righthand sides $g\in L^{2}(\Omega)$ and ${\bf f}\in L^{2}(\Omega)^{d}$.

There are several variational formulations for (1.1). The standard conforming finite element method tries to approximate $u$ in the finite element subspace of its natural space $H^{1}(\Omega)$, see [27, 7, 9]. If a good approximation of $\sigma$ in its natural space $H({\rm div})$ is sought, then one can use the mixed finite element approximation based on a dual mixed formulation, see for example, [6]. In order to approximate both $u$ and $\sigma$ in their intrinsic spaces $H^{1}(\Omega)$ and $H({\rm div})$, respectively, one natural choice is the least-squares finite element method (LSFEM).

The least-squares variational principle and the corresponding LSFEM based on a first-order system reformulation have been widely used in numerical solutions of partial differential equations; see, for example [17, 19, 36, 5, 20, 18, 13, 39, 38, 42]. Compared to the standard variational formulation and the related finite element methods, the first-order system LSFEMs have several known advantages, such as the discrete problem is stable without the inf-sup condition of the discrete spaces and mesh size restriction; the first-order system is physically meaningful, thus important physical qualities such as displacement, flux, and stress can be approximated in their intrinsic spaces; and the least-squares functional itself is a good built-in a posteriori error estimator. In Section 4 of [46], as an example of an elliptic equation with an $H^{-1}$-righthand side, the first-order system LSFEM is studied for (1.1).

On the other hand, when applying LSFEMs to (1.1), we also face an important issue: robustness. For a coefficient-dependent problem ($A$ for (1.1)), the robustness of both the a priori and a posteriori estimates, i.e., genetic constants that appeared in the estimates are independent of the coefficients, is of crucial importance. For the a priori error estimate, it is well known that the model problem (1.1) may only have $H^{1+s}$ regularity, with possibly very small $s>0$, see for example, Kellogg [37]. In [14], it is noted that the a priori analysis using the global regularity constant $s$ is not satisfactory since most of the time, the solution only has a low regularity at the elements attached to the singularity but can be very smooth away from the singularity. Thus, we should study the a priori analysis under the local regularity assumption. The a priori error estimate using local regularity is also the base for adaptive finite element methods to achieve an equal-discretization error distribution. For the a posteriori error estimate, we want the a posteriori error estimator with the efficiency constant and the reliability constant to be independent of the parameter of the problem. Unfortunately, both a priori and a posteriori error estimates of the LSFEM applying to (1.1) with discontinuous coefficients are not robust; see our discussion in Sections 5.2 and 6.3. Beside the above mentioned robustness issue, we often need extra regularity for the right-hand side $g$ in the standard $L^{2}$-based LSFEM since all the errors of the LSFEM are measured in a combined norm. On the other hand, the standard mixed finite element method for (1.1) does not require this, see discussions in [45].

Besides the full bonafide LSFEM, another idea to apply the least-squares philosophy is the so-called Galerkin-Least-Squares (GaLS) method. The GaLS method is a method combining the least-squares and Galerkin methods. Some least-squares terms are added to the original variational formulation to enhance the stability in GaLS methods. We consider a special GaLS method, the augmented mixed finite element method, in this paper. The central idea of the augmented mixed method is adding consistent least-squares terms to the original mixed formulation to guarantee coercivity or stability. The first augmented mixed finite element method is introduced in [40]. Earlier contributions of GaLS methods can be found in [33, 34]. The group of Gatica made many important and seminal contributions on developing the augmented mixed finite element method to many problems, see for example [35, 3, 32, 25, 28, 1].

As a partial least-squares method, the augmented mixed finite element method shares many properties with the classic LSFEM. It is also stable and approximates the physical quantities in their native spaces. Moreover, as we will see later in the paper, the a posteriori error estimator of the augmented mixed finite element method is a least-squares error estimator. Furthermore, since we only partially use the least-squares principle, the system is more flexible and we can show the robustness of both a priori and a posteriori error estimates.

This paper proposes two versions of robust augmented mixed finite element methods with two choices of least-squares terms based on the constitutive equation. The first is a simple $L^{2}$ least-squares term and the second is a mesh-weighted least-squares term. We show the robust and local optimal a priori error estimates for both methods. For the mesh-weighted augmented mixed finite element method, we also show that optimal error can be achieved without requiring high regularity of the righthand side. For both methods, we derive robust a posteriori error estimates. In fact, the error estimators of both versions of augmented mixed finite element methods are least-squares error estimators with corresponding least-squares functionals. As comparisons, we discuss the robustness of two corresponding LSFEMs. For the $L^{2}$-based LSFEM, we show that the a priori and a posteriori error estimates are robust only under a very special condition. In general, they are not robust. For the mesh-weighted LSFEM, the results are much worse than those of augmented mixed finite element methods since its coercivity constant depends on the mesh size.

The robust (but not local optimal) a priori and robust residual type a posteriori error estimate for the energy norm was obtained for the conforming FEM [4]. Robust and local optimal a priori error estimates have also been derived for mixed FEM in [45] and nonconforming FEM and discontinuous Galerkin FEM in [14] without a restrictive assumption on the distribution of the coefficients. We also present a detailed discussion of the robust and local optimal a priori error estimate for the conforming FEM in [14]. For recovery-based error estimators, robust a posteriori error estimates are obtained by us in [22, 23, 21]. Robust equilibrated error estimators were developed by us in [24, 11, 12]. Robust residual-type of error estimates without a restrictive assumption on the distribution of the coefficients is developed for nonconforming and DG approximations in [14].

In the original paper [40], both $u$ and $\sigma$ are approximated by continuous finite elements. As mentioned in [22, 14], the flux $\sigma$ is not continuous for discontinuous coefficients. Thus the continuous finite element is not a good candidate for the approximation. In [10], a mixed discontinuous Galerkin method is developed using the stabilization in [40]. In [29], several different stabilization formulations are proposed, with one formulation being our first augmented mixed method, although the authors still emphasized using standard conforming finite elements. In [2], $H({\rm div})$ and $H^{1}$-conforming finite elements are used. In [2], an a posteriori error estimator is proposed for the first augmented mixed formulation. The error estimator in [2] is actually a least-squares error estimator. The robustness of a priori and a posteriori error estimates are not discussed in all previous papers.

The paper is organized as follows. Section 2 describes notations, the function spaces, and local interpolation results. The generalized Darcy problem and the augmented formulations are presented in Section 3, including the robust Cea’s lemma. In Section 4, we discuss simplified assumptions on the coefficient $A$ and the quasi-monotonicity assumption and robust quasi-interpolation based on this assumption. In Sections 5 and 6, we present robust a priori and a posteriori error analyses for the first and second augmented mixed formulations, respectively. Connections and comparisons with the corresponding LSFEMs are also discussed in Sections 5 and 6. Numerical experiments are presented in Section 7 to verify the findings in the paper. We make some concluding remarks in Section 8.

We use the standard notations and definitions for the Sobolev spaces $H^{s}(\Omega)$ for $s\geq 0$. The standard associated inner product is denoted by $(\cdot,\,\cdot)_{s,\Omega}$, and its norm is denoted by $\|\cdot\|_{s,\Omega}$. The notation $|\cdot|_{s,\Omega}$ is used for the semi-norm. (We suppress the superscript $d$ because the dependence on dimension will be clear by context. We also omit the subscript $\Omega$ from the inner product and norm designation when there is no risk of confusion.) For $s=0$, $H^{s}(\Omega)$ coincides with $L^{2}(\Omega)$. The symbols $\nabla\cdot$ and $\nabla$ stand for the divergence and gradient operators, respectively. Set $H^{1}_{D}(\Omega):=\{v\in H^{1}(\Omega)\,:\,v=0\,\,\mbox{on }\Gamma_{D}\}$. We use the standard $H({\rm div};\Omega)$ space equipped with the norm $\|\mbox{\boldmath$\tau$}\|_{H({\rm div};\,\Omega)}=\left(\|\mbox{\boldmath$\tau$}\|^{2}_{0,\Omega}+\|\nabla\cdot\mbox{\boldmath$\tau$}\|^{2}_{0,\Omega}\right)^{\frac{1}{2}}.$ Denote its subspace by $H_{N}({\rm div};\Omega)=\{\mbox{\boldmath$\tau$}\in H({\rm div};\Omega)\,:\,\mbox{\boldmath$\tau$}\cdot{\bf n}|_{\Gamma_{N}}=0\},$ where ${\bf n}$ is the unit vectors outward normal to the boundary $\partial\Omega$. For simplicity, we use the following notation:

(2.1)

$${\bf X}:=H_{N}({\rm div};\,\Omega)\times H^{1}_{D}(\Omega).$$

Let ${\mathcal{T}}=\{K\}$ be a triangulation of $\Omega$ using simplicial elements. The mesh ${\mathcal{T}}$ is assumed to be regular. Let $P_{k}(K)$ for $K\in{\mathcal{T}}$ be the space of polynomials of degree $k$ on an element $K$. Denote the standard linear and quadratic conforming finite element spaces by

$$S_{1,D}=\{v\in H_{D}^{1}(\Omega):v|_{K}\in P_{1}(K),\leavevmode\nobreak\ \forall\leavevmode\nobreak\ K\in{\mathcal{T}}\}\quad\mbox{and}\quad S_{2,D}=\{v\in H_{D}^{1}(\Omega):v|_{T}\in P_{2}(K),\leavevmode\nobreak\ \forall\leavevmode\nobreak\ K\in{\mathcal{T}}\},$$

respectively.

Denote the local lowest-order Raviart-Thomas (RT) [43] on element $K\in{\mathcal{T}}$ by $RT_{0}(K)=P_{0}(K)^{d}+{\bf x}\,P_{0}(K)$. Then the lowest-order $H({\rm div};\,\Omega)$ conforming RT space with zero trace on $\Gamma_{N}$ is defined by

$$RT_{0,N}=\{\mbox{\boldmath$\tau$}\in H_{N}({\rm div};\Omega):\mbox{\boldmath$\tau$}|_{K}\in RT_{0}(K)\,\,\,\,\forall\,\,K\in{\mathcal{T}}\}.$$

Similarly, the lowest-order $H({\rm div};\,\Omega)$-conforming Brezzi-Douglas-Marini (BDM) space with zero trace on $\Gamma_{N}$ is defined by

$$BDM_{1,N}=\{\mbox{\boldmath$\tau$}\in H_{N}({\rm div};\Omega):\mbox{\boldmath$\tau$}|_{K}\in P_{1}(K)^{d},\leavevmode\nobreak\ \forall\leavevmode\nobreak\ K\in{\mathcal{T}}\}.$$

We discuss some local approximation results of the standard $S_{1,D}$ and $RT_{0,N}$ spaces. By Sobolev’s embedding theorem, $H^{1+s}(\Omega)$, with $s>0$ for two dimensions and $s>1/2$ for three dimensions, is embedded in $C^{0}(\Omega)$. Thus, we can define the nodal interpolation $I^{nodal}_{h}$ of a function $v\in H^{1+s}(\Omega)$ with $I^{nodal}_{h}v\in S_{1,D}$ and $I^{nodal}_{h}v(z)=v(z)$ for a vertex $z\in{\mathcal{N}}$. It is important to notice that the nodal interpolation is completely element-wisely defined. We have the following local interpolation estimate for the linear nodal interpolation $I^{nodal}_{h}$ with local regularity $0<s_{K}\leq 1$ in two dimensions and $1/2<s_{K}\leq 1$ in there dimensions [30, 14]:

(2.2)

$$\|\nabla(v-I^{nodal}_{h}v)\|_{0,K}\leq Ch_{K}^{s_{K}}|\nabla v|_{s_{K},K}.$$

Assume that $\mbox{\boldmath$\tau$}\in L^{r}(\Omega)^{d}\cap H({\rm div};\Omega)$, and locally $\mbox{\boldmath$\tau$}\in H^{s_{K}}(K)$ with the local regularity $1/2<s_{K}\leq 1$. Let $I^{rt}_{h}$ be the canonical RT interpolation from $L^{r}(\Omega)^{d}\cap H_{N}({\rm div};\Omega)$ to $RT_{0,N}$. Then the following local interpolation estimates hold for local regularity $1/2<s_{K}\leq 1$ with the constant $C_{rt}$ being unbounded as $s_{K}\downarrow 1/2$ (see Chapter 16 of [31]):

(2.3)

$$\|\mbox{\boldmath$\tau$}-I^{rt}_{h}\mbox{\boldmath$\tau$}\|_{0,K}\leq C_{rt}h_{K}^{s_{K}}|\mbox{\boldmath$\tau$}|_{s_{K},K}\quad\forall K\in{\mathcal{T}}.$$

Due to the commutative property of the standard RT interpolation, if we further assuming that $\nabla\cdot\mbox{\boldmath$\tau$}|_{K}\in H^{t_{K}}(K)$, $0<t_{K}\leq 1$, then

(2.4)

$$\|\nabla\cdot(\mbox{\boldmath$\tau$}-I^{rt}_{h}\mbox{\boldmath$\tau$})\|_{0,K}\leq Ch_{K}^{t_{K}}|\nabla\cdot\mbox{\boldmath$\tau$}|_{t_{K},K}\quad\forall K\in{\mathcal{T}}.$$

For $v\in H^{3}(\Omega)$, the following local interpolation result in standard for nodal interpolation $I_{h}$ in $S_{2,D}$,

(2.5)

$$\|\nabla(v-I_{h}v)\|_{0,K}\leq Ch_{K}^{2}|v|_{3,K}.$$

Also, we have the following standard interpolation result for the $BDM_{1,N}$, assuming that $I^{bdm}_{h}$ is the standard $BDM_{1}$ interpolation,

(2.6)

$$\|\mbox{\boldmath$\tau$}-I^{bdm}_{h}\mbox{\boldmath$\tau$}\|_{0,K}\leq C_{bdm}h_{K}^{2}|\mbox{\boldmath$\tau$}|_{2,K}\quad\forall K\in{\mathcal{T}}.$$

If we further assuming that $\nabla\cdot\mbox{\boldmath$\tau$}|_{K}\in H^{t_{K}}(K)$, $0<t_{K}\leq 1$, then

(2.7)

$$\|\nabla\cdot(\mbox{\boldmath$\tau$}-I^{bdm}_{h}\mbox{\boldmath$\tau$})\|_{0,K}\leq Ch_{K}^{t_{K}}|\nabla\cdot\mbox{\boldmath$\tau$}|_{t_{K},K}\quad\forall K\in{\mathcal{T}}.$$

The following mesh-dependent notation is used in the paper: for $r\geq 0$,

(2.8)

$$\|h^{r}v\|_{0}=\left(\sum_{K\in{\mathcal{T}}}h_{K}^{2r}\|v\|_{0,K}^{2}\right)^{1/2}\quad\mbox{and}\quad(h^{r}v,w)=\sum_{K\in{\mathcal{T}}}h_{K}^{r}(v,w)_{K},\quad\forall v,w\in L^{2}(\Omega).$$

For the solution $u\in H_{D}^{1}(\Omega)$, the flux $\mbox{\boldmath$\sigma$}=A{\bf f}-A\nabla u$ is in $L^{2}(\Omega)^{d}$, and $\nabla\cdot\mbox{\boldmath$\sigma$}=g$ is in $L^{2}(\Omega)$. Thus, $\mbox{\boldmath$\sigma$}\in H_{N}({\rm div}\;\Omega)$. Multiplying an arbitrary $v\in H_{D}^{1}(\Omega)$ on both sides of the first equation of (1.1), taking integration by parts, and replacing $\sigma$ by $A{\bf f}-A\nabla u$, we have the standard variational problem of (1.1): find $u\in H^{1}_{D}(\Omega)$, such that

(3.1)

$$(A\nabla u,\nabla v)=(g,v)+(A{\bf f},\nabla v)\quad\forall v\in H_{D}^{1}(\Omega).$$

We have two mixed formulations.

In the augmented mixed method, we want to approximate $u$ in its natural space $H_{D}^{1}(\Omega)$ and keep the method stable without restricting the choice of finite element subspaces. Testing $\mbox{\boldmath$\tau$}\in H_{N}({\rm div};\Omega)$ for the first equation of (1.1) and $v\in H_{D}^{1}(\Omega)$ for the second equation of (1.1), we have

(3.4)

$$(A^{-1}\mbox{\boldmath$\sigma$},\mbox{\boldmath$\tau$})+(\nabla u,\mbox{\boldmath$\tau$})+(\nabla\cdot\mbox{\boldmath$\sigma$},v)=({\bf f},\mbox{\boldmath$\tau$})+(g,v),\leavevmode\nobreak\ \forall(\mbox{\boldmath$\tau$},v)\in{\bf X}.$$

We add two consistent least-squares terms to (3.4):

(3.5)

$$-\kappa(A^{-1}\mbox{\boldmath$\sigma$}+\nabla u,\mbox{\boldmath$\tau$}-A\nabla v)=-\kappa({\bf f},\mbox{\boldmath$\tau$}-A\nabla v)\leavevmode\nobreak\ \mbox{ from the constitutive equation},$$

and

(3.6)

$$\mu(\alpha^{-1}\nabla\cdot\mbox{\boldmath$\sigma$},\nabla\cdot\mbox{\boldmath$\tau$})=\mu(\alpha^{-1}g,\nabla\cdot\mbox{\boldmath$\tau$})\leavevmode\nobreak\ \mbox{ from the equilibrium equation},$$

where

(3.7)

$$\alpha(x)={\rm trace}(A(x))/d.$$

Note that ${\rm trace}(A)$ equals to the sum of its eigenvalues. Since it is assumed that $A$ is symmetric and positive definite, the scalar function $\alpha(x)$ is between the minimum and maximum of eigenvalues of $A(x)$ for all $x\in\Omega$.

Then we obtain the following problem: find $(\mbox{\boldmath$\sigma$},u)\in{\bf X}$, such that

(3.8)

$$B((\mbox{\boldmath$\sigma$},u),(\mbox{\boldmath$\tau$},v))=({\bf f},\mbox{\boldmath$\tau$})+(g,v)-\kappa({\bf f},\mbox{\boldmath$\tau$}-A\nabla v)+\mu(\alpha^{-1}g,\nabla\cdot\mbox{\boldmath$\tau$})\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf X}$$

with the bilinear form $B$ defined as follows, for $(\mbox{\boldmath$\chi$},w)\in{\bf X}$ and $(\mbox{\boldmath$\tau$},v)\in{\bf X}$,

$$B((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v)):=(A^{-1}\mbox{\boldmath$\chi$},\mbox{\boldmath$\tau$})+(\nabla w,\mbox{\boldmath$\tau$})+(\nabla\cdot\mbox{\boldmath$\chi$},v)-\kappa(A^{-1}\mbox{\boldmath$\chi$}+\nabla w,\mbox{\boldmath$\tau$}-A\nabla v)+\mu(\alpha^{-1}\nabla\cdot\mbox{\boldmath$\chi$},\nabla\cdot\mbox{\boldmath$\tau$}).$$

Let $(\mbox{\boldmath$\sigma$},u)=(\mbox{\boldmath$\tau$},v)$ in $B((\mbox{\boldmath$\sigma$},u),(\mbox{\boldmath$\tau$},v))$, and use the fact that

(3.9)

$$(\nabla v,\mbox{\boldmath$\tau$})+(\nabla\cdot\mbox{\boldmath$\tau$},v)=0\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf X},$$

We get

$$B((\mbox{\boldmath$\tau$},v),(\mbox{\boldmath$\tau$},v))=(1-\kappa)\|A^{-1/2}\mbox{\boldmath$\tau$}\|^{2}_{0}+\mu\|\alpha^{-1/2}\nabla\cdot\mbox{\boldmath$\tau$}\|^{2}_{0}+\kappa\|A^{1/2}\nabla u\|_{0}^{2}.$$

To have the coercivity, $1-\kappa$ should be positive. For convenience, we let $\kappa=1/2$. Then by (3.9), (3.8) can be written as: find $(\mbox{\boldmath$\sigma$},u)\in{\bf X}$, such that

(3.10)

$$B_{\theta}((\mbox{\boldmath$\sigma$},u),(\mbox{\boldmath$\tau$},v))=F_{\theta}(\mbox{\boldmath$\tau$},v)\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf X},$$

where, for $(\mbox{\boldmath$\chi$},w)\in{\bf X}$ and $(\mbox{\boldmath$\tau$},v)\in{\bf X}$, the bilinear form $B_{\theta}$ and the linear form $F_{\theta}$ are defined as follows:

(3.11)		$\displaystyle B_{\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))$	$\displaystyle=$	$\displaystyle(A^{-1}\mbox{\boldmath$\chi$},\mbox{\boldmath$\tau$})+(A\nabla w,\nabla v)+(\nabla w,\mbox{\boldmath$\tau$})-(\mbox{\boldmath$\chi$},\nabla v)+(\theta\alpha^{-1}\nabla\cdot\mbox{\boldmath$\chi$},\nabla\cdot\mbox{\boldmath$\tau$})$
(3.12)			$\displaystyle=$	$\displaystyle(A^{-1}\mbox{\boldmath$\chi$}+\nabla w,\mbox{\boldmath$\tau$}+A\nabla v)-2(\mbox{\boldmath$\chi$},\nabla v)+(\theta\alpha^{-1}\nabla\cdot\mbox{\boldmath$\chi$},\nabla\cdot\mbox{\boldmath$\tau$}),$
(3.13)		$\displaystyle F_{\theta}(\mbox{\boldmath$\tau$},v)$	$\displaystyle=$	$\displaystyle({\bf f},\mbox{\boldmath$\tau$}+A\nabla v)+2(g,v)+(\theta\alpha^{-1}g,\nabla\cdot\mbox{\boldmath$\tau$}).$

Note that $B_{\theta}$ is not symmetric. We will give an equivalent symmetric version in Section 3.2. We will discuss two choices of $\theta$, $\theta=1$ and a mesh dependent $\theta$, such that $\theta|_{K}=h^{2}_{K}$ for all $K\in{\mathcal{T}}$.

The formulas (3.11) and (3.12) are two equivalent ways to write the bilinear form.

In the remaining part of the paper, for simplicity of the presentation, we assume the following assumption on $A$: