A two-grid Adaptive Finite Element Method for the Dirichlet Boundary Control Problem Governed by Stokes Equation

Thirupathi Gudi Department of Mathematics, Indian Institute of Science, Bangalore - 560012, India [email protected] and Ramesh Ch. Sau Department of Mathematics, Chinese University of Hong Kong, Hong Kong [email protected]

Abstract.

In this article, we derive a posteriori error estimates for the Dirichlet boundary control problem governed by Stokes equation. An energy-based method has been deployed to solve the Dirichlet boundary control problem. We employ an inf-sup stable finite element discretization scheme by using $\mathbf{P}_{1}$ elements(in the fine mesh) for the velocity and control variable and $P_{0}$ elements(in the coarse mesh) for the pressure variable. We derive an a posteriori error estimator for the state, adjoint state, and control error. The control error estimator generalizes the standard residual type estimator of the unconstrained Dirichlet boundary control problems, by additional terms at the contact boundary addressing the non-linearity. We prove the reliability and efficiency of the estimator. Theoretical results are illustrated by some numerical experiments.

Key words and phrases:

PDE-constrained optimization; Dirichlet boundary control problem; Finite element method; Error bounds; Stokes equation

1991 Mathematics Subject Classification:

65N30; 65N15; 65N12; 65K10

1. Introduction

The study of the optimal control problems governed by partial differential equations has been a major research area in applied mathematics and its allied areas. The optimal control problem consists of finding an optimal control variable that minimizes a cost functional subject to a partial differential equation satisfied by the optimal control and an optimal state. There are many results on the finite element analysis of optimal control problems, see for example [21, 5, 24, 4]. In an optimal control problem, the control can act on the system through either a boundary condition or through an interior force. In the latter case, the control is distributed and in the former case, the control is said to be the boundary control. The choice of the boundary condition leads to several types of boundary controls e.g., Dirichlet, Neumann or Robin boundary control. We refer to [21, 28, 24, 9, 12] for the work related to distributed control and to [5, 4, 9, 12] for the work on Neumann boundary control problem.

The Dirichlet boundary control problems are important in application areas, but the problems can be difficult to analyze mathematically because in some cases the control itself solves some PDE and in other cases the physical domain has a non-smooth boundary. The study of Dirichlet control problems posed on polygonal domains can be traced back to [6], where a control constrained problem governed by a semilinear elliptic equation posed in a convex polygonal domain is studied. There are various approaches proposed in the literature. One of the popular approach is to seek control from the $L^{2}(\Gamma)$ -space. In this case, the state equation has to be understood in a ultra weak sense, since the Dirichlet boundary condition is only in $L^{2}(\Gamma)$ . This ultra weak formulation is easy to implement and usually results in optimal controls with low regularity. Especially, when the problem is posed on a convex polygonal domain, the control $y$ vanishes on the corners and is thus continuous. This is because, it is determined by the normal derivative of the adjoint state, whereas in a nonconvex polygonal domain the control may have singularity around the corner point, for more details one can see [1]. An other approach, as in [7], is the Robin boundary penalization which transforms the Dirichlet control problem into a Robin boundary control problem.

One other popular approach is to find controls from the energy space, i.e., $H^{1/2}(\Gamma),$ we refer [30] for this approach. In [30], the Steklov-Poincaré operator was used to define the cost functional with the help of a harmonic extension of the given boundary data. The Steklov-Poincaré operator transforms the Dirichlet data into a Neumann data by using harmonic extension of the Dirichlet data; but this type of abstract operator may cause some difficulties in numerical implementation. A boundary element method for this problem is proposed and analyzed in [29].

Given any $y$ in the Sobolev space $H^{1/2}(\Gamma)$ , we can construct a function $u_{y}$ in the Sobolev space $H^{1}(\Omega)$ that is harmonic in $\Omega$ and agrees with $y$ on the boundary $\Gamma$ . This function $u_{y}$ is called a harmonic extension of $y$ . The norm of $y$ in $H^{1/2}(\Gamma)$ can be equivalently written as the norm of the gradient of $u_{y}$ in $L^{2}(\Omega)$ , i.e., $|y|_{H^{1/2}(\Gamma)}\equiv\left\|\nabla u_{y}\right\|_{0,\Omega}.$ This relation suggests that we can penalize the control $y$ in the $H^{1}(\Omega)$ space by adding the term $\left\|\nabla y\right\|_{0,\Omega}^{2}$ to the cost functional. Therefore, we consider the following optimization problem:

\text{min}~{}J(u,y)=\frac{1}{2}\left\|u-u_{d}\right\|_{0,\Omega}^{2}+\frac{\lambda}{2}\left\|\nabla y\right\|_{0,\Omega}^{2}.

The paper [10] proposed a novel method for the Dirichlet boundary control problem using the above cost functional. The method is based on the energy space, where the control belongs to the Sobolev space $H^{1}(\Omega)$ . This avoids the use of the Steklov-Poincaré operator and makes the method computationally efficient. The paper [10] only considered the unconstrained case, while the constrained case was analyzed in [19]. The paper [26] presented a similar method based on the energy space. The paper [33] obtained a sharp convergence rate for the energy space method.

The literature on Stokes Dirichlet control problem considers two types of control. One is the tangential control, i.e., the control acts only in the tangential direction of the boundary (see [17]). The paper [17] used a hybridized discontinuous Galerkin (HDG) method to solve a tangential Dirichlet boundary control problem with an $L^{2}$ penalty on the boundary control, without any constraints on the control. The other is the zero flux control, which means the control has no normal component to the boundary (i.e., $\int_{\partial\Omega}\mathbf{y}\cdot\mathbf{n}=0$ ) [18]. The paper [18] studied two different boundary control regularization terms in the cost functionals: the $L^{2}$ norm and the energy space seminorm. The zero flux condition is a natural consequence of the incompressibility condition and the Dirichlet boundary condition in the PDE. Therefore, the authors either chose the tangential control as the first option or imposed the zero flux condition as a constraint in the space as the second option. Many papers on Navier-Stokes Dirichlet control problem also used either tangential control or zero flux control, for example, see [15, 22]. The zero flux condition on the control affects the regularity of the control as discussed in [18]. To address this issue, we introduced the Stokes equation with mixed boundary conditions and the control acts only on the Dirichlet boundary. Our control is more general and has both tangential and normal components. We also added constraints on the control. As a result, the optimal control satisfies a simplified Signorini problem.

In this article, we propose, analyze, and test a new a posteriori error estimator for the control error. In order to derive and analyze the error estimator for the control variable, we adopt the framework presented in [27] for Signorini problem, because the control satisfies the Signorini problem. The discrete problem consists of a discrete variational inequality for the approximate control variable and the estimator is designed for controlling its energy error. The estimator reduces to the standard residual estimator for elliptic problem, if no contact occurs. The contributions by the estimator addressing the nonlinearity are related to the contact stresses, the complementarity condition. We prove reliability and efficiency of the estimator and ensuring the equivalence with the error up to oscillation terms. A key ingredient of this approach is the so-called Galerkin functional. It is a modification of the residual with respect to the corresponding linear problem with the help of a suitable approximation of the Lagrange multiplier and thus, may be seen as the residual of a linear auxiliary problem. The correction by the Lagrange multiplier is crucial for sharpness of the upper bounds in the actual contact regions. The theoretical results are corroborated by a variety of numerical tests

The rest of the article is organized as follows. In Section 2, we formulate the Dirichlet boundary control problem with pointwise control constraints. Therein, we discuss the well-posedness of the model problem and present the optimality system. In Section 3, we define the discrete control problem and present the discrete optimality system. In Section 4, we derive a posteriori error estimates with the help of some reconstruction solution. Section 5 is devoted to the numerical experiments.

2. Continuous Problem

We proceed over the precise formulation of the optimization problem in brief in this section. We need the following definitions and notations before we can begin the analysis:

2.1. Notation

Let $\Omega\subset\mathbb{R}^{2},$ be a bounded polygonal domain, with boundary $\partial\Omega$ consists of three non-overlapping open subsets $\Gamma_{D},\;\Gamma_{C}$ and $\Gamma_{N}$ with $\partial\Omega=\Gamma_{C}\cup\bar{\Gamma}_{D}\cup\Gamma_{N}$ . The one-dimensional measure of $\Gamma_{C}$ is positive. We denote any function and any space in bold notation can be understood in the vector form e.g., $\mathbf{x}:=(x_{1},x_{2}),$ $\mathbf{L}^{2}(\Omega):=[L^{2}(\Omega)]^{2}$ and $\mathbf{H}^{1}(\Omega):=[H^{1}(\Omega)]^{2}.$ The norm and inner product on those spaces are defined component wise. Here and throughout, the $\mathbf{L}^{2}(\Omega)$ norm is denoted by $\|\cdot\|_{0,\Omega}$ and $\|\cdot\|_{k,\Omega}(k>0)$ denotes the standard norm on the Sobolev space $\mathbf{H}^{k}(\Omega)$ , see for example [11]. The trace of a vector valued function $\mathbf{x}\in\mathbf{H}^{1}(\Omega)$ is defined to be $\bm{\gamma}_{0}(\mathbf{x}):=(\gamma_{0}(x_{1}),\gamma_{0}(x_{2})),$ where $\gamma_{0}:H^{1}(\Omega)\rightarrow L^{2}(\Gamma)$ is the trace operator. Let $\mathbf{x}$ and $\mathbf{y}$ are two functions, we say that $\mathbf{x}\leq\mathbf{y}$ iff $x_{1}\leq y_{1}$ and $x_{2}\leq y_{2}$ almost everywhere in $\Omega.$

2.2. Dirichlet Control Problem

We consider the following constrained Dirichlet boundary control problem(in energy form [10, Section 2]) governed by the Stokes equation

\text{min}~{}J(\mathbf{u},\mathbf{y})=\frac{1}{2}\left\|\mathbf{u}-\mathbf{u}_{d}\right\|_{0,\Omega}^{2}+\frac{\rho}{2}\left\|\nabla\mathbf{y}\right\|_{0,\Omega}^{2}

(2.1)

subject to,

\begin{split}-\Delta\mathbf{u}+\nabla{p}&=\mathbf{f}\quad\text{in}\;\Omega,\\ \nabla\cdot{\mathbf{u}}&=0\quad\text{in}\;\Omega,\\ \bf{u}&=\mathbf{y}\quad\text{on}\;\Gamma_{C},\\ \bf{u}&=\mathbf{0}\quad\text{on}\;\Gamma_{D},\\ \frac{\partial\mathbf{u}}{\partial\mathbf{n}}-p\mathbf{n}&=\mathbf{0}\;\;\text{on}\;\;\Gamma_{N},\end{split}

(2.2)

with the control $\mathbf{y}$ comes from the following constrained set

\displaystyle\mathbf{Q}_{ad}:=\{\mathbf{x}\in\mathbf{Q}:\mathbf{y}_{a}\leq\bm{\gamma}_{0}(\mathbf{x})\leq\mathbf{y}_{b}\text{ a.e. on }\Gamma_{C}\}.

The interior force $\mathbf{f}\in\mathbf{L}^{2}(\Omega),$ the regularization parameter $\rho>0,$ and $\mathbf{u}_{d}\in\mathbf{L}^{2}(\Omega)$ and the space $\mathbf{Q}:=\{\mathbf{x}\in\mathbf{H}^{1}(\Omega):\bm{\gamma}_{0}(\mathbf{x})=\mathbf{0}\text{ on }\Gamma_{D}\cup\Gamma_{N}\}.$ The constant vectors $\mathbf{y}_{a}=(y^{1}_{a},y^{2}_{a})$ , and $\mathbf{y}_{b}=(y^{1}_{b},y^{2}_{b})\in\mathbb{R}^{2}$ satisfying $y^{1}_{a}<y^{2}_{a}$ and $y^{1}_{b}<y^{2}_{b}$ . Furthermore whenever $\Gamma_{D}$ is non empty, we assume for consistency that $y^{1}_{a},y^{1}_{b}\leq 0$ and $y^{2}_{a},y^{2}_{b}\geq 0$ so that, the admissible set $\mathbf{Q}_{ad}$ is nonempty.

A proof of the existence of the unique solution of the control problem (2.1) can be found in [20, Theorem 2.2]. The following proposition states the first-order optimality system, a details can be found in [20, Proposition 2.3].

Proposition 2.1.

There exists a unique solution $(\mathbf{u},p,\mathbf{y})\in\mathbf{H}_{D}^{1}(\Omega)\times L^{2}(\Omega)\times\mathbf{Q}_{ad}$ for the Dirichlet control problem $\eqref{min:j}$ and there exists an adjoint state $(\bm{\phi},r)\in\mathbf{V}\times L^{2}(\Omega)$ satisfying


$\displaystyle\mathbf{u}$	$\displaystyle=\mathbf{w}+\mathbf{y},\quad\mathbf{w}\in\mathbf{V},$	(2.3a)
$\displaystyle a(\mathbf{w},\mathbf{z})+b(\mathbf{z},p)$	$\displaystyle=(\mathbf{f},\mathbf{z})-a(\mathbf{y},\mathbf{z})\;\;\;{\rm for~{}all}\;\mathbf{z}\in\mathbf{V},$	(2.3b)
$\displaystyle b(\mathbf{u},q)$	$\displaystyle=0\;\quad{\rm for~{}all}\;q\in L^{2}(\Omega),$	(2.3c)
$\displaystyle a(\mathbf{z},\bm{\phi})-b(\mathbf{z},r)$	$\displaystyle=(\mathbf{u-u_{d}},\mathbf{z})\;\;\;{\rm for~{}all}\;\mathbf{z}\in\mathbf{V},$	(2.3d)
$\displaystyle b(\bm{\phi},q)$	$\displaystyle=0\;\quad{\rm for~{}all}\;q\in L^{2}(\Omega),$	(2.3e)
$\displaystyle\rho\,a(\mathbf{y},\mathbf{x}-\mathbf{y})$	$\displaystyle\geq a(\mathbf{x}-\mathbf{y},\bm{\phi})-b(\mathbf{x}-\mathbf{y},r)-(\mathbf{u-u_{d}},\mathbf{x}-\mathbf{y})\quad\forall\mathbf{x}\in\mathbf{Q}_{ad},$	(2.3f)

where $a(\mathbf{w},\mathbf{z})=\int_{\Omega}\nabla{\mathbf{w}}:\nabla{\mathbf{z}}{\rm~{}dx}$ , $b(\mathbf{z},p)=-\int_{\Omega}p\nabla\cdot{\mathbf{z}}{\rm~{}dx},$ and the matrix product $A:B:=\sum_{i,j=1}^{n}a_{ij}b_{ij}$ when $A=(a_{ij})_{1\leq i,j\leq n}$ and $B=(b_{ij})_{1\leq i,j\leq n}.$

Remark 2.2.

It is not hard to show from the equation (2.3f), that the optimal control $\mathbf{y}$ is the weak solution of the following simplified Signorini problem:


	$\displaystyle-\rho\Delta\mathbf{y}$	$\displaystyle=\mathbf{0}\quad\text{in}\quad\Omega,$
	$\displaystyle\mathbf{y}$	$\displaystyle=\mathbf{0}\quad\text{on}\quad\Gamma_{D}\cup\Gamma_{N},$
	$\displaystyle\mathbf{y}_{a}\leq\bm{\gamma}_{0}(\mathbf{y})$	$\displaystyle\leq\mathbf{y}_{b}\;\text{ a.e. on }\Gamma_{C}$
further, the following holds for almost every $x\in\Gamma_{C}$ :

	$\displaystyle\text{if}\;\mathbf{y}_{a}<\mathbf{y}(x)<\mathbf{y}_{b}\quad\text{then}\quad\big{(}\hat{\sigma}(\mathbf{y})\big{)}(x)$	$\displaystyle=\mathbf{0},$
	$\displaystyle\text{if}\;\mathbf{y}_{a}\leq\mathbf{y}(x)<\mathbf{y}_{b}\quad\text{then}\quad\big{(}\hat{\sigma}(\mathbf{y})\big{)}(x)$	$\displaystyle\geq\mathbf{0},$
	$\displaystyle\text{if}\;\mathbf{y}_{a}<\mathbf{y}(x)\leq\mathbf{y}_{b}\quad\text{then}\quad\big{(}\hat{\sigma}(\mathbf{y})\big{)}(x)$	$\displaystyle\leq\mathbf{0},$

where $\hat{\sigma}(\mathbf{y}):=\rho\frac{\partial\mathbf{y}}{\partial\mathbf{n}}-\frac{\partial\bm{\phi}}{\partial\mathbf{n}}-r\mathbf{n}.$

3. Discrete Problem

In this section, we discuss the discrete control problem before this we need to define some notations. Let $\mathcal{T}_{H}$ be a shape-regular triangulation of the domain $\Omega$ into triangles such that $\cup_{T\in\mathcal{T}_{H}}T=\bar{\Omega}$ see [2, 11]. Also let $\mathcal{T}_{h}$ be a refinement of $\mathcal{T}_{H}$ by connecting all the midpoints of $\mathcal{T}_{H}.$ Denote the set of all interior edges of $\mathcal{T}_{h}$ by $\mathcal{E}_{h}^{i}.$ The set of all Dirichlet, Neumann and Contact boundary edges of $\mathcal{T}_{h}$ are denoted by $\mathcal{E}_{h}^{b,D},$ $\mathcal{E}_{h}^{b,N}$ and $\mathcal{E}_{h}^{b,C}$ respectively and define $\mathcal{E}_{h}=\mathcal{E}_{h}^{i}\cup\mathcal{E}_{h}^{b,D}\cup\mathcal{E}_{h}^{b,N}\cup\mathcal{E}_{h}^{b,C}.$ A typical triangle is denoted by $T$ and its diameter by $h_{T}$ . Set $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ . The length of any edge $e\in\mathcal{E}_{h}$ will be denoted by $h_{e}$ . Let $\mathcal{V}_{h}$ denote the set of all the vertices of the triangles in $\mathcal{T}_{h}$ . The set of vertices on $\overline{\Gamma}_{D},$ $\Gamma_{N}$ and $\Gamma_{C}$ are denoted by $\mathcal{V}_{h}^{D},$ $\mathcal{V}_{h}^{N}$ and $\mathcal{V}_{h}^{C}.$ Also, in the problem setting, we require the jump definitions of scalar, vector valued functions and tensors. Let us define a broken Sobolev space

