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

Yuxiang Liang and Shun Zhang  Department of Mathematics, City University of
Hong Kong, Kowloon Tong, Hong Kong, China<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Date: August 28, 2024)

###### Abstract.

In this paper, for the generalized Darcy problem (an elliptic equation with
discontinuous coefficients), we study a special partial Least-Squares
(Galerkin-least-squares) method, known as the augmented mixed finite element
method, and its relationship to the standard least-squares finite element
method (LSFEM). Two versions of augmented mixed finite element methods are
proposed in the paper with robust a priori and a posteriori error estimates.
Augmented mixed finite element methods and the standard LSFEM uses the same a
posteriori error estimator: the evaluations of numerical solutions at the
corresponding least-squares functionals. As partial least-squares methods, the
augmented mixed finite element methods are more flexible than the original
LSFEMs. As comparisons, we discuss the mild non-robustness of a priori and a
posteriori error estimates of the original LSFEMs. A special case that the
$L^{2}$-based LSFEM is robust is also presented for the first time. Extensive
numerical experiments are presented to verify our findings.

###### Key words and phrases:

augmented mixed finite element method; least-squares finite element method;
Galerkin-Least-Squares method; robust a priori and a posteriori analysis

This work was supported in part by Research Grants Council of the Hong Kong
SAR, China, under the GRF Grant Project No. CityU 11316222

##  1\. Introduction

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.

##  2\. Preliminaries

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).$

##  3\. The generalized Darcy problem and the augmented mixed formulations

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.

Dual mixed method: Find $(\mbox{\boldmath$\sigma$},u)\in H_{N}({\rm
div},\Omega)\times L^{2}(\Omega)$ such that

(3.2)
$\begin{split}(A^{-1}\mbox{\boldmath$\sigma$},\mbox{\boldmath$\tau$})-(u,\nabla\cdot\mbox{\boldmath$\tau$})&=({\bf
f},\tau),\quad\forall\mbox{\boldmath$\tau$}\in H_{N}({\rm div};\Omega),\\\
(\nabla\cdot\mbox{\boldmath$\sigma$},v)&=(g,v),\quad\forall v\in
L^{2}(\Omega).\end{split}$

Primal mixed method: Find $(\mbox{\boldmath$\sigma$},u)\in
L^{2}(\Omega)^{d}\times H_{D}^{1}(\Omega)$ such that

(3.3)
$\begin{split}(A^{-1}\mbox{\boldmath$\sigma$},\mbox{\boldmath$\tau$})+(\nabla
u,\mbox{\boldmath$\tau$})&=({\bf
f},\mbox{\boldmath$\tau$}),\quad\forall\mbox{\boldmath$\tau$}\in
L^{2}(\Omega)^{d}\\\ -(\mbox{\boldmath$\sigma$},\nabla v)&=(g,v),\quad\forall
v\in H_{D}^{1}(\Omega).\end{split}$

The finite element approximations of the mixed methods require inf-sup stable
finite element pairs.

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}}$.

######  Remark 3.1.

In fact, we have a third case that $\theta=0$. However, unlike the first two
cases, this case has no corresponding least-squares formulation. We also
cannot associate its a posteriori error estimator with a corresponding least-
squares error estimator. Thus, we will not discuss this choice in the paper.
Some discussions of the case $\theta=0$ can be found in [40, 29].

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

###  3.1. Some analysis for the augmented mixed formulations

Define

(3.14)
$|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{\theta}:=(\|A^{1/2}\nabla
v\|^{2}_{0}+\|A^{-1/2}\mbox{\boldmath$\tau$}\|^{2}_{0}+\|\sqrt{\theta/\alpha}\nabla\cdot\mbox{\boldmath$\tau$}\|^{2}_{0})^{1/2}\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

By the definition of the bilinear forms $B_{\theta}$ in (3.11), we immediately
have the coercivity:

(3.15)
$B_{\theta}((\mbox{\boldmath$\tau$},v),(\mbox{\boldmath$\tau$},v))=|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{\theta}^{2}\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

It is also easy to derive the continuity of $B_{\theta}$:

(3.16) $\displaystyle
B_{\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))$
$\displaystyle\leq$
$\displaystyle\|A^{-1/2}\mbox{\boldmath$\chi$}\|_{0}\|A^{-1/2}\mbox{\boldmath$\tau$}\|_{0}+\|A^{1/2}\nabla
w\|_{0}\|A^{1/2}\nabla v\|_{0}+\|A^{1/2}\nabla
w\|_{0}\|A^{-1/2}\mbox{\boldmath$\tau$}\|_{0}$
$\displaystyle+\|A^{-1/2}\mbox{\boldmath$\chi$}\|_{0}\|A^{1/2}\nabla
v\|_{0}+\|\sqrt{\theta/\alpha}\nabla\cdot\mbox{\boldmath$\chi$}\|_{0}\|\sqrt{\theta/\alpha}\nabla\cdot\mbox{\boldmath$\tau$}\|_{0}$
$\displaystyle=$
$\displaystyle(\|A^{-1/2}\mbox{\boldmath$\chi$}\|_{0}+\|A^{1/2}\nabla
w\|_{0})(\|A^{-1/2}\mbox{\boldmath$\tau$}\|_{0}+\|A^{1/2}\nabla
v\|_{0})+\|\sqrt{\theta/\alpha}\nabla\cdot\mbox{\boldmath$\chi$}\|_{0}\|\sqrt{\theta/\alpha}\nabla\cdot\mbox{\boldmath$\tau$}\|_{0}$
$\displaystyle\leq$ $\displaystyle
2|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{\theta},\quad\forall(\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v)\in{\bf
X}$

With the coercivity and continuity of the bilinear form $B_{\theta}$, by the
Lax-Milgram Lemma, (3.10) has a unique solution
$(\mbox{\boldmath$\sigma$},u)\in{\bf X}$.

Let $\Sigma_{h,N}\subset H_{N}({\rm div};\Omega)$ and $V_{h,D}\subset
H^{1}_{D}(\Omega)$ be two finite dimensional subspaces, then we have the
following discrete problem: find
$(\mbox{\boldmath$\sigma$}_{h},u_{h})\in\Sigma_{h,N}\times V_{h,D}$ such that

(3.17)
$B_{\theta}((\mbox{\boldmath$\sigma$}_{h},u_{h}),(\mbox{\boldmath$\tau$}_{h},v_{h}))=F_{\theta}(\mbox{\boldmath$\tau$}_{h},v_{h}),\quad\forall(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D}.$

Since $\Sigma_{h,N}\times V_{h,D}\subset{\bf X}$, we have the well-posedness
of the discrete problems (3.17). Also, the following Galerkin orthogonality is
true:

(3.18)
$B_{\theta}((\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h}),(\mbox{\boldmath$\tau$}_{h},v_{h}))=0,\quad\forall(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D}.$

The following Cea’s lemma type of best approximation property is also true.

######  Theorem 3.2.

Assume that $(\mbox{\boldmath$\sigma$}_{h},u_{h})$ is the solution of problem
(3.17), and $(\mbox{\boldmath$\sigma$},u)$ is the solution of the problem
(1.1). The following best-approximation result is true:

(3.19)
$\displaystyle|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h})|\\!|\\!|_{\theta}$
$\displaystyle\leq$ $\displaystyle
2\inf_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D}}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h})|\\!|\\!|_{\theta}.$

###### Proof.

Let $(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times V_{h,D}$. Using
the coercivity and the continuity of the bilinear form, the Galerkin
orthogonality, and the Cauchy-Schwarz inequality, we have

$\displaystyle|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h})|\\!|\\!|^{2}_{\theta}$
$\displaystyle=$ $\displaystyle
B_{\theta}((\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h}),(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h}))=B_{\theta}((\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h}),(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h}))$
$\displaystyle\leq$ $\displaystyle
2|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h})|\\!|\\!|_{\theta}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h})|\\!|\\!|_{\theta},$

which implies (3.19). ∎

To derive the a posteriori error estimate, we need the following lemma.

######  Lemma 3.3.

(Error representation) Let $(\mbox{\boldmath$\sigma$},u)$ be the solution of
(1.1), $(\mbox{\boldmath$\sigma$}_{h},u_{h})$ be the solution of (3.17), and
$v_{h}\in V_{h,D}$ be an arbitrary function in the discrete space $V_{h,D}$.
We have the following error representation with ${\bf
E}=\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}$ and $e=u-u_{h}$:

(3.20) $|\\!|\\!|({\bf E},e)|\\!|\\!|^{2}_{\theta}=({\bf
f}-A^{-1}\mbox{\boldmath$\sigma$}_{h}-\nabla u_{h},{\bf
E}+A\nabla(e-v_{h}))+(\theta\alpha^{-1}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{h}),\nabla\cdot{\bf
E})+2(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{h},e-v_{h}).$

###### Proof.

Let $\tilde{e}=e-v_{h}$. By the fact ${\bf E}\in H_{N}({\rm div};\Omega)$ and
$e\in H_{D}^{1}(\Omega)$, the coercivity (3.15), the Galerkin orthogonality

(3.21) $B_{\theta}(({\bf E},e),(0,v_{h}))=0,\quad\forall v_{h}\in V_{h,D},$

the definitions of $B_{\theta}$ (3.11) and $F_{\theta}$ (3.13), and the
integration by parts, we get

$\displaystyle|\\!|\\!|({\bf E},e)|\\!|\\!|^{2}_{\theta}$ $\displaystyle=$
$\displaystyle B_{\theta}(({\bf E},e),({\bf E},e))=B_{\theta}(({\bf
E},e),({\bf E},\tilde{e}))=F_{\theta}({\bf
E},\tilde{e})-B_{\theta}((\mbox{\boldmath$\sigma$}_{h},u_{h}),({\bf
E},\tilde{e}))$ $\displaystyle=$ $\displaystyle({\bf
f}-A^{-1}\mbox{\boldmath$\sigma$}_{h}-\nabla u_{h},{\bf
E}+A\nabla\tilde{e})+(\theta\alpha^{-1}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{h}),\nabla\cdot{\bf
E})+2(g,\tilde{e})+2(\mbox{\boldmath$\sigma$}_{h},\nabla\tilde{e})$
$\displaystyle=$ $\displaystyle({\bf
f}-A^{-1}\mbox{\boldmath$\sigma$}_{h}-\nabla u_{h},{\bf
E}+A\nabla\tilde{e})+(\theta\alpha^{-1}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{h}),\nabla\cdot{\bf
E})+2(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{h},\tilde{e}).$

The lemma is then proved. ∎

###  3.2. Symmetric formulations

The formulation (3.10) is non-symmetric. In many situations, for example,
eigenvalues problems or developing efficient linear solvers, symmetric
formulations are always preferred. Also, we can always associate a Ritz-
minimization variational principle to a symmetric problem. Luckily, the method
(3.10) is equivalent to a symmetric GLS formulation by adding least-squares
residuals

$-\displaystyle\frac{1}{2}(A^{-1}\mbox{\boldmath$\sigma$}+\nabla
u,\mbox{\boldmath$\tau$}+A\nabla v)=-\displaystyle\frac{1}{2}({\bf
f},\mbox{\boldmath$\tau$}+A\nabla
v)\quad\mbox{and}\quad\displaystyle\frac{1}{2}(\alpha^{-1}\nabla\cdot\mbox{\boldmath$\sigma$},\nabla\cdot\mbox{\boldmath$\tau$})=\displaystyle\frac{1}{2}(\alpha^{-1}g,\nabla\cdot\mbox{\boldmath$\tau$})$

to the following symmetric saddle point mixed formulation: find
$(\mbox{\boldmath$\sigma$},u)\in{\bf X}$, such that

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

We have the following symmetric formulation: find
$(\mbox{\boldmath$\sigma$},u)\in{\bf X}$, such that

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

with the forms are defined for $(\mbox{\boldmath$\chi$},w)\in{\bf X}$ and
$(\mbox{\boldmath$\tau$},v)\in{\bf X}$,

