Duality-Based Error Control for the Signorini Problem
Abstract.
In this paper we study the a posteriori bounds for a conforming piecewise linear finite element approximation of the Signorini problem. We prove new rigorous a posteriori estimates of residual type in \(\operatorname{L}^{p}\), for \(p \in (4,\infty )\) in two spatial dimensions. This new analysis treats the positive and negative parts of the discretization error separately, requiring a novel sign- and bound-preserving interpolant, which is shown to have optimal approximation properties. The estimates rely on the sharp dual stability results on the problem in \(\operatorname{W}^{2,(4 - \varepsilon )/3}\) for any \(\varepsilon \ll 1\). We summarize extensive numerical experiments aimed at testing the robustness of the estimator to validate the theory.
Keywords
MSC codes
1. Introduction.
Consider the open and convex set \({\Omega } \subset \mathbb R^2\) with a polygonal boundary \(\partial{\Omega }\), and let \(\boldsymbol n\) denote the outward pointing normal to \(\partial{\Omega }\). This paper focuses on a variational inequality that emerges from the study of the following partial differential equation (PDE):accompanied by the boundary conditionsThe specified conditions dictate that, for any segment \(I \subset \partial{\Omega }\) where \(u\gt 0\), it follows that \(\nabla u \cdot{\boldsymbol n} = 0\) over \(I\). This problem, along with its variants, finds applications in diverse areas such as elasto-plasticity, fluid flow in porous media [14], and finance [46]. Its transformation into a variational inequality has facilitated extensive analysis, notably in [31], which established the existence and uniqueness of solutions across a broad class of problems.
\begin{align} - \Delta u+u=f \text{ in }{\Omega }, \end{align}
(1.1)
\begin{align} u \geq 0, \quad \nabla u \cdot{\boldsymbol n} \geq 0, \quad u \nabla u \cdot{\boldsymbol n} = 0 \text{ on } \partial{\Omega }. \end{align}
(1.2)
Significant progress has been made in deriving a priori error estimates for finite element approximations to variational inequalities, achieving optimal order in energy norms [12, 22, 33, 36], with notable advancements also reported in [26] and [34] for estimates in \(\operatorname{L}^{\infty }({\Omega })\). Despite these developments, results in norms other than the energy norm remain relatively rare. For instance, [32] offers an a priori bound in \(\operatorname{L}^{2}({\Omega })\), albeit in a single spatial dimension. The literature also documents bounds on trace norms along the Signorini boundary [38], with more recent studies, such as [16], extending bounds to several low order global norms.
The theory of a posteriori error estimation for elliptic variational inequalities is significantly less well developed than that for elliptic boundary value problems. Nonetheless, there have been several works on the Signorini problem as well as the related problem where the obstacle is interior to the domain rather than on the boundary. In [27], hierarchical estimates are used under the assumption that quadratic finite elements give better approximation than linears. Error estimates for variational inequalities with nontrivial boundary data are given in [30], estimates for an alternative form using Lagrange multipliers are given in [10, 45], with results for the discontinuous Galerkin method given in [44].
Rigorous a posteriori error bounds of residual type were derived for the variational inequality of the first kind for the Signorini problem in [15] for piecewise linear finite elements, where the authors give an upper bound in \(\operatorname{H}^{1}({\Omega })\). Several other works have considered a posteriori control in energy norms; see [25]. Pointwise error control was established in [35]. Recent regularity results for the problem in which constraints only apply at the boundary [16] show that better rates can be achieved in this case, paving the way for a posteriori error estimates in lower order norms.
We remark that error estimates have been derived for this problem using the dual-weighted residual methodology (see, for example, [41]) to estimate the error in a target quantity. In particular, estimates follow for the \(\operatorname{L}^{2}({\Omega })\)-error; however, the necessary stability of the dual problem enters as an assumption.
A priori and a posteriori error estimates in low order norms are established using duality arguments, and the resulting error bounds in, for example, \(\operatorname{L}^{2}({\Omega })\) are standard in the literature for boundary value problems. The Aubin–Nitsche duality argument allows a priori and a posteriori bounds in nonenergy norms; see [1] and the references therein. Application of this methodology has proved to be a challenge in the case of variational inequalities since the dual problem can lack the necessary regularity. Indeed, given smooth data, solutions of the Signorini problem are smooth away from the boundary but suffer from singular points at the boundaries of the contact set, that is, the set of points at which the boundary changes type. In this work, using the ideas of an a priori argument recently given in [16], the approach taken in [32] in which different dual problems are considered for the positive and negative parts of the error is adapted for the a posteriori setting and combined with a two-sided approximation that takes the place of the usual Lagrange interpolant. The novelty in this work is that we are able to prove rigorous a posteriori error control in \(\operatorname{L}^{p}({\Omega })\) for \(p \in (4,\infty )\). The results in [16] suggest that the exponent 4 is a critical value in the sense that piecewise linear finite element approximations achieve optimal a priori orders of convergence in \(\operatorname{W}^{1,4-\varepsilon }({\Omega })\) for any \(\varepsilon \ll 1\) but become increasingly suboptimal in \(\operatorname{W}^{1,p}({\Omega })\) for \(p\gt 4\). Increasingly suboptimal convergence is also observed in \(\operatorname{L}^{p}({\Omega })\) for \(p\gt 4\). We observe the same critical exponent for our a posteriori bounds; however, the assumptions we require to derive them are weaker than the a priori setting. Furthermore, we are able to regain optimal convergence rates through adaptive refinement techniques.
The main aim of this paper to understand the theoretical framework in which duality arguments can be applied to obtain a posteriori error bounds in nonenergy norms. To that end, we give a brief discussion of the underlying assumptions and those used elsewhere in the literature. We defer technical discussions to the main body of the article but give a brief overview here to put our results in context.
1.
A regularity assumption on the structure of the contact set of the solution to the Signorini problem (known as Condition \((\mathbf{A})\)) is required to prove regularity of the dual problem, and indeed such an assumption is commonplace in the literature [6, 7, 8, 16, 24, 28]. Condition \((\mathbf{A})\) excludes the possibility of a fractal contact set. Technical details will be given in Remark 2.4 (Condition \((\mathbf{A})\)), but we note that this assumption has been proven to hold in some cases, with the general case still an open problem. All original results in this article depend on Condition \((\mathbf{A})\), and it is sufficient to prove \(h\)-robust a posteriori error control in \(\operatorname{L}^{p}({\Omega })\) for \(p \in (4,\infty )\) on the positive part of the error (Theorem 3.12).
2.
We utilize a dual problem to derive a bound for the negative part of the error that depends upon the discrete solution. We make an assumption on the structure of its contact set, denoted Condition \((\mathbf a_h)\), which, to our knowledge, cannot be proven a priori, but can be verified a posteriori once the discrete solution is obtained. Using this, we have a posteriori error control on the negative part of the error for any discrete solution which satisfies Condition \((\mathbf a_h)\) (Theorem 3.17). The bound has an \(h\)-dependent regularity constant which we are able to compute a posteriori (see Remark 3.16), and in numerical experiments remains small.
3.
A strictly stronger assumption, known as Condition \((\mathbf{A}_h)\), which was used in [16] to obtain a priori error bounds in the \(\operatorname{L}^{4}({\Omega })\)-norm using duality methods, implies that the \(h\)-dependent constant mentioned above is uniformly bounded as \(h \to 0\). Under this stronger assumption, the bound for the negative part of the error would therefore be \(h\)-robust.
The rest of the paper is set out as follows: In section 2 we introduce notation, formalize the model problem, and describe the finite element approximation. In section 3 we state the key estimates and main results. In section 4 we discuss bound-preserving interpolation from above and below. We introduce a new interpolant which is shown to preserve bilateral bounds and to have optimal approximation properties. The performance of our error estimates is then tested numerically in section 5.
2. Problem setup.
This section contains a summary of the problem at hand, as well as a brief review of the current state of the art for regularity of the variational inequality and its finite element approximation. We note that the article [16] considers the problem on the unit square. In the current paper, to make use of their results, we will do the same, but we remark that arguments can be extended to convex polygonal domains.
In what follows, \({\Omega }\subseteq \mathbb R^2\) is assumed to be the unit square. Let \(\omega\) be a measurable subset of the domain \(\Omega\), and let \(\operatorname{L}^{p}(\omega )\) denote the Lebesgue spaces of \(p\)th power Lebesgue integrable functions over \(\omega\) with corresponding norm \({\left \|\cdot \right \|}_{\operatorname{L}^{p}(\omega )}\). We denote \({\left \langle \cdot,\cdot \right \rangle }_\omega\) as the \(\operatorname{L}^{2}\) inner product over \(\omega\) and use the convention of dropping the subscript if \(\omega ={\Omega }\). Finally, we let \(\operatorname{W}^{k,p}(\omega )\) be the Sobolev space whose first \(k\) weak derivatives are \(\operatorname{L}^{p}(\omega )\) and let \(\operatorname{H}^{k}(\omega ):=\operatorname{W}^{k,2}(\omega )\) with corresponding norms \({\left \|\cdot \right \|}_{\operatorname{W}^{k,p}(\omega )}\) and \({\left \|\cdot \right \|}_{\operatorname{H}^{k}(\omega )}\), respectively. For a real-valued function \(g\), we define \(g_+ = \max (g,0)\) and \(g_- = \min (g,0)\).
We will also henceforth consistently use the notation that \(p\) and \(q\) denote Hölder conjugates satisfyingwith \(p \in (4,\infty )\) and \(q \in (1, \frac 43)\).
\begin{align} \frac 1p + \frac 1q=1, \end{align}
(2.1)
We define the convex set of test and trial functionsthat is appropriate to define the weak form of (1.1)–(1.2). It reads as follows: let \(f \in \operatorname{L}^{s}({\Omega })\) for \(s\geq 2\), and find \(u \in{\mathcal K}\) such that
\begin{align}{\mathcal K} = \{v \in H^1({\Omega }) : v \geq 0 \text{ on } \partial{\Omega } \} \subset \operatorname{H}^{1}({\Omega }) \end{align}
(2.2)
\begin{align} a(u, v-u) :={\left \langle \nabla u,\nabla v - \nabla u\right \rangle } +{\left \langle u,v-u\right \rangle } \geq{\left \langle f,v - u\right \rangle } =: l(v-u) \quad \forall v \in{\mathcal K}. \end{align}
(2.3)
This leads us to state some well-known properties of the problem (1.1)–(1.2) as well as some new regularity results which require stronger assumptions on the problem data.
Given \(f \in \operatorname{L}^{2}({\Omega })\) there exists a unique solution \(u \in H^2({\Omega })\) to problem (2.3). Moreover,
\begin{align}{\left \|u\right \|}_{\operatorname{H}^{2}({\Omega })} \leq C{\left \|f\right \|}_{\operatorname{L}^{2}({\Omega })}. \end{align}
(2.4)
There are variants of the Signorini problem for which this regularity result cannot be applied. For example, if the boundary also contains portions where Neumann and/or Dirichlet conditions are applied, then the boundary conditions cannot be expressed in the correct form to apply Proposition 2.1. In this case we are not guaranteed \(u \in \operatorname{H}^{2}({\Omega })\). Such problems do arise in practical applications. Indeed they are studied in [9, 4].
Let \(u\) be the unique solution of problem (2.3). Since we have \(u \in H^2({\Omega })\) from Proposition 2.1, \(u\) has a continuous representative on \(\overline{\Omega }\) which we may use to define the contact setwith \(\mathring{\mathcal{A}}\) the relative interior of this set.
\begin{align*} \mathcal {A} := \{ x \in \partial \Omega \mid u({\boldsymbol x}) = 0\}, \end{align*}
Stronger regularity for the solution of (2.3) follows from an assumption referred to as Condition \((\mathbf{A})\) in [16], namely that the contact set, \(\mathcal{A}\), consists of finitely many connected components and isolated points, and in particular the set of critical points (i.e., points contained in the boundary of the contact set) does not contain accumulation points.
This means that the relative interior, \(\mathring{\mathcal{A}}\), is a union of open intervals in \(\partial{\Omega }\). This is a common assumption in works relating to problems with a Signorini boundary, such as [6, 7, 8, 16, 24, 28], with a related assumption on the free boundary in a problem with an interior obstacle given in [12, 13]. We also mention [20], where a priori bounds in \(\operatorname{H}^{1}({\Omega })\) are proved without this assumption; however, their methods are not applicable to the analysis in this work.
Optimal error estimates in certain nonenergy norms can be obtained for the Signorini problem if an appropriate dual problem can be defined which has sufficient regularity. Condition \((\mathbf{A})\) is a sufficient condition under which such a dual problem can be defined. It is not required for existence and uniqueness of the continuous problem, but is rather a regularity assumption. If Condition \((\mathbf{A})\) is not satisfied, optimal a priori estimates in the energy norm will still hold [20], but bounds in low order norms may be rendered suboptimal.
In order to violate this assumption, the solution to the Signorini problem would have to have infinitely many critical points. The boundary of the contact set would therefore contain accumulation points, meaning that approaching such a point, the solution along the boundary would switch from contact to noncontact infinitely many times and over arbitrarily small distances; see Figure 1.
We finally remark that this condition was recently shown to hold in [2] in polygonal domains in some cases, although the verification in the general setting remains a challenging open problem.
Suppose that \(u\) solves (2.3), and that Condition \((\mathbf{A})\) holds. If \(f \in \operatorname{L}^{s}({\Omega })\) for \(2 \lt s \lt 4\), then \(u \in \operatorname{W}^{2,s} ({\Omega })\). If \(f\in \operatorname{L}^{s}({\Omega })\) for \(4 \lt s \lt \infty\), then \(u\) admits a decomposition \(u = u_R + u_S\) with \(u_R \in \operatorname{W}^{2,s}({\Omega })\) and \(u_S \in \operatorname{H}^{5/2-\epsilon }({\Omega })\).
Finite element discretization. We consider \(\mathscr{T}{}\) to be a conforming triangulation of \(\Omega\); namely, \(\mathscr{T}{}\) is a finite family of sets such thatWe will make the assumption that the triangulation is shape-regular. Further, we define \(h:{\Omega }\to \mathbb R\) to be the piecewise constant meshsize function of \(\mathscr{T}{}\) given by
1.
\(K\in \mathscr{T}{}\) implies \(K\) is an open simplex or box,
2.
for any \(K,J\in \mathscr{T}{}\) we have that \(\overline{K}\cap \overline{J}\) is a full lower-dimensional simplex (i.e., it is either \(\emptyset\), a vertex, an edge, or the whole of \(\overline{K}\) and \(\overline{J}\)) of both \(\overline{K}\) and \(\overline{J}\), and
3.
\(\bigcup_{K\in \mathscr{T}{}}\overline{K}=\overline{{\Omega }}\).
\begin{align} h({\boldsymbol x}):=\max_{\overline{K}\ni{\boldsymbol x}}\text{diam} (K). \end{align}
(2.5)
We let \(\mathscr{E}{}\) be the skeleton (set of common interfaces) of the triangulation \(\mathscr{T}{}\) and say \(e\in \mathscr{E}\) if \(e\) is on the interior of \(\Omega\) and \(e\in \partial{\Omega }\) if \(e\) lies on the boundary \(\partial{\Omega }\) and set \(h_e\) to be the diameter of \(e\). We let \(\mathcal R^1(K)\) denote a space of piecewise linear (respectively, bilinear) polynomials over a triangle (respectively, quadrilateral), that is, \(\mathcal R^1(K) \in \{\mathbb P^1(K)\), \(\mathbb Q^1(K)\}\), and introduce the finite element spaceto be the usual space of continuous piecewise polynomial functions. We define the element patch of \(K \in \mathscr{T}{}\) asWe also define a discrete analogue of \(\mathcal K\) as
\begin{align} \mathbb V := \{\phi \in \operatorname{H}^{1}({\Omega }) : \phi |_K \in \mathcal R^1(K)\} \end{align}
(2.6)
\begin{align*} \hat K := \{J \in \mathscr {T}{} | \overline J \cap \overline K \neq \emptyset \}. \end{align*}
\begin{align} {\mathcal K}_h := \{v \in \mathbb V : v \geq 0 \text{ on } \partial{\Omega } \} ={\mathcal K} \cap \mathbb V \subset \operatorname{H}^{1}({\Omega }) . \end{align}
(2.7)
We are now in a position to state the finite element approximation of (1.1)–(1.2) which is to find \(U\in{\mathcal K}_h\) such thatNote that we restrict our attention to either piecewise linear finite element spaces on triangular meshes or piecewise bilinear spaces on quadrilateral meshes. For the sake of this work, we do not extend to higher degree polynomial approximation and consider only lowest order.
\begin{align} a(U, \Phi - U) \geq l(\Phi - U) \quad \forall \Phi \in{\mathcal K}_h. \end{align}
(2.8)
For \(f\in \operatorname{L}^{2}({\Omega })\) and all \(h \gt 0\) there exists a unique solution \(U\) of (2.8).
Under slightly stronger assumptions, namely that \(f\) is essentially bounded and a technical assumption on the discrete solution \(U\), which we detail below, one can prove a priori estimates in \(\operatorname{L}^{4}({\Omega })\) [16]. We note that this result is included for completeness, and that we do not make these additional assumptions.
The discrete contact set is defined analogously to the contact set. With \(U\) the solution to (2.8), we define
\begin{align*} \mathcal {A}_U := \{ x \in \partial \Omega \mid U({\boldsymbol x}) = 0\}. \end{align*}
We say that Condition \((\mathbf{A}_h)\) holds if there exist points \(d_i\) lying on \(\partial{\Omega }\) and numbers \(\delta_i \gt 0\) such that the intersections of the open balls \(B_{\delta_i}(d_i)\) and the boundary have nonzero distance from each other, and from the corners of the domain, such that the following holds for all sufficiently small \(h \gt 0\).
1.
The sets \(B_{\delta_i}(d_i) \cap \partial{\Omega }\) cover the boundary of the discrete contact set, and are such that each \(B_{\delta_i}(d_i) \cap \partial{\Omega }\) contains precisely one element of the boundary of the discrete contact set.
2.
Every connected component of \(\mathcal{A}_U\) has a nonempty interior.
We now summarize the range of a priori results and the main ideas used to derive them in [16] to illustrate the interplay of primal and dual regularity.
Let \(u\) solve (2.3) for some \(f \in \operatorname{L}^{p}({\Omega })\) with \(p\in (4, \infty )\). Let \(R : \operatorname{H}^{1}({\Omega }) \to \mathbb V\) denote the Ritz projector, and let \(U \in \mathbb V\) be the solution of (2.8). Then for all \(\varepsilon \in (0,1)\) there exists \(C\gt 0\) such that the superconvergence resultholds. The authors of [16] additionally assume that \(f \in \operatorname{L}^{\infty }({\Omega })\). With this assumption, using inverse estimates and choosing arbitrarily large \(p\), the result (2.9) provides optimal control of \(R(u) - U\) in the \(\operatorname{W}^{1,s}({\Omega })\)-norm for \(s\lt 4\), i.e.,Using an error bound in \(\operatorname{W}^{1,\infty }({\Omega })\), which the authors also prove, a bound in \(\operatorname{W}^{1,4}({\Omega })\) is obtained from (2.10) of optimal order up to arbitrarily small \(\varepsilon\), namelyThis result is combined with classical results on approximability of the Ritz projector to obtainThe a priori result (2.12) is a key ingredient in an Aubin–Nitsche duality argument which uses a Hölder inequality to bound in the strongest norm possible for an appropriately defined dual solution, which does not have \(\operatorname{H}^{2}({\Omega })\) regularity. The bound (2.12) is able to compensate for this and deliver optimal rates (up to an arbitrarily small positive number \(\varepsilon\)):If we additionally assume that Condition \((\mathbf{A}_h)\) holds, then we also haveand therefore achieve full a priori error control in \(\operatorname{L}^{4}({\Omega })\).
\begin{align} {\left \|R(u) - U\right \|}_{\operatorname{H}^{1}({\Omega })} \leq C h^{\tfrac 32 - \tfrac 2p -\varepsilon } \end{align}
(2.9)
\begin{align} {\left \|R(u) - U\right \|}_{\operatorname{W}^{1,s}({\Omega })} \leq C h. \end{align}
(2.10)
\begin{align}{\left \|R(u) - U\right \|}_{\operatorname{W}^{1,4}({\Omega })} \leq C h^{1-\varepsilon }. \end{align}
(2.11)
\begin{align} {\left \|u - U\right \|}_{\operatorname{W}^{1,4}({\Omega })} \leq C h^{1-\varepsilon }. \end{align}
(2.12)
\begin{align}{\left \|(u-U)_+\right \|}_{\operatorname{L}^{4}(\Omega )} \leq Ch^{2 - \varepsilon }. \end{align}
(2.13)
\begin{align}{\left \|(u-U)_-\right \|}_{\operatorname{L}^{4}(\Omega )} \leq Ch^{2 - \varepsilon }, \end{align}
(2.14)
Note that convergence rates in \(\operatorname{W}^{1,s}({\Omega })\) for \(s\gt 4\) are suboptimal, and that while the argument above can be modified to include stronger norms at a suboptimal rate, this suboptimality carries through to the \(\operatorname{L}^{p}({\Omega })\) estimate, i.e.,
\begin{align} {\left \|u-U\right \|}_{\operatorname{L}^{p}(\Omega )} \leq Ch^{\tfrac 32 + \tfrac 2p - \varepsilon }. \end{align}
(2.15)
3. A posteriori error control.
In this section we state and prove the main result of this work, a rigorous a posteriori bound in \(\operatorname{L}^{s}({\Omega })\) for \(s \gt 4\) for the problem (1.1). To begin we state some technical results.
Let \(v\in \operatorname{W}^{1,s}(K)\) with \(s\in (1,\infty )\). Then, there exists a constant \(C_{\textit{tr}}\) depending only on \(s\) and \(K\) such that
\begin{align}{\left \|v\right \|}_{\operatorname{L}^{s}(\partial K)} \leq C_{\textit{tr}}{\!\left ({h^{-\tfrac 1 s}{\left \|v\right \|}_{\operatorname{L}^{s}(K)} + h^{1-\tfrac 1 s}{\left \|\nabla v\right \|}_{\operatorname{L}^{s}(K)} }\right )}. \end{align}
(3.1)
We define jump operators as
\begin{align*}{\unicode{x27E6} v\unicode{x27E7}} ={v}|_{K_1}{\boldsymbol n}_{K_1} + {v}|_{K_2}{\boldsymbol n}_{K_2},\qquad {\unicode{x27E6}{\boldsymbol v}\unicode{x27E7}} ={{{\boldsymbol v}|_{K_1}}} \cdot{\boldsymbol n}_{K_1} +{{{\boldsymbol v}|_{K_2}}} \cdot{\boldsymbol n}_{K_2}.\end{align*}
Let \(\mathcal I\) denote the piecewise linear Lagrange interpolator over \(\mathscr{T}{}\). The following estimate holds for all elements \(K \in \mathscr{T}{}\), faces \(e\in \mathscr{E}\), and all functions \(v \in \operatorname{W}^{2,s}(K)\) with \(s\gt 1\):
\begin{align*} \max{\!\left ({ h^{-2}{\left \|v -{\mathcal I} v\right \|}_{\operatorname{L}^{s}(K)}, h^{-2+1/s}{\left \|v -{\mathcal I} v\right \|}_{\operatorname{L}^{s}(e)} }\right )} &\leq C{\left \|v\right \|}_{\operatorname{W}^{2,s}(\hat K)}. \end{align*}
3.1. Dual variational inequality.
Motivated by the approach used in [32] to derive a priori \(\operatorname{L}^{2}({\Omega })\) estimates, we define the following set:For \(p \in (4, \infty )\) consider the problem of finding \(z \in \mathcal{M}\) such that
\begin{align} \mathcal{M}= \{v \in \operatorname{H}^{1}({\Omega })\mid v \leq 0 \text{ on } \mathring{\mathcal{A}}\}. \end{align}
(3.2)
\begin{align} a(z, v-z) \geq{\left \langle \max (0,u-U)^{p-1},v-z\right \rangle } \quad \forall v \in \mathcal{M}. \end{align}
(3.3)
Observe that \(z + u - U\) lies in \(\mathcal{M}\), since by the definition of \(\mathring{\mathcal{A}}\), \(u\) vanishes there and \(U \geq 0\) on \(\partial{\Omega }\) by definition of \({\mathcal K}_h\), (2.7).
To derive a posteriori error bounds, we require results on the stability and regularity of the dual problem (3.3), which we will now quantify.
Before presenting the necessary results, we provide the following definitions regarding the corners of the domain and the boundary conditions.
We let \(\{b_1,\ldots,b_{N_p}\}\) be the set of critical points that make up the boundary of \(\mathring{\mathcal{A}}\) and the vertices of \(\partial{\Omega }\). Note that under the assumption that Condition \((\mathbf{A})\) holds, this is a finite set consisting of the vertices of \(\partial{\Omega }\) and the points at which there is a change in boundary condition. We let \(\omega_k\) denote the angle at the boundary at \(b_k\), that is, \(\omega_k = \pi\) at a change of boundary conditions or \(\omega_k = \tfrac{\pi }{2}\) at the corners of the domain, since the domain is the unit square. We will denote by \(r_k, \theta_k\) the polar coordinates centered at point \(b_k\), by \(\zeta_k\) a smooth cut-off function identically equal to 1 in a neighborhood of \(b_k\), and by \(\phi_k = \phi_k(\theta_k)\) a trigonometric function.
The functions described in Definition 3.5 are the key tools in quantifying the behavior of the solution around the points \(b_k\). The properties of the dual problem that we will need are summarized in the following lemma.
Let \(u\) solve (2.3) for some \(f \in \operatorname{L}^{2}({\Omega })\) and satisfy Condition \((\mathbf{A})\), and let \(U \in \mathbb V\) solve (2.8). Then the dual problem (3.3) is uniquely solvable in \(\operatorname{H}^{1}({\Omega })\). In addition, we have \(z \geq 0\) in \(\Omega\) and \(z=0\) on \(\mathcal{A}\), and for any \(p\in (4,\infty )\) and its conjugate \(q = \tfrac{p}{p-1}\) there exists \(C\gt 0\) independent of \(h\) such that the stability boundholds. Furthermore, we have the following regularity result in \(\operatorname{W}^{2,q}({\Omega })\):
\begin{align} {\left \|z\right \|}_{\operatorname{H}^{1}({\Omega })} \leq C{\left \|\max (0, u-U)^{p-1}\right \|}_{\operatorname{L}^{q}({\Omega })} \end{align}
(3.4)
\begin{align} {\left \|z\right \|}_{\operatorname{W}^{2,q}({\Omega })} \leq C{\left \|\max (0, u-U)^{p-1}\right \|}_{\operatorname{L}^{q}({\Omega })}. \end{align}
(3.5)
We first claim that \(a(z_+, z_-) = 0\). To see this, note that by definition we have \(z=z_+ + z_-\), andIt is clear that \(z_+ z_- \equiv 0\) since \(z\) cannot be both positive and negative. Further, we can write \({\Omega } = \{x \in{\Omega } : z \lt 0 \} \cup \{x \in{\Omega } : z \geq 0 \}\). Note that in each set, precisely at least one of \(z_+\) and \(z_-\) is identically the zero function, and therefore if the set has nonzero measure, its weak gradient must be zero there. Contributions to the integral (3.6) must therefore be zero and the claim is proved.
\begin{align} a(z_+, z_-) = \int_{{\Omega }}\nabla z_+ \cdot \nabla z_- + z_+ z_- \mathop{}\!\textrm{d} x. \end{align}
(3.6)
Since \(z \in \mathcal{M}\), \(z \leq 0\) on \(\mathring{\mathcal{A}}\), its positive part must be zero there, and we may take \(v = z_+\) as test function in (3.3), which yieldsIn other words, since \(\max (0,u-U)^{p-1}\) is nonnegative and \(z_-\) is nonpositive,so that \(z_- \equiv 0\), which proves that \(z\geq 0\). This in turn gives \(z=0\) on \(\mathring{\mathcal{A}}\).
\begin{align} -\|z_-\|^2_{\operatorname{H}^{1}({\Omega })} = a(z, -z_-) \geq{\left \langle \max (0,u-U)^{p-1},-z_-\right \rangle }. \end{align}
(3.7)
\begin{align} 0 \leq \|z_-\|^2_{\operatorname{H}^{1}({\Omega })} \leq{\left \langle \max (0,u-U)^{p-1},z_-\right \rangle } \leq 0, \end{align}
(3.8)
The bound (3.4) follows after noting that we can take \(v=0\) in (3.3) to seeand we can take \(v=2z\) in (3.3), yieldingTherefore,where we have used the embedding of \(\operatorname{H}^{1}({\Omega })\) in \(\operatorname{L}^{p}({\Omega })\) with embedding constant \(C_{\text{Sob}}\). The result now follows after dividing by \({\left \|z\right \|}_{\operatorname{H}^{1}({\Omega })}\).
\begin{align} a(z, -z) \geq{\left \langle \max (0,u-U)^{p-1},-z\right \rangle }, \end{align}
(3.9)
\begin{align} a(z, z) \geq{\left \langle \max (0,u-U)^{p-1},z\right \rangle }. \end{align}
(3.10)
\begin{align*}{\left \|z\right \|}^2_{\operatorname{H}^{1}({\Omega })} &={\left \langle \max (0,u-U)^{p-1},z\right \rangle } \\ &\leq{\left \|\max (0,u-U)^{p-1}\right \|}_{\operatorname{L}^{q}({\Omega })}{\left \|z\right \|}_{\operatorname{L}^{p}({\Omega })} \\ &\leq C_{\text{Sob}}{\left \|\max (0,u-U)^{p-1}\right \|}_{\operatorname{L}^{q}({\Omega })}{\left \|z\right \|}_{\operatorname{H}^{1}({\Omega })}, \end{align*}
Now, \(z\) is the weak solution of the boundary value problemwith zero Dirichlet boundary condition on \(\mathcal{A}\) and zero Neumann condition on the remainder of the boundary. The problem therefore fits into the framework of [23, Theorem 4.4.3.7], which allows us to quantify the regularity of \(z\) to the following extent. With the notation of Definition 3.5, there exist unique real numbers \(c_{k,m}\) such thatwhere \(q^{\prime }\) is the Hölder conjugate of \((4-\varepsilon )/ (1 - \varepsilon )\) and \(\lambda_{k,m}\) are eigenvalues of a Laplacian operator which depends upon \(\omega_k\) (see [23, section 4.4] for complete enumeration of the \(\lambda_{k, m}\) in all cases). In our case where \(\Omega\) is a convex polygonal domain, the regularity of \(z\) is determined by the term in (3.12) with the lowest power of \(r\), which is \(r^{1/ 2} \zeta_k \phi_k\). The singularity of type \(r^{1/ 2}\) occurs at points where the boundary condition changes (compare with [2, section 3] on this point).
\begin{align} - \Delta z+z = \max (0,u-U)^{p-1}, \end{align}
(3.11)
\begin{align} z - \sum_{\substack{1 \leq k \leq N_p \\ -2/ q^{\prime } \lt \lambda_{k,m} \lt 0 \\\lambda_{k,m} \neq -1}} c_{k,m} r^{-\lambda_{k,m}} \zeta_k \phi_k \in \operatorname{W}^{2,(4-\varepsilon )/ (1 - \varepsilon )}({\Omega }), \end{align}
(3.12)
Now, we observe that \(r^{1/ 2} \zeta_k \phi_k \in \operatorname{W}^{k,s}({\Omega })\) if and only if \(k - \tfrac 2 s \lt \tfrac{1}{2}\). This gives a limit on the regularity of the second derivatives of \(z\): \(z\in \operatorname{W}^{2,s}({\Omega })\) only if \(s \lt \tfrac 4 3\).
The proof of the stability bound in \(\operatorname{W}^{2,q}({\Omega })\) is a technical argument that we will not reproduce here, but it can be found in full in [16, Lemma 5.2]. □
It is important to note that the limit on regularity for the dual problem is imposed by the assumption of convexity of the domain and the nature of the boundary conditions, not the regularity of the problem data \(f\). Indeed, the nonsingular part of \(z\), that is, the expression in (3.12), possesses higher regularity depending on the problem data, but the regularity of the singular components of \(u\) is limited by the geometry. For certain boundary conditions including the Signorini condition (1.2), \(\operatorname{H}^{2}({\Omega })\)-regularity follows from Proposition 2.1 which implies that the singular part of the solution vanishes. For the mixed boundary condition, this result does not apply, and the regularity is therefore constrained by the least regular of the components in (3.12).
Proposition 3.6 also motivates a posteriori error estimation in \(\operatorname{L}^{p}({\Omega })\) for \(p \gt 4\). We see that we can only apply a duality approach for \(q \lt \tfrac 4 3\), with the Hölder conjugate of \(\tfrac 4 3\) being 4. Thus, \(\operatorname{L}^{4 + \varepsilon }({\Omega })\) for any \(\varepsilon \ll 1\) is the weakest space in which estimates of the traditional residual form can be established. In particular, error estimates in \(\operatorname{L}^{2}({\Omega })\) are not possible in this framework.
In light of the structure of the contact set, described in Remark 2.4, it can be seen (see [16, Thm. 2.3]) that \(u\) solves the boundary value problemHere, \(\Gamma_i \subset{\Omega }\) are disjoint open line segments such that \(\bigcup_{i=1}^{N+M} \bar{\Gamma }_i = \partial{\Omega }\), and such that there is a set \(R\) of measure zero such thatThe claim now follows immediately after choosing \(w\) as test function in (3.14) and integrating by parts. □
\begin{align} \begin{split} - \Delta u+u &= f \quad \text{in} \quad{\Omega },\\ u &= 0 \quad \text{on} \quad \Gamma_i,\ i=1,\ldots,N,\\ \nabla u \cdot{\boldsymbol n} &= 0 \quad \text{on} \quad \Gamma_i,\ i=N+1,\ldots,N+M. \end{split} \end{align}
(3.14)
\begin{align*} \mathcal {A} = \bigcup_{i=1}^N \bar {\Gamma }_i \cup R. \end{align*}
3.2. Interpolation operator for bilateral approximation.
We now discuss our requirements for a finite element interpolation operator, the construction of which is postponed to section 4.
We make the assumption that there exists an interpolation operator \({\hat \Pi } : \operatorname{W}^{2,q}({\Omega }) \to \mathbb V\) such that for any \(w\in \operatorname{W}^{2,q}({\Omega })\), \(q \in (1, \tfrac 43)\), with \(w \geq 0\),and for any \(K\in \mathscr{T}{}\)
\begin{align} 0 \leq{\hat \Pi }(w) \leq w, \end{align}
(3.15)
\begin{align}{\left \|w -{\hat \Pi }(w)\right \|}_{\operatorname{L}^{q}(K)} \leq C_b h^2{\left \|w\right \|}_{\operatorname{W}^{2,q}({\hat{{K}}})}. \end{align}
(3.16)
This interpolant is a crucial part of our a posteriori error analysis, as demonstrated by the following lemma, which allows us to introduce an interpolant into the error inequality for \((u-U)_+\). This is done by using the bilateral bounds (3.15) to determine the sign of \(a{ ({{\hat \Pi }(z), u- U} )}\).
By (3.13) and since \(0 \leq{\hat \Pi }(z) \leq z=0\) on \(\mathring{\mathcal{A}}\), we have immediately thatSince \({\hat \Pi }(z) \geq 0\), we may take \(\Phi = U +{\hat \Pi }(z)\) in (2.8), givingThe result then follows after combining (3.18) and (3.19). □
\begin{align} a{\left ({u, \,{\hat \Pi }(z)}\right )} ={\left \langle f,{\hat \Pi }(z)\right \rangle }. \end{align}
(3.18)
\begin{align} a{\left ({U, \,-{\hat \Pi }(z)}\right )} \leq -{\left \langle f,{\hat \Pi }(z)\right \rangle } . \end{align}
(3.19)
3.3. A posteriori error bounds.
We are now ready to state and prove the main result of this article, a rigorous a posteriori bound in \(\operatorname{L}^{p}({\Omega })\) for the problem (1.1)– (1.2). We prove this in two parts, bounding separately the quantities \({\left \|(u-U)_+\right \|}_{\operatorname{L}^{p}({\Omega })}\) and \({\left \|(u-U)_-\right \|}_{\operatorname{L}^{p}({\Omega })}\).
Let \(u\) solve problem (2.3) and \(U\) be the finite element approximation. Then, for \(p\in \left (4,\infty \right )\), \(q = \tfrac{p}{p-1}\) its conjugate, and \(f\in \operatorname{L}^{p}({\Omega })\), we havewhere
\begin{align} {\left \|(u-U)_+\right \|}_{\operatorname{L}^{p}({\Omega })}^p \leq \eta (U, f, h) := C \sum_{K\in \mathscr{T}{}} \left ( \eta_K^p + \frac 1 2 \eta_J^p \right ), \end{align}
(3.20)
\begin{align} \begin{split} \eta_K &= h^2{\left \|-\Delta U+U - f\right \|}_{\operatorname{L}^{p}(K)}, \\ \eta_J &= h^{2-1/q}{\left \|{\unicode{x27E6} \nabla U\unicode{x27E7}}\right \|}_{\operatorname{L}^{p}(\partial K)}. \end{split} \end{align}
(3.21)
We may select \(v = z+u - U\) in (3.3); this function is shown to belong to \(\mathcal{M}\) in Remark 3.4. We therefore haveWe may use Lemma 3.11 to introduce \({\hat \Pi }(z)\) as follows:Since \(z\) and \({\hat \Pi }(z)\) both have zero trace on \(\mathcal{A}\), we can use Proposition 3.9 to arrive atThe rest of the argument follows standard a posteriori techniques; to begin,Splitting the integral elementwise and making use of the Cauchy–Schwarz inequality, we see thatInvoking the optimal approximation result of Assumption 3.10 (proven in Lemma 4.4), we obtainSimilarly, we bound \(R_2\):To bound this term we use split \(z -{\hat \Pi }(z) = (z - \mathcal{I}z) + (\mathcal{I}z -{\hat \Pi }(z))\) and use the triangle inequality. We may then use approximation properties of the Lagrange interpolant on element boundaries, trace estimates (see Proposition 3.1), inverse estimates, and Assumption 3.10 to finally arrive atwhere \(C_{\text{inv}}\) is the constant from the inverse inequality. Collecting (3.26)–(3.28), we have
\begin{align} {\left \|(u - U)_+\right \|}^p_{\operatorname{L}^{p}({\Omega })} \leq a(z, u - U). \end{align}
(3.22)
\begin{align} {\left \|(u - U)_+\right \|}^p_{\operatorname{L}^{p}({\Omega })} \leq a{\left ({z-{\hat \Pi }(z), u - U}\right )} . \end{align}
(3.23)
\begin{align*} {\left \|(u - U)_+\right \|}^p_{\operatorname {L}^{p}({\Omega })} \leq l{\!\left ({z - {\hat \Pi }(z)}\right )} - a{\left ({z-{\hat \Pi }(z), U}\right )} . \end{align*}
\begin{align} \begin{split} l{\!\left ({z-{\hat \Pi }(z)}\right )} \!-\! a{\left ({U,z \!-\!{\hat \Pi }(z)}\right )} &= \int_{\Omega } f{\!\left ({z -{\hat \Pi }(z)}\right )} \!-\! \nabla U \cdot{\!\left ({\nabla z \!-\! \nabla{\hat \Pi }(z)}\right )} - U{\!\left ({z -{\hat \Pi }(z)}\right )} \mathop{}\!\textrm{d} x \\ &= \int_{\Omega }{\!\left ({f + \Delta U - U}\right )}{\!\left ({z -{\hat \Pi }(z)}\right )} \mathop{}\!\textrm{d} x \\ &\quad - \int_{\mathscr{E}\cup \partial{\Omega }}{\unicode{x27E6} \nabla U\unicode{x27E7}}{\!\left ({z -{\hat \Pi }(z)}\right )} \mathop{}\!\textrm{d} S \\ &=: R_1 + R_2. \end{split} \end{align}
(3.24)
\begin{align} \begin{split} R_1 &= \int_{\Omega }{\!\left ({f + \Delta U - U}\right )}{\!\left ({z -{\hat \Pi }(z)}\right )} \mathop{}\!\textrm{d} x \\ &\leq \sum_{K\in \mathscr{T}{}}{\left \|h^2{\!\left ({f + \Delta U - U}\right )}\right \|}_{\operatorname{L}^{p}(K)}{\left \|h^{-2}{\!\left ({z -{\hat \Pi }(z)}\right )}\right \|}_{\operatorname{L}^{q}(K)}. \end{split} \end{align}
(3.25)
\begin{align} R_1 \leq C_b \sum_{K\in \mathscr{T}{}}{\left \|h^2{\!\left ({f + \Delta U - U}\right )}\right \|}_{\operatorname{L}^{p}(K)} |z|_{\operatorname{W}^{2,q}(\hat K)}. \end{align}
(3.26)
\begin{align} \begin{split} R_2 &= -\sum_{e\in \mathscr{E} \cup \partial{\Omega }} \int_e{\unicode{x27E6} \nabla U\unicode{x27E7}}{\!\left ({z -{\hat \Pi }(z)}\right )} \\ &\leq \frac{1}{2} \sum_{K\in \mathscr{T}{}}{\!\left ({ \sum_{e\in \partial K}{\left \|h^{2-1/q}{\unicode{x27E6} \nabla U\unicode{x27E7}}\right \|}_{\operatorname{L}^{p}(e)}{\left \|h^{-2+1/q}{\!\left ({z -{\hat \Pi }(z)}\right )}\right \|}_{\operatorname{L}^{q}(e)} }\right )}. \end{split} \end{align}
(3.27)
\begin{align} R_2 \leq C_b C_{\text{tr}}(1 + C_{\text{inv}}) \frac{1}{2} \sum_{K\in \mathscr{T}{}}{\left \|h^{2-1/q}{\unicode{x27E6} \nabla U\unicode{x27E7}}\right \|}_{\operatorname{L}^{p}(\partial K)}{\left \|z\right \|}_{\operatorname{W}^{2,q}(\hat{K})}, \end{align}
(3.28)
\begin{align}{\left \|(u - U)_+\right \|}^p_{\operatorname{L}^{p}({\Omega })} \leq l{\!\left ({z-{\hat \Pi }(z)}\right )} - a{\left ({U,z -{\hat \Pi }(z)}\right )} \leq C \sum_{K\in \mathscr{T}{}}{\!\left ({\eta_{K} + \frac{1}{2}\eta_{J}}\right )}|z|_{\operatorname{W}^{2,q}(\hat{K})}. \end{align}
(3.29)
Hence, the result follows from a discrete Cauchy–Schwarz inequality and the regularity bound on \(z\) given in (3.5), after noting that, since \(1/ p+1/ q=1\), □
\begin{align*} {\left \|(u-U)^{p-1}_+\right \|}_{\operatorname {L}^{q}({\Omega })} = {\left \|(u-U)_+\right \|}^{p-1}_{\operatorname {L}^{p}({\Omega })}. \\[-42pt]\end{align*}
We now prove that the negative part of the error also satisfies a bound of the form (3.20). Once again taking an approach analogous to that in [32], we define a second dual problem as follows: let \(\bar{\mathcal{M}} = \{v \mid v \geq 0 \text{ on } \mathcal{A}_{U} \}\). Find \(\bar{z} \in \bar{\mathcal{M}}\) such that
\begin{align} a (\bar{z}, v - \bar{z}) \geq{\left \langle -\max (0, U - u)^{p-1},v-\bar{z}\right \rangle } \quad \forall v \in \bar{\mathcal{M}}. \end{align}
(3.30)
Key properties follow in a manner analogous to Proposition 3.6, and are summarized in Proposition 3.15. This result relies on a mild assumption, a posteriori verifiable from the discrete solution \(U\), which we detail below.
For the finite element approximation \(U\), we say Condition \((\mathbf a_h)\) holds if there exist finitely many points \(d_i\) lying on \(\partial{\Omega }\) and numbers \(\delta_i \gt 0\) such that the intersections of the open balls \(B_{\delta_i}(d_i)\) and the boundary have nonzero distance from each other, and from the corners of the domain, and such that the following hold:
1.
The sets \(B_{\delta_i}(d_i) \cap \partial{\Omega }\) cover the boundary of the discrete contact set, and are such that each \(B_{\delta_i}(d_i) \cap \partial{\Omega }\) contains precisely one element of the boundary of the discrete contact set.
2.
Every connected component of \(\mathcal{A}_U\) has a nonempty interior.
We now have three related conditions: Condition \((\mathbf{A})\), Condition \((\mathbf{A}_h)\), and Condition \((\mathbf a_h)\) The latter two conditions are not discrete analogues of the first; indeed any finite element function may only have a contact set with finitely many connected components. Rather they place some restrictions on the boundaries of discrete contact sets. We note that in the continuous case, the contact set is allowed to contain singletons, whereas this is not allowed for discrete solutions. In addition, Condition \((\mathbf{A}_h)\) assumes topological stability of the discrete contact set as \(h \to 0\).
Notice that the Condition \((\mathbf{a}_h)\) is strictly weaker than its a priori sibling Condition \((\mathbf{A}_h)\) where in addition, the topology of the discrete contact set is assumed to be stable in the limit \(h \to 0\). Condition \((\mathbf{a}_h)\) ensures that the discrete contact set does not include isolated singletons, and that the boundary of the discrete contact set does not include any of the corners of the domain. This condition can be verified a posteriori from the discrete solution \(U\) by examination of the boundary nodal values. For example, since \(U\) is piecewise linear on the boundary, a boundary degree of freedom is a singleton component of \(\mathcal{A}_U\) if and only if it is zero and its neighboring boundary degrees of freedom are both nonzero. Likewise, a degree of freedom at a corner of the domain is in the boundary of \(\mathcal{A}_U\) if it is zero and precisely one of its neighbors is nonzero. Thus, one can iterate over boundary degrees of freedom and if neither of these situations occurs, then Condition \((\mathbf{a}_h)\) holds for \(U\).
We remark that this condition is only required to prove bounds on the negative part of the error. This is due to the different convex sets in which solutions are sought for the two dual problems. The set \(\mathcal{M}\) is defined using the contact set of \(u\), and is therefore \(h\)-independent, whereas the definition of \(\bar{\mathcal{M}}\) involves the discrete contact set.
Let \(u\) solve (2.3) for some \(f \in \operatorname{L}^{2}({\Omega })\), and let \(U \in \mathbb V\) solve (2.8), with \(U\) satisfying Condition \((\mathbf{a}_h)\). Then the dual problem (3.30) is uniquely solvable in \(\operatorname{H}^{1}({\Omega })\). In addition we have \(\bar{z} \leq 0\) in \(\Omega\) and \(\bar{z}=0\) on \(\mathcal{A}_{U}\). For \(p\in (4,\infty )\) and its conjugate \(q = \tfrac{p}{p-1}\) there exists \(C\gt 0\) such that the stability boundholds. In addition, we have the following regularity result:
\begin{align} {\left \|\bar{z}\right \|}_{\operatorname{H}^{1}({\Omega })} \leq C{\left \|\max (0, U-u)^{p-1}\right \|}_{\operatorname{L}^{q}({\Omega })} \end{align}
(3.31)
\begin{align} {\left \|\bar{z}\right \|}_{\operatorname{W}^{2,q}({\Omega })} \leq C{\left \|\max (0, U-u)^{p-1}\right \|}_{\operatorname{L}^{q}({\Omega })}. \end{align}
(3.32)
A closer examination of the proof of Proposition 3.15, given in [16, Lemma 5.6], reveals that the constant in (3.32) hides a dependence on the number of boundary points of the contact set of the discrete solution. To better quantify this regularity constant, we give some details of the proof. For some \(h\), under the given assumptions there are finitely many, \(N_h\), say, critical points (i.e., the boundary points of \(\mathcal{A}_U\)). If \(\psi_i\) is a smooth cut-off function that is identically 1 in a neighborhood of the \(i\)th critical point for \(i=1,\ldots,d\) and \(\psi_0=1-\sum_{i=1}^{N_h} \psi_i\), such that their supports are disjoint, then \(\bar z\) admits a decompositionwhere for each \(i\) we have defined \(z_i = \psi_i \bar{z}\). Each locally supported function is shown to satisfy a regularity boundand thereforewhere \(C\) is the maximum of the \(C_i\). The quality of the upper bound therefore depends upon the behavior of \(N_h\) as \(h \to 0\). Typically it is observed in numerical experiments that the number and location of critical points converge to their true values (see, for example, [16, 38] and section 5). In [38, section 7], the distance between a discrete critical point and its true location was always bounded from above by the mesh size.
\begin{align*} { \bar {z} = z_0 + z_1 + \cdots + z_{N_h},} \end{align*}
\begin{align*} {{\left \|z_i\right \|}_{\operatorname {W}^{2,q}({\Omega })} \leq C_i {\left \|\max (0,U - u)^{p-1}\right \|}_{\operatorname {L}^{q}({\Omega })}}, \end{align*}
\begin{align*} { {\left \|\bar {z}\right \|}_{\operatorname {W}^{2,q}({\Omega })} \leq \sum_i C_i {\left \|\max (0,U - u)^{p-1}\right \|}_{\operatorname {L}^{q}({\Omega })} \leq (N_h+1) C {\left \|\max (0,U - u)^{p-1}\right \|}_{\operatorname {L}^{q}({\Omega })},} \end{align*}
If \(N_h\) changes significantly under mesh refinement, as noted in Remark 3.14, once the finite element approximation is known we are able to determine the number of critical points by examining the boundary degrees of freedom (see Remark 3.14) so that if \(N_h\) does increase, its expected increase in size can be quantified. If \(N_h\) were to increase without limit, this would imply that discrete solutions can oscillate between contact and noncontact over arbitrarily small distances, behavior that is not present at the continuous level in light of the assumed Condition \((\mathbf{A})\). The stronger Condition \((\mathbf{A}_h)\) assumes that this number remains uniformly bounded as \(h \to 0\). Under this assumption there exists a single constant such that the results of Proposition 3.15 are true for all sufficiently small \(h\).
Since \(\bar z \leq 0\), we may now once again use the two-sided approximation. Let \(\check{\Pi }(\bar{z}) = -{\hat \Pi }(-\bar z)\). Then by Assumption 3.10 (proven in section 4), \(\check{\Pi }(\bar{z})\) is a finite element function such that \(\bar{z} \leq \check{\Pi }(\bar{z}) \leq 0\) with required local approximation properties.
Let \(u\) solve problem (2.3), and let \(U\) be the finite element approximation, with \(U\) satisfying Condition \((\mathbf{a}_h)\). Then, for \(p\in \left (4,\infty \right )\), \(q = \tfrac{p}{p-1}\) its conjugate, and \(f\in \operatorname{L}^{p}({\Omega })\) we havewhere
\begin{align} {\left \|(u-U)_-\right \|}_{\operatorname{L}^{p}({\Omega })}^p \leq \eta (U, f, h) := C \sum_{K\in \mathscr{T}{}} \left ( \eta_K^p + \frac 1 2 \eta_J^p \right ), \end{align}
(3.33)
\begin{align} \begin{split} \eta_K &= h^2{\left \|-\Delta U+U - f\right \|}_{\operatorname{L}^{p}(K)}, \\ \eta_J &= h^{2-1/q}{\left \|{\unicode{x27E6} \nabla U\unicode{x27E7}}\right \|}_{\operatorname{L}^{p}(\partial K)}. \end{split} \end{align}
(3.34)
We may immediately take \(v = \bar{z} + u - U\) in the dual problem since \(U=0\) and \(u \geq 0\) on \(\mathcal{A}_{U}\). We therefore have
\begin{align} \int_{{\Omega }} \max (0, U-u)^p ={\left \|(u-U)_-\right \|}^p_{\operatorname{L}^{p}({\Omega })} \leq a(u - U, \bar{z} ). \end{align}
(3.35)
Since we know that \(U=0\) on \(\mathcal{A}_{U}\) and is nonnegative on \(\partial{\Omega }\), and since we also have that \(\bar{z} = 0\) on \(\mathcal{A}_{U}\) and is nonpositive on \(\Omega\), we can choose sufficiently small \(s\gt 0\) so that \(U \pm s \check{\Pi }(\bar{z}) \in{\mathcal K}_h\). Taking \(\Phi = U \pm s \check{\Pi }(\bar{z})\) as test functions in the discrete problem (2.8) givesand therefore we must haveChoosing \(v = u - \check{\Pi }(\bar{z}) \in{\mathcal K}\) in (2.3) we also haveTogether with (3.35), we have shown
\begin{align} \begin{split} a{\left ({U, \check{\Pi }(\bar{z})}\right )} &\geq{\left \langle f,\check{\Pi }(\bar{z})\right \rangle }, \\ a{\left ({U, -\check{\Pi }(\bar{z})}\right )} &\geq -{\left \langle f,\check{\Pi }(\bar{z})\right \rangle }, \end{split} \end{align}
(3.36)
\begin{align} a{\left ({U, \check{\Pi }(\bar{z})}\right )} ={\left \langle f,\check{\Pi }(\bar{z})\right \rangle }. \end{align}
(3.37)
\begin{align} a{\left ({u, -\check{\Pi }(\bar{z})}\right )} \geq -{\left \langle f,\check{\Pi }(\bar{z})\right \rangle }. \end{align}
(3.38)
\begin{align} {\left \|(u-U)_-\right \|}^p_{\operatorname{L}^{p}({\Omega })} \leq a{\left ({u - U, \bar{z} - \check{\Pi }(\bar{z})}\right )}. \end{align}
(3.39)
To introduce the problem data, \(f\), of the primal problem, we choose \(v=u + \check{\Pi }(\bar{z}) - \bar{z}\) as test function in (2.3). Since \(\check{\Pi }(\bar{z})\) lies between the (nonpositive) \(\bar{z}\) and zero, this function lies in \(\mathcal K\), and we obtainor equivalently,Inserting this in (3.39), we now haveThe proof now proceeds exactly as in Theorem 3.12. □
\begin{align} a{\left ({u,\check{\Pi }(\bar{z}) - \bar{z}}\right )} \geq{\left \langle f,\check{\Pi }(\bar{z}) - \bar{z}\right \rangle }, \end{align}
(3.40)
\begin{align} a{\left ({u, \bar{z} - \check{\Pi }(\bar{z})}\right )} \leq{\left \langle f,\bar{z} - \check{\Pi }(\bar{z})\right \rangle }. \end{align}
(3.41)
\begin{align}{\left \|(u-U)_-\right \|}^p_{\operatorname{L}^{p}({\Omega })} \leq{\left \langle f,\bar{z} - \check{\Pi }(\bar{z})\right \rangle } - a{\left ({U, \bar{z} - \check{\Pi }(\bar{z})}\right )}. \end{align}
(3.42)
4. Two-sided approximation.
In this section we develop the bound-preserving interpolant used to prove the a posteriori error bounds in section 3. We construct an interpolant that is bounded from above by a particular function in section 4.2. Then, in section 4.3, we introduce a new interpolant that satisfies bilateral bounds and prove that it enjoys the same optimal approximation properties as the Lagrange interpolant.
Let \(\{ x_i \}_{i=1}^N\) denote the vertices of the triangulation \(\mathscr{T}{}\), and let \(\phi_i \in \mathbb V\) be the \(i\)th canonical basis function, with \(\phi_i(x_j) = \delta_{ij}\) for \(i,j=1, \dots, N\). LetWe also recall that \(\mathcal I\) denotes the (piecewise linear) nodal Lagrange interpolation operator. Let \(q \in (1, \tfrac 43)\), and assume we have a function \(0 \leq z\in \operatorname{W}^{2,q}({\Omega })\). Then the aim of this section is to examine the question of existence of an interpolation operator \({\hat \Pi } : \operatorname{W}^{2,q}({\Omega }) \to \mathbb V\) satisfying both
\begin{align}{\hat x}_j := \bigcup \{ K\in \mathscr{T}{} : \text{supp}(\phi_j) \cap K \neq \emptyset \}. \end{align}
(4.1)
\begin{align} \begin{split} 0 \leq{\hat \Pi }(z) &\leq z \ \text{ and } \\{\left \|z -{\hat \Pi }(z)\right \|}_{\operatorname{L}^{q}({\Omega })} &\leq C h^2{\left \|z\right \|}_{\operatorname{W}^{2,q}({\Omega })}. \end{split} \end{align}
(4.2)
In other words, we wish to show that there exists a two-sided approximation that also has the optimal approximation properties enjoyed by the Lagrange interpolant (see Lemma 3.3). This form of approximation theory arises in many areas of finite element analysis. In particular, the optimal approximation of nonsmooth functions is a nontrivial task relying on local averaging, addressed, for example, in [17, 37].
4.1. From below.
In this section we discuss positivity-preserving interpolation. We note that since we work in two spatial dimensions, and since \(q\in (1,\tfrac 43)\), point values are well defined in \(\operatorname{W}^{2,q}({\Omega })\). In this case, one can simply consider the Lagrange interpolant, which has the properties(see Figure 2(a)), as well as optimal approximation locally and globally:
\begin{align}{\mathcal I} z \geq 0\ \text{ if }\ z \geq 0 \end{align}
(4.3)
\begin{align}{\left \|z -{\mathcal I} z\right \|}_{\operatorname{L}^{q}({\Omega })} \leq C h^2{\left \|z\right \|}_{\operatorname{W}^{2,q}({\Omega })}. \end{align}
(4.4)
If point values do not exist, the approximation is considerably more involved, although it has been tackled in [15] where an interpolant is constructed through local mean-values of the function. It was shown in [15] that the constructed interpolant is \(\operatorname{L}^{q}\)-stable, second order accurate, and linear. See also [42, 43] for related ideas using a nonlinear interpolation operator [18].
4.2. From above.
In this section we examine the design of interpolants that are bounded from above by the function to be interpolated. This means that we seek an interpolation operator \(\Pi\) that satisfies \({\Pi }(z) \leq z\). In this case we can see the Lagrange interpolant does not satisfy the requirements; see Figure 2(a). Indeed, any strictly convex function will see the Lagrange interpolant violate this bound. It can, however, be modified. Indeed, consider the function
\begin{align} {\Pi }(z)(x) = \sum_{i=1}^N{\!\left ({{\mathcal I} z(x_i) - \max_{y\in{\hat x}_i} ({\mathcal I}(z)(y) - z(y))}\right )}\phi_i(x) =:{\mathcal I}(z)(x) - R(x). \end{align}
(4.5)
For \(0 \leq z \in \operatorname{W}^{2,q}({\Omega })\) the approximation \({\Pi }(z)\) defined in (4.5) satisfies
\begin{align}{\Pi }(z) \leq z \textit{ in }{\Omega }. \end{align}
(4.6)
Note that by definition of \(R\) we have □
\begin{align}{\Pi }(z) ={\mathcal I}(z) - R \leq{\mathcal I}(z) - ({\mathcal I}(z) - z) = z. \\[-36pt]\nonumber\end{align}
(4.7)
This is a nonlinear interpolant since for general \(u,v\in \operatorname{W}^{2,q}({\Omega })\)This means we cannot directly apply Bramble–Hilbert techniques to obtain optimal approximation under minimal regularity. Instead we must take a different approach.
\begin{align}{\Pi }(u+v) \neq{\Pi }(u) +{\Pi }(v). \end{align}
(4.8)
Suppose \(z\in \operatorname{W}^{2,q}({\Omega })\) for \(q\in (1,\tfrac 43)\). Then there exists \(C \gt 0\) independent of \(h\) such that for all elements \(K\in \mathscr{T}{}\) the approximation \({\Pi }(z)\) defined through (4.5) satisfies
\begin{align}{\left \|z -{\Pi }(z)\right \|}_{\operatorname{L}^{q}(K)} \leq C h^2{\left \|z\right \|}_{\operatorname{W}^{2,q}(\hat{K})}. \end{align}
(4.9)
To begin, we noteWe have from Lemma 3.3 thatTo control the second term, notice that since \(R|_K \in \mathcal{R}^1(K)\), for a \(K\in \mathscr{T}{}\) we can write a lumped \(\operatorname{L}^{q}\)-norm to see that (cf. [21, Proposition 12.5])
\begin{align} {\left \|z -{\Pi }(z)\right \|}_{\operatorname{L}^{q}(K)} \leq{\left \|z -{\mathcal I}(z)\right \|}_{\operatorname{L}^{q}(K)} +{\left \|R\right \|}_{\operatorname{L}^{q}(K)}. \end{align}
(4.10)
\begin{align} {\left \|z -{\mathcal I}(z)\right \|}_{\operatorname{L}^{q}(K)} \leq C h^2{\left \|z\right \|}_{\operatorname{W}^{2,q}(K)}. \end{align}
(4.11)
\begin{align} \begin{split}{\left \|R\right \|}_{\operatorname{L}^{q}(K)}^q &\leq C h^2 \sum_{x_i \in K}{\left \|R(x_i)\right \|}^q \\ &\leq C{ h^2} \sum_{x_i \in K}{\left \|{\mathcal I}(z) - z\right \|}_{\operatorname{L}^{\infty }({\hat x}_i)}^q. \end{split} \end{align}
(4.12)
Now, a further property of the Lagrange interpolant (see [11, Corollary 4.4.7]) isHence substituting into (4.12) we havewhich, together with (4.10) and (4.11), completes the proof. □
\begin{align}{\left \|z -{\mathcal I}(z)\right \|}_{\operatorname{L}^{\infty }(K)} \leq C h^{2-\tfrac 2 q}{\left \|z\right \|}_{\operatorname{W}^{2,q}(K)}. \end{align}
(4.13)
\begin{align} \begin{split}{\left \|R\right \|}_{\operatorname{L}^{q}(K)}^q \leq C h^{2q}{\left \|z\right \|}_{\operatorname{W}^{2,q}(\hat K)}^q, \end{split} \end{align}
(4.14)
Note the interpolant we construct would not be well defined in three spatial dimensions, \(n=3\), for the range of values of \(q\) considered here, as the construction requires point values which, in turn, requires \(z\in \operatorname{W}^{2,s}({\Omega })\) for \(2s \gt n\).
4.3. Bilateral approximation.
There is much less in the literature on the design of two-sided approximations. The only arguments date back to Mosco and Strang (see [33, 40, 39]), who consider piecewise linear approximations for dimension \(n=1,2,3\) and \(q=2\). Unfortunately these arguments do not appear to extend to the case at hand.
The difficulty in this problem can be seen by examining Figure 2(b) around where \(z\sim 0\). To keep a bilateral constraint on the approximation, any appropriate function is squeezed and the only mechanism is to force it to zero on a surrounding patch.
4.4. Optimal approximation over \({\boldsymbol{\mathcal{R}}}^{\boldsymbol{1}}\boldsymbol{(K)}\).
We modify \(\Pi\) (defined in (4.5)) in the following way at the degrees of freedom:as illustrated in Figure 2(c).
\begin{align}{\hat \Pi }(z)(x_i) = \begin{cases}{\Pi }(z)(x_i)\ \text{ if }\ {\Pi }(z) \geq 0\ \text{ on }\ {\hat x}_i, \\ 0\ \text{ otherwise,} \end{cases} \end{align}
(4.15)
Let \(0 \leq z\in \operatorname{W}^{2,q}({\Omega })\), \(q\in (1,\tfrac 43)\), and defineThen \({\hat \Pi }(z) \in{\mathcal Z}_h\). Further, there exists a constant \(C_b \gt 0\) such that for all elements \(K\in \mathscr{T}{}\)
\begin{align}{\mathcal Z}_h := \{ z_h \in \mathbb V : 0 \leq z_h \leq z \}. \end{align}
(4.16)
\begin{align}{\left \|z -{\hat \Pi }(z)\right \|}_{\operatorname{L}^{q}(K)} \leq C_{b} h^2{\left \|z\right \|}_{\operatorname{W}^{2,q}(\hat K)}. \end{align}
(4.17)
To begin, notice that due to the regularity of \(z\) we can find a constant \(C\) independent of \(h\) such thatwhereControl of the first term is given in Theorem 4.2. For the second, notice that \(S\) only has support when \({\mathcal I}(z) \gt z\) in a vicinity of when \(z\) vanishes. This can only happen if \(z\) is locally convex. Further, since \({\mathcal I}(z) \gt 0\) whenever \(z\gt 0\), we have thatand henceconcluding the proof. □
\begin{align}{\left \|z -{\hat \Pi }(z)\right \|}_{\operatorname{L}^{q}(K)} \leq{\left \|z -{\Pi }(z)\right \|}_{\operatorname{L}^{q}(K)} + C{\left \|S\right \|}_{\operatorname{L}^{q}(K)}, \end{align}
(4.18)
\begin{align} S = \max{\!\left ({-{\Pi }(z), 0}\right )}. \end{align}
(4.19)
\begin{align}{\left \|\max{\!\left ({-{\Pi }(z),0}\right )}\right \|} \leq{\left \|R\right \|} \end{align}
(4.20)
\begin{align}{\left \|S\right \|}_{\operatorname{L}^{q}(K)} \leq{\left \|R\right \|}_{\operatorname{L}^{q}(K)} \leq C h^2{\left \|z\right \|}_{\operatorname{W}^{2,q}(\hat K)}, \end{align}
(4.21)
Two-sided bounds are conjectured in [33] for higher dimension, at least with \(q=2\), although their results only seem to hold with dimension \(n\leq 3\). The conjecture we make is that optimal bilateral approximations in \(\operatorname{L}^{q}({\Omega })\) are only possible when \(2-\tfrac n q \gt 0\); hence we believe our operator is a constructive example of that studied in [33], i.e., is valid for \(n=3\) and \(q=2\).
5. Numerical examples.
In this section, we present numerical results to demonstrate the effectiveness of the error estimate and adaptive routine against an exact solution within the framework of sections 2 and 3, that is, convex \(\Omega \subseteq \mathbb{R}^2\) with polygonal boundary. We then present some more challenging situations in which the theory of the preceding sections may fail, but in which the error estimate may still prove useful. In particular, we test the estimate on a nonconvex domain in \(\mathbb{R}^2\) with a re-entrant corner as well as a three-dimensional example.
All simulations presented here are conducted using deal.II, an open source C++ software library providing tools for adaptive finite element computations [3]. We note here that deal.II uses quadrilateral meshes. Since all coarse meshes are uniform meshes consisting of squares, and refinement is by quadrisection, regularity of the mesh does not degrade under heavy refinement. At each stage of the iterative process used to solve the variational inequality, a preconditioned conjugate gradient method is used to solve the algebraic problem. As the mesh is refined, hanging nodes are created, and are constrained so that the resulting numerical solution is continuous.
5.1. Example 1.
For our first example, we choose the exact solution used by the authors of [16]. Let \(r, \theta\) denote polar coordinates centered at \((0.5, 0)\), that is,We define the functionand define \(\psi\) to be a ninth order spline defined by endpoint values and gradients as follows (see also [16, section 6]). We impose the following values to determine \(\psi\):and \(\psi (s)=0\) for any \(s \lt 0\) or \(s \geq 0.45\). Then if we choose \(f\) accordingly, \((r, \theta ) \mapsto 10 \psi (r) \tilde{u}(r, \theta )\) solves (2.3). We also point out that this exact solution satisfies Condition \((\mathbf{A})\), since by the definition its contact set isthat is, a single connected component consisting of an open interval in \(\partial{\Omega }\). All numerical approximations satisfy Condition \((\mathbf a_h)\). An example of the verification (see Remark 3.14) of this condition is shown in Figure 5, where the discrete contact set for Example 1 is shown for varying levels of mesh refinement.
\begin{align} r(x,y) = ((x-0.5)^2 + y^2)^{1/ 2}, \end{align}
(5.1)
\begin{align} \theta (x,y) = \arccos \left (\frac{x - 0.5}{r}\right ). \end{align}
(5.2)
\begin{align} \tilde{u}(r, \theta ) : = - r^{3/ 2} \sin (\tfrac{3}{2} \theta ), \end{align}
(5.3)
\begin{align} \psi (0)=1, \,\, \psi^{\prime }(0) = \cdots = \psi^{\prime \prime \prime \prime }(0) = \psi (0.45) = \psi^{\prime }(0.45) = \cdots = \psi^{\prime \prime \prime \prime }(0.45)=0, \end{align}
(5.4)
\begin{align*} {\partial {\Omega } \backslash \{(x,0) : 0.05\leq x \leq 0.5\},} \end{align*}
The problem is initialized on a coarse uniform mesh with \(h=1/ 8\). For the adaptive algorithm we use Dörfler marking (see [19]) with refinement parameter \(\beta=0.9\) and no coarsening. A numerical approximation to problem (2.3) with datum \(f\) chosen so that \(10 \psi \tilde{u}\) is the exact solution is shown in Figure 3(a), represented on the final mesh \(\mathcal{T}^{15}\) produced by the adaptive algorithm, consisting of around 80,000 degrees of freedom. The progression of the \(\operatorname{L}^{4}({\Omega })\)-error, \({\left \|u - U\right \|}_{\operatorname{L}^{4}({\Omega })}\), and error estimate \(\eta^{1/ 4}\) as defined in equation (3.20) under adaptive mesh refinement is shown in Figure 3(b). For a given number of degrees of freedom, the adaptive algorithm reduces the error and appears to give a slightly better convergence rate than uniformly refining the mesh in terms of degrees of freedom. We observe that the overestimation factor (or effectivity index as it is often called) of the error estimate is approximately 20. As expected from the theory, this factor is approximately constant once the asymptotic regime is reached. Comparing with Figure 7(a), we see that asymptotically, error is lower under adaptive refinement than uniform refinement in terms of degrees of freedom, although both are optimal.
A selection of adaptive meshes are shown in Figures 3(c)–3(f). These meshes show refinement around the main features of the solution. In particular, the mesh is heavily refined around large gradients, along with significant refinement where the boundary conditions change type and constraints are active.
Finally, we demonstrate the ability of adaptive mesh refinement, driven by our error estimate, to compensate for suboptimal convergence in \(\operatorname{L}^{p}({\Omega })\) for large \(p\), observed in [16, Table 2]. We set \(p=32\) and adapt the mesh using Dörfler marking with refinement parameter \(\beta=0.8\) and no coarsening. In Figure 4(b), we see that the error under uniform mesh refinement is asymptotically suboptimal, but that optimal rates in terms of number of degrees of freedom are achieved with the adaptive scheme. Spatial distribution of the error estimate is shown in Figure 4(a).
5.2. Example 2: Re-entrant corner.
To test the estimate in the presence of a geometric singularity, we introduce a re-entrant corner to the domain. In this example, data \(f\) is chosen to ensure that both boundary constraints are active. The problem data is selected to try to force the solution to be close towhere \(r = (x^2 + y^2)^{\frac 12}\) and \(b\) can be varied to force different behaviors of the solution. For this example we make the choice \(b=0.91\) and set \(f = \Delta w+w\).
\begin{align} w := \sin (2 \pi (r - b)^2) - 0.5, \end{align}
(5.5)
A numerical solution to this problem is shown in Figure 6(a). Under uniform mesh refinement, the error estimate converges to zero at a suboptimal rate. Optimality is restored using an adaptive routine utilizing Dörfler marking with refinement fraction of 0.8 and no coarsening. A sample of the meshes produced is given in Figures 6(b)–6(e). The behavior of the error estimate under uniform and adaptive refinement is shown in Figure 8. The slopes of the error-reduction curves reveal suboptimal convergence in the uniform case which is somewhat recovered by adaptive mesh refinement, but to different degrees. This time we see most refinement around gradients and features in the solution. Interestingly, there is little refinement around the re-entrant corner where the solution is expected to lose regularity. Instead, the error estimate prioritizes resolution of the solution near where the boundary constraints are active. This simulation was performed with the Dörfler marking criterion with refinement parameter 0.8 and no coarsening. The algorithm was initialized on a coarse uniform mesh with \(h=1/ 16\).
5.3. Example 3, \(\boldsymbol{d=3}\).
Working in three dimensions, let \((r,\theta, \phi )\) denote spherical polar coordinates centered at \((0.5,0,0)\). We note that the functionsatisfies appropriate constraints on the boundary and so solves (1.1), (1.2) for appropriately defined \(f\). Contours of (5.6) are shown in Figure 9. Again we use Dörfler marking with refinement parameter of 0.8 and no coarsening.
\begin{align} (r,\theta, \phi ) \mapsto - 10 \psi (r) r^{3/ 2} \sin \left ( \frac{3}{2} \varphi \right ) \sin \left (\frac{3}{4} \theta \right ) \end{align}
(5.6)
Adaptive meshes from several stages of the adaptive algorithm are displayed in Figure 10 (note that these are slices of a three-dimensional mesh). We observe that mesh resolution is greater close to the boundary where the Signorini constraints are active (in this case the plane defined by \(y=0\)) and provides good resolution of the contact set. Refinement is also present around key features of the solution, as observed for Example 1 in Figures 3(a) and 3(c)–3(f).
6. Conclusions and discussion.
In this article, rigorous error estimates were derived for the scalar Signorini problem in nonenergy norms, specifically \(\operatorname{L}^{p}({\Omega })\) with \(p\in (4,\infty )\), using a duality argument. The a posteriori error estimates were benchmarked against known exact solutions and shown to have the same asymptotic rates as the a priori analysis, and to be sharp to a satisfactory degree (overestimation factor of approximately 20). Adaptive mesh refinement was shown to decrease error for a given number of degrees of freedom compared to uniform meshes, particularly for a challenging three-dimensional problem. Furthermore, since asymptotic convergence is suboptimal in \(\operatorname{L}^{p}({\Omega })\) for large \(p\) we demonstrated that adaptive schemes are able to regain optimal rates in terms of the number of degrees of freedom.
The error estimate was tested on more challenging test cases, both in situations where the assumptions of the main theorem were satisfied, and in a case where they were not (three spatial dimensions, re-entrant corners). In the three-dimensional case where the theory does not hold (see section 5.3 and Figure 8(b)), we remark that the error estimate is too optimistic in the sense that the error estimate decreases at the expected rate under uniform refinement, whereas the error does not (see Figure 7(b)). This means that although adaptivity did produce superior error reduction in this case, it can only be used as an indicator, not an estimator.
This paper presents many avenues for further investigation which include the analysis in higher spatial dimensions, as already indicated. Further work is underway on nontrivial obstacles as well as the generalization to higher polynomial degrees. While we do not believe the bilateral interpolation result can be true globally, which is critical to the analysis, we note the recent work [5], where the authors give an alternative nodal interpretation that may shed some light on the quest for higher order methods.
References
1.
M. Ainsworth and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Pure Appl. Math. (N. Y.), John Wiley & Sons, 2011.
2.
T. Apel and S. Nicaise, Regularity of the solution of the scalar Signorini problem in polygonal domains, Results Math., 75 (2020), 75.
3.
D. Arndt, W. Bangerth, T. C. Clevenger, D. Davydov, M. Fehling, D. Garcia-Sanchez, G. Harper, T. Heister, L. Heltai, M. Kronbichler, et al., The deal.II library, Version 9.1, J. Numer. Math., 27 (2019), pp. 203–213.
4.
B. Ashby, C. Bortolozo, A. Lukyanov, and T. Pryer, Adaptive modelling of variably saturated seepage problems, Quart. J. Mech. Appl. Math., 74 (2021), pp. 55–81.
5.
G. R. Barrenechea, E. Georgoulis, T. Pryer, and A. Veeser, A nodally bound-preserving finite element method, IMA J. Numer. Anal., (2023).
6.
F. B. Belgacem, P. Hild, and P. Laborde, Extension of the mortar finite element method to a variational inequality modeling unilateral contact, Math. Models Methods Appl. Sci., 9 (1999), pp. 287–303.
7.
Z. Belhachmi and F. Belgacem, Quadratic finite element approximation of the Signorini problem, Math. Comp., 72 (2003), pp. 83–104.
8.
F. Ben Belgacem, Numerical simulation of some variational inequalities arisen from unilateral contact problems by the finite element methods, SIAM J. Numer. Anal., 37 (2000), pp. 1198–1216, https://doi.org/10.1137/S0036142998347966.
9.
H. Blum and F. Suttmeier, An adaptive finite element discretisation for a simplified Signorini problem, Calcolo, 37 (2000), pp. 65–77.
10.
D. Braess, A posteriori error estimators for obstacle problems—another look, Numer. Math., 101 (2005), pp. 415–421.
11.
S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer Science & Business Media, 2008.
12.
F. Brezzi, W. W. Hager, and P. Raviart, Error estimates for the finite element solution of variational inequalities: Part I. Primal theory, Numer. Math., 28 (1977), pp. 431–443.
13.
F. Brezzi, W. W. Hager, and P. Raviart, Error estimates for the finite element solution of variational inequalities: Part II. Mixed methods, Numer. Math., 31 (1978), pp. 1–16.
14.
F. Brezzi and G. Sacchi, A finite element approximation of a variational inequality related to hydraulics, Calcolo, 13 (1976), pp. 257–273.
15.
Z. Chen and R. H. Nochetto, Residual type a posteriori error estimates for elliptic obstacle problems, Numer. Math., 84 (2000), pp. 527–548.
16.
C. Christof and C. Haubner, Finite element error estimates in non-energy norms for the two-dimensional scalar Signorini problem, Numer. Math., 145 (2020), pp. 513–551.
17.
P. Clément, Approximation by finite element functions using local regularization, Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge Anal. Numér., 9 (1975), pp. 77–84.
18.
R. A. DeVore, Nonlinear approximation, Acta Numer., 7 (1998), pp. 51–150.
19.
W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal., 33 (1996), pp. 1106–1124, https://doi.org/10.1137/0733054.
20.
G. Drouet and P. Hild, Optimal convergence for discrete variational inequalities modelling Signorini contact in 2D and 3D without additional assumptions on the unknown contact set, SIAM J. Numer. Anal., 53 (2015), pp. 1488–1507, https://doi.org/10.1137/140980697.
21.
A. Ern and J. Guermond, Finite Elements I: Approximation and Interpolation, Springer, 2021.
22.
R. S. Falk, Error estimates for the approximation of a class of variational inequalities, Math. Comput., 28 (1974), pp. 963–971.
23.
P. Grisvard, Elliptic Problems in Nonsmooth Domains, Classics Appl. Math. 69, SIAM, Philadelphia, 2011, https://doi.org/10.1137/1.9781611972030.
24.
J. Haslinger, I. Hlaváček, and J. Nečas, Numerical methods for unilateral problems in solid mechanics, in Handbook of Numerical Analysis, Vol. IV, North-Holland, Amsterdam, 1996, pp. 313–485.
25.
P. Hild and S. Nicaise, A posteriori error estimations of residual type for Signorini’s problem, Numer. Math., 101 (2005), pp. 523–549.
26.
P. Hild and Y. Renard, An improved a priori error analysis for finite element approximations of Signorini’s problem, SIAM J. Numer. Anal., 50 (2012), pp. 2400–2419, https://doi.org/10.1137/110857593.
27.
R. H. W. Hoppe and R. Kornhuber, Adaptive multilevel methods for obstacle problems, SIAM J. Numer. Anal., 31 (1994), pp. 301–323, https://doi.org/10.1137/0731016.
28.
S. Hüeber and B. I. Wohlmuth, An optimal a priori error estimate for nonlinear multibody contact problems, SIAM J. Numer. Anal., 43 (2005), pp. 156–173, https://doi.org/10.1137/S0036142903436678.
29.
D. Kinderlehrer and G. Stampacchia, An Introduction to Variational Inequalities and Their Applications, SIAM, Philadelphia, 2000, https://doi.org/10.1137/1.9780898719451.
30.
R. Krause, A. Veeser, and M. Walloth, An efficient and reliable residual-type a posteriori error estimator for the Signorini problem, Numer. Math., 130 (2015), pp. 151–197.
31.
J.-L. Lions and G. Stampacchia, Variational inequalities, Comm. Pure Appl. Math., 20 (1967), pp. 493–519.
32.
U. Mosco, Error estimates for some variational inequalities, in Mathematical Aspects of Finite Element Methods, Springer, 1977, pp. 224–236.
33.
U. Mosco and G. Strang, One-sided approximation and variational inequalities, Bull. Amer. Math. Soc., 80 (1974), pp. 308–312.
34.
J. Nitsche, \(L^{\infty}\)-convergence of finite element approximations, in Mathematical Aspects of Finite Element Methods, Springer, 1977, pp. 261–274.
35.
R. H. Nochetto, K. G. Siebert, and A. Veeser, Pointwise a posteriori error control for elliptic obstacle problems, Numer. Math., 95 (2003), pp. 163–195.
36.
F. Scarpini and M. A. Vivaldi, Error estimates for the approximation of some unilateral problems, RAIRO. Anal. Numér., 11 (1977), pp. 197–208.
37.
L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp., 54 (1990), pp. 483–493.
38.
O. Steinbach, B. Wohlmuth, and L. Wunderlich, Trace and flux a priori error estimates in finite-element approximations of Signorni-type problems, IMA J. Numer. Anal., 36 (2015), pp. 1072–1095.
39.
G. Strang, The dimension of piecewise polynomial spaces, and one-sided approximation, in Conference on the Numerical Solution of Differential Equations, Springer, 1974, pp. 144–152.
40.
G. Strang, One-sided approximation and plate bending, in Computing Methods in Applied Sciences and Engineering, Part 1, Springer-Verlag, 1974, pp. 140–155.
41.
F.-T. Suttmeier, Numerical Solution of Variational Inequalities by Adaptive Finite Elements, Springer, 2008.
42.
A. Veeser, Approximating gradients with continuous piecewise polynomial functions, Found. Comput. Math., 16 (2016), pp. 723–750.
43.
A. Veeser, Positivity preserving gradient approximation with linear finite elements, Comput. Methods Appl. Math., 19 (2019), pp. 295–310.
44.
M. Walloth, A reliable, efficient and localized error estimator for a discontinuous Galerkin method for the Signorini problem, Appl. Numer. Math., 135 (2019), pp. 276–296.
45.
A. Weiss and B. I. Wohlmuth, A posteriori error estimator for obstacle problems, SIAM J. Sci. Comput., 32 (2010), pp. 2627–2658, https://doi.org/10.1137/090773921.
46.
C. Zhang, Adaptive Finite Element Methods for Variational Inequalities: Theory and Applications in Finance, Ph.D. thesis, 2007.
Information & Authors
Information
Published In
Copyright
© 2024 SIAM. Published by SIAM under the terms of the Creative Commons 4.0 license. This work is licensed under a Creative Commons Attribution 4.0 International License.
History
Submitted: 14 November 2022
Accepted: 29 March 2024
Published online: 23 July 2024
Keywords
MSC codes
Authors
Funding Information
Leverhulme Trust: RPG-2021-238
Funding: The first author was supported through a Ph.D. scholarship awarded by the EPSRC Centre for Doctoral Training in the Mathematics of Planet Earth at Imperial College London and the University of Reading EP/L016613/1. The second author received partial support from the EPSRC programme grant EP/W026899/1. Both authors were supported by the Leverhulme RPG-2021-238.
Metrics & Citations
Metrics
Citations
If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.
Cited By
There are no citations for this item
View Options
View options
- Access via your Institution
- Questions about how to access this content? Contact SIAM at [email protected].