H^{1}(\Omega,\mathcal{T}_{h})=\{v\in L^{2}(\Omega):v|_{T}\in H^{1}(T)\;\;{\rm for~{}all~{}}T\in\mathcal{T}_{h}\}.

For any $e\in\mathcal{E}^{i}_{h}$ , there are two triangles $T_{+}$ and $T_{-}$ such that $e=\partial T_{+}\cap\partial T_{-}$ . Let $\mathbf{n}_{+}$ be the unit normal of $e$ pointing from $T_{+}$ to $T_{-}$ and let $\mathbf{n}_{-}=-\mathbf{n}_{+}$ (cf. Fig.3.1). For any $v\in H^{1}(\Omega,\mathcal{T}_{h})$ , we define the jump of $v$ on an edge $e$ by $[\hskip-1.5pt[v]\hskip-1.5pt]=v_{+}\mathbf{n}_{+}+v_{-}\mathbf{n}_{-}$ where $v_{\pm}=v|_{T_{\pm}}.$

Figure 3.1. Here

T_{-}

and

T_{+}

are the two neighboring triangles that share the edge

e=\partial T_{-}\cap\partial T_{+}

with initial node

A

and end node

B

and unit normal

\mathbf{n}_{e}

. The orientation of

\mathbf{n}_{e}=\mathbf{n}_{-}=-\mathbf{n}_{+}

equals the outer normal of

T_{-}

, and hence, points into

T_{+}

For $\mathbf{v}\in[H^{1}(\Omega,\mathcal{T}_{h})]^{2}$ we define the jump of $\mathbf{v}$ on $e\in\mathcal{E}^{i}_{h}.$ by $[\hskip-1.5pt[\mathbf{v}]\hskip-1.5pt]=\mathbf{v}_{+}\cdot\mathbf{n}_{+}+\mathbf{v}_{-}\cdot\mathbf{n}_{-}.$ Similarly, for tensors $\tau\in[H^{1}(\Omega,\mathcal{T}_{h})]^{2\times 2}$ , the jump on $e\in\mathcal{E}^{i}_{h}$ are defined by $[\hskip-1.5pt[\tau]\hskip-1.5pt]=\tau_{+}\mathbf{n}_{+}+\tau_{-}\mathbf{n}_{-}.$ For notational convenience, we also define the jump on the boundary faces $e\in\mathcal{E}_{h}^{b}$ by modifying them appropriately. We use the definition of jump by understanding that $v_{-}=0$ (similarly, $\mathbf{v}_{-}=0$ and $\tau_{-}=0$ ).

Define the discrete space for velocity $\mathbf{V}_{h}\subset\mathbf{V}$ by

\mathbf{V}_{h}:=\{\mathbf{v}_{h}\in\mathbf{V}:\mathbf{v}_{h}|_{T}\in\mathbf{P}_{1}(T)\;\;\forall T\in\mathcal{T}_{h}\},

and the discrete space for pressure is

M_{H}:=\{p_{H}\in L^{2}(\Omega):p_{H}|_{T}\in P_{0}(T)\;\;\forall T\in\mathcal{T}_{H}\},

and the discrete control space $\mathbf{Q}_{h}\subset\mathbf{Q}$ by

\mathbf{Q}_{h}:=\{\mathbf{x}_{h}\in\mathbf{Q}:\mathbf{x}_{h}|_{T}\in\mathbf{P}_{1}(T),\;\;\forall T\in\mathcal{T}_{h}\},

where $\mathbf{P}_{1}(T)$ is the space of polynomials of degree less than or equal to one on the triangle $T$ . The discrete admissible set of controls is defined by

\mathbf{Q}^{h}_{ad}:=\{\mathbf{x}_{h}\in\mathbf{Q}_{h}:\mathbf{y}_{a}\leq\mathbf{x}_{h}(z)\leq\mathbf{y}_{b}\text{ for all }z\in\mathcal{V}_{h}^{C}\}.

It is easy to check that $\mathbf{Q}_{ad}^{h}\subset\mathbf{Q}_{ad}.$ Throughout the article, we assume that $C$ denotes a generic positive constant that is independent of the mesh parameter $h$ . A proof of the following proposition on the existence and uniqueness of the solution of discrete problem can be found in [20, Proposition 3.1].

Proposition 3.1 (Discrete Optimality System).

There exists unique $((\mathbf{w}_{h},p_{H}),(\bm{\phi}_{h},r_{H}),\mathbf{y}_{h})\in\big{(}\mathbf{V}_{h}\times M_{H}\big{)}\times\big{(}\mathbf{V}_{h}\times M_{H}\big{)}\times\mathbf{Q}^{h}_{ad}$ satisfying the following:


$\displaystyle\mathbf{u}_{h}$	$\displaystyle=\mathbf{w}_{h}+\mathbf{y}_{h},\quad\mathbf{w}_{h}\in\mathbf{V}_{h},$
$\displaystyle a(\mathbf{w}_{h},\mathbf{z}_{h})+b(\mathbf{z}_{h},p_{H})$	$\displaystyle=(\mathbf{f},\mathbf{z}_{h})-a(\mathbf{y}_{h},\mathbf{z}_{h})\;\;\;{\rm for~{}all}\;\mathbf{z}_{h}\in\mathbf{V}_{h},$	(3.1a)
$\displaystyle b(\mathbf{u}_{h},q_{H})$	$\displaystyle=0\;\quad{\rm for~{}all}\;q_{H}\in M_{H},$	(3.1b)
$\displaystyle a(\mathbf{z}_{h},\bm{\phi}_{h})-b(\mathbf{z}_{h},r_{H})$	$\displaystyle=(\mathbf{u}_{h}-\mathbf{u_{d}},\mathbf{z}_{h})\;\;\;{\rm for~{}all}\;\mathbf{z}_{h}\in\mathbf{V}_{h},$	(3.1c)
$\displaystyle b(\bm{\phi}_{h},q_{H})$	$\displaystyle=0\;\quad{\rm for~{}all}\;q_{H}\in M_{H},$	(3.1d)
$\displaystyle\rho\,a(\mathbf{y}_{h},\mathbf{x}_{h}-\mathbf{y}_{h})\geq$	$\displaystyle a(\mathbf{x}_{h}-\mathbf{y}_{h},\bm{\phi}_{h})-b(\mathbf{x}_{h}-\mathbf{y}_{h},r_{H})$
	$\displaystyle-(\mathbf{u}_{h}-\mathbf{u_{d}},\mathbf{x}_{h}-\mathbf{y}_{h})\;\quad{\rm for~{}all}\;\mathbf{x}_{h}\in\mathbf{Q}^{h}_{ad}.$	(3.1e)

4. A posteriori Error Analysis

This section is devoted to a posteriori error analysis. Define reconstructions $\mathbf{R}\mathbf{w}\in\mathbf{V}$ , $R_{0}p\in L^{2}(\Omega),$ $\bar{\mathbf{R}}\bm{\phi}\in\mathbf{V}$ and $\bar{R}_{0}r\in L^{2}(\Omega),$ $\mathbf{R}\mathbf{y}\in\mathbf{Q}_{ad}$ by


$\displaystyle\mathbf{R}\mathbf{u}=$	$\displaystyle\mathbf{R}\mathbf{w}+\mathbf{y}_{h}$
$\displaystyle a(\mathbf{R}\mathbf{w},\mathbf{v})+b(\mathbf{v},R_{0}p)=$	$\displaystyle{\langle\mathbf{f},\mathbf{z}\rangle}-a(\mathbf{y}_{h},\mathbf{v})\;~{}~{}\hskip 31.2982pt{\rm for~{}all}\;\mathbf{v}\in\mathbf{V},$	(4.1a)
$\displaystyle b(\mathbf{R}\mathbf{w},q)=$	$\displaystyle-b(\mathbf{y}_{h},q)\;~{}~{}\hskip 59.75095pt{\rm for~{}all}\;q\in L^{2}(\Omega),$	(4.1b)
$\displaystyle a(\mathbf{z},\bar{\mathbf{R}}\bm{\phi})-b(\mathbf{z},\bar{R_{0}}r)=$	$\displaystyle{\langle\mathbf{u}_{h}-\mathbf{u}_{d},\mathbf{z}\rangle}_{W}\;~{}~{}\hskip 45.52458pt{\rm for~{}all}\;\mathbf{z}\in\mathbf{V},$	(4.1c)
$\displaystyle b(\bar{\mathbf{R}}\bm{\phi},q)=$	$\displaystyle 0\;~{}~{}\hskip 110.96556pt{\rm for~{}all}\;q\in L^{2}(\Omega).$	(4.1d)
$\displaystyle\rho a(\mathbf{R}\mathbf{y},\mathbf{x}-\mathbf{R}\mathbf{y})\geq$	$\displaystyle a(\mathbf{x}-\mathbf{R}\mathbf{y},\bm{\phi}_{h})-b(\mathbf{x}-\mathbf{R}\mathbf{y},r_{H})$	(4.1e)
	$\displaystyle-(\mathbf{u}_{h}-\mathbf{u}_{d},\mathbf{x}-\mathbf{R}\mathbf{y})\;~{}~{}\hskip 14.22636pt{\rm for~{}all}\;\mathbf{x}\in\mathbf{Q}_{ad}.$

The well-posedness of the above system (4.1) follows from the facts that the right-hand side of (4.1a) is a bounded linear functional on $\mathbf{V}$ , the bilinear forms $a$ and $b$ are continuous, $a$ is elliptic and $b$ is inf-sup stable, and hence the system (4.1a)-(4.1b) has a unique solution [16, pp. 81] . Similarly, the system (4.1c)-(4.1d) is well-posed.

Subtraction of (4.1) from (2.3) yields,


$\displaystyle a(\mathbf{w}-\mathbf{R}\mathbf{w},\mathbf{z})+b(\mathbf{z},p-R_{0}p)$	$\displaystyle=-a(\mathbf{y}-\mathbf{y}_{h},\mathbf{z})\;\;~{}~{}{\rm for~{}all}\;\mathbf{z}\in\mathbf{V},$	(4.2a)
$\displaystyle b(\mathbf{u}-\mathbf{R}\mathbf{u},q)$	$\displaystyle=0\;~{}~{}\hskip 56.9055pt{\rm for~{}all}\;q\in L^{2}(\Omega),$	(4.2b)
$\displaystyle a(\mathbf{z},\bm{\phi}-\bar{\mathbf{R}}\bm{\phi})+b(\mathbf{z},\bar{R_{0}}r-r)$	$\displaystyle=(\mathbf{u}-\mathbf{u}_{h},\mathbf{z})\;~{}~{}\hskip 8.5359pt{\rm for~{}all}\;\mathbf{z}\in\mathbf{V},$	(4.2c)
$\displaystyle b(\bm{\phi}-\bar{\mathbf{R}}\bm{\phi},q)$	$\displaystyle=0\;~{}~{}\hskip 68.28644pt{\rm for~{}all}\;q\in L^{2}(\Omega).$	(4.2d)

Theorem 4.1 (Energy error estimate of control and $L^{2}$ -estimate of velocity).

There holds,

	$\displaystyle\rho\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y})\right\\|_{0,\Omega}+\left\\|\mathbf{R}\mathbf{u}-\mathbf{u}\right\\|_{0,\Omega}\leq$	$\displaystyle C\big{(}\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\\|_{0,\Omega}$
		$\displaystyle+\left\\|\nabla(\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})\right\\|_{0,\Omega}+\left\\|r_{H}-\bar{R}_{0}r\right\\|_{0,\Omega}\big{)}.$

Proof.

Selecting $\mathbf{x=y}$ in (4.1e), $\mathbf{x=Ry}$ in (2.3f) and adding the result, we obtain

	$\displaystyle\rho a(\mathbf{R}\mathbf{y}-\mathbf{y},\mathbf{y}-\mathbf{R}\mathbf{y})\geq$	$\displaystyle a(\mathbf{y}-\mathbf{R}\mathbf{y},\bm{\phi}_{h}-\bm{\phi})-b(\mathbf{y}-\mathbf{R}\mathbf{y},r_{H}-r)-(\mathbf{u}_{h}-\mathbf{u},\mathbf{y}-\mathbf{R}\mathbf{y})$
	$\displaystyle\geq$	$\displaystyle a(\mathbf{y}-\mathbf{R}\mathbf{y},\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})+a(\mathbf{y}-\mathbf{y}_{h},\bar{\mathbf{R}}\bm{\phi}-\bm{\phi})+a(\mathbf{y}_{h}-\mathbf{R}\mathbf{y},\bar{\mathbf{R}}\bm{\phi}-\bm{\phi})$
		$\displaystyle-b(\mathbf{y}-\mathbf{R}\mathbf{y},r_{H}-r)-(\mathbf{u}_{h}-\mathbf{u},\mathbf{y}-\mathbf{R}\mathbf{y})$
	$\displaystyle\geq$	$\displaystyle a(\mathbf{y}-\mathbf{R}\mathbf{y},\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})+a(\mathbf{y}-\mathbf{y}_{h},\bar{\mathbf{R}}\bm{\phi}-\bm{\phi})+a(\mathbf{y}_{h}-\mathbf{R}\mathbf{y},\bar{\mathbf{R}}\bm{\phi}-\bm{\phi})$
		$\displaystyle-b(\mathbf{y}-\mathbf{R}\mathbf{y},r_{H}-r)-(\mathbf{u}_{h}-\mathbf{u},\mathbf{y}-\mathbf{R}\mathbf{y})$
	$\displaystyle\geq$	$\displaystyle a(\mathbf{y}-\mathbf{R}\mathbf{y},\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})-b(\mathbf{y}-\mathbf{y}_{h},r-\bar{R}_{0}r)+a(\mathbf{y}_{h}-\mathbf{R}\mathbf{y},\bar{\mathbf{R}}\bm{\phi}-\bm{\phi})$
		$\displaystyle-b(\mathbf{y}-\mathbf{R}\mathbf{y},r_{H}-r)-(\mathbf{u}-\mathbf{u}_{h},\mathbf{y}+\mathbf{w}-\mathbf{R}\mathbf{y}-\mathbf{R}\mathbf{w})$
	$\displaystyle\geq$	$\displaystyle a(\mathbf{y}-\mathbf{R}\mathbf{y},\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})-b(\mathbf{y}-\mathbf{y}_{h},r-\bar{R}_{0}r)+a(\mathbf{y}_{h}-\mathbf{R}\mathbf{y},\bar{\mathbf{R}}\bm{\phi}-\bm{\phi})$
		$\displaystyle-b(\mathbf{y}-\mathbf{R}\mathbf{y},r_{H}-r)+\left\\|\mathbf{u}-\mathbf{R}\mathbf{u}\right\\|_{0,\Omega}^{2}+(\mathbf{R}\mathbf{u}-\mathbf{u}_{h},\mathbf{u}-\mathbf{R}\mathbf{u})$
		$\displaystyle+(\mathbf{u}-\mathbf{u}_{h},\mathbf{y}_{h}-\mathbf{R}\mathbf{y}).$

Now we have

	$\displaystyle\rho\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y})\right\\|_{0,\Omega}^{2}+\left\\|\mathbf{u}-\mathbf{R}\mathbf{u}\right\\|_{0,\Omega}^{2}\leq$	$\displaystyle-a(\mathbf{y}-\mathbf{R}\mathbf{y},\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})+b(\mathbf{y}-\mathbf{y}_{h},r-\bar{R}_{0}r)$
		$\displaystyle-a(\mathbf{y}_{h}-\mathbf{R}\mathbf{y},\bar{\mathbf{R}}\bm{\phi}-\bm{\phi})+b(\mathbf{y}-\mathbf{R}\mathbf{y},r_{H}-r)$
		$\displaystyle-(\mathbf{R}\mathbf{u}-\mathbf{u}_{h},\mathbf{u}-\mathbf{R}\mathbf{u})-(\mathbf{u}-\mathbf{u}_{h},\mathbf{y}_{h}-\mathbf{R}\mathbf{y})$
	$\displaystyle\leq$	$\displaystyle-a(\mathbf{y}-\mathbf{R}\mathbf{y},\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})+b(\mathbf{R}\mathbf{y}-\mathbf{y}_{h},r-\bar{R}_{0}r)$
		$\displaystyle+b(\mathbf{y}-\mathbf{R}\mathbf{y},r_{H}-\bar{R}_{0}r)-a(\mathbf{y}_{h}-\mathbf{R}\mathbf{y},\bar{\mathbf{R}}\bm{\phi}-\bm{\phi})$
		$\displaystyle-(\mathbf{R}\mathbf{u}-\mathbf{u}_{h},\mathbf{u}-\mathbf{R}\mathbf{u})-(\mathbf{R}\mathbf{u}-\mathbf{u}_{h},\mathbf{y}_{h}-\mathbf{R}\mathbf{y})$
		$\displaystyle+(\mathbf{u}-\mathbf{R}\mathbf{u},\mathbf{y}_{h}-\mathbf{R}\mathbf{y}).$

Applying Cauchy-Schwarz inequality and Young’s inequality we arrive at the desired estimate. ∎

Theorem 4.2 (Energy error estimate of velocity).

There holds,

	$\displaystyle\left\\|\nabla(\mathbf{u}-\mathbf{R}\mathbf{u})\right\\|_{0,\Omega}\leq$	$\displaystyle C\big{(}\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})\right\\|_{0,\Omega}$
		$\displaystyle+\left\\|r_{H}-\bar{R}_{0}r\right\\|_{0,\Omega}\big{)}.$

Proof.

The splitting $\mathbf{u=w+y}$ and $\mathbf{Ru=Rw+y_{h}}$ yields