(3.23) $\displaystyle
B_{sym,\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))$
$\displaystyle:=$
$\displaystyle(A^{-1}\mbox{\boldmath$\chi$},\mbox{\boldmath$\tau$})+(\nabla
w,\mbox{\boldmath$\tau$})+(\mbox{\boldmath$\chi$},\nabla v)-(A\nabla w,\nabla
v)+(\theta\alpha^{-1}\nabla\cdot\mbox{\boldmath$\chi$},\nabla\cdot\mbox{\boldmath$\tau$}).$
(3.24) $\displaystyle F_{sym,\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$}).$

Let $\Sigma_{h,N}\subset H_{N}({\rm div};\Omega)$ and $V_{h,D}\subset
H^{1}_{D}(\Omega)$ be two finite dimensional subspaces, then we have the
following discrete problem corresponding to (3.22): find
$(\mbox{\boldmath$\sigma$}_{h},u_{h})\in\Sigma_{h,N}\times V_{h,D}$ such that

(3.25)
$B_{sym,\theta}((\mbox{\boldmath$\sigma$}_{h},u_{h}),(\mbox{\boldmath$\tau$}_{h},v_{h}))=F_{sum,\theta}(\mbox{\boldmath$\tau$}_{h},v_{h}),\quad\forall(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D}.$

It is easy to see that (3.10) and (3.22) are equivalent since replacing the
test function $v$ by $-v$ in one formulation leads to the other formulation.
By doing so, we know that (3.10) and (3.22) and their corresponding finite
element formulations (3.17) and (3.25) produce identical solutions. Thus all
the analysis of the non-symmetric formulations can be applied to the symmetric
versions.

Alternatively, we can establish the inf-sup stability of the symmetric
formulation directly.

######  Lemma 3.4.

The following robust inf-sup stability with the stability constant being $1$
holds:

(3.26) $\inf_{(\mbox{\boldmath$\chi$},w)\in{\bf
X}}\sup_{(\mbox{\boldmath$\tau$},v)\in{\bf
X}}\displaystyle\frac{B_{sym,\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))}{|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{\theta}}=\inf_{(\mbox{\boldmath$\tau$},v)\in{\bf
X}}\sup_{(\mbox{\boldmath$\chi$},w)\in{\bf
X}}\displaystyle\frac{B_{sym,\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))}{|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{\theta}}\geq
1.$

###### Proof.

Due to the fact $B_{sym,\theta}$ is symmetric, we only need to show one inf-
sup condition in (3.26). Let
$(\mbox{\boldmath$\tau$},v)=(\mbox{\boldmath$\chi$},-w)$ and using the facts
that
$B_{\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\chi$},w)=B_{sym,\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\chi$},-w))$
and (3.15) then

$\displaystyle\sup_{(\mbox{\boldmath$\tau$},v)\in{\bf
X}}\displaystyle\frac{B_{sym,\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))}{|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{\theta}}$
$\displaystyle\geq$
$\displaystyle\displaystyle\frac{B_{sym,\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\chi$},-w))}{|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}}=\displaystyle\frac{B_{\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\chi$},w))}{|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}}=\displaystyle\frac{|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}^{2}}{|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}}=|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}.$

∎

Similarly, we have the following robust discrete inf-sup stability with ${\bf
X}_{h}=\Sigma_{h,N}\times V_{h,D}$:

(3.27) $\inf_{(\mbox{\boldmath$\chi$}_{h},w_{h})\in{\bf
X}_{h}}\sup_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in{\bf
X}_{h}}\displaystyle\frac{B_{sym,\theta}((\mbox{\boldmath$\chi$}_{h},w_{h}),(\mbox{\boldmath$\tau$}_{h},v_{h}))}{|\\!|\\!|(\mbox{\boldmath$\chi$}_{h},w_{h})|\\!|\\!|_{\theta}|\\!|\\!|(\mbox{\boldmath$\tau$}_{h},v_{h})|\\!|\\!|_{\theta}}=\inf_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in{\bf
X}_{h}}\sup_{(\mbox{\boldmath$\chi$}_{h},w_{h})\in{\bf
X}_{h}}\displaystyle\frac{B_{sym,\theta}((\mbox{\boldmath$\chi$}_{h},w_{h}),(\mbox{\boldmath$\tau$}_{h},v_{h}))}{|\\!|\\!|(\mbox{\boldmath$\chi$}_{h},w_{h})|\\!|\\!|_{\theta}|\\!|\\!|(\mbox{\boldmath$\tau$}_{h},v_{h})|\\!|\\!|_{\theta}}\geq
1.$

Also, it is easy to show the continuity of $B_{sym,\theta}$ as in (3.16):

(3.28)
$B_{sym,\theta}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))\leq
2|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{\theta}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{\theta},\quad\forall(\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

Using the theorem in [44], we immediately have following robust best
approximation theorem, which is identical to Theorem 3.2.

######  Theorem 3.5.

Assume that $(\mbox{\boldmath$\sigma$}_{h},u_{h})$ is the solution of problem
(3.25), and $(\mbox{\boldmath$\sigma$},u)$ is the solution of the problem
(1.1). The following robust best-approximation result is true:

(3.29)
$\displaystyle|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h})|\\!|\\!|_{\theta}$
$\displaystyle\leq$ $\displaystyle
2\inf_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D}}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h})|\\!|\\!|_{\theta}.$

######  Remark 3.6.

The formulation (3.22) (with $\theta=1$ and $0$) can be found in Section 4.1
(Symmetric stabilizations in $H({\rm div};\Omega)\times H^{1}(\Omega)$) of
[29]. The paper [29] mainly wanted to use continuous finite elements to
approximate both $\sigma$ and $u$. It did mention the possible usage of
$H({\rm div})$-conforming finite elements. Robust a priori error estimate with
respect to the coefficient matrix $A$ is not sought in [29]. A posteriori
error estimate is not discussed in [29].

##  4\. Assumptions on the coefficient matrix $A$ and Robust Interpolations

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

######  Assumption 4.1.

Piecewise constant assumption on $A$. We assume that $A=\alpha(x)I$ where
$\alpha(x)$ is positive and piecewise constant function in $\Omega$ with
possible large jumps across subdomain boundaries (interfaces):

$\alpha(x)=\alpha_{i}>0\mbox{ in }\Omega_{i}$

for $i=1,\cdots,n$. Here, $\\{\Omega_{i}\\}_{i=1}^{n}$ is a partition of the
domain $\Omega$ with $\Omega_{i}$ being an open polygonal domain.

######  Remark 4.2.

For the more general cases of $A$, the analysis in this paper will still be
valid. However, the genetic constants appeared in the paper will depend on the
ratio $\lambda_{\max,K}/\lambda_{\min,K}$, for all $K\in{\mathcal{T}}$, where
$\lambda_{\max,K}$ and $\lambda_{\min,K}$ are the respective maximal and
minimal eigenvalues of $A_{K}:=A|_{K}$. See discussion in [4, 15].

We then discuss the quasi-monotonicity assumption and robust quasi-
interpolation based on this assumption.

######  Assumption 4.3.

Quasi-monotonicity assumption (QMA). Assume that any two different subdomains
$\overline{\Omega}_{i}$ and $\overline{\Omega}_{i}$, which share at least one
point, have a connected path passing from $\overline{\Omega}_{i}$ to
$\overline{\Omega}_{j}$ through adjacent subdomains such that the diffusion
coefficient $\alpha(x)$ is monotone along this path.

It is also common to use Clément-type interpolation operators (see, e.g., [4,
41]) for establishing the reliability bound of a posteriori error estimators.
Following [4], one can define the interpolation operator
$I_{rcl}:L^{2}(\Omega)\rightarrow S^{1}_{D}$ (see [22] for more details) so
that the following estimates are true under Assumption 4.3 (QMA):

(4.1)
$\|\alpha_{K}^{\frac{1}{2}}(v-I_{rcl}v)\|_{0,K}+h_{K}\|\alpha_{K}^{\frac{1}{2}}\nabla(v-I_{rcl}v)\|_{0,K}\leq
C\,h_{K}\|\alpha^{\frac{1}{2}}\nabla v\|_{0,\Delta_{K}}\quad\forall
K\in{\mathcal{T}},v\in H^{1}_{D}(\Omega),$

where $\Delta_{K}$ is the union of all elements that share at least one vertex
with K,

In general, we do not have the following robust Poincaré-Friedrichs inequality

(4.2) $\|\alpha^{1/2}v\|_{0}\leq C\|\alpha^{1/2}\nabla v\|_{0},\quad v\in
H^{1}_{D}(\Omega),$

where $C>0$ is independent of $\alpha$.

For the following special case, (4.2) does holds.

######  Lemma 4.4.

Assume that each $\Omega_{i}$ in the Assumption 4.1 has a part of the
Dirichlet boundary condition with a positive measure, i.e.,

$\mbox{measure}\\{\Gamma_{D}^{i}\\}>0,\mbox{ where
}\Gamma_{D}^{i}=\partial\Omega_{i}\cap\Gamma_{D}\quad\mbox{for }i=1,\cdots n,$

then the robust Poincaré inequality (4.2) is true with $C>0$ being independent
of $\alpha$.

The proof is strghtforward since $v_{i}=v|_{\Omega_{i}}\in\\{w\in
H^{1}(\Omega_{i}),w|_{\Gamma_{D}^{i}}=0\\}$, and $\alpha$ in $\Omega_{i}$ is
$\alpha_{i}$, a single constant, we have

$\|\alpha_{i}^{1/2}v_{i}\|_{0,\Omega_{i}}\leq C\|\alpha_{i}^{1/2}\nabla
v_{i}\|_{0,\Omega_{i}},$

with $C>0$ independent of $\alpha_{i}$. Summing up all the subdomains we get
the robust Poincaré-Friedrichs inequality for this special case.

##  5\. The first augmented mixed formulation: $\theta=1$

In this section, we consider that case that $\theta=1$. For simplicity, we
consider the finite element approximation in $RT_{0,N}\times
S_{1,D}\subset{\bf X}$. The discrete problem is: find
$(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})\in RT_{0,N}\times S_{1,D}$ such that

(5.1)
$B_{1}((\mbox{\boldmath$\sigma$}_{1,h},u_{1,h}),(\mbox{\boldmath$\tau$}_{h},v_{h}))=F_{1}(\mbox{\boldmath$\tau$}_{h},v_{h}),\quad\forall(\mbox{\boldmath$\tau$}_{h},v_{h})\in
RT_{0,N}\times S_{1,D},$

where $B_{1}$ and $F_{1}$ are the corresponding forms (3.11) and (3.13) with
$\theta=1$. Based on the discussions in Section 3.1, we have the well-
posedness of discrete problem (5.1). Let $|\\!|\\!|(\cdot,\cdot)|\\!|\\!|_{1}$
be the norm defined in (3.14) with $\theta=1$. We have the following locally
robust and optimal a priori error estimate.

######  Theorem 5.1.

Let $(\mbox{\boldmath$\sigma$},u)$ be the solution of (1.1) and
$(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$ be the solution of problem (5.1),
respectively. We have the following a priori estimate:

(5.2)
$\displaystyle|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|_{1}$
$\displaystyle\leq$ $\displaystyle
2\inf_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in RT_{0,N}\times
S_{1,D}}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h})|\\!|\\!|_{1}.$

Under Assumption 4.1 on the coefficients, if we further assume that $u|_{K}\in
H^{1+s_{K}}(K)$, $(\nabla u-{\bf f})|_{K}\in H^{q_{K}}(K)^{d}$, $g|_{K}\in
H^{t_{K}}(K)$ for $K\in{\mathcal{T}}$, where the local regularity indexes
$s_{K}$, $q_{K}$, and $t_{K}$ satisfy the following assumptions: $0<s_{K}\leq
1$ in two dimensions and $1/2<s_{K}\leq 1$ in there dimensions, $1/2<q_{K}\leq
1$ with the constant $C_{rt}>0$ being unbounded as $q_{K}\downarrow 1/2$, and
$0<t_{K}\leq 1$, then the following local robust and local optimal a priori
error estimate holds: there exists a constant $C$ independent of $\alpha$ and
the mesh-size, such that

(5.3)
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|_{1}\leq
C\sum_{K\in{\mathcal{T}}}\alpha_{K}^{1/2}\left(h_{K}^{s_{K}}|\nabla
u|_{s_{K},K}+C_{rt}h_{K}^{q_{K}}|\nabla u-{\bf
f}|_{q_{K},K}+\alpha_{K}^{-1}h_{K}^{t_{K}}|g|_{t_{K},K}\right).$

###### Proof.

The result (5.2) is from (3.19). The result of (5.3) can be derived from
regularity assumptions and interpolation results (2.2), (2.3), and (2.7). ∎

######  Remark 5.2.

In theorem 5.1, when ${\bf f}=0$, then $s_{K}=q_{K}$ for each
$K\in{\mathcal{T}}$. When ${\bf f}\neq{\bf 0}$, the local regularity of
$\sigma$ and $\nabla u$ can be different, and the regularity of $g$ (which is
$\nabla\cdot\mbox{\boldmath$\sigma$}$) is also independent of that of $u$.

###  5.1. A least-squares a posteriori error estimator for the first augmented
mixed formulation

We discuss a posteriori error estimator for the first augmented mixed
formulation in this subsection.

Define an indicator on each $K\in{\mathcal{T}}$:

$\displaystyle\eta_{1,K}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$
$\displaystyle=$
$\displaystyle\Big{\\{}\|\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{1,h})\|^{2}_{0,K}+\|\alpha^{1/2}({\bf
f}-\nabla
u_{1,h}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{1,h})\|^{2}_{0,K}\Big{\\}}^{1/2}.$

Define the corresponding global a posteriori error estimator:

(5.4)
$\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})=\Big{\\{}\|\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{1,h})\|^{2}_{0}+\|\alpha^{1/2}({\bf
f}-\nabla
u_{1,h}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{1,h})\|^{2}_{0}\Big{\\}}^{1/2}.$

Note that the error estimator
$\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$ is actually a least-squares
error estimator, see discussion below in section 5.2.

######  Theorem 5.3.

Let $(\mbox{\boldmath$\sigma$},u)$ be the solution of (1.1) and
$(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$ be the solution of problem (5.1),
repectively. Assume that Assumption 4.1 on the coefficients and Assumption 4.3
(QMA) are true, then there exists positive constants $C$ independent of
$\alpha$ and the mesh-size such that the following reliability holds:

(5.5)
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|_{1}\leq
C\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h}).$

###### Proof.

Let ${\bf E}_{1}=\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h}$ and
$e_{1}=u-u_{1,h}$. In (3.20), let $\theta=1$ and $v_{h}=I_{rcl}e_{1}$, we have

$|\\!|\\!|({\bf E}_{1},e_{1})|\\!|\\!|^{2}_{1}=({\bf
f}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{1,h}-\nabla u_{1,h},{\bf
E}_{1}+\alpha\nabla(e_{1}-I_{rcl}e_{1}))+(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{1,h},\alpha^{-1}\nabla\cdot{\bf
E}_{1}+2(e_{1}-I_{rcl}e_{1})).$

Applying Cauchy-Schwarz and triangle inequalities, we get

$\displaystyle|\\!|\\!|({\bf E}_{1},e_{1})|\\!|\\!|^{2}_{1}$
$\displaystyle\leq$ $\displaystyle\|\alpha^{1/2}({\bf
f}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{1,h}-\nabla
u_{1,h})\|_{0}(\|\alpha^{-1/2}{\bf
E}_{1}\|_{0}+\|\alpha^{1/2}\nabla(e_{1}-I_{rcl}e_{1})\|_{0})$
$\displaystyle+\|\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{1,h})\|_{0}(\|\alpha^{-1/2}\nabla\cdot{\bf
E}_{1}\|_{0}+2\|\alpha^{1/2}(e_{1}-I_{rcl}e_{1})\|_{0}.$

By the robust Clément interpolation result (4.1) under Assumption 4.3 and the
fact that the mesh-size is bounded, we have the following robust results,

$\|\alpha^{1/2}\nabla(e_{1}-I_{rcl}e_{1})\|_{0}\leq C\|\alpha^{1/2}\nabla
e_{1}\|_{0}\quad\mbox{and}\quad\|\alpha^{1/2}(e_{1}-I_{rcl}e_{1})\|_{0}\leq
C\|\alpha^{1/2}h\nabla e_{1}\|_{0}\leq C\|\alpha^{1/2}\nabla e_{1}\|_{0}.$

Substitute these two robust results into (5.1), we get

$\displaystyle|\\!|\\!|({\bf E}_{1},e_{1})|\\!|\\!|^{2}_{1}$
$\displaystyle\leq$ $\displaystyle C(\|\alpha^{1/2}({\bf
f}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{1,h}-\nabla
u_{1,h})\|_{0}+\|\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{1,h})\|_{0})$
$\displaystyle(\|\alpha^{-1/2}{\bf E}_{1}\|_{0}+\|\alpha^{-1/2}\nabla\cdot{\bf
E}_{1}\|_{0}+\|\alpha^{1/2}\nabla e_{1}\|_{0})$ $\displaystyle\leq$
$\displaystyle C\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})|\\!|\\!|({\bf
E}_{1},e_{1})|\\!|\\!|_{1}.$

The robust reliability result (5.5) is proved. ∎

The efficiency of the proposed error indicator is the same as the standard
least-squares a posteriori error estimator.

######  Theorem 5.4.

Let $(\mbox{\boldmath$\sigma$},u)$ be the solution of (1.1) and
$(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$ be the solution of problem (5.1),
respectively. Assume that Assumption 4.1 on the coefficients is true, we have
the following efficiency:

(5.7)
$\eta_{1,K}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})\leq\sqrt{2}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|_{1,K}\quad\forall
K\in{\mathcal{T}}.$

###### Proof.

By the triangle inequality and the first-order system (1.1), for any
$K\in{\mathcal{T}}$, we have

$\displaystyle\eta_{1,K}^{2}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$
$\displaystyle=$
$\displaystyle\|\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{1,h})\|^{2}_{0,K}+\|\alpha^{1/2}{\bf
f}-\alpha^{-1/2}\mbox{\boldmath$\sigma$}_{1,h}-\alpha^{1/2}\nabla
u_{1,h}\|^{2}_{0,K}$ $\displaystyle=$
$\displaystyle\|\alpha^{-1/2}(\nabla\cdot\mbox{\boldmath$\sigma$}-\nabla\cdot\mbox{\boldmath$\sigma$}_{1,h})\|^{2}_{0,K}+\|\alpha^{-1/2}(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h})+\alpha^{1/2}\nabla(u-u_{1,h})\|^{2}_{0,K}$
$\displaystyle\leq$
$\displaystyle\|\alpha^{-1/2}\nabla\cdot(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h})\|^{2}_{0,K}+2(\|\alpha^{-1/2}(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h})\|^{2}_{0,K}+\|\alpha^{1/2}\nabla(u-u_{1,h})\|^{2}_{0,K}\big{)}$
$\displaystyle=$ $\displaystyle
2|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|^{2}_{1,K}.$

The theorem is proved. ∎

######  Remark 5.5.

Note that Assumption 4.3 (QMA) is not required for the robustness of the
efficiency bound.

###  5.2. Comparison with the $L^{2}$-based LSFEM

For the first-order system (1.1), the $L^{2}$-based least-squares functional
is

(5.8) $J(\mbox{\boldmath$\tau$},v;{\bf f},g):=\|\alpha^{1/2}\nabla
v+\alpha^{-1/2}\mbox{\boldmath$\tau$}-\alpha^{1/2}{\bf
f}\|_{0}^{2}+\|\alpha^{-1/2}(\nabla\cdot\mbox{\boldmath$\tau$}-g)\|_{0}^{2},\quad(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

Then the $L^{2}$-based least-squares minimization problem is: find
$(\mbox{\boldmath$\sigma$},u)\in{\bf X}$, such that

(5.9) $J(\mbox{\boldmath$\sigma$},u;{\bf
f},g)=\inf_{(\mbox{\boldmath$\tau$},v)\in{\bf
X}}J(\mbox{\boldmath$\tau$},v;{\bf f},g).$

Equivalently, it can be written in a weak form as: find
$(\mbox{\boldmath$\sigma$},u)\in{\bf X}$, such that

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

where the least-squares bilinear form $b$ is defined as follows:

(5.11)
$b((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))=(\alpha^{-1}\mbox{\boldmath$\chi$}+\nabla
w,\mbox{\boldmath$\tau$}+\alpha\nabla
v)+(\alpha^{-1}\nabla\cdot\mbox{\boldmath$\chi$},\nabla\cdot\mbox{\boldmath$\tau$}),\quad\forall(\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

The lowest-order least-squares finite element method (LSFEM) is: find
$(\mbox{\boldmath$\sigma$}_{h}^{ls},u_{h}^{ls})\in RT_{0,N}\times S_{1,D}$,
such that

(5.12) $J(\mbox{\boldmath$\sigma$}^{ls}_{h},u^{ls}_{h};{\bf
f},g)=\inf_{(\mbox{\boldmath$\tau$},v)\in RT_{0,N}\times
S_{1,D}}J(\mbox{\boldmath$\tau$},v;{\bf f},g).$

Or, equivalently, find $(\mbox{\boldmath$\sigma$}_{h}^{ls},u_{h}^{ls})\in
RT_{0,N}\times S_{1,D}$, such that,

(5.13)
$b((\mbox{\boldmath$\sigma$}_{h}^{ls},u_{h}^{ls}),(\mbox{\boldmath$\tau$},v))=({\bf
f},\mbox{\boldmath$\tau$}+\alpha\nabla
v)+(\alpha^{-1}g,\nabla\cdot\mbox{\boldmath$\tau$}),\quad\forall(\mbox{\boldmath$\tau$},v)\in
RT_{0,N}\times S_{1,D}.$

The key ingredient to establish a priori and a posteriori error estimates of
the LSFEM (5.12) or (5.13) is the following norm equivalence:

(5.14) $C_{coe}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{1}^{2}\leq
J(\mbox{\boldmath$\tau$},v;0,0)\leq
C_{con}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{1}^{2},\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf
X},$

for some positive constants $C_{coe}$ and $C_{con}$. The first inequality of
(5.14) is the coercivity of the least-squares bilinear form (5.10):

(5.15) $b((\mbox{\boldmath$\tau$},v),(\mbox{\boldmath$\tau$},v))\geq
C_{coe}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{1}^{2}\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

The second inequality of (5.14) is equivalent to the continuity of the least-
squares bilinear form (5.10) since it is symmetric:

(5.16) $b((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))\leq
C_{con}|\\!|\\!|(\mbox{\boldmath$\chi$},w)|\\!|\\!|_{1}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{1},\quad\forall(\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

With the coercivity and continuity, we can easily get the a priori error
estimate of the LSFEM (5.13):

(5.17)
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{ls},u-u_{h}^{ls})|\\!|\\!|_{1}\leq\displaystyle\frac{C_{con}}{C_{coe}}\inf_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in
RT_{0,N}\times
S_{1,D}}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h})|\\!|\\!|_{1}.$

By the triangle inequality, it is easy to prove that continuity constant
$C_{con}$ is a constant independent of $\alpha$ (actually, $C_{con}=2$ for
this case, see arguments of the proof in Theorem 5.4). On the other hand, the
coercivity constant $C_{coe}$ usually depends $\alpha$. The reason is that the
proof of the coercivity of the least-squares bilinear form (5.10) requires the
Poincaré inequality, which is usually not robust with respect to $\alpha$. For
the special case discussed in Lemma 4.4, we have the following result:

######  Theorem 5.6.

Assume that each $\Omega_{i}$ in the Assumption 4.1 has a part of the
Dirichlet boundary condition with a positive measure, then the coercivity
constant $C_{coe}$ the norm equivalence (5.14) and the coercivity (5.15) is
independent of $\alpha$.

###### Proof.

For a $v\in H^{1}_{D}(\Omega)$, let $\tau$ be an arbitrary vector in
$H_{N}({\rm div};\Omega)$, then