\displaystyle\left\|\nabla(\mathbf{u}-\mathbf{R}\mathbf{u})\right\|_{0,\Omega}\leq\left\|\nabla(\mathbf{w}-\mathbf{R}\mathbf{w})\right\|_{0,\Omega}+\left\|\nabla(\mathbf{y}-\mathbf{y}_{h})\right\|_{0,\Omega}.

(4.3)

Now we need to estimate the term $\left\|\nabla(\mathbf{w}-\mathbf{R}\mathbf{w})\right\|_{0,\Omega}.$ A selection $\mathbf{z}=\mathbf{w}-\mathbf{R}\mathbf{w}$ in (4.2a) provides

\displaystyle\left\|\nabla(\mathbf{w}-\mathbf{R}\mathbf{w})\right\|_{0,\Omega}^{2}+b(\mathbf{w}-\mathbf{R}\mathbf{w},p-R_{0}p)=-a(\mathbf{y}-\mathbf{y}_{h},\mathbf{w}-\mathbf{R}\mathbf{w}).

(4.4)

Using the fact that $b(\mathbf{w}-\mathbf{R}\mathbf{w},p-R_{0}p)=-b(\mathbf{y}-\mathbf{y}_{h},p-R_{0}p)$ and appplying Cauchy-Schwarz and Young’s inequality in (4.4), we obtain

\displaystyle(1-\frac{\epsilon}{2})\left\|\nabla(\mathbf{w}-\mathbf{R}\mathbf{w})\right\|_{0,\Omega}^{2}\leq\frac{1}{\epsilon}\left\|\nabla(\mathbf{y}-\mathbf{y}_{h})\right\|_{0,\Omega}^{2}+\frac{\epsilon}{2}\left\|p-R_{0}p\right\|_{0,\Omega}^{2}.

(4.5)

The estimate of $\left\|p-R_{0}p\right\|_{0,\Omega}$ is in the following:

$\displaystyle\left\\|p-R_{0}p\right\\|_{0,\Omega}$	$\displaystyle\leq\sup_{\mathbf{v}\in\mathbf{V}}\frac{b(\mathbf{v},p-R_{0}p)}{\left\\|v\right\\|_{\mathbf{V}}}$
	$\displaystyle\leq\sup_{\mathbf{v}\in\mathbf{V}}\frac{-a(\mathbf{y}-\mathbf{y_{h}},\mathbf{v})-a(\mathbf{w}-\mathbf{R}\mathbf{w},\mathbf{v})}{\left\\|v\right\\|_{\mathbf{V}}}$
	$\displaystyle\leq\left\\|\nabla(\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{w}-\mathbf{R}\mathbf{w})\right\\|_{0,\Omega}.$	(4.6)

Using (4) in (4.5), we obtain $\left\|\nabla(\mathbf{w}-\mathbf{R}\mathbf{w})\right\|_{0,\Omega}\leq C\left\|\nabla(\mathbf{y}-\mathbf{y}_{h})\right\|_{0,\Omega}.$ Hence from (4.3) we get

\displaystyle\left\|\nabla(\mathbf{u}-\mathbf{R}\mathbf{u})\right\|_{0,\Omega}

\displaystyle\leq C\left\|\nabla(\mathbf{y}-\mathbf{y}_{h})\right\|_{0,\Omega}

(4.7)

Introducing the reconstruction $\mathbf{R}\mathbf{y}$ in (4.7) we obtain

\displaystyle\left\|\nabla(\mathbf{u}-\mathbf{R}\mathbf{u})\right\|_{0,\Omega}

\displaystyle\leq C\big{(}\left\|\nabla(\mathbf{y}-\mathbf{R}\mathbf{y})\right\|_{0,\Omega}+\left\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\|_{0,\Omega}\big{)}.

(4.8)

Substituting the estimate of $\left\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y})\right\|_{0,\Omega}$ (from Theorem 4.1) in (4.8), we obtain the desired estimate. ∎

Theorem 4.3 (Energy error estimate of adjoint velocity).

There holds,

	$\displaystyle\left\\|\nabla(\bm{\phi}-\mathbf{R}\bm{\phi})\right\\|_{0,\Omega}\leq$	$\displaystyle C\big{(}\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})\right\\|_{0,\Omega}$
		$\displaystyle+\left\\|r_{H}-\bar{R}_{0}r\right\\|_{0,\Omega}\big{)}.$

Proof.

Selecting $\mathbf{z}=\bm{\phi}-\mathbf{R}\bm{\phi}$ in (4.2c) we have the following:

\displaystyle\left\|\nabla(\bm{\phi}-\bar{\mathbf{R}}\bm{\phi})\right\|_{0,\Omega}^{2}-b(\bm{\phi}-\bar{\mathbf{R}}\bm{\phi},r-\bar{R}_{0}r)=(\mathbf{u}-\mathbf{u}_{h},\bm{\phi}-\bar{\mathbf{R}}\bm{\phi}).

(4.9)

Using the fact that $b(\bm{\phi}-\bar{\mathbf{R}}\bm{\phi},r-\bar{R}_{0}r)=0$ and applying Cauchy-Schwarz inequality in (4.9), we obtain

$\displaystyle\left\\|\nabla(\bm{\phi}-\bar{\mathbf{R}}\bm{\phi})\right\\|_{0,\Omega}$	$\displaystyle\leq\left\\|\mathbf{u}-\mathbf{u}_{h}\right\\|_{0,\Omega}$
	$\displaystyle\leq\left\\|\mathbf{u}-\mathbf{Ru}\right\\|_{0,\Omega}+\left\\|\mathbf{Ru}-\mathbf{u}_{h}\right\\|_{0,\Omega}$
	$\displaystyle\leq\left\\|\mathbf{u}-\mathbf{Ru}\right\\|_{0,\Omega}+\left\\|\mathbf{Rw}-\mathbf{w}_{h}\right\\|_{0,\Omega}$
	$\displaystyle\leq\left\\|\mathbf{u}-\mathbf{Ru}\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{Rw}-\mathbf{w}_{h})\right\\|_{0,\Omega}.$	(4.10)

In the above, we have used the fact that $\mathbf{u}_{h}=\mathbf{w}_{h}+\mathbf{y}_{h}$ , $\mathbf{Ru}=\mathbf{Rw}+\mathbf{y}_{h}$ and Poincaré inequality. Using the estimates of $\left\|\mathbf{u}-\mathbf{Ru}\right\|_{0,\Omega}$ from Theorem 4.1 in (4.10) we obtain the desired result. ∎

In the next theorem we state the estimate of pressure and adjoint pressure the proof can be easily derived by using the inf-sup condition.

Theorem 4.4 (Error estimate of pressure and adjoint pressure).

There holds,

	$\displaystyle\left\\|p-R_{0}p\right\\|_{0,\Omega}\leq$	$\displaystyle C\big{(}\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})\right\\|_{0,\Omega}$
		$\displaystyle+\left\\|r_{H}-\bar{R}_{0}r\right\\|_{0,\Omega}\big{)},$

and

	$\displaystyle\left\\|r-\bar{R}_{0}r\right\\|_{0,\Omega}\leq$	$\displaystyle C\big{(}\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})\right\\|_{0,\Omega}$
		$\displaystyle+\left\\|r_{H}-\bar{R}_{0}r\right\\|_{0,\Omega}\big{)}.$

Combining all the above results, Theorem 4.1-4.4, we obtain the following theorem:

Theorem 4.5.

There holds

	$\displaystyle\left\\|\nabla(\mathbf{u}-\mathbf{u}_{h})\right\\|_{0,\Omega}+$	$\displaystyle\left\\|p-p_{H}\right\\|_{0,\Omega}+\left\\|\nabla(\bm{\phi}-\bm{\phi}_{h})\right\\|_{0,\Omega}+\left\\|r-r_{H}\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}\leq C\big{(}\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}$
		$\displaystyle+\left\\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\\|_{0,\Omega}+\left\\|p_{H}-R_{0}p\right\\|_{0,\Omega}+\left\\|\nabla(\bm{\phi}_{h}-\bar{\mathbf{R}}\bm{\phi})\right\\|_{0,\Omega}+\left\\|r_{H}-\bar{R}_{0}r\right\\|_{0,\Omega}\big{)}.$

Before going to derive the a posteriori error estimator we need to define some preliminary definitions which are given in the following: Let $p\in\mathcal{V}_{h}$ be a node and $\omega_{p}$ be the node patch of $p$ and define $h_{p}=\text{diam}\;\omega_{p}.$ Denote $\gamma_{p}$ be the union of all sides of $\bar{\omega}_{p}$ and union of interior sides of $\bar{\omega}_{p}$ are denoted by $\gamma_{p,I}.$ Given any $p\in\mathcal{V}_{h}^{C}$ we subdivide the intersections between $\partial\Omega$ and $\partial\omega_{p}$ in the three following sets:

	$\displaystyle\gamma_{p,C}$	$\displaystyle:=\Gamma_{C}\cap\partial\omega_{p}$
	$\displaystyle\gamma_{p,D}$	$\displaystyle:=\Gamma_{D}\cap\partial\omega_{p}$
	$\displaystyle\gamma_{p,N}$	$\displaystyle:=\Gamma_{N}\cap\partial\omega_{p}.$

We define the discrete contact stress by $\hat{\bm{\sigma}}(\mathbf{y}_{h}):=\rho\frac{\partial\mathbf{y}_{h}}{\partial n}-\frac{\partial\bm{\phi}_{h}}{\partial n}-r_{H}\mathbf{n},$ clearly it is a vector quantity so the sign of this quantity would be in the sense of componentwise. We classify the nodes on $\Gamma_{C}$ as follows.

We classify the actual contact nodes $p\in\mathcal{V}_{h}^{C}$ with $\mathbf{y}_{h}(p)=\mathbf{y}_{a}(p)(=\mathbf{y}_{a})$ in two different categories. At so-called full-contact nodes $p\in\mathcal{V}^{fC}_{h,a}$ if the discrete solution satisfies $\mathbf{y}_{h}=\mathbf{y}_{a},$ $\hat{\bm{\sigma}}(\mathbf{y}_{h}):=(\hat{\sigma}_{1}(\mathbf{y}_{h}),\hat{\sigma}_{2}(\mathbf{y}_{h}))\geq\mathbf{0}$ on $\gamma_{p,C},$ which means that the conditions of actual contact are satisfied. The remaining actual contact nodes with $\mathbf{y}_{h}(p)=\mathbf{y}_{a}(p)$ are called semi-contact nodes and the set is denoted by $\mathcal{V}^{sC}_{h,a}.$

Similarly, we classify the actual contact nodes $p\in\mathcal{V}_{h}^{C}$ with $\mathbf{y}_{h}(p)=\mathbf{y}_{b}(p)(=\mathbf{y}_{b})$ in two different categories. At so-called full-contact nodes $p\in\mathcal{V}^{fC}_{h,b}$ if the discrete solution satisfies $\mathbf{y}_{h}=\mathbf{y}_{b},$ $\hat{\bm{\sigma}}(\mathbf{y}_{h})\leq\mathbf{0}$ on $\gamma_{p,C},$ which means that the conditions of actual contact are satisfied. The remaining actual contact nodes with $\mathbf{y}_{h}(p)=\mathbf{y}_{b}(p)$ are called semi-contact nodes and the set is denoted by $\mathcal{V}^{sC}_{h,b}.$

Also, we define $\mathcal{V}_{*}^{C}=\{p\in\mathcal{V}^{C}_{h}:\mathbf{y}_{a}<\mathbf{y}_{h}(p)<\mathbf{y}_{b}\}.$ It is clear that

\mathcal{V}^{C}_{h}=\mathcal{V}^{fC}_{h,a}\cup\mathcal{V}^{fC}_{h,b}\cup\mathcal{V}^{sC}_{h,a}\cup\mathcal{V}^{sC}_{h,b}\cup\mathcal{V}_{*}^{C}.

We define the residual $\mathbf{R}_{h}$ by the following:

\displaystyle(\mathbf{R}_{h},\mathbf{x})_{-1,1}=a(\mathbf{x},\bm{\phi}_{h})-b(\mathbf{x},r_{H})-(\mathbf{u}_{h}-\mathbf{u}_{d},\mathbf{x})-\rho a(\mathbf{y}_{h},\mathbf{x})\quad\text{for~{}all}\;\mathbf{x}\in\mathbf{Q}.

Also, for the a posteriori estimator we need to define the the Lagrange multiplier $\bm{\lambda}$ by the following:

\displaystyle(\bm{\lambda},\mathbf{x})_{-1,1}=a(\mathbf{x},\bm{\phi}_{h})-b(\mathbf{x},r_{H})-(\mathbf{u}_{h}-\mathbf{u}_{d},\mathbf{x})-\rho a(\mathbf{Ry},\mathbf{x})\quad\text{for~{}all}\;\mathbf{x}\in\mathbf{Q}.

Here, $(\cdot,\cdot)_{-1,1}$ denotes the duality pairing between $\mathbf{Q}$ and its dual. Clearly, $(\bm{\lambda},\mathbf{x}-\mathbf{Ry})_{-1,1}\leq 0$ for all $\mathbf{x}\in\mathbf{Q}_{ad}.$ Define the discrete Lagrange multiplier $\bm{\lambda}_{h}$ by the following:

\displaystyle(\bm{\lambda}_{h},\mathbf{x}_{h})_{-1,1}=a(\mathbf{x}_{h},\bm{\phi}_{h})-b(\mathbf{x}_{h},r_{H})-(\mathbf{u}_{h}-\mathbf{u}_{d},\mathbf{x}_{h})-\rho a(\mathbf{y}_{h},\mathbf{x}_{h})\quad\text{for~{}all}\;\mathbf{x}_{h}\in\mathbf{Q}_{h}.

It is clear that, $(\bm{\lambda}_{h},\mathbf{x}_{h})_{-1,1}=(\mathbf{R}_{h},\mathbf{x}_{h})_{-1,1}$ for all $\mathbf{x}_{h}\in\mathbf{Q}^{h}_{ad}$ and

\displaystyle(\bm{\lambda}_{h},\mathbf{x}_{h}-\mathbf{y}_{h})_{-1,1}\leq 0\quad\text{for~{}all}\;\mathbf{x}_{h}\in\mathbf{Q}^{h}_{ad}.

In order to investigate $\bm{\lambda}_{h}$ further, we use the partition of unity and integration by parts

	$\displaystyle(\bm{\lambda}_{h},\mathbf{x}_{h})_{-1,1}=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}(\bm{\lambda}_{h},x_{h,i}(p)\phi_{p}\mathbf{e}_{i})_{-1,1}$
	$\displaystyle=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}\int_{\omega_{p}}(-\Delta\bm{\phi}_{h}-\nabla r_{H}-(\mathbf{u}_{h}-\mathbf{u}_{d})+\rho\Delta\mathbf{y}_{h})\cdot x_{h,i}(p)\phi_{p}\mathbf{e}_{i}$
		$\displaystyle+\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}\int_{\gamma_{p,I}}[\hskip-1.5pt[\nabla\bm{\phi}_{h}+r_{H}\mathbb{I}-\rho\nabla\mathbf{y}_{h}]\hskip-1.5pt]\cdot x_{h,i}(p)\phi_{p}\mathbf{e}_{i}$
		$\displaystyle-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}^{C}_{h}}\int_{\gamma_{p,C}}\hat{\bm{\sigma}}(\mathbf{y}_{h})\cdot x_{h,i}(p)\phi_{p}\mathbf{e}_{i},$

where $\{\mathbf{e}_{1},\mathbf{e}_{2}\}$ standard basis of $\mathbb{R}^{2},$ and $\phi_{p}$ is the Lagrange basis at node $p$ . This motivates the representation $(\bm{\lambda}_{h},\mathbf{x}_{h})_{-1,1}=\sum_{i=1}^{2}(\lambda_{h,i},x_{h,i})_{-1,1}.$

Now for the a posteriori estimator we replace the residual $\mathbf{R}_{h}$ by a Galerkin functional, whose abstract definition is given by

\displaystyle(\mathbf{G}_{h},\mathbf{x})_{-1,1}:=\rho a(\mathbf{Ry}-\mathbf{y}_{h},\mathbf{x})+(\bm{\lambda}-\tilde{\bm{\lambda}}_{h},\mathbf{x})_{-1,1}

(4.11)

for all $\mathbf{x}\in\mathbf{Q},$ where $\tilde{\bm{\lambda}}_{h}$ is an approximation of $\bm{\lambda}_{h}$ and it depends on the discrete solution, data and reflects the properties of $\bm{\lambda},$ details can be found in [14]. We call it quasi-discrete contact force density.

By using the partition of unity we have the following definition of $\tilde{\bm{\lambda}}_{h},$ for all $\bm{\psi}\in\mathbf{Q}$ :

\displaystyle(\tilde{\bm{\lambda}}_{h},\bm{\psi})_{-1,1}:=\sum_{i=1}^{2}(\tilde{\lambda}_{h,i},\psi_{i})_{-1,1}=\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}}(\tilde{\lambda}^{p}_{h,i},\psi_{i}\phi_{p})_{-1,1}

(4.12)

and it adjusts the local contributions so that, on one hand, the Galerkin functional is prepared for the derivation of an upper bound, and on the other hand, tries to maximize the cancellation within $\mathbf{G}_{h}.$

For semi-contact nodes $p\in\mathcal{V}^{sC}_{h}:=\mathcal{V}^{sC}_{h,a}\cup\mathcal{V}^{sC}_{h,b},$ we define

\displaystyle(\tilde{\lambda}^{p}_{h,i},\psi_{i}\phi_{p})_{-1,1}:=(R_{h,i},\phi_{p})_{-1,1}m_{p}(\psi_{i})=s_{p,i}m_{p}(\psi_{i})\int_{\gamma_{p,C}}\phi_{p},