$\displaystyle\|\alpha^{1/2}\nabla v\|_{0}^{2}$ $\displaystyle=$
$\displaystyle(\alpha\nabla v,\nabla v)=(\alpha\nabla
v+\mbox{\boldmath$\tau$},\nabla v)-(\mbox{\boldmath$\tau$},\nabla
v)=(\alpha\nabla v+\mbox{\boldmath$\tau$},\nabla
v)+(\nabla\cdot\mbox{\boldmath$\tau$},v)$ $\displaystyle\leq$
$\displaystyle\|\alpha^{1/2}\nabla
v+\alpha^{-1/2}\mbox{\boldmath$\tau$}\|_{0}\|\alpha^{1/2}\nabla
v\|_{0}+\|\alpha^{-1/2}\nabla\cdot\mbox{\boldmath$\tau$}\|_{0}\|\alpha^{1/2}v\|_{0}.$

By lemma 4.4, for the special setting of the theorem, a robust Poincaré
inequality $\|\alpha^{1/2}v\|_{0}\leq C\|\alpha^{1/2}\nabla v\|_{0}$ for $v\in
H^{1}_{D}(\Omega)$ is true. Thus, we have

(5.18) $\|\alpha^{1/2}\nabla v\|_{0}\leq\|\alpha^{1/2}\nabla
v+\alpha^{-1/2}\mbox{\boldmath$\tau$}\|_{0}+C\|\alpha^{-1/2}\nabla\cdot\mbox{\boldmath$\tau$}\|_{0},\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf
X},$

with the constant $C$ independent of $\alpha$. It is also simple to see that,
for $(\mbox{\boldmath$\tau$},v)\in{\bf X}$,

(5.19) $\|\alpha^{-1/2}\mbox{\boldmath$\tau$}\|_{0}\leq\|\alpha^{1/2}\nabla
v+\alpha^{-1/2}\mbox{\boldmath$\tau$}\|_{0}+\|\alpha^{1/2}\nabla v\|_{0}\leq
2\|\alpha^{1/2}\nabla
v+\alpha^{-1/2}\mbox{\boldmath$\tau$}\|_{0}+C\|\alpha^{-1/2}\nabla\cdot\mbox{\boldmath$\tau$}\|_{0}.$

Thus, we prove that robust coercivity. ∎

From the above proof, we also confirm that the weight $\alpha^{-1/2}$ in
$\|\alpha^{-1/2}\nabla\cdot\mbox{\boldmath$\tau$}\|_{0}$ in (5.8) is the right
choice.

From (5.17), except in the special cases where the robust version of Poincaré
inequality holds, the a priori error estimate of the LSFEM (5.12), (5.13) in
general is not robust with respect to $\alpha$.

Let $(\mbox{\boldmath$\sigma$}_{a},u_{a})\in{\bf X}$, we can define the
following least-squares based a posteriori error estimator:

(5.20)
$\eta_{ls}(\mbox{\boldmath$\sigma$}_{a},u_{a})=\Big{(}\|\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{a})\|^{2}_{0}+\|\alpha^{1/2}({\bf
f}-\nabla
u_{a}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{a})\|^{2}_{0}\Big{)}^{1/2}=J(\mbox{\boldmath$\sigma$}_{a},u_{a};{\bf
f},g)^{1/2}.$

Let ${\bf E}_{a}=\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{a}$ and
$e_{a}=u-u_{a}$. Using the facts that ${\bf f}=\nabla
u+\alpha^{-1}\mbox{\boldmath$\sigma$}$ and
$g=\nabla\cdot\mbox{\boldmath$\sigma$}$ from (1.1), we have the following
identity:

$\displaystyle J(\mbox{\boldmath$\sigma$}_{a},u_{a};{\bf f},g)$
$\displaystyle=$
$\displaystyle\|\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{a})\|^{2}_{0}+\|\alpha^{1/2}({\bf
f}-\nabla u_{a}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{a})\|^{2}_{0}$
$\displaystyle=$
$\displaystyle\|\alpha^{-1/2}(\nabla\cdot\mbox{\boldmath$\sigma$}-\nabla\cdot\mbox{\boldmath$\sigma$}_{a})\|^{2}_{0}+\|\alpha^{1/2}(\nabla
u+\alpha^{-1}\mbox{\boldmath$\sigma$}-\nabla
u_{a}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{a})\|^{2}_{0}=J({\bf
E}_{a},e_{a};0,0).$

By (5.14), the following reality and efficiency bounds are true,

(5.21) $C_{coe}|\\!|\\!|({\bf E}_{a},e_{a})|\\!|\\!|_{1}^{2}\leq J({\bf
E}_{a},e_{a};0,0)=J(\mbox{\boldmath$\sigma$}_{a},u_{a};{\bf f},g)\leq
C_{con}|\\!|\\!|({\bf E}_{a},e_{a})|\\!|\\!|_{1}^{2}.$

An important fact of (5.21) is that
$(\mbox{\boldmath$\sigma$}_{a},u_{a})\in{\bf X}$ does not need to be the
numerical solution of the LSFEM problem (5.12). In fact, the pair can be any
functions in ${\bf X}$. Let $(\mbox{\boldmath$\sigma$}_{h}^{ls},u_{h}^{ls})\in
RT_{0,N}\times S_{1,D}\subset{\bf X}$ be the numerical solution of the LSFEM
problem (5.13), we immediately have the reliability and efficiency of the
least-squares error estimator for the LSFEM approximation (5.13),

(5.22)
$C_{coe}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{ls},u-u_{h}^{ls})|\\!|\\!|_{1}^{2}\leq\eta_{ls}(\mbox{\boldmath$\sigma$}_{h}^{ls},u_{h}^{ls})^{2}=J(\mbox{\boldmath$\sigma$}_{h}^{ls},u_{h}^{ls};{\bf
f},g)\leq
C_{con}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{ls},u-u_{h}^{ls})|\\!|\\!|_{1}^{2}.$

Of course, since $C_{coe}$ depends on $\alpha$, the a posteriori error
estimator $\eta_{ls}(\mbox{\boldmath$\sigma$}_{h}^{ls},u_{h}^{ls})$ is not
robust for the LSFEM approximation (5.13).

If the robustness of the estimator is not our goal, we can have the non-robust
reliability and robust efficiency of the a posteriori error estimator
$\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$ of the first augmented
mixed method (5.4) by using the fact that
$(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})\in RT_{0,N}\times S_{1,D}\subset{\bf
X}$ and (5.21):

(5.23) $C_{coe}|\\!|\\!|({\bf
E}_{1},e_{1})|\\!|\\!|_{1}^{2}\leq\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})^{2}=\eta_{ls}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})^{2}\leq
C_{con}|\\!|\\!|({\bf E}_{1},e_{1})|\\!|\\!|_{1}^{2}.$

This result is weaker than that of Theorem 5.3.

######  Remark 5.7.

It is interesting to see that for the first augmented mixed method (5.1) and
the LSFEM (5.12) or (5.13), their a posteriori error estimators are the same.
However, one is robust, and the other one is not. The subtle difference is
that the numerical solution in the a posteriori error estimator (5.4) is
obtained by the robust first augmented mixed method (5.1), and the Galerkin
orthogonality (3.21) is used in the error representation Lemma 3.3. On the
other hand, the reliability and efficiency of a general least-squares a
posteriori error estimator do not require that the approximations are the
numerical solutions of the corresponding LSFEM or augmented mixed method.

######  Remark 5.8.

We take a comparison of the least-squares problem (5.10) and the augmented
problem (3.10) with $\theta=1$ with the bilinear form in the form of (3.12).
The first augmented mixed method can be viewed as adding a consistent term,

$-2(\mbox{\boldmath$\sigma$},\nabla
v)=2(\nabla\cdot\mbox{\boldmath$\sigma$},v)=(g,v)\quad\forall v\in
H^{1}_{D}(\Omega)$

to the least-squares problem (5.10). The extra term makes the new formulation
lose the least-squares energy minimization principle, but
$-2(\mbox{\boldmath$\tau$},\nabla v)$ cancels the cross-term in
$\|A^{-1/2}\mbox{\boldmath$\tau$}+A^{1/2}\nabla v\|_{0}^{2}$, thus the new
augmented mixed method is robust in the energy norms.

######  Remark 5.9.

We also need to mention that the non-robustness of the $L^{2}$-LSFEM is mild.
Take a close look at the proof of the robustness in Theorem 5.6, a robust
Poincaré inequality for any $v\in H^{1}_{D}(\Omega)$ is needed. In general,
the robust Poincaré inequality is not true. But, for the robust error analysis
of LSFEM (5.22), we only need the fact

(5.24) $\|\alpha^{1/2}(u-u^{ls}_{h})\|_{0}\leq
C\|\alpha^{1/2}\nabla(u-u^{ls}_{h})\|_{0},$

for a constant $C>0$ independent of $\alpha$. Assuming that we have some
$L^{2}$-error bound $\|u-u^{ls}_{h}\|_{0}\leq
Ch^{r}\|\nabla(u-u^{ls}_{h})\|_{0}$ for some regularity $r>0$ and $h$ being
the maximum size of the mesh, see discussions in [16], then
$\|\alpha^{1/2}(u-u^{ls}_{h})\|_{0}\leq
C\|\alpha^{1/2}\nabla(u-u^{ls}_{h})\|_{0}$ for $C$ independent of $\alpha$ is
possible with a small enough $h$. The situation of a non-uniform mesh and a
solution with a low regularity is less clear, but realizing that (5.24) for
the error $u-u^{ls}_{h}$ is what we needed is helpful to explain the mildness
of non-robustness of the LSFEM.

######  Remark 5.10.

In [2], the same error estimator is proposed for the first augmented mixed
method. In the proof of [2], the same technique as the coercivity proof of the
LSFEM (5.15) is used. Since the Galerkin orthogonality (3.21) is not used, the
analysis in [2] is not robust.

##  6\. The second augmented mixed formulation: a mesh-weighted version

###  6.1. The second augmented mixed formulation

In this formulation, we choose $\theta$ to be a piecewisely defined function
such that $\theta|_{K}=h^{2}_{K}$, for $K\in{\mathcal{T}}$.

We consider the finite element approximation in two pairs: $RT_{0,N}\times
S_{1,D}\subset{\bf X}$ and $BDM_{1,N}\times S_{2,D}\subset{\bf X}$. The
notation $\Sigma_{h,N}\times V_{h,D}$ is used to represent these two choices.
The discrete problem is: find
$(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})\in\Sigma_{h,N}\times V_{h,D}$ such
that

(6.1)
$B_{2}((\mbox{\boldmath$\sigma$}_{2,h},u_{2,h}),(\mbox{\boldmath$\tau$}_{h},v_{h}))=F_{2}(\mbox{\boldmath$\tau$}_{h},v_{h}),\quad\forall(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D},$

where $B_{2}$ and $F_{2}$ are the corresponding forms (3.11) and (3.13) with
$\theta|_{K}=h_{K}^{2}$ for any $K\in{\mathcal{T}}$. Based on the discussions
in Section 3.1, we have the well-posedness of discrete problem (6.1). Let
$|\\!|\\!|(\cdot,\cdot)|\\!|\\!|_{2}$ be the norm defined in (3.14) with
$\theta|_{K}=h^{2}_{K}$, for $K\in{\mathcal{T}}$. We also have the following
locally robust and optimal a priori error estimate.

######  Theorem 6.1.

Let $(\mbox{\boldmath$\sigma$},u)$ be the solution of (1.1) and
$(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$ be the solution of problem (6.1),
respectively. The following best approximation result is true:

(6.2)
$\displaystyle|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$
$\displaystyle\leq$ $\displaystyle
2\inf_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D}}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h})|\\!|\\!|_{2}.$