where $s_{p,i}:=\frac{(\lambda_{h,i},\phi_{p})_{-1,1}}{\int_{\gamma_{p,C}}\phi_{p}},$ $i=1,2.$ The constant $s_{p,i}$ is the nodal value of the discrete contact force density obtained by lumping the boundary mass matrix and $m_{p}(\psi_{i})$ is defined below. Sign of $s_{p,i}$ is the following:
i) For $p\in\Gamma_{C}$ with $y_{h,i}(p)=y^{i}_{a}$ then $s_{p,i}\leq 0.$
ii) For $p\in\Gamma_{C}$ with $y_{h,i}(p)=y^{i}_{b}$ then $s_{p,i}\geq 0.$
iii) For $p\in\mathcal{V}_{*}^{C}$ then $s_{p,i}=0.$

For full-contact nodes $p\in\mathcal{V}^{fC}_{h}:=\mathcal{V}^{fC}_{h,a}\cup\mathcal{V}^{fC}_{h,b},$ we define

	$\displaystyle(\tilde{\lambda}^{p}_{h,i},\psi_{i}\phi_{p})_{-1,1}:=$	$\displaystyle(\tilde{R}_{h,i},\phi_{p})_{-1,1}m_{p}(\psi_{i})-\int_{\gamma_{p,C}}\hat{\sigma}_{i}(\mathbf{y}_{h})\psi_{i}\phi_{p}$
	$\displaystyle=$	$\displaystyle s_{p,i}m_{p}(\psi_{i})\int_{\gamma_{p,C}}\phi_{p}-\int_{\gamma_{p,C}}\hat{\sigma}_{i}(\mathbf{y}_{h})(\psi_{i}-m_{p}(\psi_{i}))\phi_{p},$

where,

\displaystyle(\tilde{R}_{h,i},\theta)_{-1,1}:=(R_{h,i},\theta)_{-1,1}+\int_{\Gamma_{C}}\hat{\sigma}_{i}(\mathbf{y}_{h})\theta.

We need to specify the choices of $m_{p}(\psi_{i})$ for semi- and full-contact nodes.

For full-contact nodes $p\in\mathcal{V}^{fC}_{h,b}$ we use

\displaystyle m_{p}(\psi_{i})=\max_{e\subset\gamma_{p,C}}\frac{\int_{e}\psi_{i}\phi_{p}}{\int_{e}\phi_{p}}

(4.13)

and for full-contact nodes $p\in\mathcal{V}^{fC}_{h,a}$ we use

\displaystyle m_{p}(\psi_{i})=\min_{e\subset\gamma_{p,C}}\frac{\int_{e}\psi_{i}\phi_{p}}{\int_{e}\phi_{p}}.

(4.14)

This choice is important for the derivation of the upper bound, see for instance (LABEL:fourth_term) and (LABEL:sixth_term).

For semi-contact nodes $p\in\mathcal{V}^{sC}_{h,a}\cup\mathcal{V}^{sC}_{h,b},$ we set

\displaystyle m_{p}(\psi_{i})=\frac{\int_{\tilde{\gamma}_{p,C}}\psi_{i}\phi_{p}}{\int_{\tilde{\gamma}_{p,C}}\phi_{p}},

(4.15)

where $\tilde{\gamma}_{p,C}$ is a strict subset of $\gamma_{p,C}.$ Now we are ready to derive the a posteriori error estimators.

Theorem 4.6 (A posteriori error estimator).

It holds,

	$\displaystyle\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\Omega}+$	$\displaystyle\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\Omega}+\lVert p-p_{H}\rVert_{0,\Omega}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\Omega}$
		$\displaystyle+\lVert r-r_{H}\rVert_{0,\Omega}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\Omega}\lesssim\big{(}\eta_{(\mathbf{u},p)}+\eta_{(\bm{\phi},r)}+\eta_{\mathbf{y}}\big{)},$

where the estimators are defined as

	$\displaystyle\eta_{(\mathbf{u},p)}^{2}=$	$\displaystyle\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\lVert\mathbf{f}+\Delta\mathbf{u}_{h}-\nabla p_{H}\rVert^{2}_{0,T}+\sum_{e\in\mathcal{E}_{h}^{i}}\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[p_{H}\mathbb{I}-\nabla{\mathbf{u}_{h}}]\hskip-1.5pt]\rVert_{0,e}^{2}$
		$\displaystyle+\sum_{e\in\mathcal{E}_{h}^{b,N}}\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[p_{H}\mathbb{I}-\nabla{\mathbf{u}_{h}}]\hskip-1.5pt]\rVert_{0,e}^{2}+\sum_{T\in\mathcal{T}_{h}}\left\\|\nabla\cdot\mathbf{u}_{h}\right\\|_{0,T}^{2}$

and

	$\displaystyle\eta_{(\bm{\phi},r)}^{2}=$	$\displaystyle\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\lVert\Delta\bm{\phi}_{h}+\nabla r_{H}+\mathbf{u}_{h}-\mathbf{u}_{d}\rVert^{2}_{0,T}+\sum_{e\in\mathcal{E}_{h}^{i}}\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[r_{H}\mathbb{I}+\nabla{\bm{\phi}_{h}}]\hskip-1.5pt]\rVert_{0,e}^{2}$
		$\displaystyle+\sum_{e\in\mathcal{E}_{h}^{b,N}}\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[r_{H}\mathbb{I}+\nabla{\bm{\phi}_{h}}]\hskip-1.5pt]\rVert_{0,e}^{2}+\sum_{T\in\mathcal{T}_{h}}\left\\|\nabla\cdot\bm{\phi}_{h}\right\\|_{0,T}^{2}$

	$\displaystyle\eta_{\mathbf{y}}^{2}=$	$\displaystyle\sum_{p\in\mathcal{V}_{h}}\big{(}\underbrace{h_{p}^{2}\left\\|\bm{\mathcal{R}}(\mathbf{y}_{h})\right\\|^{2}_{0,\omega_{p}}}_{\eta_{\mathbf{y}}^{a}}+\underbrace{h_{p}\left\\|\bm{\mathcal{J}}^{I}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,I}}^{2}}_{\eta_{\mathbf{y}}^{b}}\big{)}+\sum_{p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}}\underbrace{h_{p}\left\\|\hat{\bm{\sigma}}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}^{2}}_{\eta_{\mathbf{y}}^{c}}$
		$\displaystyle+\sum_{i=1}^{2}\big{[}\sum_{p\in\mathcal{V}^{sC}_{h,a}}\underbrace{s_{p,i}d^{a}_{p,i}}_{\eta_{\mathbf{y}}^{d}}+\sum_{p\in\mathcal{V}^{sC}_{h,b}}\underbrace{s_{p,i}d^{b}_{p,i}}_{\eta_{\mathbf{y}}^{e}}\big{]},$

where, $\bm{\mathcal{R}}(\mathbf{y}_{h}):=(\mathcal{R}_{1}(\mathbf{y}_{h}),\mathcal{R}_{2}(\mathbf{y}_{h}))=-\Delta\bm{\phi}_{h}-\nabla r_{H}-(\mathbf{u}_{h}-\mathbf{u}_{d})+\rho\Delta\mathbf{y}_{h},$ $\bm{\mathcal{J}}^{I}(\mathbf{y}_{h}):=(\mathcal{J}^{I}_{1}(\mathbf{y}_{h}),\mathcal{J}^{I}_{2}(\mathbf{y}_{h}))=[\hskip-1.5pt[\nabla\bm{\phi}_{h}+r_{H}\mathbb{I}-\rho\nabla\mathbf{y}_{h}]\hskip-1.5pt],$ $d^{a}_{p,i}=\int_{\tilde{\gamma}_{p,C}}(y^{i}_{a}-y_{h,i})\phi_{p}.$ and $d^{b}_{p,i}=\int_{\tilde{\gamma}_{p,C}}(y^{i}_{b}-y_{h,i})\phi_{p},$ $\mathbb{I}$ be the $2\times 2$ identity matrix. and the set $\tilde{\gamma}_{p,C}$ is a strict subset of $\gamma_{p,C}.$

Proof.

From Theorem 4.5, we have

$\displaystyle\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\Omega}+$	$\displaystyle\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\Omega}+\lVert p-p_{H}\rVert_{0,\Omega}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\Omega}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\Omega}$
	$\displaystyle+\lVert r-r_{H}\rVert_{0,\Omega}\lesssim\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\Omega}+\left\\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\\|_{0,\Omega}$
	$\displaystyle+\left\\|p_{H}-R_{0}p\right\\|_{0,\Omega}+\left\\|\nabla(\bar{\mathbf{R}}\bm{\phi}-\bm{\phi}_{h})\right\\|_{0,\Omega}+\left\\|r_{H}-\bar{R}_{0}r\right\\|_{0,\Omega}.$	(4.16)

The a posteriori error analysis in [25, Theorem 3.1] and [3, Section 7] gives the following error estimators of state and adjoint state variables:

	$\displaystyle\left\\|\nabla(\mathbf{R}\mathbf{u}-\mathbf{u}_{h})\right\\|^{2}_{0,\Omega}+\left\\|p_{H}-R_{0}p\right\\|^{2}_{0,\Omega}\lesssim$	$\displaystyle\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\lVert\mathbf{f}+\Delta\mathbf{u}_{h}-\nabla p_{H}\rVert^{2}_{0,T}+\sum_{e\in\mathcal{E}_{h}^{i}}\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[p_{H}\mathbb{I}-\nabla{\mathbf{u}_{h}}]\hskip-1.5pt]\rVert_{0,e}^{2}$
		$\displaystyle+\sum_{e\in\mathcal{E}_{h}^{b,N}}\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[p_{H}\mathbb{I}-\nabla{\mathbf{u}_{h}}]\hskip-1.5pt]\rVert_{0,e}^{2}+\sum_{T\in\mathcal{T}_{h}}\left\\|\nabla\cdot\mathbf{u}_{h}\right\\|_{0,T}^{2}$		(4.17)

	$\displaystyle\left\\|\nabla(\bar{\mathbf{R}}\bm{\phi}-\bm{\phi}_{h})\right\\|^{2}_{0,\Omega}+\left\\|r_{H}-\bar{R}_{0}r\right\\|^{2}_{0,\Omega}\lesssim$	$\displaystyle\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\lVert\Delta\bm{\phi}_{h}+\nabla r_{H}+\mathbf{u}_{h}-\mathbf{u}_{d}\rVert^{2}_{0,T}+\sum_{e\in\mathcal{E}_{h}^{i}}\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[r_{H}\mathbb{I}+\nabla{\bm{\phi}_{h}}]\hskip-1.5pt]\rVert_{0,e}^{2}$
		$\displaystyle+\sum_{e\in\mathcal{E}_{h}^{b,N}}\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[r_{H}\mathbb{I}+\nabla{\bm{\phi}_{h}}]\hskip-1.5pt]\rVert_{0,e}^{2}+\sum_{T\in\mathcal{T}_{h}}\left\\|\nabla\cdot\bm{\phi}_{h}\right\\|_{0,T}^{2}.$		(4.18)

We denote the right hand side of (4) by $\eta^{2}_{(\mathbf{u},p)}$ and the right hand side of (4) by $\eta^{2}_{(\bm{\phi},r)}.$ The splitting $\mathbf{Ru}=\mathbf{Rw}+\mathbf{y}_{h}$ and $\mathbf{u}_{h}=\mathbf{w}_{h}+\mathbf{y}_{h},$ yields $\left\|\nabla(\mathbf{R}\mathbf{u}-\mathbf{u}_{h})\right\|_{0,\Omega}=\left\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\|_{0,\Omega}.$ Thus we have estimates for all the last four terms of (4). Therefore we only need to estimate $\left\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\|_{0,\Omega}+\left\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\|_{-1,\Omega}.$ Estimator for the control variable needs some special care. To find an upper bound of the term $\left\|\nabla(\mathbf{Ry}-\mathbf{y}_{h})\right\|_{0,\Omega},$ we choose $\mathbf{x}=\mathbf{Ry}-\mathbf{y}_{h}$ in (4.11)

$\displaystyle\rho\left\\|\nabla(\mathbf{Ry}-\mathbf{y}_{h})\right\\|_{0,\Omega}^{2}$	$\displaystyle=\rho\;a(\mathbf{Ry}-\mathbf{y}_{h},\mathbf{Ry}-\mathbf{y}_{h})$
	$\displaystyle=\mathbf{G}_{h}(\mathbf{Ry}-\mathbf{y}_{h})-(\bm{\lambda}-\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}$
	$\displaystyle\leq\left\\|\mathbf{G}_{h}\right\\|_{-1}\left\\|\nabla(\mathbf{Ry}-\mathbf{y}_{h})\right\\|_{0,\Omega}-(\bm{\lambda}-\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}.$	(4.19)

Applying the Young’s inequality in (4.19), we get the following estimates:

\displaystyle\frac{\rho}{2}\left\|\nabla(\mathbf{Ry}-\mathbf{y}_{h})\right\|_{0,\Omega}^{2}\leq

\displaystyle\frac{1}{2\rho}\left\|\mathbf{G}_{h}\right\|^{2}_{-1}-(\bm{\lambda}-\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}.

(4.20)

Also, using (4.11) one can easily derive

\displaystyle\left\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\|_{-1,\Omega}^{2}\leq(2+\frac{2}{\rho})\left\|\mathbf{G}_{h}\right\|^{2}_{-1}-4(\bm{\lambda}-\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}.

(4.21)

Adding (4.20) and (4.21), we get

\displaystyle\left\|\nabla(\mathbf{Ry}-\mathbf{y}_{h})\right\|_{0,\Omega}^{2}+\left\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\|_{-1,\Omega}^{2}\leq(2+\frac{2}{\rho}+\frac{1}{\rho^{2}})\left\|\mathbf{G}_{h}\right\|^{2}_{-1}-(4+\frac{2}{\rho})(\bm{\lambda}-\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}.

(4.22)

From the equation (4.22) it is clear that the error in control and contact force densities are bounded by the dual norm of the Galerkin functional $\left\|\mathbf{G}_{h}\right\|_{-1}$ and the duality pairing between the contact force densities and the controls.

First, we estimate the term $\left\|\mathbf{G}_{h}\right\|_{-1}.$ Using the definition of $\mathbf{G}_{h}$ in (4.11), with $\bm{\psi}\in\mathbf{Q}$ , we obtain

$\displaystyle(\mathbf{G}_{h},\bm{\psi})_{-1,1}=$	$\displaystyle\rho a(\mathbf{Ry}-\mathbf{y}_{h},\bm{\psi})+(\bm{\lambda}-\tilde{\bm{\lambda}}_{h},\bm{\psi})_{-1,1}$
$\displaystyle=$	$\displaystyle\rho a(\mathbf{Ry}-\mathbf{y}_{h},\bm{\psi})+a(\bm{\psi},\bm{\phi}_{h})-b(\bm{\psi},r_{H})-(\mathbf{u}_{h}-\mathbf{u}_{d},\bm{\psi})$
	$\displaystyle-\rho a(\mathbf{Ry},\bm{\psi})-(\tilde{\bm{\lambda}}_{h},\bm{\psi})_{-1,1}$
$\displaystyle=$	$\displaystyle a(\bm{\psi},\bm{\phi}_{h})-b(\bm{\psi},r_{H})-(\mathbf{u}_{h}-\mathbf{u}_{d},\bm{\psi})-\rho a(\mathbf{y}_{h},\bm{\psi})-(\tilde{\bm{\lambda}}_{h},\bm{\psi})_{-1,1}$
$\displaystyle=$	$\displaystyle(\mathbf{R}_{h},\bm{\psi})_{-1,1}-(\tilde{\bm{\lambda}}_{h},\bm{\psi})_{-1,1}$
$\displaystyle=$	$\displaystyle\sum_{i=1}^{2}(R_{h,i},\psi_{i})_{-1,1}-\sum_{i=1}^{2}(\tilde{\lambda}_{h,i},\psi_{i})_{-1,1}$
$\displaystyle=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}(R_{h,i},\psi_{i}\phi_{p})_{-1,1}-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}}(\tilde{\lambda}^{p}_{h,i},\psi_{i}\phi_{p})_{-1,1}$
$\displaystyle=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}\setminus\mathcal{V}_{h}^{C}}(R_{h,i},(\psi_{i}-m_{p}(\psi_{i}))\phi_{p})_{-1,1}$
	$\displaystyle+\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}}(R_{h,i},(\psi_{i}-m_{p}(\psi_{i}))\phi_{p})_{-1,1}-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}}(\tilde{\lambda}^{p}_{h,i},\psi_{i}\phi_{p})_{-1,1}$
$\displaystyle=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}\setminus\mathcal{V}_{h}^{C}}(R_{h,i},(\psi_{i}-m_{p}(\psi_{i}))\phi_{p})_{-1,1}$
	$\displaystyle+\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}}(\tilde{R}_{h,i},(\psi_{i}-m_{p}(\psi_{i}))\phi_{p})_{-1,1}$
	$\displaystyle-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}}\int_{\Gamma_{C}}\hat{\sigma}_{i}(\mathbf{y}_{h})(\psi_{i}-m_{p}(\psi_{i}))\phi_{p}-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}}(\tilde{\lambda}^{p}_{h,i},\psi_{i}\phi_{p})_{-1,1}$
$\displaystyle=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}\setminus\mathcal{V}_{h}^{C}}(R_{h,i},(\psi_{i}-m_{p}(\psi_{i}))\phi_{p})_{-1,1}$
	$\displaystyle+\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}}(\tilde{R}_{h,i},(\psi_{i}-m_{p}(\psi_{i}))\phi_{p})_{-1,1}$
	$\displaystyle-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}^{C}\setminus\mathcal{V}^{fC}_{h}}\int_{\Gamma_{C}}\hat{\sigma}_{i}(\mathbf{y}_{h})(\psi_{i}-m_{p}(\psi_{i}))\phi_{p}.$	(4.23)