For the $RT_{0,N}\times S_{1,D}$ approximation, under Assumption 4.1 on the
coefficients, if we further assume that $u|_{K}\in H^{1+s_{K}}(K)$, $(\nabla
u-{\bf f})|_{K}\in H^{q_{K}}(K)^{d}$, $g|_{K}\in H^{t_{K}}(K)$ for
$K\in{\mathcal{T}}$, where the local regularity indexes $s_{K}$, $q_{K}$, and
$t_{K}$ satisfy the following assumptions: $0<s_{K}\leq 1$ in two dimensions
and $1/2<s_{K}\leq 1$ in three dimensions, $1/2<q_{K}\leq 1$ with the constant
$C_{rt}>0$ being unbounded as $q_{K}\downarrow 1/2$, and $0<t_{K}\leq 1$, then
the following local robust and local optimal a priori error estimate holds:
there exists a constant $C$ independent of $\alpha$ and the mesh-size, such
that

(6.3)
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}\leq
C\sum_{K\in{\mathcal{T}}}\alpha_{K}^{1/2}\left(h_{K}^{s_{K}}|\nabla
u|_{s_{K},K}+C_{rt}h_{K}^{q_{K}}|\nabla u-{\bf
f}|_{q_{K},K}+\alpha_{K}^{-1}h_{K}^{1+t_{K}}|g|_{t_{K},K}\right).$

For the $BDM_{1,N}\times S_{2,D}$ approximation, under Assumption 4.1 on the
coefficients, if we further assume that $u\in H^{3}(\Omega)$, $(\nabla u-{\bf
f})|_{K}\in H^{2}(K)^{d}$, and $g|_{K}\in H^{1}(K)$ for $K\in{\mathcal{T}}$,
then the following local robust and local optimal a priori error estimate
holds: there exists a constant $C$ independent of $\alpha$ and the mesh-size,
such that

(6.4)
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}\leq
C\sum_{K\in{\mathcal{T}}}\alpha_{K}^{1/2}h_{K}^{2}\left(|\nabla
u|_{2,K}+|\nabla u-{\bf f}|_{2,K}+\alpha_{K}^{-1}|g|_{1,K}\right).$

The proof of the theorem is almost identical to that of Theorem 5.1 with some
necessary changes due to the extra factor $h$.

######  Remark 6.2.

We give some explanations that why the mesh-weighted second augmented mixed
formulation is suggested. In the standard mixed formulation (3.2) for Darcy’s
equation, the divergence of the flux $\nabla\cdot\mbox{\boldmath$\sigma$}=g$
is a known quantity. For the standard dual mixed formulation (3.2), we can
easily derive the robust best approximation in the $L^{2}$-norm of the flux
alone without invoking the approximation of the divergence of $\sigma$ and the
smoothness of $g$, see for example, Theorems 2 and 3 of [45]. On the contrary,
for the first augmented mixed formulation (5.1), its a priori error estimate
(5.3) is done for the combined norm $|\\!|\\!|\cdot|\\!|\\!|_{1}$. Thus, the
error of the flux measured in the $L^{2}$-norm is influenced by the
approximation of the divergence of $\sigma$, which depends on the regularity
of $g$ on the element $K\in{\mathcal{T}}$. For example, for the case ${\bf
f}=0$ and $u\in H^{2}(\Omega)$, we get $s_{K}=q_{K}=1$ for all
$K\in{\mathcal{T}}$ in (5.3). However, if the regularity of $g$ is low ($<1$),
then the error of the flux measured in the $L^{2}$-norm is dominated by the
bad approximation of the divergence of $\sigma$, which is worse than the
standard mixed formulation. Such sub-optimal result also appears in the LSFEM
(5.13) since all the terms are also coupled in (5.13).

On the other hand, for the second augmented mixed method, the convergence
order of $\|\alpha^{1/2}\nabla(u-u_{2,h})\|_{0}$ and
$\|\alpha^{-1/2}\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h})\|_{0}$
can still be optimal, even if the regularity of $g$ is low. For example, for
the case ${\bf f}=0$ and $u\in H^{2}(\Omega)$, as long as the regularity of
$g|_{K}$ on each element $K\in{\mathcal{T}}$, $t_{K}>0$, we can still get
order $h$ convergence in (6.3) for the $RT_{0,N}\times S_{1,D}$ approximation.
For the $BDM_{1,N}\times S_{2,D}$ approximation, assuming that $u\in
H^{3}(\Omega)$ and $(\nabla u-{\bf f})|_{K}\in H^{3}(K)^{d}$, we only requires
$g|_{K}\in H^{1}(K)$ to get the optimal convergence.

###  6.2. A posteriori error analysis

Define the global a posteriori error estimator:

(6.5)
$\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})=\Big{\\{}\|\alpha^{-1/2}h(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{2,h})\|^{2}_{0}+\|\alpha^{1/2}({\bf
f}-\nabla
u_{2,h}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{2,h})\|^{2}_{0}\Big{\\}}^{1/2},$

and its local indicator

$\displaystyle\eta_{2,K}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$
$\displaystyle=$
$\displaystyle\Big{\\{}\|\alpha^{-1/2}h_{K}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{2,h})\|^{2}_{0,K}+\|\alpha^{1/2}({\bf
f}-\nabla
u_{2,h}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{2,h})\|^{2}_{0,K}\Big{\\}}^{1/2}.$

######  Theorem 6.3.

Let $(\mbox{\boldmath$\sigma$},u)$ be the solution of (1.1) and
$(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$ be the solution of problem (6.1),
respectively. Assume that Assumption 4.1 on the coefficients and Assumption
4.3 (QMA) are true, then there exists positive constants $C$ independent of
$\alpha$ and the mesh-size such that the following reliability holds:

(6.6)
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}\leq
C\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h}).$

###### Proof.

Let ${\bf E}_{2}=\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h}$ and
$e_{2}=u-u_{2,h}$. In (3.20), let $\theta|_{K}=h^{2}_{K}$ and
$v_{h}=I_{rcl}e_{2}$, we have

$|\\!|\\!|({\bf E}_{2},e_{2})|\\!|\\!|^{2}_{2}=({\bf
f}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{2,h}-\nabla u_{2,h},{\bf
E}_{2}+\alpha\nabla(e_{2}-I_{rcl}e_{2}))+(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{2,h},\alpha^{-1}h^{2}\nabla\cdot{\bf
E}_{2}+2(e_{2}-I_{rcl}e_{2})).$

Applying Cauchy-Schwarz and triangle inequalities, we get

$\displaystyle|\\!|\\!|({\bf E}_{2},e_{2})|\\!|\\!|^{2}_{2}$
$\displaystyle\leq$ $\displaystyle\|\alpha^{1/2}({\bf
f}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{2,h}-\nabla
u_{2,h})\|_{0}(\|\alpha^{-1/2}{\bf
E}_{2}\|_{0}+\|\alpha^{1/2}\nabla(e_{2}-I_{rcl}e_{2})\|_{0})$
$\displaystyle+\|h\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{2,h})\|_{0}(\|h\alpha^{-1/2}\nabla\cdot{\bf
E}_{2}\|_{0}+2\|h^{-1}\alpha^{1/2}(e_{2}-I_{rcl}e_{2})\|_{0}.$

By the robust Clément interpolation result (4.1) under Assumption 4.3. we have
the following robust results,

$\|\alpha^{1/2}\nabla(e_{2}-I_{rcl}e_{2})\|_{0}\leq C\|\alpha^{1/2}\nabla
e_{2}\|_{0}\quad\mbox{and}\quad\|\alpha^{1/2}(e_{2}-I_{rcl}e_{2})\|_{0}\leq
C\|h\alpha^{1/2}\nabla e_{2}\|_{0}.$

Substitute these two robust results into (6.2), we get

$\displaystyle|\\!|\\!|({\bf E}_{2},e_{2})|\\!|\\!|^{2}_{2}$
$\displaystyle\leq$ $\displaystyle C(\|\alpha^{1/2}({\bf
f}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{2,h}-\nabla
u_{2,h})\|_{0}+\|h\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{2,h})\|_{0})$
$\displaystyle(\|\alpha^{-1/2}{\bf
E}_{2}\|_{0}+\|h\alpha^{-1/2}\nabla\cdot{\bf E}_{2}\|_{0}+\|\alpha^{1/2}\nabla
e_{2}\|_{0})$ $\displaystyle\leq$ $\displaystyle
C\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})|\\!|\\!|({\bf
E}_{2},e_{2})|\\!|\\!|_{2}.$

The robust reliability result (6.6) is proved. ∎

Then using similar techniques of the proof of Theorem 5.4, we give the
efficiency for the error estimators or indicators of
$\eta_{2,K}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$.

######  Theorem 6.4.

Let $(\mbox{\boldmath$\sigma$},u)$ be the solution of (1.1) and
$(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$ be the solution of problem (6.1),
respectively.. Assume that Assumption 4.1 on the coefficients is true, we have
the following efficiency:

(6.8)
$\eta_{2,K}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})\leq\sqrt{2}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2,K},\quad\forall
K\in{\mathcal{T}}.$

###  6.3. Comparison with the mesh-weighted LSFEM

We have a corresponding mesh-weighted least-squares method for the second
augmented mixed formulation. Define the mesh-weighted least-squares functional
as

(6.9) $J_{h}(\mbox{\boldmath$\tau$},v;{\bf f},g):=\|\alpha^{1/2}\nabla
v+\alpha^{-1/2}\mbox{\boldmath$\tau$}-\alpha^{1/2}{\bf
f}\|_{0}^{2}+\|h\alpha^{-1/2}(\nabla\cdot\mbox{\boldmath$\tau$}-g)\|_{0}^{2},\quad(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

Then the mesh-weighted $L^{2}$-based least-squares minimization problem is:
find $(\mbox{\boldmath$\sigma$},u)\in{\bf X}$, such that

(6.10) $J_{h}(\mbox{\boldmath$\sigma$},u;{\bf
f},g)=\inf_{(\mbox{\boldmath$\tau$},v)\in{\bf
X}}J_{h}(\mbox{\boldmath$\tau$},v;{\bf f},g).$

Equivalently, it can be written in a weak form as: find
$(\mbox{\boldmath$\sigma$},u)\in{\bf X}$, such that

(6.11) $b_{h}((\mbox{\boldmath$\sigma$},u),(\mbox{\boldmath$\tau$},v))=({\bf
f},\mbox{\boldmath$\tau$}+\alpha\nabla
v)+(h^{2}\alpha^{-1}g,\nabla\cdot\mbox{\boldmath$\tau$}),\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf
X},$

where the mesh-weighted least-squares bilinear form $b_{h}$ is defined as
follows:

(6.12)
$b_{h}((\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v))=(\mbox{\boldmath$\chi$}+\nabla
w,\mbox{\boldmath$\tau$}+\alpha\nabla
v)+(h^{2}\alpha^{-1}\nabla\cdot\mbox{\boldmath$\chi$},\nabla\cdot\mbox{\boldmath$\tau$}),\quad\forall(\mbox{\boldmath$\chi$},w),(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

This kind of least-squares method is called the weighted $L^{2}$-discrete-
least-squares principle; see Section 5.6.1 of [5].

Using the same approximation as the second augmented mixed formulation, the
mesh-weighted least-squares finite element method is: find
$(\mbox{\boldmath$\sigma$}_{h}^{hls},u_{h}^{hls})\in\Sigma_{h,N}\times
V_{h,D}$, such that

(6.13) $J_{h}(\mbox{\boldmath$\sigma$}^{hls}_{h},u^{hls}_{h};{\bf
f},g)=\inf_{(\mbox{\boldmath$\tau$},v)\in\Sigma_{h,N}\times
V_{h,D}}J_{h}(\mbox{\boldmath$\tau$},v;{\bf f},g).$

Or, equivalently, find
$(\mbox{\boldmath$\sigma$}_{h}^{hls},u_{h}^{hls})\in\Sigma_{h,N}\times
V_{h,D}$, such that,

(6.14)
$b_{2}(\mbox{\boldmath$\sigma$}_{h}^{hls},u_{h}^{hls},(\mbox{\boldmath$\tau$},v))=(\alpha^{-1}{\bf
f},\mbox{\boldmath$\tau$}+\alpha\nabla
v)+(h^{2}\alpha^{-1}g,\nabla\cdot\mbox{\boldmath$\tau$}),\quad\forall(\mbox{\boldmath$\tau$},v)\in\Sigma_{h,N}\times
V_{h,D}.$

The mathematical theory of the mesh-weighted LSFEM (6.13) or (6.14) is much
less satisfactory. Most importantly, we only have the following quasi-norm
equivalence with $C_{1}$ and $C_{2}$ independent of the mesh size:

(6.15)
$C_{1}h_{min}^{2}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{2}^{2}\leq
J_{h}(\mbox{\boldmath$\tau$},v;0,0)\leq
C_{2}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{2}^{2},\quad\forall(\mbox{\boldmath$\tau$},v)\in{\bf
X},$

where $h_{min}=\min\\{h_{K},K\in{\mathcal{T}}\\}$. The coercivity can be
easily derived from (5.14) and the definition of the norm. With only (6.15)
available, we can not expect mesh-independent a priori and a posteriori error
estimates in the standard norm $|\\!|\\!|\cdot|\\!|\\!|_{2}$.

Another way to establish the analysis for the mesh-weighted LSFEM (6.14) is to
adopt the non-standard least-squares norm. Define

(6.16)
$|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{hls}=J_{h}^{1/2}(\mbox{\boldmath$\tau$},v;0,0)=(\|\alpha^{1/2}\nabla
v+\alpha^{-1/2}\mbox{\boldmath$\tau$}\|_{0}^{2}+\|h\alpha^{-1/2}\nabla\cdot\mbox{\boldmath$\tau$}\|_{0}^{2})^{1/2},\quad(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

It is easy to see that

$|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{hls}\leq\sqrt{2}|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{2}.$

Then we can establish the following best approximation a priori error estimate
for the least-squares induced norm
$|\\!|\\!|(\mbox{\boldmath$\tau$},v)|\\!|\\!|_{hls}$:

(6.17)
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{hls},u-u_{h}^{hls})|\\!|\\!|_{hls}=\inf_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D}}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h})|\\!|\\!|_{hls}\leq\sqrt{2}\inf_{(\mbox{\boldmath$\tau$}_{h},v_{h})\in\Sigma_{h,N}\times
V_{h,D}}|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\tau$}_{h},u-v_{h})|\\!|\\!|_{2}.$

Then we can get similar convergence results as in Theorem 6.1.

Let $(\mbox{\boldmath$\sigma$}_{a},u_{a})\in{\bf X}$, we can define the
following mesh-weighted least-squares a posteriori error estimator:

(6.18)
$\eta_{hls}(\mbox{\boldmath$\sigma$}_{a},u_{a}):=\Big{(}\|h\alpha^{-1/2}(g-\nabla\cdot\mbox{\boldmath$\sigma$}_{a})\|^{2}_{0}+\|\alpha^{1/2}({\bf
f}-\nabla
u_{a}-\alpha^{-1}\mbox{\boldmath$\sigma$}_{a})\|^{2}_{0}\Big{)}^{1/2}=J_{h}(\mbox{\boldmath$\sigma$}_{a},u_{a};{\bf
f},g)^{1/2}.$

Let ${\bf E}_{a}=\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{a}$ and
$e_{a}=u-u_{a}$. We have the following identity:

(6.19)
$\displaystyle\eta_{hls}^{2}(\mbox{\boldmath$\sigma$}_{a},u_{a})=J_{h}(\mbox{\boldmath$\sigma$}_{a},u_{a};{\bf
f},g)=J_{h}({\bf
E}_{a},e_{a};0,0)=|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{a},u-u_{a})|\\!|\\!|_{hls}^{2}.$

Thus, if we are satisfied with the non-traditional mesh-weighted least-squares
norm $|\\!|\\!|\cdot|\\!|\\!|_{hls}$, then we can use
$\eta_{hls}^{2}(\mbox{\boldmath$\sigma$}_{h}^{hls},u_{h}^{hls})$ as the a
posteriori error estimator for the method (6.14).

In summary, we have robust and local optimal a priori error estimates for the
mesh-weighted augmented mixed formulation with less regularity requirement on
$g$. The mesh-weighted least-squares a posteriori error estimator is a robust
error estimator for the method. This is an improvement for the mesh-weighted
least-squares formulation where neither the a priori nor a posteriori error
estimates are mesh-size and coefficient robust with respect to the standard
norms.

######  Remark 6.5.

The mesh-weighted least-squares formulation can be viewed as a practical
version of the $H^{-1}$-least-squares method, see [8]. Define the $H^{-1}$
least-squares functional as

(6.20) $J_{-1}(\mbox{\boldmath$\tau$},v;{\bf f},g):=\|\alpha^{1/2}\nabla
v+\alpha^{-1/2}\mbox{\boldmath$\tau$}-\alpha^{1/2}{\bf
f}\|_{0}^{2}+\|\alpha^{-1/2}(\nabla\cdot\mbox{\boldmath$\tau$}-g)\|_{-1}^{2},\quad(\mbox{\boldmath$\tau$},v)\in{\bf
X}.$

A multilevel preconditioner of the Laplace operator is proposed to replace the
$H^{-1}$ norm in [8]. To avoid the multilevel computation, a very coarse
replacement is to replace the $H^{-1}$-norm with the mesh-weighted norm. Of
course, we can only get a mesh-dependent norm-equivalence (6.15).

##  7\. Numerical Experiments

In this section, we present serval numerical experiments to verify our
findings in previous sections. Our main test problem is the interface problem
with discontinuous coefficients from [37, 26, 22]. Let $\Omega=(-1,1)^{2}$ and
let $\widetilde{u}(r,\theta)=r^{\gamma}\mu(\theta)$ in polar coordinates with

(7.1)
$\mu(\theta)=\left\\{\begin{array}[]{ll}\cos((\pi/2-\phi)\gamma)\cdot((\theta-\pi/2+\rho)\gamma)&\mbox{
if }0\leq\theta\leq\pi/2,\\\
\cos(\rho\gamma)\cdot\cos((\theta-\pi+\phi)\gamma)&\mbox{ if
}\pi/2\leq\theta\leq\pi,\\\
\cos(\rho\phi)\cdot\cos((\theta-\pi-\rho)\gamma)&\mbox{ if }\pi\leq\theta\leq
3\pi/2,\\\ \cos((\pi/2-\rho)\gamma)\cdot((\theta-3\pi/2-\phi)\gamma)&\mbox{ if
}3\pi/2\leq\theta\leq 2\pi.\end{array}\right.$

The coefficient $\alpha$ is

$\alpha(x)=\left\\{\begin{array}[]{ll}R&\mbox{ in }(0,1)^{2}\cup(-1,0)^{2},\\\
1&\mbox{ in }\Omega\backslash([0,1]^{2}\cup[-1,0]^{2}).\end{array}\right.$

We choose the numbers $\gamma$, $\rho$, $\phi$, and $R$ such that
$\nabla\cdot(\alpha\nabla\widetilde{u})=0$ in $\Omega$. The following
nonlinear relations of $\gamma$, $\rho$, $\phi$, and $R$ can be found in [26]:

(7.2)
$\left\\{\begin{array}[]{l}R=-\tan((\pi/2-\phi)\gamma)\cdot\cot(\rho\gamma),\\\
1/R=-\tan(\rho\gamma)\cdot\cot(\phi\gamma),\\\
R=-\tan(\phi\gamma)\cdot\cot((\pi/2-\rho)\gamma),\\\ 0<\gamma<2,\\\
\max(0,\pi\gamma-\pi)<2\gamma\rho<\min(\pi\gamma,\pi),\\\
\max(0,\pi-\pi\gamma)<-2\gamma\phi<\min(\pi,2\pi-\pi\gamma).\end{array}\right.$

The function $\widetilde{u}(r,\theta)\in H^{1+\gamma-\epsilon}(\Omega),$ for
any $\epsilon>0$. The regularity index $\gamma$ is less than $1$, The function
$\widetilde{u}$ is singular at the origin. We give serval examples of the
coefficients $\alpha$ and the corresponding numbers $\gamma$, $\rho$, and
$\phi$ in Table 1. These numbers can be computed by solving (7.2) using
Newton’s method. Some of them can also be found in [26].

Table 1. Numbers of $\phi,\gamma$ and $R$ with $\rho=\pi/4.$ | $\phi$ | $\gamma$ | $R$
---|---|---|---
$\mbox{Data}1$ | $-2.3561944901923448$ | $0.50$ | $5.82842712474619$
$\mbox{Data}2$ | $-7.06858347058882$ | $0.20$ | $39.8634581884533$
$\mbox{Data}3$ | $-9.68657734859297$ | $0.15$ | $71.3848801304590$
$\mbox{Data}4$ | $-14.92256510455152$ | $0.10$ | $161.447638797588$

Let $u:=\widetilde{u}+u_{0}$ with
$u_{0}(x,y)=\left\\{\begin{array}[]{ll}x+1,&x\leq 0,\\\
1&x>0.\end{array}\right.$ Then $u$ is the solution of the following
generalized Darcy’s problem

(7.3)
$\left\\{\begin{array}[]{lllll}\nabla\cdot\mbox{\boldmath$\sigma$}&=&0&\mbox{in
}\Omega,\\\\[2.84526pt] \alpha\nabla u+\mbox{\boldmath$\sigma$}&=&\alpha\nabla
u_{0}&\mbox{in }\Omega.\end{array}\right.$

This will be our main problem to do numerical tests.

We propose the following criteria to test the robustness of the methods with
respect to $\alpha$. Let $(\mbox{\boldmath$\sigma$}_{h},u_{h})$ be the
numerical solutions computed by the numerical methods (augmented mixed methods
and LSFEMs) discussed in the paper; we want to compare the ratio of the error
in the energy norm and its corresponding a posteriori error estimator. For the
first kind augmented mixed methods and the $L^{2}$-LSFEM, we measure the error
in
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h})|\\!|\\!|_{1}$
(the definition in (3.14) with $\theta=1$) and the a posteriori error
estimator is the one based on the $L^{2}$-LS-functional
$\eta_{ls}(\mbox{\boldmath$\sigma$}_{h},u_{h})=J^{1/2}(\mbox{\boldmath$\sigma$}_{h},u_{h};0,0)$.
For the second kind augmented mixed methods and mesh-weighted LSFEM, we
measure the error in
$|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h})|\\!|\\!|_{2}$
(the definition in (3.14) with $\theta=h_{K}^{2}$) and the a posteriori error
estimator is the one based on the mesh-weighted LS-functional
$\eta_{hls}(\mbox{\boldmath$\sigma$}_{h},u_{h})=J^{1/2}_{h}(\mbox{\boldmath$\sigma$}_{h},u_{h};0,0)$.

If the a posteriori error estimator is robust, then for different $\alpha$,
the ratio or the so-called effectivity index

(7.4) $\mbox{eff-
index}=\displaystyle\frac{|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h})|\\!|\\!|}{\eta}$

should be a constant.

The standard adaptive finite element method is based on the following loop:
SOLVE, ESTIMATE, MARK, and REFINE. We use the Dörfler’s bulk marking strategy.
The relative error is computed by $\mbox{rel-
err}=|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h},u-u_{h})|\\!|\\!|/|\\!|\\!|(\mbox{\boldmath$\sigma$},u)|\\!|\\!|$.

We do not seek numerical examples to check the robust a priori error
estimates.

###  7.1. Convergence tests for adaptive augmented mixed methods with pure
Dirichlet boundary conditions

In this subsection, we use Data4 in Table 1. The main purpose is to show the
convergence history and the final adaptive mesh for the augmented mixed
methods. We choose the pure Dirichlet boundary condition. The parameter in
Dörfler’s bulk marking strategy is $0.3$ and the stopping criteria is
$\mbox{rel-err}\leq 0.010$.