Here, we set $m_{p}(\psi_{i})=0$ for Dirichlet and Neumann nodes. We exploited $(R_{h,i},\phi_{p})_{-1,1}=0$ for all non-contact nodes and we inserted the definition of $\tilde{\bm{\lambda}}_{h}.$ Inserting definition of $R_{h,i}$ and $\tilde{R}_{h,i}$ in (4.23), we get

$\displaystyle(\mathbf{G}_{h},\bm{\psi})_{-1,1}=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}\int_{\omega_{p}}\mathcal{R}_{i}(\mathbf{y}_{h})(\psi_{i}-m_{p}(\psi_{i}))\phi_{p}$
	$\displaystyle+\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}\int_{\gamma_{p,I}}\mathcal{J}^{I}_{i}(\mathbf{y}_{h})(\psi_{i}-m_{p}(\psi_{i}))\phi_{p}$
	$\displaystyle-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}}\int_{\Gamma_{C}}\hat{\sigma}_{i}(\mathbf{y}_{h})(\psi_{i}-m_{p}(\psi_{i}))\phi_{p}.$	(4.24)

where $\mathcal{R}_{i}(\mathbf{y}_{h})=-\Delta\phi_{h,i}-\frac{\partial}{\partial x_{i}}r_{H}-(u_{h,i}-u_{d,i})+\rho\Delta y_{h,i}$ and $\mathcal{J}^{I}_{i}(\mathbf{y}_{h})=[\hskip-1.5pt[\nabla\phi_{h,i}+\begin{pmatrix}r_{H}\\ 0\end{pmatrix}^{\prime}-\rho\nabla\mathbf{y}_{h}]\hskip-1.5pt].$ For all nodes $p\in\mathcal{V}_{h}\setminus(\mathcal{V}^{D}_{h}\cup\mathcal{V}^{C}_{h}\cup\mathcal{V}^{N}_{h})$ we choose the constants

\displaystyle m_{p}(\psi_{i})=\frac{\int_{\omega_{p}}\psi_{i}\phi_{p}}{\int_{\omega_{P}}\phi_{p}}.

(4.25)

The mean value (4.25) satisfies the following $L^{2}-$ approximation properties:

	$\displaystyle\left\\|\psi_{i}-m_{p}(\psi_{i})\right\\|_{0,\omega_{p}}$	$\displaystyle\leq Ch_{p}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}},$
	$\displaystyle\left\\|\psi_{i}-m_{p}(\psi_{i})\right\\|_{0,\gamma_{p}}$	$\displaystyle\leq Ch_{p}^{\frac{1}{2}}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}}$

can be found in [31, Lemma 3.1, Proposition 4.2]. For Dirichlet and Neumann nodes we have at least one edge $e\subset\gamma_{p,C}$ , where the test function $\psi_{i}$ is zero, therefore we can deduce $\left\|\psi_{i}\right\|_{0,\omega_{p}}\leq Ch_{p}\left\|\nabla\psi_{i}\right\|_{0,\omega_{p}}$ directly from the Poincaré-Friedrichs inequality. The above $L^{2}-$ approximation properties hold also for the constants $m_{p}(\psi_{i})$ defined in (4.14), (4.13) and (4.15) see [31, Lemma 3.1, Proposition 4.2]. Using the above estimates and applying the Cauchy-Schwarz inequality in (4) we arrive at

$\displaystyle(\mathbf{G}_{h},\bm{\psi})_{-1,1}=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}h_{p}\left\\|\mathcal{R}_{i}(\mathbf{y}_{h})\right\\|_{0,\omega_{p}}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}}$
	$\displaystyle+\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}h^{\frac{1}{2}}_{p}\left\\|\mathcal{J}^{I}_{i}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,I}}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}}$
	$\displaystyle-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}}h_{p}^{\frac{1}{2}}\left\\|\hat{\sigma}_{i}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}}$
	$\displaystyle\lesssim\left(\sum_{p\in\mathcal{V}_{h}}\big{(}h_{p}^{2}\left\\|\mathcal{R}(\mathbf{y}_{h})\right\\|^{2}_{0,\omega_{p}}+h_{p}\left\\|\mathcal{J}^{I}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,I}}^{2}\big{)}\right.$
	$\displaystyle\left.+\sum_{p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}}h_{p}\left\\|\hat{\bm{\sigma}}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}^{2}\right)^{\frac{1}{2}}\left\\|\nabla\bm{\psi}\right\\|.$	(4.26)

Thus from (4), we have

\displaystyle\left\|\mathbf{G}_{h}\right\|_{-1}\lesssim\left(\sum_{p\in\mathcal{V}_{h}}\big{(}h_{p}^{2}\left\|\bm{\mathcal{R}}(\mathbf{y}_{h})\right\|^{2}_{0,\omega_{p}}+h_{p}\left\|\bm{\mathcal{J}}^{I}(\mathbf{y}_{h})\right\|_{0,\gamma_{p,I}}^{2}\big{)}+\sum_{p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}}h_{p}\left\|\hat{\bm{\sigma}}(\mathbf{y}_{h})\right\|_{0,\gamma_{p,C}}^{2}\right)^{\frac{1}{2}}

(4.27)

Now, we need to find the upper bound of the term $(\tilde{\bm{\lambda}}_{h}-\bm{\lambda},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}.$ We can write

\displaystyle(\tilde{\bm{\lambda}}_{h}-\bm{\lambda},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}=(\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}+(\bm{\lambda},\mathbf{y}_{h}-\mathbf{Ry})_{-1,1}.

(4.28)

From the equation (4.1e), we have $(\bm{\lambda},\mathbf{x}-\mathbf{Ry})_{-1,1}\leq 0$ for all $x\in\mathbf{Q}_{ad},$ taking $\mathbf{x}=\mathbf{y}_{h}$ we have $(\bm{\lambda},\mathbf{y}_{h}-\mathbf{Ry})_{-1,1}\leq 0.$ Thus it is enough to estimate the term $(\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}.$ From the definition of the quasi discrete contact force density (4.12)-(4.15), we have

$\displaystyle(\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}^{C}_{h}}(\tilde{\lambda}_{h,i}^{p},(R_{i}\mathbf{y}-y_{h,i})\phi_{p})_{-1,1}$
$\displaystyle=$	$\displaystyle\sum_{i=1}^{2}\left[\sum_{p\in\mathcal{V}^{sC}_{h,a}}s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\int_{\gamma_{p,C}}\phi_{p}+\sum_{p\in\mathcal{V}^{sC}_{h,b}}s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\int_{\gamma_{p,C}}\phi_{p}\right.$
	$\displaystyle\left.+\sum_{p\in\mathcal{V}^{fC}_{h,a}}s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\int_{\gamma_{p,C}}\phi_{p}+\sum_{p\in\mathcal{V}^{fC}_{h,b}}s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\int_{\gamma_{p,C}}\phi_{p}\right.$
	$\displaystyle\left.-\sum_{p\in\mathcal{V}^{fC}_{h,a}}\int_{\gamma_{p,C}}\hat{\sigma}_{i}(\mathbf{y}_{h})\big{(}R_{i}\mathbf{y}-y_{h,i}-m_{p}(R_{i}\mathbf{y}-y_{h,i})\big{)}\phi_{p}\right.$
	$\displaystyle\left.-\sum_{p\in\mathcal{V}^{fC}_{h,b}}\int_{\gamma_{p,C}}\hat{\sigma}_{i}(\mathbf{y}_{h})\big{(}R_{i}\mathbf{y}-y_{h,i}-m_{p}(R_{i}\mathbf{y}-y_{h,i})\big{)}\phi_{p}\right].$	(4.29)

We need to bound each term in the right hand side of (4).
First term:

	$\displaystyle\sum_{p\in\mathcal{V}^{sC}_{h,a}}s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\int_{\gamma_{p,C}}\phi_{p}$	$\displaystyle=\sum_{p\in\mathcal{V}^{sC}_{h,a}}s_{p,i}\int_{\gamma_{p,C}}\phi_{p}\frac{\int_{\tilde{\gamma}_{p,C}}(R_{i}\mathbf{y}-y_{h,i})\phi_{p}}{\int_{\tilde{\gamma}_{p,C}}\phi_{p}}$
		$\displaystyle=\sum_{p\in\mathcal{V}^{sC}_{h,a}}s_{p,i}\int_{\gamma_{p,C}}\phi_{p}\frac{\int_{\tilde{\gamma}_{p,C}}(R_{i}\mathbf{y}-y^{i}_{a})\phi_{p}}{\int_{\tilde{\gamma}_{p,C}}\phi_{p}}$
		$\displaystyle+\sum_{p\in\mathcal{V}^{sC}_{h,a}}s_{p,i}\int_{\gamma_{p,C}}\phi_{p}\frac{\int_{\tilde{\gamma}_{p,C}}(y^{i}_{a}-y_{h,i})\phi_{p}}{\int_{\tilde{\gamma}_{p,C}}\phi_{p}}$
		$\displaystyle\leq\sum_{p\in\mathcal{V}^{sC}_{h,a}}s_{p,i}\int_{\gamma_{p,C}}\phi_{p}\frac{\int_{\tilde{\gamma}_{p,C}}(y^{i}_{a}-y_{h,i})\phi_{p}}{\int_{\tilde{\gamma}_{p,C}}\phi_{p}}=\sum_{p\in\mathcal{V}^{sC}_{h,a}}s_{p,i}d^{a}_{p,i},$

where $d^{a}_{p,i}=\int_{\tilde{\gamma}_{p,C}}(y^{i}_{a}-y_{h,i})\phi_{p}.$ In the above we exploit $s_{p,i}\leq 0,$ $R_{i}\mathbf{y}-y^{i}_{a}\geq 0$ and $\frac{\int_{\gamma_{p,C}}\phi_{p}}{\int_{\tilde{\gamma}_{p,C}}\phi_{p}}$ is a constant independent of $h_{p}$ if $\tilde{\gamma}_{p,C}$ is always a fixed fraction of ${\gamma}_{p,C}.$

Second term: A similar arguments from the first term we have

\displaystyle\sum_{p\in\mathcal{V}^{sC}_{h,b}}s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\int_{\gamma_{p,C}}\phi_{p}

\displaystyle\leq\sum_{p\in\mathcal{V}^{sC}_{h,b}}s_{p,i}d^{b}_{p,i},

where $d^{b}_{p,i}=\int_{\tilde{\gamma}_{p,C}}(y^{i}_{b}-y_{h,i})\phi_{p}.$

Third term: For full contact nodes $p\in\mathcal{V}^{fC}_{h,a},$ we have $R_{i}\mathbf{y}\geq y^{i}_{a}=y_{h,i}$ on $\gamma_{p,C}$ which implies $R_{i}\mathbf{y}-y_{h,i}\geq 0$ on $\gamma_{p,C}$ and therefore $m_{p}(R_{i}\mathbf{y}-y_{h,i})\geq 0.$ As further $s_{p,i}\leq 0$ we have $s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\leq 0.$ Hence,

\displaystyle\sum_{p\in\mathcal{V}^{fC}_{h,a}}s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\int_{\gamma_{p,C}}\phi_{p}\leq 0.

Fourth term: Similar to the third term for $p\in\mathcal{V}^{fC}_{h,b}$ we have $s_{p,i}\geq 0$ and $R_{i}\mathbf{y}-y_{h,i}\leq 0$ and hence $m_{p}(R_{i}\mathbf{y}-y_{h,i})\leq 0.$ Therefore

\displaystyle\sum_{p\in\mathcal{V}^{fC}_{h,b}}s_{p,i}m_{p}(R_{i}\mathbf{y}-y_{h,i})\int_{\gamma_{p,C}}\phi_{p}\leq 0.

Fifth term:

	$\displaystyle-$	$\displaystyle\int_{\gamma_{p,C}}\hat{\sigma}_{i}(\mathbf{y}_{h})\big{(}R_{i}\mathbf{y}-y_{h,i}-m_{p}(R_{i}\mathbf{y}-y_{h,i})\big{)}\phi_{p}$
		$\displaystyle=\sum_{e\subset\gamma_{p,C}}-\hat{\sigma}_{i}(\mathbf{y}_{h})\|_{e}\int_{e}\big{(}R_{i}\mathbf{y}-y_{h,i}-m_{p}(R_{i}\mathbf{y}-y_{h,i})\big{)}\phi_{p}.$

Here we have used the fact that $\hat{\sigma}_{i}(\mathbf{y}_{h})|_{e}$ is piecewise constant on each edge. Since $p\in\mathcal{V}^{fC}_{h,a}$ we have $-\hat{\sigma}_{i}(\mathbf{y}_{h})|_{e}\leq 0$ and using the definition of $m_{p}$ from (4.14) we get

\int_{e}\big{(}R_{i}\mathbf{y}-y_{h,i}-m_{p}(R_{i}\mathbf{y}-y_{h,i})\big{)}\phi_{p}\geq 0.

Hence,

\displaystyle-\int_{\gamma_{p,C}}\hat{\sigma}_{i}(\mathbf{y}_{h})\big{(}R_{i}\mathbf{y}-y_{h,i}-m_{p}(R_{i}\mathbf{y}-y_{h,i})\big{)}\phi_{p}\leq 0.

Sixth term: Similar to the fourth term for $p\in\mathcal{V}^{fC}_{h,b}$ we have $-\hat{\sigma}_{i}(\mathbf{y}_{h})|_{e}\geq 0$ and using the definition of $m_{p}$ from (4.13) we get

\int_{e}\big{(}R_{i}\mathbf{y}-y_{h,i}-m_{p}(R_{i}\mathbf{y}-y_{h,i})\big{)}\phi_{p}\leq 0.

Hence,

\displaystyle-\int_{\gamma_{p,C}}\hat{\sigma}_{i}(\mathbf{y}_{h})\big{(}R_{i}\mathbf{y}-y_{h,i}-m_{p}(R_{i}\mathbf{y}-y_{h,i})\big{)}\phi_{p}\leq 0.

Using the all the above estimates from First term to Sixth term, in the right hand side of the equation (4), we obtain

\displaystyle(\tilde{\bm{\lambda}}_{h},\mathbf{Ry}-\mathbf{y}_{h})_{-1,1}\leq

\displaystyle\sum_{i=1}^{2}\big{[}\sum_{p\in\mathcal{V}^{sC}_{h,a}}s_{p,i}d^{a}_{p,i}+\sum_{p\in\mathcal{V}^{sC}_{h,b}}s_{p,i}d^{b}_{p,i}\big{]}.

(4.30)

Substituting (4.27), (4.28), and (4.30) in the right hand side of (4.22), we obtain the following upper bound of control error:

	$\displaystyle\left\\|\nabla(\mathbf{Ry}-\mathbf{y}_{h})\right\\|_{0,\Omega}^{2}+$	$\displaystyle\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\Omega}^{2}\lesssim\sum_{p\in\mathcal{V}_{h}}\big{(}\underbrace{h_{p}^{2}\left\\|\bm{\mathcal{R}}(\mathbf{y}_{h})\right\\|^{2}_{0,\omega_{p}}}_{\eta_{\mathbf{y}}^{a}}+\underbrace{h_{p}\left\\|\bm{\mathcal{J}}^{I}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,I}}^{2}}_{\eta_{\mathbf{y}}^{b}}\big{)}$
		$\displaystyle+\sum_{p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}}\underbrace{h_{p}\left\\|\hat{\bm{\sigma}}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}^{2}}_{\eta_{\mathbf{y}}^{c}}+\sum_{i=1}^{2}\big{[}\sum_{p\in\mathcal{V}^{sC}_{h,a}}\underbrace{s_{p,i}d^{a}_{p,i}}_{\eta_{\mathbf{y}}^{d}}+\sum_{p\in\mathcal{V}^{sC}_{h,b}}\underbrace{s_{p,i}d^{b}_{p,i}}_{\eta_{\mathbf{y}}^{e}}\big{]}.$		(4.31)

Denote the right hand side of (4) to be $\eta^{2}_{\mathbf{y}}.$ Thus substituting the estimates (4), (4), and (4) in (4), we prove Theorem 4.6. ∎

Remark 4.7.

The term $\eta_{\mathbf{y}}^{e}$ reminds of a complementarity condition. In fact, for a semi-contact node $s_{p,i}d_{p,i}^{b}$ would be a complementarity condition with respect to the quasi-discrete contact force density $(\tilde{\lambda}_{h,i}^{p},(y_{b}^{i}-y_{h,i})\phi_{p})_{-1,1}$ if $\tilde{\gamma}_{p,C}$ was replaced by $\gamma_{p,C}$ . Thus we refer to $\eta_{\mathbf{y}}^{d}$ and $\eta_{\mathbf{y}}^{e}$ as complementarity residual and call $\eta_{\mathbf{y}}^{c}$ contact stress residual. The contributions $\eta_{\mathbf{y}}^{c},$ $\eta_{\mathbf{y}}^{d}$ are localized to semi-contact nodes and nodes which are not actually in contact. In the unconstrained case, we have $\eta_{\mathbf{y}}^{d}=\eta_{\mathbf{y}}^{e}=0$ and $\eta_{\mathbf{y}}^{c}$ has contributions from all potential contact nodes such that $\eta_{\mathbf{y}}$ is a residual error estimator for linear elliptic boundary value problems where the potential contact boundary is replaced by a Neumann boundary with zero Neumann data.

Theorem 4.8 (Local Efficiency).

Let $\mathcal{T}_{e}$ be the set of two triangles sharing the edge $e\in\mathcal{E}_{h}^{i}.$ Then, it hold

$\displaystyle h_{T}\lVert\mathbf{f}+\Delta\mathbf{u}_{h}-\nabla p_{H}\rVert_{0,T}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,T}+\lVert p-p_{H}\rVert_{0,T}$
	$\displaystyle+osc(\mathbf{f},T)\big{)},$
$\displaystyle h_{T}\lVert\Delta\bm{\phi}_{h}+\nabla r_{H}+\mathbf{u}_{h}-\mathbf{u}_{d}\rVert_{0,T}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,T}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,T}$
	$\displaystyle+\lVert r-r_{H}\rVert_{0,T}+osc(\mathbf{u}_{d},T)\big{)},$
$\displaystyle h_{p}\left\\|\bm{\mathcal{R}}(\mathbf{y}_{h})\right\\|_{0,\omega_{p}}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}$
	$\displaystyle+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}+\lVert r-r_{H}\rVert_{0,\omega_{p}}$
	$\displaystyle+osc(\mathbf{u}_{d},\omega_{p})\big{)}$	(4.32)
$\displaystyle\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[p_{H}I-\nabla{\mathbf{u}_{h}}]\hskip-1.5pt]\rVert_{0,e}\lesssim$	$\displaystyle\sum_{T\in\mathcal{T}_{e}}\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,T}+\lVert p-p_{H}\rVert_{0,T}$
	$\displaystyle+osc(\mathbf{f},T)\big{)},$
$\displaystyle\lVert h^{\frac{1}{2}}[\hskip-1.5pt[r_{H}I+\nabla{\bm{\phi}_{h}}]\hskip-1.5pt]\rVert_{0,e}\lesssim$	$\displaystyle\sum_{T\in\mathcal{T}_{e}}\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,T}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,T}$
	$\displaystyle+\lVert r-r_{H}\rVert_{0,T}+osc(\mathbf{u}_{d},T)\big{)},$
$\displaystyle\lVert h_{p}^{\frac{1}{2}}\bm{\mathcal{J}}^{I}(\mathbf{y}_{h})\rVert_{0,\gamma_{p,I}}\lesssim$	$\displaystyle\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,T}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}$
	$\displaystyle+\lVert r-r_{H}\rVert_{0,\omega_{p}}+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}$
	$\displaystyle+osc(\mathbf{u}_{d},\omega_{p}),$	(4.33)
$\displaystyle\left\\|\nabla\cdot\mathbf{u}_{h}\right\\|_{0,T}\lesssim$	$\displaystyle\left\\|\nabla(\mathbf{u}-\mathbf{u}_{h})\right\\|_{0,T}$
$\displaystyle\left\\|\nabla\cdot\bm{\phi}_{h}\right\\|_{0,T}\lesssim$	$\displaystyle\left\\|\nabla(\bm{\phi}-\bm{\phi}_{h})\right\\|_{0,T}.$

Further, for any Neumann boundary edge $e\in\mathcal{E}_{h}^{b,N}$ , it hold

	$\displaystyle\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[p_{H}I-\nabla{\mathbf{u}_{h}}]\hskip-1.5pt]\rVert_{0,e}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,T}+\lVert p-p_{H}\rVert_{0,T}+osc(\mathbf{f},T)\big{)},$
	$\displaystyle\lVert h_{e}^{\frac{1}{2}}[\hskip-1.5pt[r_{H}I+\nabla{\bm{\phi}_{h}}]\hskip-1.5pt]\rVert_{0,e}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,T}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,T}+\lVert r-r_{H}\rVert_{0,T}$
		$\displaystyle+osc(\mathbf{u}_{d},T)\big{)}.$

In the contact zone, for $p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}$

$\displaystyle h^{\frac{1}{2}}_{p}\lVert\hat{\bm{\sigma}}(\mathbf{y}_{h})\rVert$	${}_{0,\gamma_{p,C}}\lesssim\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}$
	$\displaystyle+\lVert r-r_{H}\rVert_{0,\omega_{p}}+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}$
	$\displaystyle+osc(\mathbf{u}_{d},\omega_{p})\big{)},$	(4.34)

for $p\in\mathcal{V}^{sC}_{h,a}$

	$\displaystyle\sum_{i=1}^{2}s_{p,i}d^{a}_{p,i}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}+\lVert r-r_{H}\rVert_{0,\omega_{p}}$
		$\displaystyle+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}+osc(\mathbf{u}_{d},\omega_{p})\big{)}^{2}$		(4.35)

for $p\in\mathcal{V}^{sC}_{h,b}$

	$\displaystyle\sum_{i=1}^{2}s_{p,i}d^{b}_{p,i}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}+\lVert r-r_{H}\rVert_{0,\omega_{p}}$
		$\displaystyle+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}+osc(\mathbf{u}_{d},\omega_{p})\big{)}^{2}$		(4.36)

where oscillation of a given function $\mathbf{f},\mathbf{u}_{d}\in\mathbf{L}^{2}(T)$ is defined by

osc(\mathbf{f},T)=h_{T}\min_{\mathbf{f}_{h}\in\mathbf{P}_{0}(T)}\left\|\mathbf{f}-\mathbf{f}_{h}\right\|_{0,T},osc(\mathbf{u}_{d},T)=h_{T}\min_{\mathbf{g}_{h}\in\mathbf{P}_{0}(T)}\left\|\mathbf{u}_{d}-\mathbf{g}_{h}\right\|_{0,T},

and similarly, we define the oscillation, $osc(\mathbf{u}_{d},\omega_{p})=h_{T}\min_{\mathbf{g}_{h}\in\mathbf{P}_{0}(\omega_{p})}\left\|\mathbf{u}_{d}-\mathbf{g}_{h}\right\|_{0,\omega_{p}}.$

Proof.

The local efficiencies in the above theorem can be deduced by the standard bubble function techniques in [32], except the terms (4), (4), and (4). First, we will prove the efficiency (4). We make use of the relation between the Galerkin functional and the quantity of interest which here is the boundary stress. It directly follows from the definition of the Galerkin functional (4.11), (2.3f) and (4.1e) that

$\displaystyle\left\\|\mathbf{G}_{h}\right\\|_{-1,\omega_{p}}$	$\displaystyle\lesssim\left\\|\nabla(\mathbf{Ry}-\mathbf{y}_{h})\right\\|_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}$
	$\displaystyle\lesssim\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}$
	$\displaystyle+\lVert r-r_{H}\rVert_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}.$	(4.37)

Let $\bar{p}\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}$ be an arbitrary but fixed node. In the following $s$ denotes a side which belongs to $\gamma_{\bar{p},C}$ . We take the corresponding side bubble function $\xi_{s}:=\prod_{p\in s}\phi_{p}.$ Test the function $\xi_{s}\mathbf{e}_{1}$ in the equation (4.23), we get

$\displaystyle\sum_{p\in\mathcal{V}_{h}^{C}\setminus\mathcal{V}^{fC}_{h}}\int_{\Gamma_{C}}\hat{\sigma}_{1}(\mathbf{y}_{h})\xi_{s}\phi_{p}=$	$\displaystyle-(\mathbf{G}_{h},\xi_{s}\mathbf{e}_{1})_{-1,1}+\sum_{p\in\mathcal{V}_{h}\setminus\mathcal{V}_{h}^{C}}(R_{h,1},\xi_{s}\phi_{p})_{-1,1}$
	$\displaystyle+\sum_{p\in\mathcal{V}_{h}^{C}}(\tilde{R}_{h,1},(\xi_{s}-m_{p}(\xi_{s}))\phi_{p})_{-1,1}+\sum_{p\in\mathcal{V}_{h}^{C}\setminus\mathcal{V}^{fC}_{h}}\int_{\Gamma_{C}}\hat{\sigma}_{1}(\mathbf{y}_{h})m_{p}(\xi_{s})\phi_{p}$
	$\displaystyle=-(\mathbf{G}_{h},\xi_{s}\mathbf{e}_{1})_{-1,1}+\sum_{p\in\mathcal{V}_{h}}\int_{\omega_{s}}\mathcal{R}_{1}(\mathbf{y}_{h})\xi_{s}\phi_{p}$
	$\displaystyle-\sum_{p\in\mathcal{V}_{h}^{sC}}s_{p,1}m_{p}(\xi_{s})\int_{\gamma_{p,C}}\phi_{p}-\sum_{p\in\mathcal{V}_{h}^{fC}}(\tilde{R}_{h,1},\phi_{p})_{-1,1}m_{p}(\xi_{s})$	(4.38)

If the side $s$ is not contained in any patch $\gamma_{p,C}$ of semi- or full-contact nodes $p$ , the two last terms are zero and we can proceed similar to the case of (4.32) and (4.33). Otherwise, in order to get rid of the last two terms, we replace $\xi_{s}$ by a suitable function $\theta_{s}$ such that $m_{p}(\theta_{s})=0$ for all semi- and full-contact nodes. The value of $m_{p}(.)$ for a semi-contact node $p$ depends on $\tilde{\gamma}_{p,C}$ which is a strict subset of $\gamma_{p,C}$ compare (4.15). If $\gamma_{p,C}$ consists of two intervals we choose the inner third of $\gamma_{p,C}$ containing $p$ as $\tilde{\gamma}_{p,C}.$

Figure 4.1. Subgrid of boundary patch

\gamma_{p,C}

For example in Fig. 4.1, the dark blue region is $\tilde{\gamma}_{p,C}$ for $p=p_{1}$ . A side $s$ has two nodes $p_{1},p_{2}$ . We denote the sides of the subgrid containing $p_{i}$ by $s_{i}$ and the middle part by $s_{M}$ , see Fig. 4.1. For the function $\theta_{s}$ we make the following ansatz

\displaystyle\theta_{s}=\sum_{i=1}^{2}a_{i}\xi_{i}+a_{M}\xi_{M},

(4.39)

where $\xi_{i}$ and $\xi_{M}$ are side bubble functions to $s_{i}$ and $s_{M}$ . The coefficients $a_{1},a_{2},a_{M}$ are determined so that $\int_{s}1=\sum_{p\in\mathcal{V}_{h}^{C}\setminus\mathcal{V}^{fC}_{h}}\int_{s}\theta_{s}\phi_{p},$ and $\int_{s_{i}}\theta_{s}\phi_{p_{i}}=0$ for $p_{i}$ full-contact or semi-contact nodes. As $\bar{p}$ is not a full-contact node there is at least one contribution in the right hand side of the first condition. Inserting the ansatz (4.39) in the aforementioned conditions, we get a solvable system of three equations with three coefficients (degrees of freedom) $a_{1},a_{2},a_{M}$ . At this point the special choice of $m_{p}(\phi)=\frac{\int_{\tilde{\gamma}_{p,C}}\phi\phi_{p}}{\int_{\tilde{\gamma}_{p,C}}\phi_{p}}$ as mean value on $\tilde{\gamma}_{p,C}$ for semi-contact nodes becomes important because the choice $m_{p}(\phi)=\frac{\int_{\gamma_{p,C}}\phi\phi_{p}}{\int_{\gamma_{p,C}}\phi_{p}}$ as mean value over the whole patch $\gamma_{p,C}$ would lead to a contradiction of the conditions. In the second condition $s_{i}$ would be replaced by $s$ and the condition $\int_{s}\theta_{s}\phi_{p_{i}}=0$ for all $p_{i}$ of the side $s$ would imply $\sum_{p\in\mathcal{V}_{h}^{C}\setminus\mathcal{V}^{fC}_{h}}\int_{s}\theta_{s}\phi_{p}=0$ such that the first condition could not be fulfilled. As we assumed that the mesh is made of simplices, $\hat{\sigma}_{1}(\mathbf{y}_{h})$ is constant on $s$ . Consequently, $m_{p}(\theta_{s})=0$ implies $m_{p}(\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s})=0$ and it follows from the first condition

\displaystyle\left\|\hat{\sigma}_{1}(\mathbf{y}_{h})\right\|_{0,s}^{2}=\sum_{p\in\mathcal{V}_{h}^{C}\setminus\mathcal{V}^{fC}_{h}}\int_{s}\hat{\sigma}_{1}(\mathbf{y}_{h})\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s}\phi_{p}.

(4.40)

Putting together (4.40), (4.38) with test function $\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s}$ instead of $\xi_{s}$ and exploiting the conditions $m_{p}(\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s})=0$ for all contact nodes we end up with

$\displaystyle\left\\|\hat{\sigma}_{1}(\mathbf{y}_{h})\right\\|_{0,s}^{2}$	$\displaystyle=\sum_{p\in\mathcal{V}_{h}^{C}\setminus\mathcal{V}^{fC}_{h}}\int_{s}\hat{\sigma}_{1}(\mathbf{y}_{h})\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s}\phi_{p}$
	$\displaystyle=-(\mathbf{G}_{h},\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s}\mathbf{e}_{1})_{-1,1}+\int_{\omega_{s}}\mathcal{R}_{1}(\mathbf{y}_{h})\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s}$
	$\displaystyle\lesssim\left\\|\mathbf{G}_{h}\right\\|_{-1,\omega_{\bar{p}}}\left\\|\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s}\right\\|_{1,\omega_{s}}+\left\\|\mathcal{R}_{1}(\mathbf{y}_{h})\right\\|_{0,\omega_{s}}\left\\|\hat{\sigma}_{1}(\mathbf{y}_{h})\theta_{s}\right\\|_{1,\omega_{s}}$
	$\displaystyle\lesssim\left\\|\mathbf{G}_{h}\right\\|_{-1,\omega_{\bar{p}}}h_{s}^{-\frac{1}{2}}\left\\|\hat{\sigma}_{1}(\mathbf{y}_{h})\right\\|_{0,s}+h_{s}^{\frac{1}{2}}\left\\|\mathcal{R}_{1}(\mathbf{y}_{h})\right\\|_{0,\omega_{s}}\left\\|\hat{\sigma}_{1}(\mathbf{y}_{h})\right\\|_{0,s}$	(4.41)

where $h_{s}:=diam(s).$ In the last line of (4.41) we used the properties of the bubble functions on the subgrid and the fact that $\tilde{\gamma}_{p,C}$ is a fixed portion of $\gamma_{p,C}$ so that $h_{s}=ch_{s_{i}}$ for a mesh-independent constant $c$ . We divide by $h_{s}^{-\frac{1}{2}}\left\|\hat{\sigma}_{1}(\mathbf{y}_{h})\right\|_{0,s}$ leading to

\displaystyle h_{s}^{\frac{1}{2}}\left\|\hat{\sigma}_{1}(\mathbf{y}_{h})\right\|_{0,s}\lesssim\left\|\mathbf{G}_{h}\right\|_{-1,\omega_{\bar{p}}}+h_{s}\left\|\bm{\mathcal{R}}(\mathbf{y}_{h})\right\|_{0,\omega_{s}}.

(4.42)

By means of the triangle inequality, the shape-regularity, $h_{s}\approx h_{p}$ and the upper bounds (4.37) and (4.32) of $\left\|\mathbf{G}_{h}\right\|_{-1,\omega_{\bar{p}}}$ and $\left\|\bm{\mathcal{R}}(\mathbf{y}_{h})\right\|_{0,\omega_{s}},$ we get

	$\displaystyle h_{s}^{\frac{1}{2}}\left\\|\hat{\sigma}_{1}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}$	$\displaystyle\lesssim\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}+\lVert r-r_{H}\rVert_{0,\omega_{p}}$
		$\displaystyle+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}+osc(\mathbf{u}_{d},\omega_{p}).$		(4.43)

Similarly, one can derive

	$\displaystyle h_{s}^{\frac{1}{2}}\left\\|\hat{\sigma}_{2}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}$	$\displaystyle\lesssim\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}+\lVert r-r_{H}\rVert_{0,\omega_{p}}$
		$\displaystyle+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}+osc(\mathbf{u}_{d},\omega_{p})$		(4.44)

Adding (4) and (4), we get the desired result (4).

Now we turn back to the terms (4) and (4). The proof of both (4) and (4) proceeds similarly so we will give a sketch of the proof of (4) here and details can be found in[27, Sec. 5.2]. To prove (4), we derive a lower bound of the local error in terms of the local contributions of $s_{p,1}d^{b}_{p,1}.$ If $s_{p,1}=0$ or $(y^{1}_{b}-y_{h,1})(q)=0$ for all neighbouring nodes of $p$ we have $s_{p,1}d^{b}_{p,1}=0$ . Therefore, we assume $s_{p,1}>0$ and $(y^{1}_{b}-y_{h,1})(q)>0$ for at least one node on $\gamma_{p,C}$ . Let $\hat{q}$ be a node which fulfills $(y^{1}_{b}-y_{h,1})(\hat{q})\geq(y^{1}_{b}-y_{h,1})(q)$ for all neighboring nodes $q$ of $p$ . Due to $s_{p,1}>0$ we have $(y^{1}_{b}-y_{h,1})(p)=0.$ As we consider boundary meshes of triangles and intervals the discrete functions are piecewise linear. Using Taylor series expansion of $(y^{1}_{b}-y_{h,1})$ about $p,$ we get

\displaystyle(y^{1}_{b}-y_{h,1})(\hat{q})=\nabla|_{\hat{e}}(y^{1}_{b}-y_{h,1})\cdot(\hat{q}-p)\lesssim h_{p}\nabla|_{\hat{e}}(y^{1}_{b}-y_{h,1})\cdot\tau

(4.45)

where $\hat{e}\subset\gamma_{p,C}$ is an edge containing the nodes $\hat{q}$ and $p$ and $\tau$ is the unit tangential vector pointing from $p$ to $\hat{q}.$ The following estimate of (4.45) follows from [27, Sec. 5.2]:

\displaystyle(y^{1}_{b}-y_{h,1})(\hat{q})

\displaystyle\lesssim h_{p}^{\frac{1}{2}}\lVert\bm{\mathcal{J}}^{I}(\mathbf{y}_{h})\rVert_{0,\gamma_{p,I}}.

Since, $(y^{1}_{b}-y_{h,1})(q)\leq(y^{1}_{b}-y_{h,1})(\hat{q})$ for all $q\in\gamma_{p,C}$ , we can conclude that

\displaystyle d^{b}_{p,1}=\int_{\tilde{\gamma}_{p,C}}(y^{1}_{b}-y_{h,1})\phi_{p}\lesssim h^{\frac{3}{2}}_{p}\lVert\bm{\mathcal{J}}^{I}(\mathbf{y}_{h})\rVert_{0,\gamma_{p,I}}.

(4.46)