In Figure 1, we present the numerical test of the first adaptive augmented
mixed method (5.1), with $RT_{0,N}\times S_{1,D}$ being the finite element
space. On the left of Figure 1, we show the reference line
$\mbox{Dofs}^{-1/2}$, the decay of the error estimator, and the error measured
in $|\\!|\\!|(\cdot,\cdot)|\\!|\\!|_{1}$. On the right of Figure 1, the final
mesh generated by $\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$ is
presented after $75$ loops of bisection. The final DOF is $4921$. The
convergence and the final mesh are both optimal, and the final mesh is similar
to the results presented in [26, 22].

Figure 1. Adaptive convergence results and final mesh for the first adaptive
augmented mixed method (5.1)

In Figure 2, we present the numerical test of the second (mesh-weighted)
adaptive augmented mixed method (6.1), with $RT_{0,N}\times S_{1,D}$ being the
finite element space. On the left of Figure 2, we show the reference line
$\mbox{Dofs}^{-1/2}$, the decay of the error estimator, and the error measured
in $|\\!|\\!|(\cdot,\cdot)|\\!|\\!|_{2}$. On the right of Figure 2, the final
mesh generated by $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$ is
presented after $86$ loops of bisection. The final DOF is $4621$. The
convergence and the final mesh are both optimal, and the final mesh is similar
to the results presented in [26, 22].

Figure 2. Adaptive convergence results and final mesh for the second (mesh-
weighted) adaptive augmented mixed method using $RT_{0,N}\times S_{1,D}$

In Figure 3, we present the numerical test of the second (mesh-weighted)
adaptive augmented mixed method (6.1), with $BDM_{1,N}\times S_{2,D}$. The
reference line in this case is $\mbox{Dofs}^{-1}$. The final DOFs is $1997$
after $49$ times of bisection.

Figure 3. Adaptive convergence results and final mesh for the second (mesh-
weighted) adaptive augmented mixed method using $BDM_{1,N}\times S_{2,D}$

###  7.2. Robustness tests for adaptive augmented mixed methods with pure
Dirichlet boundary conditions

In this subsection, we present the numerical results of the effectivity index
for adaptive augmented mixed methods with pure Dirichlet boundary conditions
for different $\alpha$ to check the robustness of the methods. The same
marking strategy and stopping criteria of the previous subsection are used.
From Table 2 to Table 4, we show the eff-indexes for different $\alpha$ for
the first and second augmented methods. As seen from these tables, the eff-
index is almost a constant for different $\alpha$; this verifies least-squares
a posteriori error estimators for both augmented mixed methods are robust.

Table 2. The first adaptive augmented mixed method for different $\alpha$, pure Dirichlet BC: eff-index, number of refinements $k$, number of elements $n$, the final $\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$, and the final error $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|_{1}$ | eff-index | $k$ | $n$ | $\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$ | $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|_{1}$
---|---|---|---|---|---
$\mbox{Data}1$ | $1.0006$ | $34$ | $15824$ | $0.0250$ | $0.0250$
$\mbox{Data}2$ | $1.0075$ | $64$ | $7216$ | $0.0621$ | $0.0616$
$\mbox{Data}3$ | $1.0179$ | $71$ | $4648$ | $0.0842$ | $0.0828$
$\mbox{Data}4$ | $1.0605$ | $75$ | $2448$ | $0.1340$ | $0.1264$
Table 3. The second (mesh-weighted) adaptive augmented mixed method using $RT_{0,N}\times S_{1,D}$ for different $\alpha$, pure Dirichlet BC: eff-index, number of refinements $k$, number of elements $n$, the final $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$, and the final error $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$ | eff-index | $k$ | $n$ | $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$ | $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$
---|---|---|---|---|---
$\mbox{Data}1$ | $0.9963$ | $35$ | $15184$ | $0.0254$ | $0.0255$
$\mbox{Data}2$ | $0.9966$ | $70$ | $7088$ | $0.0616$ | $0.0618$
$\mbox{Data}3$ | $0.9949$ | $80$ | $4524$ | $0.0823$ | $0.0827$
$\mbox{Data}4$ | $0.9847$ | $86$ | $2300$ | $0.1236$ | $0.1256$
Table 4. The second (mesh-weighted) adaptive augmented mixed method using $BDM_{1,N}\times S_{2,D}$ for different $\alpha$, pure Dirichlet BC: eff-index, number of refinements $k$, number of elements $n$, the final $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$, and the final error $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$ | eff-index | $k$ | $n$ | $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$ | $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$
---|---|---|---|---|---
$\mbox{Data}1$ | $1.0728$ | $20$ | $320$ | $0.0262$ | $0.0244$
$\mbox{Data}2$ | $1.0792$ | $43$ | $416$ | $0.0659$ | $0.0611$
$\mbox{Data}3$ | $1.1068$ | $47$ | $380$ | $0.0929$ | $0.0839$
$\mbox{Data}4$ | $1.1014$ | $49$ | $396$ | $0.1383$ | $0.1256$

###  7.3. Convergence and robustness tests for the adaptive $L^{2}$-LSFEM with
pure Dirichlet boundary conditions

In this subsection, we present the numerical results of convergence result of
the adaptive $L^{2}$-LSFEM (5.10) with pure Dirichlet boundary conditions and
the effectivity index for different $\alpha$.

In Figure 4, we present the numerical test of the adaptive $L^{2}$-LSFEM
(5.10) with pure Dirichlet boundary condition with Data4. The convergence and
the final mesh are both optimal with enough mesh grids, and the final mesh is
similar to the results presented in [26, 22].

Figure 4. Adaptive convergence results and final mesh for the $L^{2}$-LSFEM
with pure Dirichlet BC

We show the eff-indexes for different $\alpha$ for the $L^{2}$-LSFEM with pure
Dirichlet BC in Table 5. As seen from these tables, the eff-index is almost a
constant for different $\alpha$, this verifies that for this special case that
each subdomain has a non-empty Dirichlet boundary condition, the $L^{2}$-LSFEM
is robust.

Table 5. The $L^{2}$-LSFEM for different $\alpha$, pure Dirichlet BC: eff-index, number of refinements $k$, number of elements $n$, the final $\eta_{ls}(\mbox{\boldmath$\sigma$}^{ls}_{h},u^{ls}_{h})$, and the final error $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}^{ls}_{h},u-u^{ls}_{h})|\\!|\\!|_{1}$ | eff-index | $k$ | $n$ | $\eta_{ls}(\mbox{\boldmath$\sigma$}^{ls}_{h},u^{ls}_{h})$ | $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}^{ls}_{h},u-u^{ls}_{h})|\\!|\\!|_{1}$
---|---|---|---|---|---
$\mbox{Data}1$ | $1.0019$ | $31$ | $14476$ | $0.0263$ | $0.0262$
$\mbox{Data}2$ | $1.0171$ | $58$ | $7772$ | $0.0628$ | $0.0617$
$\mbox{Data}3$ | $1.0406$ | $63$ | $5400$ | $0.0871$ | $0.0837$
$\mbox{Data}4$ | $1.1216$ | $65$ | $3744$ | $0.1395$ | $0.1244$

###  7.4. Robustness tests for adaptive augmented mixed methods with mixed
boundary conditions

In this subsection, we present the numerical results of the effectivity index
of adaptive augmented mixed methods with mixed boundary conditions with
different $\alpha$ to check the robustness of the methods.

The mixed boundary are chosen as follows:
$\Gamma_{D}=\\{(x,y)\in\partial\Omega;x\in(-1,1),y=-1\\}$ and
$\Gamma_{N}=\partial\Omega\setminus\Gamma_{D}$. The boundary conditions
$g_{D}$ and $g_{N}$ are given by the true solution.

From Table 6 to Table 8, we show the eff-indexes for different $\alpha$ for
the first and second augmented methods. As seen from these tables, the eff-
index is almost a constant for different $\alpha$; this verifies least-squares
a posteriori error estimators for augmented mixed methods are robust for the
general mixed boundary condition case.

Table 6. The first adaptive augmented mixed method for different $\alpha$, mixed BC: eff-index, number of refinements $k$, number of elements $n$, the final $\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$, and the final error $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|_{1}$ | eff-index | $k$ | $n$ | $\eta_{1}(\mbox{\boldmath$\sigma$}_{1,h},u_{1,h})$ | $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{1,h},u-u_{1,h})|\\!|\\!|_{1}$
---|---|---|---|---|---
$\mbox{Data}1$ | $1.0006$ | $37$ | $41031$ | $0.0156$ | $0.0155$
$\mbox{Data}2$ | $1.0058$ | $73$ | $19970$ | $0.0375$ | $0.0373$
$\mbox{Data}3$ | $1.0138$ | $84$ | $13622$ | $0.0497$ | $0.0490$
$\mbox{Data}4$ | $1.0497$ | $93$ | $7605$ | $0.0795$ | $0.0758$
Table 7. The second (mesh-weighted) adaptive augmented mixed method using $RT_{0,N}\times S_{1,D}$ for different $\alpha$, mixed BC: eff-index, number of refinements $k$, number of elements $n$, the final $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$, and the final error $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$ | eff-index | $k$ | $n$ | $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$ | $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$
---|---|---|---|---|---
$\mbox{Data}1$ | $0.9965$ | $39$ | $45355$ | $0.0147$ | $0.0147$
$\mbox{Data}2$ | $0.9966$ | $79$ | $19103$ | $0.0374$ | $0.0375$
$\mbox{Data}3$ | $0.9962$ | $93$ | $12478$ | $0.0492$ | $0.0494$
$\mbox{Data}4$ | $0.9900$ | $108$ | $6046$ | $0.0748$ | $0.0756$
Table 8. The second (mesh-weighted) adaptive augmented mixed method using $BDM_{1,N}\times S_{2,D}$ for different $\alpha$, mixed BC: eff-index, number of refinements $k$, number of elements $n$, the final $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$, and the final error $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$ | eff-index | $k$ | $n$ | $\eta_{2}(\mbox{\boldmath$\sigma$}_{2,h},u_{2,h})$ | $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{2,h},u-u_{2,h})|\\!|\\!|_{2}$
---|---|---|---|---|---
$\mbox{Data}1$ | $1.0735$ | $23$ | $560$ | $0.0154$ | $0.0144$
$\mbox{Data}2$ | $1.0775$ | $51$ | $704$ | $0.0386$ | $0.0358$
$\mbox{Data}3$ | $1.0831$ | $59$ | $670$ | $0.0529$ | $0.0489$
$\mbox{Data}4$ | $1.0553$ | $70$ | $573$ | $0.0800$ | $0.0758$

###  7.5. Non-robustness for $L^{2}$-LSFEM with mixed boundary conditions

We present the non-robustness of the $L^{2}$-LSFEM (5.10) with mixed boundary
conditions. We use the same mixed boundary setting as the previous subsection.

As we can see from Table 9, the eff-index is not a constant. This verifies the
conclusion that the $L^{2}$-LSFEM is not robust with respect to $\alpha$. On
the other hand, as we can see, the eff-index does not change as strongly as
$\alpha$; this is explained in Remark 5.9.

Table 9. The $L^{2}$-LSFEM for different $\alpha$, mixed BC: eff-index, number of refinements $k$, number of elements $n$, the final $\eta_{ls}(\mbox{\boldmath$\sigma$}^{ls}_{h},u^{ls}_{h})$, and the final error $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}^{ls}_{h},u-u^{ls}_{h})|\\!|\\!|_{1}$ | eff-index | $k$ | $n$ | $\eta_{ls}(\mbox{\boldmath$\sigma$}^{ls}_{h},u^{ls}_{h})$ | $|\\!|\\!|(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}^{ls}_{h},u-u^{ls}_{h})|\\!|\\!|_{1}$
---|---|---|---|---|---
$\mbox{Data}1$ | $0.9972$ | $80$ | $14434$ | $0.0259$ | $0.0260$
$\mbox{Data}2$ | $0.8641$ | $91$ | $9542$ | $0.0528$ | $0.0611$
$\mbox{Data}3$ | $0.7079$ | $110$ | $8713$ | $0.0590$ | $0.0833$
$\mbox{Data}4$ | $0.4787$ | $139$ | $11754$ | $0.0596$ | $0.1244$