Now,

	$\displaystyle s_{p,1}d^{b}_{p,1}$	$\displaystyle=\frac{\int_{\omega_{p}}\mathcal{R}_{1}(\mathbf{y}_{h})\phi_{p}+\int_{\gamma_{p,I}}\mathcal{J}^{I}_{1}(\mathbf{y}_{h})\phi_{p}-\int_{\gamma_{p,C}}\hat{\sigma}_{1}(\mathbf{y}_{h})\phi_{p}}{\int_{\gamma_{p,C}}\phi_{p}}\cdot\int_{\tilde{\gamma}_{p,C}}(y^{1}_{b}-y_{h,1})\phi_{p}$
		$\displaystyle\lesssim\big{(}h_{p}\left\\|\mathcal{R}_{1}(\mathbf{y}_{h})\right\\|_{0,\omega_{p}}+h_{p}^{\frac{1}{2}}\left\\|\mathcal{J}^{I}_{1}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,I}}+h_{p}^{\frac{1}{2}}\left\\|\hat{\sigma}_{1}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}\big{)}\cdot h^{-1}_{p}\cdot\int_{\tilde{\gamma}_{p,C}}(y^{1}_{b}-y_{h,1})\phi_{p}$
		$\displaystyle\lesssim\big{(}h_{p}\left\\|\bm{\mathcal{R}}(\mathbf{y}_{h})\right\\|_{0,\omega_{p}}+h_{p}^{\frac{1}{2}}\left\\|\bm{\mathcal{J}}^{I}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,I}}+h_{p}^{\frac{1}{2}}\left\\|\hat{\bm{\sigma}}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}\big{)}\cdot h^{\frac{1}{2}}_{p}\lVert\bm{\mathcal{J}}^{I}(\mathbf{y}_{h})\rVert_{0,\gamma_{p,I}}$

Applying Cauchy-Schwarz inequality, (4), (4.32) and (4.33) we obtain

	$\displaystyle s_{p,1}d^{b}_{p,1}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}+\lVert r-r_{H}\rVert_{0,\omega_{p}}$
		$\displaystyle+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}+osc(\mathbf{u}_{d},\omega_{p})\big{)}^{2}.$

Thus, one can obtain the upper bound of $s_{p,2}d^{b}_{p,2}.$ Hence, for $p\in\mathcal{V}^{sC}_{h,b}$ we have

	$\displaystyle\sum_{i=1}^{2}s_{p,i}d^{b}_{p,i}\lesssim$	$\displaystyle\big{(}\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\omega_{p}}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\omega_{p}}+\lVert r-r_{H}\rVert_{0,\omega_{p}}$
		$\displaystyle+\rho\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\omega_{p}}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\omega_{p}}+osc(\mathbf{u}_{d},\omega_{p})\big{)}^{2}.$

Similarly, we can prove (4). ∎

5. Numerical Experiments

The aim of the given section is to numerically illustrate the theoretical results derived in Sections 3 and 4, respectively. We conduct two experiments with two model problems, one is a smooth solution on square mesh the other is a non-smooth solution on a non-convex domain. We construct the model problems with known solutions. The numerical experiments are performed on two model problems using MATLAB(version R2021a) software. For the computational simplicity, we slightly modify the cost functional $J$ , denoted by $\tilde{J}$ , by

\displaystyle\tilde{J}(\mathbf{w},\mathbf{x})=\frac{1}{2}\|\mathbf{w}-\mathbf{u}_{d}\|^{2}+\frac{\rho}{2}\|\nabla(\mathbf{x}-\mathbf{y}_{d})\|^{2}

(5.1)

subject to PDE,

\begin{split}-\Delta\mathbf{w}+\nabla{p}&=\mathbf{f}\quad\text{in}\;\Omega,\\ \nabla\cdot{\mathbf{w}}&=0\quad\text{in}\;\Omega,\\ \bf{w}&=\mathbf{x}\quad\text{on}\;\Gamma_{C},\\ \bf{w}&=\mathbf{0}\quad\text{on}\;\Gamma_{D},\end{split}

(5.2)

the set of controls is given by,

\mathbf{Q}_{ad}:=\{\mathbf{x}\in\mathbf{H}_{D}^{1}(\Omega):\mathbf{y}_{a}\leq\bm{\gamma}_{0}(\mathbf{x})\leq\mathbf{y}_{b}\text{ a.e. on }\Gamma_{C}\},

where the space $\mathbf{H}^{1}_{D}(\Omega)$ consists of $\mathbf{H}^{1}(\Omega)$ functions with vanishing trace on $\Gamma_{D}.$ The function $\mathbf{y}_{d}$ is given, and $\partial\Omega=\Gamma_{C}\cup\bar{\Gamma}_{D}$ . Then the minimization problem reads: Find $(\mathbf{u},\mathbf{y})\in\mathbf{H}_{D}^{1}(\Omega)\times\mathbf{Q}_{ad}$ satisfies (5.2) such that

\displaystyle\tilde{J}(\mathbf{u},\mathbf{y})=\min_{(\mathbf{v},\mathbf{x})\in\mathbf{H}_{D}^{1}(\Omega)\times\mathbf{Q}_{ad}}\tilde{J}(\mathbf{v},\mathbf{x}).

The corresponding discrete optimality system is given by

\begin{split}\mathbf{u}_{h}&=\mathbf{w}_{h}+\mathbf{y}_{h},\quad\mathbf{w}_{h}\in\mathbf{V}_{h},\\ a(\mathbf{w}_{h},\mathbf{z}_{h})+b(\mathbf{z}_{h},p_{H})&=(\mathbf{f},\mathbf{z}_{h})-a(\mathbf{y}_{h},\mathbf{z}_{h})\;\;\;{\rm for~{}all}\;\mathbf{z}_{h}\in\mathbf{V}_{h},\\ b(\mathbf{u}_{h},q_{H})&=0\;\quad{\rm for~{}all}\;q_{H}\in M_{H},\\ a(\mathbf{z}_{h},\bm{\phi}_{h})-b(\mathbf{z}_{h},r_{H})&=(\mathbf{u}_{h}-\mathbf{u_{d}},\mathbf{z}_{h})\;\;\;{\rm for~{}all}\;\mathbf{z}_{h}\in\mathbf{V}_{h},\\ b(\bm{\phi}_{h},q_{H})&=0\;\quad{\rm for~{}all}\;q_{H}\in M_{H},\\ \rho\,a(\mathbf{y}_{h},\mathbf{x}_{h}-\mathbf{y}_{h})\geq&a(\mathbf{x}_{h}-\mathbf{y}_{h},\bm{\phi}_{h})-b(\mathbf{x}_{h}-\mathbf{y}_{h},r_{H})\\ &-(\mathbf{u}_{h}-\mathbf{u_{d}},\mathbf{x}_{h}-\mathbf{y}_{h})+\rho a(\mathbf{y}_{d},\mathbf{x}_{h}-\mathbf{y}_{h})\;\quad{\rm for~{}all}\;\mathbf{x}_{h}\in\mathbf{Q}^{h}_{ad},\end{split}

(5.3)

where, $\mathbf{V}_{h}:=\{\mathbf{v}_{h}\in\mathbf{H}_{0}^{1}(\Omega):\mathbf{v}_{h}|_{T}\in\mathbf{P}_{1}(T)\;\;\forall T\in\mathcal{T}_{h}\}$ and $\mathbf{Q}^{h}_{ad}=\mathbf{Q}_{h}\cap\mathbf{Q}_{ad}.$ The set $\mathbf{Q}_{h}$ is defined by $\mathbf{Q}_{h}=\{\mathbf{x}_{h}\in\mathbf{H}^{1}_{D}(\Omega):\mathbf{x}_{h}|_{T}\in\mathbf{P}_{1}(T),\;\;\forall T\in\mathcal{T}_{h}\}.$ The discrete control space $M_{H}=\{p_{H}\in L^{2}(\Omega):p_{H}|_{T}\in P_{0}(T),\;\forall T\in\mathcal{T}_{H}\}$ . We solve the above discrete optimality (5.3) system using primal-dual active set strategy [23, Section 2].

To illustrate the primal-dual active set strategy algorithm, let us define some notations. Let the dimension of $\mathbf{V}_{h}$ and $\mathbf{Q}_{h}$ be denoted by $2n$ and $2m$ , respectively. Also, let the dimension of discrete pressure space $M_{H}$ to be $\kappa.$ Let $\mathcal{V}_{h}^{D}$ and $\mathcal{V}_{h}^{C}$ denote the set of vertices on $\overline{\Gamma}_{D}$ and the set of vertices interior to $\Gamma_{C},$ in the fine mesh $\mathcal{T}_{h}$ , respectively. The active and inactive sets corresponding to the bilateral constraints are

$\displaystyle\mathcal{A}^{k}_{b}$	$\displaystyle:=$	$\displaystyle\{i\in\mathcal{V}_{h}^{C}\;\|\;\bm{\mu}_{i}^{k}+(\mathbf{y}_{h}^{k}-\mathbf{y}_{b})_{i}>0\},$
$\displaystyle\mathcal{A}^{k}_{a}$	$\displaystyle:=$	$\displaystyle\{i\in\mathcal{V}_{h}^{C}\;\|\;\bm{\mu}_{i}^{k}+(\mathbf{y}_{h}^{k}-\mathbf{y}_{a})_{i}<0\},$
$\displaystyle\mathcal{I}^{k}$	$\displaystyle:=$	$\displaystyle\{i\in\mathcal{V}_{h}^{C}\;\|\;\bm{\mu}_{i}^{k}+(\mathbf{y}_{h}^{k}-\mathbf{y}_{b})_{i}\leq 0\leq\bm{\mu}_{i}^{k}+(\mathbf{y}_{h}^{k}-\mathbf{y}_{a})_{i}\},$

where $\mathbf{y}_{h}^{k}$ is the $k^{th}$ iterate of $\mathbf{y}_{h},$ and $\bm{\mu}\in\mathbb{R}^{2m}$ is the Lagrange multiplier. We denote by $I$ , an identity matrix of size $2m\times 2m$ . Moreover, $\mathcal{V}^{I}_{h}$ denotes the set of all interior nodes in $\mathcal{T}_{h}$ . Let $\bm{\phi}_{i}$ are basis functions of $V_{h}$ and $Q_{h}.$ Also, let $\chi_{i}$ are basis function of $M_{H}.$ We define the following matrices and vectors:

	$\displaystyle A_{i,j}$	$\displaystyle:=$	$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla\bm{\phi}_{i}:\nabla\bm{\phi}_{j},\quad M_{i,j}:=\sum_{K\in\mathcal{T}_{h}}\int_{K}\bm{\phi}_{i}\cdot\bm{\phi}_{j},\quad B_{j,k}:=\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla\cdot\bm{\phi}_{j}\;\chi_{k},$
	$\displaystyle F^{1}_{j}$	$\displaystyle:=$	$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}\mathbf{f}\cdot\bm{\phi}_{j},\quad F^{2}_{j}:=\sum_{K\in\mathcal{T}_{h}}\int_{K}\mathbf{u}_{d}\cdot\bm{\phi}_{j},\quad L=\rho A+M.$

Now the primal-dual active set algorithm for the Dirichlet boundary control problem (5.1) reads as:

Step 1.:

Initialize $\mathbf{y}^{0}_{h}$ $\bm{\mu}^{0}$ and set $k=0$ .

Step 2.:

Set the active and inactive sets ( $\mathcal{A}^{k}_{a},\mathcal{A}^{k}_{b},\mathcal{I}^{k}$ ).

Step 3.:

Solve

$\displaystyle A_{2n\times 2n}\left[\mathbf{w}_{h}^{k+1}\right]_{2n\times 1}+A_{2n\times 2m}\left[\mathbf{y}_{h}^{k+1}\right]_{2m\times 1}+B_{2n\times\kappa}[p^{k+1}_{H}]_{\kappa\times 1}$	$\displaystyle=$	$\displaystyle[F^{1}]_{2n\times 1},$
$\displaystyle B^{T}_{\kappa\times 2n}[\mathbf{w}_{h}^{k+1}]_{2n\times 1}+B^{T}_{\kappa\times 2m}[\mathbf{y}_{h}^{k+1}]_{2m\times 1}$	$\displaystyle=$	$\displaystyle[\mathbf{0}]_{\kappa\times 1}$
$\displaystyle A_{2n\times 2n}\left[\bm{\phi}_{h}^{k+1}\right]_{2n\times 1}-B_{2n\times\kappa}[r^{k+1}_{H}]_{\kappa\times 1}-M_{2n\times 2n}\left[\mathbf{w}_{h}^{k+1}\right]_{2n\times 1}-M_{2n\times 2m}\left[\mathbf{y}_{h}^{k+1}\right]_{2m\times 1}$	$\displaystyle=$	$\displaystyle-[F^{2}]_{2n\times 1},$
$\displaystyle B^{T}_{\kappa\times 2n}[\bm{\phi}^{k+1}]_{2n\times 1}$	$\displaystyle=$	$\displaystyle[\mathbf{0}]_{\kappa\times 1}$
$\displaystyle L_{2m\times 2m}\left[\mathbf{y}_{h}^{k+1}\right]_{2m\times 1}-A_{2m\times 2n}\left[\bm{\phi}_{h}^{k+1}\right]_{2n\times 1}+B_{2m\times\kappa}[r^{k+1}_{H}]_{\kappa\times 1}+M_{2m\times 2n}\left[\mathbf{w}_{h}^{k+1}\right]_{2n\times 1}+$
$\displaystyle I_{2m\times 2m}\left[\bm{\mu}^{k+1}\right]_{2m\times 1}=[F^{2}]_{2m\times 1}$	$\displaystyle+$	$\displaystyle\rho A_{2m\times 2m}\left[q_{d}\right]_{2m\times 1},$
$\displaystyle\left[\mathbf{y}_{h}^{k+1}\right]_{2m\times 1}$	$\displaystyle=$	$\displaystyle\mathbf{y}_{a}\;\;\mbox{on}\;\;\mathcal{A}^{k}_{a}$
$\displaystyle\left[\mathbf{y}_{h}^{k+1}\right]_{2m\times 1}$	$\displaystyle=$	$\displaystyle\mathbf{y}_{b}\;\;\mbox{on}\;\;\mathcal{A}^{k}_{b},$
$\displaystyle\left[\bm{\mu}^{k+1}\right]_{2m\times 1}$	$\displaystyle=$	$\displaystyle\mathbf{0}\;\;\mbox{on}\;\;\mathcal{I}^{k}\cup\mathcal{V}^{I}_{h}.$

Step 4.:

Stop using the criterion $\mathcal{A}^{k+1}=\mathcal{A}^{k}$ and $\mathcal{I}^{k+1}=\mathcal{I}^{k}$ or $\left\|\mathbf{y}_{h}^{k+1}-\mathbf{y}_{h}^{k}\right\|_{L^{\infty}(\Omega)}<\epsilon$ for $\epsilon>0$ , or set $k=k+1$ and return to Step 2.

In each of our model problems, we compute the error and estimator, which are defined as follows:

\displaystyle\textbf{Error}:=\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\Omega}+\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\Omega}+\lVert p-p_{H}\rVert_{0,\Omega}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\Omega}+\lVert r-r_{H}\rVert_{0,\Omega},

(5.5)

and,

\displaystyle\textbf{Estimator}:=\eta_{(\mathbf{u},p)}+\eta_{(\bm{\phi},r)}+\eta_{\mathbf{y}},

(5.6)

where $\eta_{(\mathbf{u},p)},$ $\eta_{(\bm{\phi},r)}$ and $\eta_{\mathbf{y}}$ are defined in Theorem 4.6.

For the adaptive algorithm, we use the following paradigm:

\displaystyle{\it Solve\rightarrow Estimate\rightarrow Mark\rightarrow Refine}.

First, we compute the discrete solutions ( $\mathbf{u}_{h},p_{H},\bm{\phi}_{h},r_{H},\mathbf{y}_{h}$ ) using the above-described primal-dual active set algorithm. Then in the second step using the discrete solution, we compute the error estimator (Estimator = $\eta_{(\mathbf{u},p)}+\eta_{(\bm{\phi},r)}+\eta_{\mathbf{y}}$ ) over each element. We use the Dörlfer marking technique [13] with bulk parameter $\theta=0.3$ for the mark step. Then the marked elements are refined using the newest vertex bisection algorithm [8] to obtain a new mesh and the algorithm is repeated. The convergence rate for Estimator is defined as follows:

\displaystyle\text{rate of convergence}(\ell):=\frac{\log(\textbf{Estimator}_{\ell+1}/\textbf{Estimator}_{\ell})}{\log(N_{\ell}/N_{\ell+1})},

for $\ell:=1,2,3,\cdots$ , where $\textbf{Estimator}_{\ell}$ and $N_{\ell}$ denotes the estimator and number of degrees of freedom at $\ell$ th level respectively. Similarly, one can define the rate of convergence for Error.

Example 5.1.

In this example we consider the optimal control problem (5.1) with the computational domain $\Omega=(0,1)^{2}$ and $\Gamma_{D}=(0,1)\times\{0\},$ $\Gamma_{C}=\partial\Omega\backslash\Gamma_{D}$ . We choose the constants $\rho=10^{-2},$ $\mathbf{y}_{a}=(-4,-2),$ and $\mathbf{y}_{b}=(2,2.5).$ The state and adjoint state variables are given by

\displaystyle{\bf u}={\bf y}=\left(\begin{array}[]{c}-\exp(x)(y\cos(y)+\sin(y))\\ \exp(x)y\sin(y)\end{array}\right),~{}~{}~{}p=\sin(2\pi x)\sin(2\pi y),

and

\displaystyle\bm{\phi}=\left(\begin{array}[]{c}(\sin(\pi x))^{2}\sin(\pi y)\cos(\pi y)\\ -(\sin(\pi y))^{2}\sin(\pi x)\cos(\pi x)\end{array}\right),~{}~{}~{}r=\sin(2\pi x)\sin(2\pi y).