In Figure 5, we can also see that even for a fixed $\alpha$, the eff-index is
not a constant when mesh is refined, which means it is not even robust with
respect to the mesh-size. In contract, in Figure 4, the eff-index is almost a
constant for a fixed $\alpha$ for the special case of $L^{2}$-LSFEM where the
robust Poincaré inequality holds.

Figure 5. Adaptive convergence results for the $L^{2}$-LSFEM with mixed
boundary conditions with Data3.

###  7.6. Numerical tests for mesh-weighted LSFEMs

For the mesh-weighted LSFEM (6.14), there are no rigorous a priori and a
posteriori error estimates with respect to standard norms. We show the
convergence result of adaptive convergence history for the problem with pure
Dirichlet boundary conditions with Data4 using $RT_{0,N}\times S_{1,D}$ with
$39$ times of refinement in Figure 6. The error measured in norm
$|\\!|\\!|(\cdot,\cdot)|\\!|\\!|_{2}$ (the blue line in Figure 6) does not
decay as the a posteriori error estimator do. Also, the final mesh it obtained
is skewed, which is a non-optimal case as discussed in [26, 22]. This suggests
that this method may not be a good choice for this class of problem.

Figure 6. Adaptive convergence results and final mesh for the mesh-weighted
LSFEM with pure Dirichlet boundary conditions.

##  8\. Final Comments

In this paper, for the generalized Darcy problem, we study a special Galerkin-
Least-Squares method, the augmented mixed finite element method, and its
relationship to the standard least-squares finite element method (LSFEM). One
of the paper’s main contributions is to connect the augmented mixed finite
element method and the bonafide LSFEM and discuss their shared properties and
differences. Both methods share some good properties: both methods are based
on the physically meaningful first-order system. Thus, important physical
qualities can be approximated in their intrinsic spaces; both methods are
coercive and stable, and their finite element discrete problems are coercive
and stable as long as the discrete spaces are subspaces of the abstract spaces
of the variational problems. Thus, no inf-sup condition of the discrete spaces
and mesh size restriction are needed; both methods can use the least-squares
functional as the build-in a posteriori error estimator. On the other hand,
the augmented mixed methods and the LSFEMs have their advantages. For the
augmented mixed finite element methods, we show that the a priori and a
posteriori error estimates are robust with respect to the coefficients of the
problem. In contrast, we discuss the non-robustness of standard least-squares
finite element methods. With the flexibility of being a partial least-squares
method, it is possible that the augmented mixed finite element method can have
better numerical properties than the bona-fide least-squares method. However,
being a partial least-squares method, the augmented mixed method also loses
one main property of the bonafide least-squares method, which may be very
useful for minimization-based methods: the minimization of the least-squares
energy, even though we can always associate a Ritz-minimization variational
principle to the symmetric version of the augmented mixed method.

## References

  * [1] Javier A. Almonacid, Gabriel N. Gatica, and Ricardo Ruiz-Baier.  Ultra-weak symmetry of stress for augmented mixed finite element formulations in continuum mechanics.  Calcolo, 57(2), 2020.
  * [2] Tomas P. Barrios, J. Manuel Cascon, and Maria Gonzalezlez.  A posteriori error analysis of an augmented mixed finite element method for darcy flow.  Computer Methods in Applied Mechanics and Engineering, 283:909–922, 2015.
  * [3] Tomas P. Barrios, Gabriel N. Gatica, Maria Gonzalez, and Norbert Heuer.  A residual based a posteriori error estimator for an augmented mixed finite element method in linear elasticity.  ESAIM: Mathematical Modelling and Numerical Analysis, 40(5):843–869, 2006.
  * [4] Christine Bernardi and Rüdiger Verfürth.  Adaptive finite element methods for elliptic equations with non-smooth coefficients.  Numer. Math., 85:579–608, 2000.
  * [5] Pavel B. Bochev and Max D Gunzburger.  Least-Squares Finite Element Methods.  Applied Mathematical Sciences, 166. Springer, 2009.
  * [6] Daniele Boffi, Franco Brezzi, and Michel Fortin.  Mixed Finite Element Methods and Applications, volume 44 of Springer Series in Computational Mathematics.  Springer, 2013.
  * [7] Dietrich Braess.  Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics.  Cambridge University Press, 2007.
  * [8] James H. Bramble, Raytcho D. Lazarov, and Joseph E. Pasciak.  A least-squares approach based on a discrete minus one inner product for first order systems.  Mathematics of Computation, 66(219):935–055, 1997.
  * [9] Susanne Brenner and Ridgway Scott.  The Mathematical Theory of Finite Element Methods, volume 15 of Texts in Applied Mathematics.  Springer, third edition, 2008.
  * [10] Franco Brezzi, Thomas J.R. Hughes, L. Donatella Marini, and Arif Masud.  Mixed discontinuous galerkin methods for darcy flow.  Journal of Scientific Computing, 84:119–145, 2005.
  * [11] Difeng Cai, Zhiqiang Cai, and Shun Zhang.  Robust equilibrated a posteriori error estimator for higher order finite element approximations to diffusion problems.  Numer. Math., 144(1):1–21, 2020.
  * [12] Difeng Cai, Zhiqiang Cai, and Shun Zhang.  Robust equilibrated error estimator for diffusion problems: Mixed finite elements in two dimensions.  Journal of Scientific Computing, 83(1), 2020.
  * [13] Zhiqiang Cai, Rob Falgout, and Shun Zhang.  Div first-order system LL* (FOSLL*) least-squares for second-order elliptic partial differential equations.  SIAM J. Numer. Anal., 53(1):405–420, 2015.
  * [14] Zhiqiang Cai, Cuiyu He, and Shun Zhang.  Discontinuous finite element methods for interface problems: Robust a priori and a posteriori error estimates.  SIAM J. Numer. Anal., 55:400–418, 2017.
  * [15] Zhiqiang Cai, Cuiyu He, and Shun Zhang.  Improved zz a posteriori error estimators for diffusion problems: Discontinuous element.  Applied Numerical Mathematics, 159:174–189, 2021.
  * [16] Zhiqiang Cai and JaEun Ku.  The $L^{2}$ norm error estimates for the div least-squares methods.  SIAM J. Numer. Anal., 44(4):1721–1734, 2006.
  * [17] Zhiqiang Cai, R. Lazarov, T. Manteuffel, and S. McCormick.  First order system least-squares for second order partial differential equations: Part I.  SIAM J. Numer. Anal., 31:1785–1799, 1994.
  * [18] Zhiqiang Cai, Barry Lee, and Ping Wang.  Least-squares methods for incompressible newtonian fluid flow: linear stationary problems,.  SIAM J. Numer. Anal., 42:843–859, 2004.
  * [19] Zhiqiang Cai, Tom Manteuffel, and Stephen F. McCormick.  First-order system least squares for second-order partial differential equations: Part ii.  SIAM J. Numer. Anal., 34(2):425–454, 1997.
  * [20] Zhiqiang Cai and Gerhard Starke.  Least-squares methods for linear elasticity.  SIAM J. Numer. Anal., 42:826–842, 2004.
  * [21] Zhiqiang Cai, Xiu Ye, and Shun Zhang.  Discontinuous galerkin finite element methods for interface problems: a priori and a posteriori error estimations.  SIAM J. Numer. Anal., 49(5):1761–1787, 2011.
  * [22] Zhiqiang Cai and Shun Zhang.  Recovery-based error estimator for interface problems: Conforming linear elements.  SIAM J. Numer. Anal., 47(3):2132–2156, 2009.
  * [23] Zhiqiang Cai and Shun Zhang.  Recovery-based error estimators for interface problems: Mixed and nonconforming finite elements.  SIAM J. Numer. Anal., 48(1):30–52, 2010.
  * [24] Zhiqiang Cai and Shun Zhang.  Robust equilibrated residual error estimator for diffusion problems: conforming elements.  SIAM J. Numer. Anal., 50:151–170, 2012.
  * [25] Jessika Camano, Gabriel N. Gatica, Ricardo Oyarzua, and Giordano Tierra.  An augmented mixed finite element method for the navier–stokes equations with variable viscosity.  SIAM J. Numer. Anal., 54(2):1069–1092, 2016.
  * [26] Zhiming Chen and Shibin Dai.  On the efficiency of adaptive finite element methods for elliptic problems with discontinuous coefficients.  SIAM J. Sci. Comput., 24(2):443–462, 2002.
  * [27] Philippe G. Ciarlet.  The Finite Element Method for Elliptic Problems, volume 40 of Classics in Applied Mathematics.  SIAM, 2002.
  * [28] Eligio Colmenares, Gabriel N. Gatica, and Ricardo Oyarzúa.  An augmented fully-mixed finite element method for the stationary boussinesq problem.  Calcolo, 54:167–205, 2017.
  * [29] Maicon R. Correa and Abimael F D Loula.  Unconditionally stable mixed finite element methods for darcy flow.  Comput. Methods Appl. Mech. Engrg., 197:1525–1540, 2008.
  * [30] T. Dupont and R. Scott.  Polynomial approximation of functions in sobolev spaces.  Math. Comp., 34:441–463, 1980.
  * [31] Alexandre Ern and Jean-Luc Guermond.  Finite Elements I: Approximation and Interpolation, volume 72 of Texts in Applied Mathematics.  Springer, 2021.
  * [32] Leonardo E. Figuero, Gabriel N. Gatica, and Antonio Marquez.  Augmented mixed finite element methods for the stationary stokes equations.  SIAM J. Sci. Comput., 31(2):1082–1119, 2008.
  * [33] L. P. Franca.  New Mixed Finite Element Methods.  PhD thesis, Stanford University, 1987.
  * [34] L. P. Franca and T.J.R. Hughes.  Two classes of finite element methods.  Computer Methods in Applied Mechanics and Engineering, 69:89–129, 1988.
  * [35] Gabriel N. Gatica.  Analysis of a new augmented mixed finite element method for linear elasticity allowing rt0-p1-p0 approximation.  Math. Model. Numer. Anal., 40(1):1–28, 2006.
  * [36] Bo-nan Jiang.  The Least-Squares Finite Element Method Theory and Applications in Computational Fluid Dynamics and Electromagnetics.  Scientifc Computation. Springer, 1998.
  * [37] R. Bruce kellogg.  On the poisson equation with intersecting interfaces on the poisson equation with intersecting interfaces.  Applicable Analysis, 4(2):101–129, 1974.
  * [38] Qunjie Liu and Shun Zhang.  Adaptive flux-only least-squares finite element methods for linear transport equations.  Journal of Scientific Computing, 84:26, 2020.
  * [39] Qunjie Liu and Shun Zhang.  Adaptive least-squares finite element methods for linear transport equations based on an H(div) flux reformulation.  Comput. Methods Appl. Mech. Engrg., 366:113041, 2020.
  * [40] Arif Masud and Thomas J.R. Hughes.  A stabilized mixed finite element method for darcy flow.  Computer Methods in Applied Mechanics and Engineering, 191:4341–4370, 2002.
  * [41] M. Petzoldt.  A posteriori error estimators for elliptic equations with discontinuous coefficients.  Adv. Comput. Math., 16:47–75, 2002.
  * [42] Weifeng Qiu and Shun Zhang.  Adaptive first-order system least-squares finite element methods for second order elliptic equations in non-divergence form.  SIAM J. Numer. Anal., 58(6):3286–3308, 2020.
  * [43] P. A. Raviart and J. M. Thomas.  A mixed finite element method for second order elliptic problems.  In I. Galligani and E. Magenes, editors, Mathematical Aspects of the Finite Element Method, volume 606 of Lectures Notes in Mathematics,. Springer, 1977.
  * [44] Jinchao Xu and Ludmil Zikatanov.  Some observations on Babuška and Brezzi theories.  Numer. Math., 94:195–202, 2003.
  * [45] Shun Zhang.  Robust and local optimal a priori error estimates for interface problems with low regularity: Mixed finite element approximations.  Journal of Scientific Computing, 84(40), 2020.
  * [46] Shun Zhang.  A simple proof of coerciveness of first-order system least-squares methods for general second-order elliptic pdes.  Computers & Mathematics with Applications, pages 98–104, 2023\.