We choose $\mathbf{u}$ and $\bm{\phi}$ such that $\nabla\cdot{\mathbf{u}}=\nabla\cdot{\bm{\phi}}=0\quad\text{in}\;\Omega$ and $\bm{\phi}=\mathbf{0}~{}~{}\text{on}~{}~{}\partial\Omega.$ The data of the problem are chosen such that $\mathbf{f}=-\Delta\mathbf{u}+\nabla{p},~{}\mathbf{u}_{d}=\mathbf{u}+\Delta\bm{\phi}+\nabla{r}~{}\text{and}~{}\mathbf{y}_{d}=\mathbf{y}$ .

We have used the above-described primal-dual active set algorithm to solve the optimal control problem. Figure 5.1(A) and Figure 5.1(B) show the coarse and refine meshes respectively. Figure 5.1(C) shows the convergence of error and estimator in terms of the number of degrees of freedom. It is clear from Figure 5.1(C), that the Estimator and the Error show the optimal rate of convergence. Here the optimal rate of convergence means the rate of convergence is $0.5$ with respect to the number of degrees of freedom( $N$ ).

Refer to caption — (a) Coarse mesh $\mathcal{T}_{H}$

Example 5.2.

In this example we consider the optimal control problem (5.1) with the L-shaped domain $\Omega=(-1,1)^{2}\setminus([0,1]\times[-1,0])$ , $\Gamma_{C}=\partial\Omega,$ and the exact solutions

	$\displaystyle{\bf u}={\bf y}=r^{\alpha}\left(\begin{array}[]{c}(1+\alpha)\sin(\theta)\omega(\theta)+\cos(\theta)\omega^{\prime}(\theta)\\ -(1+\alpha)\cos(\theta)\omega(\theta)+\sin(\theta)\omega^{\prime}(\theta)\end{array}\right),$
	$\displaystyle p=-r^{\alpha-1}((1+\alpha)^{2}\omega^{\prime}(\theta)+\omega{{}^{\prime\prime\prime}}(\theta))/(1-\alpha),$

where

	$\displaystyle\omega(\theta)=$	$\displaystyle 1/(1+\alpha)\sin(\alpha+1)\theta)\cos(\alpha w)-\cos((\alpha+1)\theta)$
		$\displaystyle+1/(1+\alpha)\sin(\alpha-1)\theta)\cos(\alpha\omega)-\cos((\alpha-1)\theta)$

and $\alpha=856399/1572864$ and $w=3\pi/2$ . The adjoint variables ${\bm{\phi}},~{}r$ are considered the same as in Example 5.1. The data of the problem is chosen such that $\mathbf{f}=-\Delta\mathbf{u}+\nabla{p},~{}\mathbf{u}_{d}=\mathbf{u}+\Delta\bm{\phi}+\nabla{r}~{}\text{and}~{}\mathbf{y}_{d}=\mathbf{y}$ . The constants $\rho=10^{-2},\;\mathbf{y}_{a}=(-3,-3),$ and $\mathbf{y}_{b}=(4,4).$

This problem is defined on the L-shaped domain, and the derivative of the solution $({\bf u},p)$ has a singularity at the origin. It is well known that for this problem the uniform refinements will not provide an optimal convergence rate. We have a similar observation from Figure 5.2, for uniform refinements convergence rate with respect to the number of degrees of freedom ( $N$ ) is $0.33$ Hence, one can use the adaptive algorithm to improve the convergence rate. Figure 5.3(A) and Figure 5.3(B) show the adaptive coarse and refined meshes respectively. Figure 5.3(C) shows the adaptive convergence of error and estimator in terms of the number of degrees of freedom( $N$ ). We see that the convergence rate has been improved from $0.33$ (Figure 5.2) to $0.50$ (Figure 5.3(C)). Thus the optimal convergence is achieved using the adaptive algorithm for the error in energy norm in the state and adjoint state velocity, control approximation, in $L2-$ norm of pressure and adjoint pressure variables. Hence, the optimal convergence for the a posteriori estimator(Estimator) and the total error(Error), which are defined in (5.5) and (5.6). Here, the optimal convergence means the rate of convergence is $0.5$ with respect to the number of degrees of freedom( $N$ ).

6. Conclusions

In this article, we propose, analyze, and test an a posteriori error estimator for the Dirichlet boundary control problem governed by Stokes equation. We develop an inf-sup stable finite element discretization scheme by using $\mathbf{P}_{1}$ elements(in the fine mesh) for the velocity and control variable and $P_{0}$ elements(in the coarse mesh) for the pressure variable. The optimal control satisfies a bilateral Signorini contact problem, thus the discrete optimality system consists of a discrete variational inequality for the approximate control variable. We derive and analyze the error estimator for the control variable and the estimator is designed for controlling its energy error. The estimator reduces to the standard residual estimator for elliptic problem, if no contact occurs. The contributions by the estimator addressing the nonlinearity are related to the contact stresses and the complementarity condition. We prove the reliability and efficiency of the estimator and ensure the equivalence with the error up to oscillation terms. Our numerical experiments confirm the theoretical results.

References

Apel et al. [2015] T. Apel, M. Mateos, J. Pfefferer, and A. Rösch. On the regularity of the solutions of Dirichlet optimal control problems in polygonal domains. SIAM J. Control Optim., 53(6):3620–3641, 2015.
Brenner and Scott [2008] S. C. Brenner and L. R. Scott. The mathematical theory of finite element methods. Texts in Applied Mathematics. Springer, New York, third edition, 2008.
Carstensen et al. [2012] C. Carstensen, M. Eigel, R. H. W. Hoppe, and C. Löbhard. A review of unified a posteriori finite element error control. Numer. Math. Theory Methods Appl., 5(4):509–558, 2012.
Casas and Dhamo [2012] E. Casas and V. Dhamo. Error estimates for the numerical approximation of Neumann control problems governed by a class of quasilinear elliptic equations. Comput. Optim. Appl., 52(3):719–756, 2012.
Casas and Mateos [2008] E. Casas and M. Mateos. Error estimates for the numerical approximation of Neumann control problems. Comput. Optim. Appl., 39(3):265–295, 2008.
Casas and Raymond [2006] E. Casas and J.-P. Raymond. Error estimates for the numerical approximation of Dirichlet boundary control for semilinear elliptic equations. SIAM J. Control Optim., 45(5):1586–1611, 2006.
Casas et al. [2009] E. Casas, M. Mateos, and J.-P. Raymond. Penalization of Dirichlet optimal control problems. ESAIM Control Optim. Calc. Var., 15(4):782–809, 2009.
Chen and Zhang [2010] L. Chen and C. Zhang. A coarsening algorithm on adaptive grids by newest vertex bisection and its applications. J. Comput. Math., 28(6):767–789, 2010.
Chowdhury et al. [2015] S. Chowdhury, T. Gudi, and A. K. Nandakumaran. A framework for the error analysis of discontinuous finite element methods for elliptic optimal control problems and applications to $C^{0}$ IP methods. Numer. Funct. Anal. Optim., 36(11):1388–1419, 2015.
Chowdhury et al. [2017] S. Chowdhury, T. Gudi, and A. K. Nandakumaran. Error bounds for a Dirichlet boundary control problem based on energy spaces. Math. Comp., 86(305):1103–1126, 2017.
Ciarlet [1978] P. G. Ciarlet. The finite element method for elliptic problems. North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978. Studies in Mathematics and its Applications, Vol. 4.
Dond et al. [2019] A. K. Dond, T. Gudi, and R. C. Sau. An error analysis of discontinuous finite element methods for the optimal control problems governed by Stokes equation. Numer. Funct. Anal. Optim., 40(4):421–460, 2019.
Dörfler [1996] W. Dörfler. A convergent adaptive algorithm for poisson’s equation. SIAM Journal on Numerical Analysis, 33(3):1106–1124, 1996.
Fierro and Veeser [2003] F. Fierro and A. Veeser. A posteriori error estimators for regularized total variation of characteristic functions. SIAM J. Numer. Anal., 41(6):2032–2055, 2003.
Fursikov et al. [1998] A. V. Fursikov, M. D. Gunzburger, and L. S. Hou. Optimal Dirichlet control and inhomogeneous boundary value problems for the unsteady Navier-Stokes equations. In Control and partial differential equations (Marseille-Luminy, 1997), volume 4 of ESAIM Proc., pages 97–116. Soc. Math. Appl. Indust., Paris, 1998.
Girault and Raviart [1979] V. Girault and P.-A. Raviart. Finite element approximation of the Navier-Stokes equations. Lecture Notes in Mathematics. Springer-Verlag, Berlin-New York, 1979.
Gong et al. [2020a] W. Gong, W. Hu, M. Mateos, J. R. Singler, and Y. Zhang. Analysis of a hybridizable discontinuous Galerkin scheme for the tangential control of the Stokes system. ESAIM Math. Model. Numer. Anal., 54(6):2229–2264, 2020a.
Gong et al. [2020b] W. Gong, M. Mateos, J. R. Singler, and Y. Zhang. Analysis and approximations of Dirichlet boundary control of Stokes flows in the energy space. https://arxiv.org/abs/2011.08221., 2020b.
Gudi and Sau [2020] T. Gudi and R. C. Sau. Finite element analysis of the constrained Dirichlet boundary control problem governed by the diffusion problem. ESAIM Control Optim. Calc. Var., 26:Paper No. 78, 19, 2020.
Gudi and Sau [2023] T. Gudi and R. C. Sau. A two level finite element method for Stokes constrained Dirichlet boundary control problem. Comput. Math. Appl., 129:126–135, 2023.
Günther and Hinze [2011] A. Günther and M. Hinze. Elliptic control problems with gradient constraints—variational discrete versus piecewise constant controls. Comput. Optim. Appl., 49(3):549–566, 2011.
Gunzburger et al. [1991] M. D. Gunzburger, L. S. Hou, and T. P. Svobodny. Analysis and finite element approximation of optimal control problems for the stationary Navier-Stokes equations with Dirichlet controls. RAIRO Modél. Math. Anal. Numér., 25(6):711–748, 1991.
Hintermüller et al. [2002] M. Hintermüller, K. Ito, and K. Kunisch. The primal-dual active set strategy as a semismooth newton method. SIAM Journal on Optimization, 13(3):865–888, 2002.
Hinze [2005] M. Hinze. A variational discretization concept in control constrained optimization: the linear-quadratic case. Comput. Optim. Appl., 30(1):45–61, 2005.
Houston et al. [2005] P. Houston, D. Schötzau, and T. P. Wihler. Energy norm a posteriori error estimation for mixed discontinuous Galerkin approximations of the Stokes problem. J. Sci. Comput., 22/23:347–370, 2005.
Karkulik [2020] M. Karkulik. A finite element method for elliptic Dirichlet boundary control problems. Comput. Methods Appl. Math., 20(4):827–843, 2020.
Krause et al. [2015] R. Krause, A. Veeser, and M. Walloth. An efficient and reliable residual-type a posteriori error estimator for the Signorini problem. Numer. Math., 130(1):151–197, 2015.
Meyer and Rösch [2004] C. Meyer and A. Rösch. Superconvergence properties of optimal control problems. SIAM J. Control Optim., 43(3):970–985, 2004.
Of et al. [2010] G. Of, T. X. Phan, and O. Steinbach. Boundary element methods for Dirichlet boundary control problems. Math. Methods Appl. Sci., 33(18):2187–2205, 2010.
Of et al. [2015] G. Of, T. X. Phan, and O. Steinbach. An energy space finite element approach for elliptic Dirichlet boundary control problems. Numer. Math., 129(4):723–748, 2015.
Veeser and Verfürth [2009] A. Veeser and R. Verfürth. Explicit upper bounds for dual norms of residuals. SIAM J. Numer. Anal., 47(3):2387–2405, 2009.
Verfürth [1995] R. Verfürth. A Review of A Posteriori Error Estimation and Adaptive Mesh Refinement Techniques. Wiley-Teubner, Chichester, 1995.
Winkler [2020] M. Winkler. Error estimates for variational normal derivatives and Dirichlet control problems with energy regularization. Numer. Math., 144(2):413–445, 2020.

$\displaystyle\left\\|p-R_{0}p\right\\|_{0,\Omega}$	$\displaystyle\leq\sup_{\mathbf{v}\in\mathbf{V}}\frac{b(\mathbf{v},p-R_{0}p)}{\left\\|v\right\\|_{\mathbf{V}}}$
	$\displaystyle\leq\sup_{\mathbf{v}\in\mathbf{V}}\frac{-a(\mathbf{y}-\mathbf{y_{h}},\mathbf{v})-a(\mathbf{w}-\mathbf{R}\mathbf{w},\mathbf{v})}{\left\\|v\right\\|_{\mathbf{V}}}$
	$\displaystyle\leq\left\\|\nabla(\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{w}-\mathbf{R}\mathbf{w})\right\\|_{0,\Omega}.$	(4.6)

$\displaystyle\left\\|\nabla(\bm{\phi}-\bar{\mathbf{R}}\bm{\phi})\right\\|_{0,\Omega}$	$\displaystyle\leq\left\\|\mathbf{u}-\mathbf{u}_{h}\right\\|_{0,\Omega}$
	$\displaystyle\leq\left\\|\mathbf{u}-\mathbf{Ru}\right\\|_{0,\Omega}+\left\\|\mathbf{Ru}-\mathbf{u}_{h}\right\\|_{0,\Omega}$
	$\displaystyle\leq\left\\|\mathbf{u}-\mathbf{Ru}\right\\|_{0,\Omega}+\left\\|\mathbf{Rw}-\mathbf{w}_{h}\right\\|_{0,\Omega}$
	$\displaystyle\leq\left\\|\mathbf{u}-\mathbf{Ru}\right\\|_{0,\Omega}+\left\\|\nabla(\mathbf{Rw}-\mathbf{w}_{h})\right\\|_{0,\Omega}.$	(4.10)

$\displaystyle\lVert\nabla(\mathbf{y}-\mathbf{y}_{h})\rVert_{0,\Omega}+$	$\displaystyle\lVert\nabla(\mathbf{u}-\mathbf{u}_{h})\rVert_{0,\Omega}+\lVert p-p_{H}\rVert_{0,\Omega}+\lVert\nabla(\bm{\phi}-\bm{\phi}_{h})\rVert_{0,\Omega}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\Omega}$
	$\displaystyle+\lVert r-r_{H}\rVert_{0,\Omega}\lesssim\left\\|\nabla(\mathbf{R}\mathbf{y}-\mathbf{y}_{h})\right\\|_{0,\Omega}+\left\\|\bm{\lambda}-\tilde{\bm{\lambda}}_{h}\right\\|_{-1,\Omega}+\left\\|\nabla(\mathbf{R}\mathbf{w}-\mathbf{w}_{h})\right\\|_{0,\Omega}$
	$\displaystyle+\left\\|p_{H}-R_{0}p\right\\|_{0,\Omega}+\left\\|\nabla(\bar{\mathbf{R}}\bm{\phi}-\bm{\phi}_{h})\right\\|_{0,\Omega}+\left\\|r_{H}-\bar{R}_{0}r\right\\|_{0,\Omega}.$	(4.16)

	$\displaystyle\left\\|\psi_{i}-m_{p}(\psi_{i})\right\\|_{0,\omega_{p}}$	$\displaystyle\leq Ch_{p}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}},$
	$\displaystyle\left\\|\psi_{i}-m_{p}(\psi_{i})\right\\|_{0,\gamma_{p}}$	$\displaystyle\leq Ch_{p}^{\frac{1}{2}}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}}$

$\displaystyle(\mathbf{G}_{h},\bm{\psi})_{-1,1}=$	$\displaystyle\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}h_{p}\left\\|\mathcal{R}_{i}(\mathbf{y}_{h})\right\\|_{0,\omega_{p}}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}}$
	$\displaystyle+\sum_{i=1}^{2}\sum_{p\in\mathcal{V}_{h}}h^{\frac{1}{2}}_{p}\left\\|\mathcal{J}^{I}_{i}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,I}}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}}$
	$\displaystyle-\sum_{i=1}^{2}\sum_{p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}}h_{p}^{\frac{1}{2}}\left\\|\hat{\sigma}_{i}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}\left\\|\nabla\psi_{i}\right\\|_{0,\omega_{p}}$
	$\displaystyle\lesssim\left(\sum_{p\in\mathcal{V}_{h}}\big{(}h_{p}^{2}\left\\|\mathcal{R}(\mathbf{y}_{h})\right\\|^{2}_{0,\omega_{p}}+h_{p}\left\\|\mathcal{J}^{I}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,I}}^{2}\big{)}\right.$
	$\displaystyle\left.+\sum_{p\in\mathcal{V}^{C}_{h}\setminus\mathcal{V}^{fC}_{h}}h_{p}\left\\|\hat{\bm{\sigma}}(\mathbf{y}_{h})\right\\|_{0,\gamma_{p,C}}^{2}\right)^{\frac{1}{2}}\left\\|\nabla\bm{\psi}\right\\|.$	(4.26)

A two-grid Adaptive Finite Element Method for the Dirichlet Boundary Control Problem Governed by Stokes Equation

Abstract.

Key words and phrases:

1991 Mathematics Subject Classification:

1. Introduction

2. Continuous Problem

2.1. Notation

2.2. Dirichlet Control Problem

Proposition 2.1.

Remark 2.2.

3. Discrete Problem

Proposition 3.1 (Discrete Optimality System).

4. A posteriori Error Analysis

Theorem 4.1 (Energy error estimate of control and L2L^{2}-estimate of velocity).

Proof.

Theorem 4.2 (Energy error estimate of velocity).

Proof.

Theorem 4.3 (Energy error estimate of adjoint velocity).

Proof.

Theorem 4.4 (Error estimate of pressure and adjoint pressure).

Theorem 4.5.

Theorem 4.6 (A posteriori error estimator).

Proof.

Remark 4.7.

Theorem 4.8 (Local Efficiency).

Proof.

5. Numerical Experiments

Example 5.1.

Example 5.2.

6. Conclusions

References

Theorem 4.1 (Energy error estimate of control and $L^{2}$ -estimate of velocity).