Residual-based a posteriori error analysis for symmetric mixed Arnold-Winther FEM

TL;DR

This paper develops an explicit residual-based a posteriori error estimator for the symmetric mixed finite element method in linear elasticity, enabling reliable adaptive mesh refinement and improved convergence rates.

Contribution

It introduces a new explicit residual-based a posteriori error estimator that is reliable, efficient, and easy to implement for Arnold-Winther FEM in linear elasticity.

Findings

Estimator is reliable and efficient

Adaptive refinement improves convergence rates

Higher accuracy with adaptive mesh in non-smooth cases

Abstract

This paper introduces an explicit residual-based a posteriori error analysis for the symmetric mixed finite element method in linear elasticity after Arnold-Winther with pointwise symmetric and H(div)-conforming stress approximation. Opposed to a previous publication, the residual-based a posteriori error estimator of this paper is reliable and efficient and truly explicit in that it solely depends on the symmetric stress and does neither need any additional information of some skew symmetric part of the gradient nor any efficient approximation thereof. Hence it is straightforward to implement an adaptive mesh-refining algorithm obligatory in practical computations. Numerical experiments verify the proven reliability and efficiency of the new a posteriori error estimator and illustrate the improved convergence rate in comparison to uniform mesh-refining. A higher convergence rates for piecewise affine data is observed in the L2 stress error and reproduced in non-smooth situations by the adaptive mesh-refining strategy.

Equations

$$\underset{v\in V}{\inf }\Vert {\mathbb{C}}^{-1/2}{\sigma }_h-{\mathbb{C}}^{1/2}\epsilon (v){\Vert }_{L^2(\Omega )}\lesssim \eta (\mathcal{T},{\sigma }_h).$$

$$h_E^{1/2}\Vert {\tau }_E\cdot [{\mathbf{\epsilon }}_h{]}_E{\tau }_E{\Vert }_{L^2(E)}+h_E^{3/2}\Vert {\tau }_E\cdot [{\mathrm{rot}}_{NC}{\mathbf{\epsilon }}_h{]}_E-\partial ({\nu }_E\cdot [{\mathbf{\epsilon }}_h{]}_E{\tau }_E)/\partial s{\Vert }_{L^2(E)},$$

$$\begin{array}{r}-\mathrm{div}\sigma =f\text{and}\sigma =\mathbb{C}\epsilon (u)\text{mbox}in\Omega ,\\ u=u_D\text{mbox}on{\Gamma }_D,\sigma \nu =g\text{mbox}on{\Gamma }_N.\end{array}$$

$$V:=\{v\in H^{1}(\Omega;\mathbb{R}^{2})\ \big{|}\ v|_{\Gamma_{D}}=0\}\text{ and }H(\operatorname{div},\Omega;\mathbb{S}):=\{\tau\in L^{2}(\Omega;\mathbb{S})\ \big{|}\ \operatorname{div}\tau\in L^{2}(\Omega;\mathbb{R}^{2})\},$$

$$\Sigma(g):=\Big{\{}\sigma\in H(\operatorname{div},\Omega;\mathbb{S})\ \big{|}\int_{\Gamma_{N}}\psi\cdot(\sigma\nu)\;ds=\int_{\Gamma_{N}}\psi\cdot g\;ds\text{ for all }\psi\in V\Big{\}}$$

$$\begin{array}{rl}{\int }_{\Omega }\sigma :{\mathbb{C}}^{-1}\tau dx+{\int }_{\Omega }u\cdot \mathrm{div}\tau dx& ={\int }_{{\Gamma }_D}{u}_D\cdot (\tau \nu )ds\text{for all }\tau \in {\Sigma }_0,\\ {\int }_{\Omega }v\cdot \mathrm{div}\sigma dx& =-{\int }_{\Omega }f\cdot vdx\text{for all }v\in L^2(\Omega ;{\mathbb{R}}^2).\end{array}$$

$$AW_{k}(\mathcal{T}):=\big{\{}\tau\in P_{k+2}(\mathcal{T};\mathbb{S})\cap H(\operatorname{div},\Omega;\mathbb{S})\,\big{|}\,\operatorname{div}\tau\in P_{k}(\mathcal{T};\mathbb{R}^{2})\big{\}}.$$

$$G(\mathcal{T}):=\{(\tau_{h}\nu)|_{\Gamma_{N}}\in L^{2}(\Gamma_{N};\mathbb{R}^{2})\,\big{|}\,\tau_{h}\in AW_{k}(\mathcal{T})\}$$

$$\Sigma (g_h,\mathcal{T}):=\Sigma (g)\cap A{W}_k(\mathcal{T})$$

$$\begin{array}{rl}{\int }_{\Omega }{\sigma }_h:{\mathbb{C}}^{-1}{\tau }_{h}dx+{\int }_{\Omega }{u}_h\cdot \mathrm{div}{\tau }_{h}dx& ={\int }_{{\Gamma }_D}{u}_D\cdot ({\tau }_h\nu )ds\text{for all }{\tau }_h\in \Sigma (0,\mathcal{T}),\\ {\int }_{\Omega }{v}_h\cdot \mathrm{div}{\sigma }_{h}dx& ={\int }_{\Omega }f\cdot v_{h}dx\text{ for all }{v}_h\in V_h.\end{array}$$

$$\Vert \sigma -{\sigma }_h{\Vert }_{L^2(\Omega )}$$

$$\Vert u-u_h{\Vert }_{L^2(\Omega )}$$

$$[v{]}_E:=(v{|}_{T_{+}}){|}_E-(v{|}_{T_{-}}){|}_E.$$

$${\eta }^2(\mathcal{T},{\sigma }_h):=\sum _{T\in \mathcal{T}}{h}_T^4\Vert \mathrm{rot}\mathrm{rot}{\mathbf{\epsilon }}_h{\Vert }_{L^2(T)}^2+{\mathrm{osc}}^2(f,\mathcal{T})+{\mathrm{osc}}^2(g-g_h,\mathcal{E}({\Gamma }_N))$$

$$+\sum _{E\in \mathcal{E}(\Omega )}\left(h_E\Vert {\tau }_E\cdot [{\mathbf{\epsilon }}_h{]}_E{\tau }_E{\Vert }_{L^2(E)}^2+h_E^3\Vert {\tau }_E\cdot ([{\mathrm{rot}}_{NC}{\mathbf{\epsilon }}_h{]}_E-\frac{\partial [{\mathbf{\epsilon }}_h{]}_E{\nu }_E}{\partial s}){\Vert }_{L^2(E)}^2\right)$$

$$+\sum _{E\in \mathcal{E}({\Gamma }_D)}\left(h_E\Vert {\tau }_E\cdot ({\mathbf{\epsilon }}_h{\tau }_E-\frac{\partial u_D}{\partial s}){\Vert }_{L^2(E)}^2+h_E^3\Vert {\tau }_E\cdot \mathrm{rot}{\mathbf{\epsilon }}_h-{\nu }_E\cdot (\frac{\partial {\mathbf{\epsilon }}_h{\tau }_E}{\partial s}+\frac{{\partial }^2{u}_D}{\partial s^2}){\Vert }_{L^2(E)}^2\right)$$

$$\Vert \sigma -{\sigma }_h{\Vert }_{L^2(\Omega )}\le C_{\text{rel}}\eta (\mathcal{T},{\sigma }_h).$$

$${\mathrm{osc}}^2(f,\mathcal{T}):=\sum _{T\in \mathcal{T}}{h}_T^2\Vert f-f_h{\Vert }_{L^2(T)}^2\text{ for the }{L}^2\text{ projection }{f}_h\text{ of }f\text{ onto }{P}_k(\mathcal{T};{\mathbb{R}}^2);$$

$${\mathrm{osc}}^2(g-g_h,\mathcal{E}({\Gamma }_N)):=\sum _{E\in \mathcal{E}({\Gamma }_N)}{h}_E\Vert g-g_h{\Vert }_{L^2(E)}^2.$$

$${\mathrm{osc}}_I^2(u_D,E)$$

$${\mathrm{osc}}_{II}^2(u_D,E)$$

$${\mathrm{osc}}^2(u_D,\mathcal{E}({\Gamma }_D)):=\sum _{E\in \mathcal{E}({\Gamma }_D)}\left({\mathrm{osc}}_I^2(u_D,E)+{\mathrm{osc}}_{II}^2(u_D,E)\right).$$

$$C_{\text{eff}}^{-1}\eta (\mathcal{T},{\sigma }_h)\le \Vert \sigma -{\sigma }_h{\Vert }_{L^2(\Omega )}+\mathrm{osc}(f,\mathcal{T})+\mathrm{osc}(g-g_h,\mathcal{E}({\Gamma }_N))+\mathrm{osc}(u_D,\mathcal{E}({\Gamma }_D)).$$

$$\mathrm{Curl}\Phi :=(-{\partial }_2{\Phi }_1,{\partial }_1{\Phi }_1;-{\partial }_2{\Phi }_2,{\partial }_1{\Phi }_2)=D\Phi (0,1;-1,0).$$

$$({\mathbf{\epsilon }}_h,{\mathrm{Curl}}^2\varphi {)}_{L^2(\Omega )}=({\mathrm{rot}}_{NC}{\mathrm{rot}}_{NC}{\mathbf{\epsilon }}_h,\varphi {)}_{L^2(\Omega )}$$

$$+\sum _{E\in \mathcal{E}(\Omega )}\left(({\tau }_E\cdot [{\mathbf{\epsilon }}_h{]}_E{\tau }_E,{\partial }_{{\nu }_E}\varphi {)}_{L^2(E)}-([{\mathrm{rot}}_{NC}{\mathbf{\epsilon }}_h{]}_E-\frac{\partial [{\mathbf{\epsilon }}_h{]}_E{\nu }_E}{\partial s},\varphi {\tau }_E{)}_{L^2(E)}\right)$$

$$+\sum _{E\in \mathcal{E}(\partial \Omega )}\left(({\tau }_E\cdot {\mathbf{\epsilon }}_h{\tau }_E,{\partial }_{{\nu }_E}\varphi {)}_{L^2(E)}-(\mathrm{rot}{\mathbf{\epsilon }}_h-\frac{\partial {\mathbf{\epsilon }}_h{\nu }_E}{\partial s},\varphi {\tau }_E{)}_{L^2(E)}\right).$$

$${\int }_{\Omega }\Psi :\mathrm{Curl}\varphi dx+{\int }_{\Omega }\varphi \cdot \mathrm{rot}\Psi dx={\int }_{\partial \Omega }\varphi \cdot \Psi {\tau }_{E}ds.$$

$${\tau }_E\cdot \mathrm{Curl}\phi =\partial \phi /\partial {\nu }_E\text{and}{\nu }_E\cdot \mathrm{Curl}\phi =-\partial \phi /\partial s=-\partial \phi /\partial {\tau }_E\text{on }E\in \mathcal{E}.$$

$$({\mathbf{\epsilon }}_h,{\mathrm{Curl}}^2\varphi {)}_{L^2(\Omega )}=({\mathrm{rot}}_{NC}^2{\mathbf{\epsilon }}_h,\varphi {)}_{L^2(\Omega )}$$

$$+\sum _{E\in \mathcal{E}(\Omega )}\left(([{\mathbf{\epsilon }}_h{]}_E{\tau }_E,\mathrm{Curl}\varphi {)}_{L^2(E)}-([{\mathrm{rot}}_{NC}{\mathbf{\epsilon }}_h{]}_E\cdot {\tau }_E,\varphi {)}_{L^2(E)}\right)$$

$$+\sum _{E\in \mathcal{E}(\partial \Omega )}\left(({\mathbf{\epsilon }}_h{\tau }_E,\mathrm{Curl}\varphi {)}_{L^2(E)}-(\mathrm{rot}{\mathbf{\epsilon }}_h\cdot {\tau }_E,\varphi {)}_{L^2(E)}\right).$$

$$\mathrm{Curl}\varphi =({\tau }_E\cdot \mathrm{Curl}\varphi ){\tau }_E+({\nu }_E\cdot \mathrm{Curl}\varphi ){\nu }_E=({\tau }_E\cdot \mathrm{Curl}\varphi ){\tau }_E-(\partial \varphi /\partial s){\nu }_E\text{on }E\in \mathcal{E}.$$

$$([{\mathbf{\epsilon }}_h{]}_E{\tau }_E,\mathrm{Curl}\varphi {)}_{L^2(E)}=({\tau }_E\cdot [{\mathbf{\epsilon }}_h{]}_E{\tau }_E,{\partial }_{{\nu }_E}\varphi {)}_{L^2(E)}+(\frac{\partial {\nu }_E\cdot [{\mathbf{\epsilon }}_h{]}_E{\tau }_E}{\partial s},\varphi {)}_{L^2(E)}.$$

$$(\epsilon (v),{\mathrm{Curl}}^2\varphi {)}_{L^2(\Omega )}=\sum _{E\in \mathcal{E}({\Gamma }_D)}\left((\frac{\partial u_D}{\partial s},\frac{\partial \varphi }{\partial {\nu }_E}{\tau }_E{)}_{L^2(E)}+(\frac{{\partial }^2{u}_D}{\partial s^2},\varphi {\nu }_E{)}_{L^2(E)}\right).$$

Full text

Numerical experiments verify the proven reliability and efficiency of the new a posteriori error estimator and illustrate the improved convergence rate in comparison to uniform mesh-refining. A higher convergence rates for piecewise affine data is observed in the $L^{2}$ stress error and reproduced in non-smooth situations by the adaptive mesh-refining strategy.

Keywords linear elasticity, mixed finite element method, a posteriori, explicit residual-based, error estimator, reliability, efficiency, symmetric stress finite elements, Arnold-Winther finite element

AMS subject classification 65N15, 65N30

1 Introduction

1.1 Overview

The design of a pointwise symmetric stress approximation $\sigma_{h}\in L^{2}(\Omega;\mathbb{S})$ with divergence in $L^{2}(\Omega;\mathbb{R}^{d})$, written $\sigma_{h}\in H(\operatorname{div},\Omega;\mathbb{S})$, has been a long-standing challenge [Arn02] and the first positive examples in [AW02] initiated what nowadays is called the finite element exterior calculus [AFW06]. The a posteriori error analysis of mixed finite element methods in elasticity started with [CD98] on PEERS [ABD84], where the asymmetric stress approximation $\gamma_{h}$ arises in the discretization as a Lagrange multiplier to enforce weakly the stress symmetry. This allows the treatment of the term $\mathbb{C}^{-1}\sigma_{h}+\gamma_{h}$ as an approximation of the (non symmetric) functional matrix $Du$ for the displacement field [CD98] with the arguments of [Alo96, Car97] developed for mixed finite element schemes for a Poisson model problem. Here and throughout, $\mathbb{C}$ denotes a fourth-order elasticity tensor with two Lamé constants $\lambda$ and $\mu$ and $\mathbb{C}^{-1}$ is its inverse. Mixed finite element methods appear attractive in the incompressible limit for they typically avoid the locking phenomenon [CEG11] as $\lambda\to\infty$.

For mixed finite element methods like the symmetric Arnold-Winther finite element schemes [AW02], the subtle term is the nonconforming residual: Given any piecewise polynomial $\sigma_{h}\in H(\operatorname{div},\Omega;\mathbb{S})$, compute an upper bound $\eta(\mathcal{T},\sigma_{h})$ of

[TABLE]

Despite general results in this direction [Car05, CH07, CPS16], this task had been addressed only by the computation of an approximation to the optimal $v$ with Green strain $\epsilon(v):=\operatorname{sym}Dv$ or of some skew-symmetric approximation $\gamma_{h}$ motivated from the first results in [CD98] on PEERS. In fact, any choice of a piecewise smooth and pointwise skew-symmetric $\gamma_{h}$ allows for an a posteriori error control of the symmetric stress error $\sigma-\sigma_{h}$ in [CG16]. Its efficiency, however, depends on the (unknown and uncontrolled) efficiency of the choice of $\gamma_{h}$ as an approximation to the skew-symmetric part $\gamma$ of $Du$.

This paper proposes the first reliable and efficient explicit residual-based a posteriori error estimator of the nonconforming residual with the typical contributions to $\eta(\mathcal{T},\sigma_{h})$ computed from the (known) Green strain approximation $\bm{\varepsilon}_{h}:=\mathbb{C}^{-1}\sigma_{h}$. Besides oscillations of the applied forces in the volume and along the Neumann boundary, there is a volume contribution $h_{T}^{2}\|\operatorname{rot}\operatorname{rot}\bm{\varepsilon}_{h}\|_{L^{2}(T)}$ for each triangle $T\in\mathcal{T}$ and an edge contribution with the jump $[\bm{\varepsilon}_{h}]_{E}$ across an interior edge $E$ with unit normal $\nu_{E}$, tangential unit vector $\tau_{E}$, and length $h_{E}$, namely

[TABLE]

and corresponding modification on the edges on the Dirichlet boundary with the (possibly inhomogeneous) Dirichlet data; cf. Remark 9 for some partial simplification of the last term displayed.

For the ease of this presentation, the analysis involves explicit calculations in two dimensions without any reference to the exterior calculus but with inhomogeneous Dirichlet and Neumann boundary data. The main result is reliability and efficiency to control the stress error robustly in the sense that the multiplicative generic constants hidden in the notation $\lesssim$ do neither depend on the (local or global) mesh-size nor on the parameter $\lambda>0$ but may depend on $\mu>0$ and on the shape regularity of the underlying triangulation $\mathcal{T}$ of the domain $\Omega$ into triangles through a lower bound of the minimal angle therein.

1.2 Linear elastic model problem

The elastic body $\Omega$ is a simply-connected bounded Lipschitz domain $\Omega\subset\mathbb{R}^{2}$ in the plane with a (connected) polygonal boundary $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$ split into parts. The displacement boundary $\Gamma_{D}$ is compact and of positive surface measure, while the traction boundary is the relative open complement $\Gamma_{N}=\partial\Omega\backslash\Gamma_{D}$ with outer unit normal vector $\nu$. Given $u_{D}\in H^{1}(\Omega;\mathbb{R}^{2})$, the volume force $f\in L^{2}(\Omega;\mathbb{R}^{2})$, and the applied surface traction $g\in L^{2}(\Gamma_{N};\mathbb{R}^{2})$, the linear elastic problem seeks a displacement $u\in H^{1}(\Omega;\mathbb{R}^{2})$ and a symmetric stress tensor $\sigma\in H(\operatorname{div},\Omega;\mathbb{S})$ with

[TABLE]

Throughout this paper, given the Lamé parameters $\lambda,\mu>0$ for isotropic linear elasticity, the positive definite fourth-order elasticity tensor $\mathbb{C}$ acts as $\mathbb{C}E:=2\mu\,E+\lambda\,\operatorname{tr}(E)\,1_{2\times 2}$ on any matrix $E\in\mathbb{S}$ with trace $\operatorname{tr}(E)$ and the $2\times 2$ unit matrix $1_{2\times 2}$. Note that $u_{D}$ acts in (1.1) only on $\Gamma_{D}$ and is an extension of the continuous function $u_{D}\in C(\Gamma_{D};\mathbb{R}^{2})$ also supposed to belong to the edgewise second order Sobolev space $H^{2}(\mathcal{E}(\Gamma_{D}))$ below to allow second derivatives with respect to the arc length along boundary edges.

More essential will be a discussion on the precise conditions on the Neumann data $g$ and its discrete approximation $g_{h}$ below for they belong to the essential boundary conditions in the mixed finite element method based on the dual formulation.

In addition to the set of homogeneous displacements $V$ and the aforementioned stress space $H(\operatorname{div},\Omega;\mathbb{S})$, namely,

[TABLE]

and with the exterior unit normal vector $\nu$ along $\partial\Omega$, the inhomogeneous stress space

[TABLE]

is defined with respect to the Neumann data $g\in L^{2}(\Gamma_{N})$ and, in particular, $\Sigma_{0}:=\Sigma(0)$ abbreviates the stress space with homogeneous Neumann boundary conditions.

Given data $u_{D},f,g$ as before, the dual weak formulation of (1.1) seeks $(\sigma,u)\in\Sigma(g)\times L^{2}(\Omega;\mathbb{R}^{2})$ with

[TABLE]

It is well known that the two formulations are equivalent and well posed in the sense that they allow for unique solutions in the above spaces and are actually slightly more regular according to the reduced elliptic regularity theory. The reader is refereed to textbooks on finite element methods [BS08, Bra07, BBF13] for proofs and further descriptions.

Throughout this paper, the model problem considers truly mixed boundary conditions with the hypothesis that both $\Gamma_{D}$ and $\Gamma_{N}$ have positive length. The remaining cases of a pure Neumann problem or a pure Dirichlet problem require standard modification and are immediately adopted. The presentation focusses on the case of isotropic linear elasticity with constant Lamé parameters $\lambda$ and $\mu$ for brevity and many results carry over to more general situations (cf. Remark 8 and 9 for instance).

1.3 Mixed finite element discretization

Let $\mathcal{T}$ denote a shape-regular triangulation of $\Omega$ into triangles (in the sense of Ciarlet [BS08]) with set of nodes $\mathcal{N}$, set of interior edges $\mathcal{E}(\Omega)$, set of Dirichlet edges $\mathcal{E}(\Gamma_{D})$ and set of Neumann edges $\mathcal{E}(\Gamma_{N})$. The triangulation is compatible with the boundary pieces $\Gamma_{D}$ and $\Gamma_{N}$ in that the boundary condition changes only at some vertex $\mathcal{N}$ and $\Gamma_{D}$ (resp. $\overline{\Gamma_{N}}$) is partitioned in $\mathcal{E}(\Gamma_{D})$ (resp. $\mathcal{E}(\Gamma_{N})$).

The piecewise polynomials (piecewise with respect to the triangulation $\mathcal{T}$) of total degree at most $k\in\mathbb{N}_{0}$ are denoted as $P_{k}(\mathcal{T})$, their vector- or matrix-valued versions as $P_{k}(\mathcal{T};\mathbb{R}^{2})$ or $P_{k}(\mathcal{T};\mathbb{R}^{2\times 2})$ etc. The subordinated Arnold-Winther finite element space $AW_{k}(\mathcal{T})$ of index $k\in\mathbb{N}$ [AW02] reads

[TABLE]

The Neumann boundary conditions are essential conditions and are traditionally implemented by some approximation $g_{h}$ to $g$ in the normal trace space

[TABLE]

(recall that $\nu$ is the exterior unit normal along the boundary). Given any $g_{h}\in G(\mathcal{T})$, the discrete stress approximations are sought in the non-void affine subspace

[TABLE]

of $AW_{k}(\mathcal{T})$ with test functions in the linear subspace $\Sigma(0,\mathcal{T}):=\Sigma_{0}\cap AW_{k}(\mathcal{T})$. Then there exists a unique discrete solution $\sigma_{h}\in\Sigma(g_{h},\mathcal{T})$ and $u_{h}\in V_{h}:=P_{k}(\mathcal{T};\mathbb{R}^{2})$ to

[TABLE]

The explicit design of a Fortin projection leads in [AW02] to quasi-optimal a priori error estimates for an exact solution $(\sigma,u)\in(\Sigma(g)\cap H^{k+2}(\Omega;\mathbb{S}))\times H^{k+2}(\Omega)$ to (1.1) and the approximate solution $(\sigma_{h},u_{h})$ to (1.3), namely (with the maximal mesh-size $h$)

[TABLE]

Another stable pair of different and mesh-depending norms in [CGS16] implies the $L^{2}$ best approximation of the stress error $\sigma-\sigma_{h}$ up to a generic multiplicative constant and data oscillations on $f$ under some extra condition (N) on the Neumann data approximation $g_{h}$ implied by the first and zero moment orthogonality assumption $g-g_{h}\perp P_{1}(\mathcal{E}(\Gamma_{N});\mathbb{R}^{2})$ ($\perp$ indicates orthogonality in $L^{2}(\Gamma_{N})$) met in all the numerical examples of this paper.

For simple benchmark examples with piecewise polynomial data $f$ and $g$, there is even a superconvergence phenomenon visible in numerical examples. The arguments of this paper allow a proof of fourth-order convergence of the $L^{2}$ stress error $\|\sigma-\sigma_{h}\|=o(h^{4})$ in the lowest-order Arnold-Winther method with $k=1$ for a smooth stress $\sigma\in H^{4}(\Omega;\mathbb{S})$ with $f=f_{h}\in P_{1}(\mathcal{T};\mathbb{R}^{2})$ and $g=g_{h}\in G(\mathcal{T})$. (In fact, once the data are not piecewise affine, the arising oscillation terms are only of at most third order and the aforementioned convergence estimates are sharp.)

This is stated as Theorem 17 in the appendix, because the a priori error analysis lies outside of the main focus of this work. It is surprising though that adaptive mesh-refining suggested with this paper recovers this higher convergence rate even for the inconsistent Neumann data in the Cook membrane benchmark example below.

1.4 Explicit residual-based a posteriori error estimator

The novel explicit residual-based error estimator for the discrete solution $(\sigma_{h},u_{h})$ to (1.3) depends only on the Green strain approximation $\mathbb{C}^{-1}\sigma_{h}$ and its piecewise derivatives and jumps across edges.

Given any edge $E$ of length $h_{E}$, let $\nu_{E}$ denote the unit normal vector (chosen with a fixed orientation such that it points outside along the boundary $\partial\Omega$ of $\Omega$) and let $\tau_{E}$ denote its tangential unit vector; by convention $\tau_{E}=(0,-1;1,0)\nu_{E}$ with the indicated asymmetric $2\times 2$ matrix. The tangential derivative $\tau_{E}\cdot\nabla\bullet$ along an edge (or boundary) is identified with the one-dimensional derivative $\partial\bullet/\partial s$ with respect to the arc-length parameter $s$. The jump $[v]_{E}$ of any piecewise continuous scalar, vector, or matrix $v$ across an interior edge $E=\partial T_{+}\cap\partial T_{-}$ shared by the two triangles $T_{+}$ and $T_{-}$ such that $\nu_{E}$ points outside $T_{+}$ along $E\subset\partial T_{+}$ reads

[TABLE]

The rotation acts on a vector field $\Phi$ (and row-wise on matrices) via $\operatorname{rot}\Phi:=\partial_{1}\Phi_{2}-\partial_{2}\Phi_{1}$ and $\operatorname{rot}_{NC}$ denotes its piecewise application.

Under the present notation and the throughout abbreviation $\bm{\varepsilon}_{h}:=\mathbb{C}^{-1}\sigma_{h}$, the explicit residual-based a posteriori error estimator reads

[TABLE]

for the oscillations $\operatorname{osc}(f,\mathcal{T})$ of the volume force and the oscillations $\operatorname{osc}(g-g_{h},\mathcal{E}(\Gamma_{N}))$ of the traction boundary condition, defined below.

Theorem 1** (reliability).**

There exists a mesh-size and $\lambda$ independent constant $C_{\textrm{rel}}$ (which may depend on $\mu$ and on the shape-regularity of the triangulation $\mathcal{T}$ through a global lower bound of the minimal angle therein) such that the exact (resp. discrete) stress $\sigma$ from (1.1) (resp. $\sigma_{h}$ from (1.3) with $g-g_{h}\perp P_{0}(\mathcal{E}(\Gamma_{N});\mathbb{R}^{2})$ and the error estimator (1.4) satisfy

[TABLE]

The a posteriori error estimator $\eta(\mathcal{T},\sigma_{h})$ already involves two data oscillation terms $\operatorname{osc}(f,\mathcal{T})$ and $\operatorname{osc}(g-g_{h},\mathcal{E}(\Gamma_{N}))$ defined as the square roots of the respective terms in

[TABLE]

For any edge $E$ and a degree $m\geq k+2$, let $\Pi_{m,E}:L^{2}(E)\to P_{m}(E)$ denote the $L^{2}$ projection onto polynomials of degree at most $m$. For any $E\in\mathcal{E}(\Gamma_{D})$ define the two Dirichlet data oscillation terms

[TABLE]

Their sum defines the overall Dirichlet data approximation $\operatorname{osc}(u_{D},\mathcal{E}(\Gamma_{D}))$ as the square root of

[TABLE]

The analysis of Section 3 is local and states for each of the five local residuals an upper bound related to the error in a neighborhood. The global efficiency is displayed as follows.

Theorem 2** (efficiency).**

There exists a mesh-size and $\lambda,\mu$ independent constant $C_{\textrm{eff}}$ (which may depend on the shape-regularity of the triangulation $\mathcal{T}$ through a global lower bound of the minimal angle therein) such that the exact (resp. discrete) stress $\sigma$ from (1.1) (resp. $\sigma_{h}$ from (1.3) with $g-g_{h}\perp P_{0}(\mathcal{E}(\Gamma_{N});\mathbb{R}^{2})$ and the error estimator (1.4) satisfy

[TABLE]

1.5 Outline of the paper

The remaining parts of this paper provide a mathematical proof of Theorem 1 and Theorem 2 and numerical evidence in computational experiments on the novel a posteriori error estimation and its robustness as well as on associated mesh-refining algorithms.

The proof of the reliability of Theorem 1 in Section 2 adopts arguments of [CD98, CG16] and carries out two integration by parts on each triangle plus one-dimensional integration by parts along all edges. The resulting terms are in fact locally efficient in Section 3 with little generalizations of the bubble-function methodology due to Verfürth [Ver96]. The five lemmas of Section 3 give slightly sharper results and in total imply Theorem 2.

The point in Theorem 1 and 2 is that the universal constants $C_{\textrm{rel}}$ and $C_{\textrm{eff}}$ may depend on the Lamé parameter $\mu$ but are independent of the critical Lamé parameter $\lambda$ as supported by the benchmark examples of the concluding Section 4. Adaptive mesh-refining proves to be highly effective with the novel a posteriori error estimator even for incompatible Neumann data. Four benchmark examples with the Poisson ratio $\nu=0.3$ or $0.4999$ provide numerical evidence of the robustness of the reliabile and efficient a posteriori error estimation and for the fourth-order convergence of Theorem 17.

1.6 Comments on general notation

Standard notation on Lebesgue and Sobolev spaces and norms is adopted throughout this paper and, for brevity, $\|\cdot\|:=\|\cdot\|_{L^{2}(\Omega)}$ denotes the $L^{2}$ norm. The piecewise action of a differential operator is denoted with a subindex ${NC}$, e.g., $\nabla_{NC}$ denotes the piecewise gradient $(\nabla_{NC}\bullet)|_{T}:=\nabla(\bullet|_{T})$ for all $T\in\mathcal{T}$. Sobolev functions are usually defined on open sets and the notation $W^{m,p}(T)$ (resp. $W^{m,p}(\mathcal{T})$) substitutes $W^{m,p}(\operatorname{int}(T))$ for a (compact) triangle $T$ and its interior $\operatorname{int}(T)$ (resp. $W^{m,p}(\operatorname{int}(\mathcal{T}))$) and their vector and matrix versions.

For a differentiable function $\phi$, $\operatorname{Curl}\phi:=(-\partial_{2}\phi,\partial_{1}\phi)$ is the rotated gradient; for a two-dimensional vector field $\Phi$, $\operatorname{Curl}\Phi$ is the $2\times 2$ matrix-valued rotated gradient

[TABLE]

(The signs are not uniquely determined in the literature and some care is required.)

The colon denotes the scalar product $A:B:=\sum_{\alpha,\beta=1,2}A_{\alpha,\beta}B_{\alpha,\beta}$ of $2\times 2$ matrices $A,B$. The inequality $A\lesssim B$ between two terms $A$ and $B$ abbreviates $A\leq C\,B$ with some multiplicative generic constant $C$, which is independent of the mesh-size and independent of the one Lamé parameter $\lambda\geq 0$ but may depend on the other $\mu>0$ and may depend on the shape-regularity of the underlying triangulation $\mathcal{T}$.

2 Proof of reliability

This section is devoted to the proof of Theorem 1 based on a Helmholtz decomposition on [CD98] with two parts as decomposed in Theorem 5 below. The critical part is the $L^{2}$ product of $\mathbb{C}^{-1}(\sigma-\sigma_{h})$ times the $\operatorname{Curl}$ of an unknown function $\operatorname{Curl}\beta$. The observation from [CG16] is that one can find an Argyris finite element approximation $\beta_{h}$ to $\beta\in H^{2}(\Omega)$ such that the continuous function $\phi:=\beta-\beta_{h}\in H^{2}(\Omega)$ vanishes at all vertices $\mathcal{N}$ of the triangulation. Two integration by parts on each triangle plus one-dimensional integration by parts along the edges $\mathcal{E}$ of the triangulation eventually lead to a key identity.

Lemma 3** (representation formula).**

Any function $\bm{\varepsilon}_{h}\in H^{2}(\mathcal{T};\mathbb{S})$ (i.e. $\bm{\varepsilon}_{h}$ is piecewise in $H^{2}$ with values in $\mathbb{S}$) and any $\phi\in H^{2}(\Omega)$ with $\phi(z)=0$ at all vertices $z\in\mathcal{N}$ in the regular triangulation $\mathcal{T}$ satisfy

[TABLE]

The subsequent integration by parts formula is utilized frequently throughout this paper for $\phi\in H^{1}(\Omega;\mathbb{R}^{2})$ and $\Psi\in H^{1}(\Omega;\mathbb{R}^{2\times 2})$

[TABLE]

Any differentiable (scalar) function $\varphi$, satisfies the elementary relations

[TABLE]

Proof.

Integrate by parts twice on each triangle and rearrange the remaining boundary terms to deduce (with the abbreviation $\operatorname{rot}_{NC}\operatorname{rot}_{NC}\equiv\operatorname{rot}_{NC}^{2}$)

[TABLE]

The term $([\bm{\varepsilon}_{h}]_{E}\tau_{E},\operatorname{Curl}\phi)_{L^{2}(E)}$ in the above sum is split into orthogonal components

[TABLE]

Since $\phi$ vanishes at the vertices, an integration by parts along each interior edge $E$ for the last term shows $([\bm{\varepsilon}_{h}]_{E}\tau_{E},(\partial\phi/\partial s)\nu_{E})_{L^{2}(E)}=-(\partial[\bm{\varepsilon}_{h}]_{E}\tau_{E}/\partial s,\phi\nu_{E})_{L^{2}(E)}$. This proves

[TABLE]

The same formula holds for any boundary edge $E$ when $[\bm{\varepsilon}_{h}]_{E}$ is replaced by $\bm{\varepsilon}_{h}$. The combination of the latter identities with the first displayed formula of this proof verifies the asserted representation formula. ∎

The contribution of $\epsilon(u)=\mathbb{C}^{-1}\sigma$ times the $\operatorname{Curl}^{2}\phi\in L^{2}(\Omega;\mathbb{S})$ exclusively leads to boundary terms. Throughout this paper, suppose that the Dirichlet data $u_{D}$ satisfies $u_{D}\in C(\Gamma_{D})\cap H^{2}(\mathcal{E}(\Gamma_{D})$ in the sense that $u_{D}$ is globally continuous with $u_{D}|_{E}\in H^{2}(E;\mathbb{R}^{2})$ for all $E\in\mathcal{E}(\Gamma_{D})$.

Lemma 4** (boundary terms).**

Any Sobolev function $v\in H^{1}(\Omega;\mathbb{R}^{2})$ with boundary values $u_{D}\in C(\Gamma_{D})\cap H^{2}(\mathcal{E}(\Gamma_{D})$ on $\Gamma_{D}$ and any $\phi\in H^{2}(\Omega)$ with $\phi=\partial\phi/\partial\nu=0$ along $\Gamma_{N}$ with $\phi(z)=0$ for any vertex $z$ of $\Gamma_{D}$ in its relative interior satisfy

[TABLE]

Proof.

A density argument shows that it suffices to prove this identity for smooth functions $v$ and $\phi$, when integration by parts arguments show that the left-hand side is equal to

[TABLE]

The equality follows from an orthogonal split $\operatorname{Curl}\phi=(\tau\cdot\operatorname{Curl}\phi)\tau+(\nu\cdot\operatorname{Curl}\phi)\nu$ into the normal and tangential directions of $\nu$ and $\tau$ along the boundary $\partial\Omega$ followed by an integration by parts along $\partial\Omega$ with $\phi(z)=0$ for vertices $z$ in $\Gamma_{D}$ with a jump of the normal unit vector. The substitution of the boundary conditions concludes the proof. ∎

The consequence of the previous two lemmas is a representation formula for the error times a typical function $\operatorname{Curl}^{2}\phi$. To understand why the contributions on the Neumann boundary of $\phi$ and $\nabla\phi$ disappear along $\Gamma_{N}$, some details on the Helmholtz decomposition are recalled from the literature. For this, let $\Gamma_{0},\dots,\Gamma_{J}$ denote the compact connectivity components of $\overline{\Gamma_{N}}$.

Theorem 5** (Helmholtz decomposition [CD98, Lemma 3.2]).**

Given any $\sigma-\sigma_{h}\in L^{2}(\Omega;\mathbb{S})$, there exists $\alpha\in V$, constant vectors $c_{0},\dots,c_{J}\in\mathbb{R}^{2}$ with $c_{0}=0$ and $\beta\in H^{2}(\Omega)$ with $\int_{\Omega}\beta\,dx=0$ and $\operatorname{Curl}\beta=c_{j}$ on $\Gamma_{j}\subseteq\Gamma_{N}$ for all $j=0,\dots,J$ such that

[TABLE]

The second ingredient is an approximation $\beta_{h}$ of $\beta$ from the Helmholtz decomposition in Theorem 5 based on the Argyris finite element functions $A(\mathcal{T})\subset C^{1}(\Omega)\cap P_{5}(\mathcal{T})$ [BS08, Bra07, Cia78]. The local mesh-size $h_{\mathcal{T}}\in P_{0}(\mathcal{T})$ in the triangulation $\mathcal{T}$ is defined as its diameter $h_{\mathcal{T}}|_{T}:=h_{T}$ on each triangle $T\in\mathcal{T}$.

Lemma 6** (quasiinterpolation).**

Given any $\beta$ as in Theorem 5 there exists some $\beta_{h}\in A(\mathcal{T})$ such that $\phi:=\beta-\beta_{h}\in H^{2}(\Omega)$ vanishes at any vertex $z\in\mathcal{N}$ of the triangulation, $\phi$ and its gradient $\nabla\phi$ vanish on $\Gamma_{N}$, and the local approximation and stability property holds in the sense that

[TABLE]

Proof.

This has been (partly) utilized in [CG16] and also follows from [GS02]. ∎

The combination of all aforementioned arguments leads to the following estimate as an answer to the question of Subsection 1.1 in terms of directional derivatives of $\bm{\varepsilon}_{h}:=\mathbb{C}^{-1}\sigma_{h}$. Recall the definition of $\eta(\mathcal{T},\sigma_{h})$ from (1.4).

Theorem 7** (key result).**

Let $\sigma\in H(\operatorname{div},\Omega;\mathbb{S})$ solve (1.1) and let $\sigma_{h}\in AW_{k}(\mathcal{T})$ solve (1.3). Given $\beta$ from Theorem 5 and its quasiinterpolation $\beta_{h}$ from Lemma 6, the difference $\phi:=\beta-\beta_{h}$ satisfies

[TABLE]

Proof.

Lemma 3 and Lemma 4 lead to a formula for $(\bm{\varepsilon}_{h},\operatorname{Curl}^{2}\phi)_{L^{2}(\Omega)}$, $\bm{\varepsilon}_{h}:=\mathbb{C}^{-1}\sigma_{h}$, in which all the contributions for $E\in\mathcal{E}(\Gamma_{N})$ with $\phi$ and $\nabla\phi$ vanish along $\Gamma_{N}$. The remaining formula reads

[TABLE]

The proof concludes with Cauchy-Schwarz inequalities, trace inequalities, and the approximation estimates of Lemma 6. The remaining details are nowadays standard arguments in the a posteriori error analysis of nonconforming and mixed finite element methods and hence are omitted. ∎

Before the proof of Theorem 1 concludes this section, three remarks and one lemma are in order.

Remark 8** (nonconstant coefficients).**

The main parts of the reliability analysis of this section hold for rather general material tensors $\mathbb{C}$ as long as $\bm{\varepsilon}_{h}:=\mathbb{C}^{-1}\sigma_{h}$ allows for the existence of the traces and the derivatives in the error estimator (1.4) in the respective $L^{2}$ spaces. For instance, if $\lambda$ and $\mu$ are piecewise smooth with respect to the underlying triangulation $\mathcal{T}$.

Remark 9** (constant coefficients).**

The overall assumption of constant Lamé parameters $\lambda$ and $\mu$ allows a simplification in the error estimator (1.4). It suffices to have $\mu$ globally continuous and $\mu$ and $\lambda$ piecewise smooth to guarantee

[TABLE]

(The proof utilizes the structure of $\mathbb{C}$ and $\mathbb{C}^{-1}$ with $\mathbb{C}^{-1}E=\frac{1}{2\mu}(E-\frac{\lambda}{2(\lambda+\mu)}\operatorname{tr}(E)1_{2\times 2})$ for any $E\in\mathbb{S}$ as a linear combination of the identity and some scalar multiple of the $2\times 2$ unit matrix. The terms with the identity lead to $1/(2\mu)$ times the jump $[\sigma_{h}]_{E}\nu_{E}=0$ of the $H(\operatorname{div})$ conforming stress approximations. The jump terms with the unit matrix (even with jumps of $\lambda$) are multiplied with the orthogonal unit vectors $\nu_{E}$ and $\tau_{E}$ and so vanish as well.)

Remark 10** (Related work).**

Although the work [HHX11] concerns a different problem (bending of a plate of fourth order) with a different discretisation (even nonconforming in $H(\operatorname{div})$), some technical parts of that paper are related to those of this by a rotation of the underlying coordinate system and the substitution of $\operatorname{div}\operatorname{div}$ by $\operatorname{rot}\operatorname{rot}$ etc. Another Helmholtz decomposition also allows for a discrete version and thereby enables a proof of optimal convergence of an adaptive algorithm with arguments from [CFPP14, CR16].

A technical detail related to the robustness in $\lambda\to\infty$ is a well known lemma that controls the trace of a matrix $E\in\mathbb{R}^{2\times 2}$ by its deviatoric part $\operatorname{dev}E:=E-\operatorname{tr}(E)/2\,1_{2\times 2}$ and its divergence measured in the dual $V^{*}\subset H^{-1}(\Omega;\mathbb{R}^{2})$ of $V$, namely

[TABLE]

Lemma 11** (tr-dev-div).**

Let $\Sigma_{0}$ be a closed subspace of $H(\operatorname{div},\Omega;\mathbb{R}^{2\times 2})$, which does not contain the constant tensor $1_{2\times 2}$. Then any $\tau\in\Sigma_{0}$ satisfies

[TABLE]

Proof.

There are several variants of the tr-dev-div lemma known in the literature [BF91, Proposition 7 in Section IV.3]. The version in [CD98, Theorem 4.1] explicitly displays a version with $\|\operatorname{div}\tau\|$ replacing $\|\operatorname{div}\tau\|_{-1}$. Since its proof is immediately adopted to prove the asserted version, further details are omitted. ∎

The remaining part of this section outlines why Theorem 1 follows from Theorem 7 with the arguments from [CD98, CG16]. The energy norms of $v\in V$ and $\tau\in H(\operatorname{div},\Omega;\mathbb{S})$ read

[TABLE]

The remaining residual

[TABLE]

with its dual norm

[TABLE]

allow the following error decomposition [CG16, Theorem 3.1], with the weighted $L^{2}$ norms $\|\bullet\|_{\mathbb{C}^{-1}}:=\|\mathbb{C}^{-1/2}\bullet\|_{L^{2}(\Omega)}$ and $\|\bullet\|_{\mathbb{C}}:=\|\mathbb{C}^{1/2}\bullet\|_{L^{2}(\Omega)}$,

[TABLE]

In fact, it is shown in the proof of [CG16, Theorem 3.1] that $\alpha\in V$ and $\beta\in H^{2}(\Omega)$ from the Helmholtz decomposition of the error $\sigma-\sigma_{h}$ in Theorem 5 are orthogonal with respect to the $L^{2}$ scalar product weighted with $\mathbb{C}$ and that the last term in (2.2) equals $\|\operatorname{Curl}^{2}\beta\|_{\mathbb{C}^{-1}}^{2}=(\mathbb{C}^{-1}(\sigma-\sigma_{h}),\operatorname{Curl}^{2}\phi)_{L^{2}(\Omega)}$ and hence is controlled in the key estimate of Theorem 7.

Lemma 11 applies to $\Sigma_{0}$ as the subspace of all $\tau\in H(\operatorname{div},\Omega;\mathbb{S})$ with homogeneous Neumann data $\tau\nu=0$ along $\Gamma_{N}$. Since $\tau=\operatorname{Curl}^{2}\beta$ is divergence free and since $\tau\nu=-\partial\operatorname{Curl}\beta/\partial s$ along $\Gamma_{N}$ (owing to the aforementioned elementary relations and the convention that the first $\operatorname{Curl}$ acts row-wise on $\operatorname{Curl}\beta$), where $\operatorname{Curl}\beta$ in Theorem 5 is piecewise constant, it follows that $\tau\in\Sigma_{0}$. On the other hand $1_{2\times 2}\notin\Sigma_{0}$ because $\Gamma_{N}\neq\emptyset$. Consequently, Lemma 11 implies $\|\operatorname{Curl}^{2}\beta\|\lesssim\|\operatorname{dev}\operatorname{Curl}^{2}\beta\|$. This and elementary calculations with $\mathbb{C}^{-1}$ lead to

[TABLE]

The combination with the estimate resulting from Theorem 7 proves

[TABLE]

This implies $\|\operatorname{Curl}^{2}\beta\|_{\mathbb{C}^{-1}}\lesssim\eta(\mathcal{T},\sigma_{h})$ and leads in (2.2) to

[TABLE]

The remaining term is the estimate of the dual norm $\lvert\!\lvert\!\lvert\operatorname{Res}\rvert\!\rvert\!\rvert_{*}$ of the residual which is done, e.g., in [CG16, Lemma 3.3] (under the assumption $g-g_{h}\perp P_{0}(\mathcal{E}(\Gamma_{N}))$)

[TABLE]

This and (2.3) imply

[TABLE]

For any test function $v\in V$ with $|v|_{H^{1}(\Omega)}=1$, $\int_{\Omega}(\sigma-\sigma_{h}):Dv\,dx=\operatorname{Res}(v)$ and so

[TABLE]

(In the second last step one utilizes that $2\mu\,E:E\leq E:\mathbb{C}E$ for all $E\in\mathbb{S}$.) The combination of Lemma 11 for $\tau=\sigma-\sigma_{h}$ with the previous displayed estimates concludes the proof of $\|\sigma-\sigma_{h}\|\lesssim\eta(\mathcal{T},\sigma_{h})$. There exist several appropriate choices of $\Sigma_{0}\subset H(\operatorname{div},\Omega;\mathbb{S})$ in this last step. Recall that $\Gamma_{N}$ is the union of connectivity components and so pick one edge $E_{0}$ in this polygon and consider $\Sigma_{0}:=\{\tau\in H(\operatorname{div},\Omega;\mathbb{S}):\int_{E_{0}}\tau\nu\,ds=0\}$ with $1_{2\times 2}\notin\Sigma_{0}$. This choice of $E_{0}$ and so $\Sigma_{0}$ depend only on $\Gamma_{N}$ (independent of $\mathcal{T}$). Since $g-g_{h}=(\sigma-\sigma_{h})\nu$ along $E_{0}$ has (piecewise on $\mathcal{E}(E_{0})$, whence in total) an integral mean zero, Lemma 11 indeed applies to $\tau=\sigma-\sigma_{h}\in\Sigma_{0}$. ∎

3 Local efficiency analysis

The local efficiency follows with the bubble-function technique for $C^{1}$ finite elements [Ver96, Sec 3.7]. This section focusses on a constant $\mathbb{C}$ for linear isotropic elasticity with constant Lamé parameters $\lambda$ and $\mu$ such that $\bm{\varepsilon}_{h}:=\mathbb{C}^{-1}\sigma_{h}\in P_{k+2}(\mathcal{T})$ for some $\sigma_{h}\in AW_{k}(\mathcal{T})$ is a polynomial of degree at most $k+2$. Appart from this, the Lamé parameters do not further arise in this section.

The moderate point of departure is the volume term for each triangle $T\in\mathcal{T}$ with barycentric coordinates $\lambda_{1},\lambda_{2},\lambda_{3}\in P_{1}(T)$ and their product, the cubic volume bubble function, $b_{T}:=27\,\lambda_{1}\lambda_{2}\lambda_{3}\in W^{1,\infty}_{0}(T)$ plus its square $b_{T}^{2}\in W^{2,\infty}_{0}(T)$ with $0\leq b_{T}^{2}\leq 1$, $\|b_{T}\|_{L^{2}(T)}\lesssim 1$, and $|b_{T}|_{H^{2}(T)}\lesssim h_{T}^{-2}$ etc.

Lemma 12** (efficiency of volume residual).**

Any $v\in H^{1}(T;\mathbb{R}^{2})$, $T\in\mathcal{T}$, satisfies

[TABLE]

Proof.

An inverse estimate for the polynomial $\operatorname{rot}\operatorname{rot}\bm{\varepsilon}_{h}\equiv\operatorname{rot}^{2}\bm{\varepsilon}_{h}$ implies the estimate

[TABLE]

Lemma 3 with $\phi=b_{T}\operatorname{rot}^{2}\bm{\varepsilon}_{h}$ and $(\varepsilon(v),\operatorname{Curl}^{2}\phi)_{L^{2}(T)}=0$ lead to

[TABLE]

This and the inverse estimate $\|\operatorname{Curl}^{2}(b_{T}\operatorname{rot}^{2}\bm{\varepsilon}_{h})\|_{L^{2}(T)}\lesssim h_{T}^{-2}\|b_{T}\operatorname{rot}^{2}\bm{\varepsilon}_{h}\|_{L^{2}(T)}$ imply

[TABLE]

This concludes the proof. ∎

The localization of the first edge residual involves the piecewise quadratic edge-bubble function $b_{E}$ with support $T_{+}\cup T_{-}$ for an interior edge $E=\partial T_{+}\cap\partial T_{-}$ shared by the two triangles $T_{+}$ and $T_{-}$ with edge-patch $\omega_{E}:=\operatorname{(}T_{+}\cup T_{-})$. With an appropriate scaling $b_{E}|_{T}=4\lambda_{1}\lambda_{2}$ for the two barycentric coordinates $\lambda_{1},\lambda_{2}$ on $T\in\{T_{+},T_{-}\}$ associated with the two vertices of $E$. Then $b_{E}\in W^{1,\infty}(\omega_{E})$ and $b_{E}^{2}\in W^{1,\infty}(\omega_{E})$ satisfy $0\leq b_{E}^{2}\leq b_{E}\leq 1$ and $|b_{E}|_{H^{1}(E)}\lesssim h_{E}^{-1}$ etc.

The remaining technical detail is an extension of functions on the edge $E$ to $\omega_{E}$. Throughout this sections those functions are polynomials and given $\rho_{E}\in P_{m}(E)$, their coefficients (in some fixed basis) already define an algebraic object that is a natural extension $\rho\in P_{m}(\hat{E})$ along the straight line $\hat{E}:=\operatorname{mid}(E)+\mathbb{R}\,\tau_{E}$ that extends $E$ with midpoint $\operatorname{mid}(E)$ and tangential unit vector $\tau_{E}$. This and

[TABLE]

define a linear extension operator $P_{E}:P_{m}(E)\to C^{\infty}(\mathbb{R}^{2})$ with $P_{E}(\rho_{E})=\rho_{E}$ on $E$ for any $\rho_{E}\in P_{m}(E)$, which is constant in the normal direction, $\nabla P_{E}(\rho_{E})\cdot\nu_{E}\equiv 0$. This design is different from that in [Ver96].

Lemma 13** (efficiency of first interior edge residual).**

Any $v\in H^{1}(\omega_{E};\mathbb{R}^{2})$, $E\in\mathcal{E}(\Omega)$, satisfies

[TABLE]

Proof.

Since $\tau_{E}\cdot[\bm{\varepsilon}_{h}]_{E}\tau_{E}\in P_{k+2}(E)$ is a polynomial, the above extension $P_{E}(\tau_{E}\cdot[\bm{\varepsilon}_{h}]_{E}\tau_{E})$ and the function $b\in W^{2,\infty}_{0}(\omega_{E})$ with

[TABLE]

define some function $\phi:=b\,P_{E}(\tau_{E}\cdot[\bm{\varepsilon}_{h}]_{E}\tau_{E})$. Since $b=0$ and $\nabla b_{E}\cdot\nu_{E}=b_{E}^{2}$ along $E$, the test function $\phi\in H^{2}_{0}(\omega_{E})\subset H_{0}^{2}(\Omega)$ leads in Lemma 3 to

[TABLE]

Since $\partial_{\nu_{E}}\phi=b_{E}^{2}\,\tau_{E}\cdot[\bm{\varepsilon}_{h}]_{E}\tau_{E}$ on $E$ and $\varepsilon(v)\perp\operatorname{Curl}^{2}\phi$, an inverse estimate shows

[TABLE]

At the heart of the bubble-function methodology are inverse and trace inequalities that allow for appropriate scaling properties [Ver96] under the overall assumption of shape-regularity. In the present case, one power of $h_{E}\approx h_{T_{\pm}}$ is hidden in the function $b$ and

[TABLE]

The combination with the previous estimate results in

[TABLE]

This and Lemma 12 conclude the proof. ∎

For any boundary edge $E\in\mathcal{E}(\Omega)$, the edge-bubble function $b_{E}=4\lambda_{1}\lambda_{2}\in W^{1,\infty}(\omega_{E})$ for the two barycentric coordinates $\lambda_{1},\lambda_{2}$ associated with the two vertices of $E$ and $b_{E}$ vanishes on the remaining sides $\partial\omega_{E}\setminus E$ of the aligned triangle $\overline{\omega_{E}}$. The Dirichlet data $u_{D}$ allows for some polynomial approximation $\Pi_{m,E}u_{D}\in P_{m}(E)$ of a maximal degree bounded by $m\geq k+2$; recall the definition of $\operatorname{osc}_{I}(u_{D},E)$ from (1.5).

Lemma 14** (efficiency of first boundary edge residual).**

Any function $v\in H^{1}(\omega_{E};\mathbb{R}^{2})$ with $v|_{E}=u_{D}|_{E}$ along $E\in\mathcal{E}(\Gamma_{D})$ satisfies

[TABLE]

Proof.

Since $\tau_{E}\cdot\bm{\varepsilon}_{h}\tau_{E}$ is a polynomial of degree at most $k+2\leq m$ along the exterior edge $E$, the residual $\tau_{E}\cdot(\bm{\varepsilon}_{h}\tau_{E}-\partial u_{D}/\partial s)$ is well approximated by its $L^{2}$ projection $\rho_{E}:=(\tau_{E}\cdot(\bm{\varepsilon}_{h}\tau_{E}-\Pi_{m,E}\partial u_{D}/\partial s))$ onto $P_{m}(E)$. The Pythagoras theorem based on the $L^{2}$ orthogonality reads

[TABLE]

and it remains to bound $h_{E}^{1/2}\|\rho_{E}\|_{L^{2}(E)}$ by the right-hand side of the claimed inequality. The extension $P_{E}\rho_{E}\in C^{\infty}(\mathbb{R}^{2})$ and the function $b$ from (3.1) lead to an admissible test function $\phi:=bP_{E}\rho_{E}\in W^{1,\infty}_{0}(\omega_{E})$. Two successive integration by parts as in Lemma 3 show

[TABLE]

This and Lemma 3 lead to

[TABLE]

Since $\partial_{\nu_{E}}\phi=b_{E}^{2}\rho_{E}$ along $E$ and $\rho_{E}$ is the $L^{2}$ projection of $\tau_{E}\cdot(\bm{\varepsilon}_{h}\tau_{E}-\partial u_{D}/\partial s)$, the left-hand side equals $\|b_{E}\rho_{E}\|_{L^{2}(E)}^{2}-((1-\Pi_{m,E})\partial u_{D}/\partial s,b_{E}^{2}\rho_{E})_{L^{2}(E)}$. The scaling argument which leads to (3.2) shows that the left-hand side of (3.2) is $\lesssim\|\rho_{E}\|_{L^{2}(E)}$. The combination with the previously displayed identity leads to

[TABLE]

This and Lemma 12 conclude the proof. ∎

The edge-bubble functions for the second edge residuals are defined slightly differently to ensure some vanishing normal derivative.

Lemma 15** (efficiency of second interior edge residual).**

Any $v\in H^{1}(\omega_{E};\mathbb{R}^{2})$, $E\in\mathcal{E}(\Omega)$, satisfies

[TABLE]

Proof.

There are many ways to define an edge-bubble function for this situation and one may first select a maximal open ball $B(x_{E},2r_{E})\subset\omega_{E}$ around a point $x_{E}\in E$ with maximal radius $2r_{E}$, which is entirely included in $\omega_{E}$. The characteristic function $\chi_{B(x_{E},r_{E})}$ of the smaller ball $B(x_{E},r_{E})$ may be regularized with a standard mollification $\eta_{r_{E}}$ to define the smooth function $b:=\chi_{B(x_{E},r_{E})}*\eta_{r_{E}}\in C_{c}^{\infty}(\Omega_{E})$ with values in $[0,1]$ and with $\nabla b\cdot\nu_{E}=0$ along $E$. The polynomial $\rho_{E}:=\tau_{E}\cdot([\operatorname{rot}_{NC}\bm{\varepsilon}_{h}]_{E}-\partial[\bm{\varepsilon}_{h}]_{E}/\partial s\,\nu_{E})$ and its extension $P_{E}\rho_{E}$ define the test function $\phi:=bP_{E}\rho_{E}\in C_{0}^{\infty}(\omega_{E})$ in Lemma 3. The representation formula and $(\varepsilon(v),\operatorname{Curl}^{2}\phi)_{L^{2}(\omega_{E})}=0$ lead to

[TABLE]

The inverse inequality $\|\rho_{E}\|_{L^{2}(E)}\|\leq\|b^{1/2}\rho_{E}\|_{L^{2}(E)}$, Cauchy inequalities, and the right scaling properties of $\phi$ lead to

[TABLE]

This and Lemma 12 conclude the proof. ∎

The efficiency of the last edge contribution involves the second Dirichlet data oscillation $\operatorname{osc}_{II}(u_{D},E)$ from (1.6).

Lemma 16** (efficiency of second boundary edge residual).**

Any function $v\in H^{1}(\omega_{E};\mathbb{R}^{2})$ with $v|_{E}=u_{D}|_{E}$ along $E\in\mathcal{E}(\Gamma_{D})$ satisfies

[TABLE]

Proof.

Select a maximal open ball $B(x_{E},2r_{E})\cap\Omega\subset\omega_{E}$ around a point $x_{E}\in E$ with maximal radius $2r_{E}$ such that $B(x_{E},2r_{E})\cap\omega_{E}$ is a half ball. The regularization $b:=\chi_{B(x_{E},r_{E})}*\eta_{r_{E}}\in C_{c}^{\infty}(\mathbb{R}^{2})$ of the characteristic function $\chi_{B(x_{E},r_{E})}$ attains values in $[0,1]$ and a positive integral mean $h_{E}^{-1}\int_{E}b\,ds\approx 1$ along $E$ (depending only on the shape regularity of $\mathcal{T}$); $b$ vanishes on $\partial\omega_{E}\setminus E$ and its normal derivative $\nabla b\cdot\nu=0$ vanishes along the entire boundary $\partial\omega_{E}$.

The Pythagoras theorem $\|\rho\|_{L^{2}(E)}^{2}=\|\rho_{E}\|_{L^{2}(E)}^{2}+h_{E}^{-3}\operatorname{osc}_{II}^{2}(u_{D},E)$ for the residual $\rho:=\tau_{E}\cdot\operatorname{rot}\bm{\varepsilon}_{h}-\nu_{E}\cdot(\frac{\partial\bm{\varepsilon}_{h}\tau_{E}}{\partial s}+\frac{\partial^{2}u_{D}}{\partial s^{2}})$ and its $L^{2}$ projection $\rho_{E}:=\Pi_{m,E}\rho$ onto $P_{m}(E)$ reduces the proof to the estimation of $\|\rho_{E}\|_{L^{2}(E)}$. The normal derivative of $\phi:=b\,P_{E}\rho_{E}\in C^{\infty}(\overline{\omega_{E}})$ vanishes along the boundary $\partial\omega_{E}$ and Lemma 3 shows

[TABLE]

The arguments in Lemma 4 show $(\partial^{2}u_{D}/\partial s^{2},b\rho_{E}\,\nu_{E})_{L^{2}(E)}=(\varepsilon(v),\operatorname{Curl}^{2}\phi)_{L^{2}(\omega_{E})}$. The combination of the two identities leads to a formula for $(\rho,b\rho_{E})_{L^{2}(E)}$. Since $\rho-\rho_{E}$ is controlled in $\operatorname{osc}_{II}^{2}(u_{D},E)$, this and an inverse inequality in the beginning result in

[TABLE]

The scaling properties of $\phi$ and its derivatives are as in the proof of the previous lemma. With Lemma 12 in the end, this concludes the proof. ∎

4 Numerical examples

This section is devoted to numerical experiments for four different domains to demonstrate robustness in the reliability and efficiency of the a posteriori error estimator $\eta(\mathcal{T}_{\ell},\sigma_{\ell})$. The implementation follows [CGRT08, CEG11, CG16] for $k=1$ with Lamé parameters $\lambda$ and $\mu$ from $\lambda=E\nu/((1+\nu)(1-2\nu))$ and $\mu={E}/(2(1+\nu))$ for a Young’s modulus $E=10^{5}$ and various Poisson ratios $\nu=0.3$ and $\nu=0.4999$.

4.1 Academic example

The model problem (1.1) on the unit square $\Omega=(0,1)^{2}$ with homogeneous Dirichlet boundary conditions and the right-hand side $f=(f_{1},f_{2})$,

[TABLE]

allows the smooth exact solution

[TABLE]

Note that $f$ depends only on the Lamé parameter $\mu$ and not on the critical Lamé parameter $\lambda$. For uniform mesh refinement, Figure 1 displays the robust third-order convergence of the a posteriori error estimator $\eta(\mathcal{T}_{\ell},\sigma_{\ell})$ as well as the Arnold-Winther finite element stress error. The convergence is robust in the Poisson ratio $\nu\to 1/2$ and the a posteriori error estimator proves to be reliable and efficient. In this example, the oscillations $\operatorname{osc}(f,\mathcal{T}_{\ell})$ dominate the a posteriori error estimator.

This typical observation motivates numerical examples with $f\equiv 0$ in the sequel.

4.2 Circular inclusion

The second benchmark example from the literature models a rigid circular inclusion in an infinite plate for the domain $\Omega$ with rather mixed boundary conditions indicated with mechanical symbols in Figure 2. The exact solution [KS95] to the model problem (1.1) reads (with polar coordinates $(r,\phi)$ and $\kappa=3-4\nu$, $\gamma=2\nu-1$, $a=1/4$)

[TABLE]

The approximation of the circular inclusion through a polygon is rather critical for the higher-order Arnold–Winther MFEM. In the absence of an implementation of parametric boundaries, adaptive mesh refinement is necessary for higher improvements. The adaptive algorithm of this section is the same for all examples and acts on polygons; in particular, it does not monitor the curved boundary, but whenever some edge at the curved part $\Gamma_{D}$ is refined in this example, the midpoint is a new node and projected onto $\Gamma_{D}$. The convergence history plot in Figure 3 shows a reduced convergence for uniform refinement, while adaptive refinement (of the circular boundary) leads to optimal third-order convergence.

4.3 L-shaped benchmark

Consider the rotated L-shaped domain with Dirichlet and Neumann boundary depicted in Figure 4. The exact solution reads in polar coordinates

[TABLE]

The constants are $C_{1}:=-\cos((\alpha+1)\omega)/\cos((\alpha-1)\omega)$ and $C_{2}:=2(\lambda+2\mu)/(\lambda+\mu)$, where $\alpha=0.544483736782$ is the first root of $\alpha\sin(2\omega)+\sin(2\omega\alpha)=0$ for $\omega=3\pi/4$. The volume force $f\equiv 0$ and the Neumann boundary data $g\equiv 0$ vanish, and the Dirichlet boundary conditions $u_{D}$ are extracted from the exact solution.

Figure 5 shows suboptimal convergence $\mathcal{O}(N_{\ell}^{-0.27})$, namely an expected rate $\alpha$ in terms of the maximal mesh-size, for uniform and fourth-order $L^{2}$ stress convergence for adaptive mesh-refinement.

Despite the singular solution, the adaptive algorithm recovers the higher convergence of Theorem 17 as in [CG16].

4.4 Cook membrane problem

One of the more popular benchmarks in computational mechanics is the tapered panel $\Omega$ with the vertices $A,B,C,D$ of Figure 6 clamped on the left side $\Gamma_{D}=\operatorname{conv}(D,A)$ (with $u_{D}\equiv 0$) under no volume force ($f\equiv 0$) but applied surface tractions $g=(0,1)$ along $\operatorname{conv}(B,C)$ and traction free on the remaining parts $\operatorname{conv}(A,B)$ and $\operatorname{conv}(C,D)$ along the Neumann boundary.

This example is a particular difficult one for the Arnold–Winther MFEM because of the incompatible Neumann boundary conditions on the right corners [CGRT08, CEG11, CG16]. That means, although $g$ is piecewise constant, $g$ does not belong to $G(\mathcal{T})$ for any triangulation. In the two Neumann corner vertices $B$ and $C$ we therefore strongly impose the values $\sigma_{\ell}(B)=(0.2491,0.7283;0.7283,0.6676)$ and $\sigma_{\ell}(C)=(3/20,11/20;11/20,11/60)$ for the design of $g_{\ell}\in G(\mathcal{T}_{\ell})$.

Since the exact solution is unknown, the error approximation rests on a reference solution $\tilde{\sigma}$ computed as $P_{5}(\mathcal{T})$ displacement approximation on the uniform refinement of the finest adapted triangulation.

The large pre-asymptotic range of the convergence history plot in Figure 7 illustrates the difficulties of the Arnold-Winther finite element method in case of incompatible Neumann boundary conditions according to its nodal degrees of freedom. Once the resulting and dominating boundary oscillations (caused by the necessary choice of discrete compatible Neumann conditions in $G(\mathcal{T}_{\ell})$) $\operatorname{osc}(g-g_{\ell},\mathcal{E}_{\ell}(\Gamma_{N}))$ are resolved through adaptive mesh-refining, even the fourth-order $L^{2}$ stress convergence is visible in a long asymptotic range in (the approximated error and) the equivalent error estimator.

This example underlines that adaptive mesh-refining is unavoidable in computational mechanics with optimal rates and a large saving in computational time and memory compared to naive uniform mesh-refining.

Acknowledgements

The work has been written, while the three authors enjoyed the hospitality of the Hausdorff Research Institute of Mathematics in Bonn, Germany, during the Hausdorff Trimester Program Multiscale Problems: Algorithms, Numerical Analysis and Computation. The research of the first author (CC) has been supported by the Deutsche Forschungsgemeinschaft in the Priority Program 1748 ‘Reliable simulation techniques in solid mechanics. Development of non-standard discretization methods, mechanical and mathematical analysis’ under the project ‘Foundation and application of generalized mixed FEM towards nonlinear problems in solid mechanics’ (CA 151/22-1). The second author (DG) has been supported by the Deutsche Forschungsgemeinschaft (DFG) through CRC 1173. The third author (JG) has been funded by the Austrian Science Fund (FWF) through the project P 29197-N32.

Appendix: Fourth-order convergence of the stress in $L^{2}$

This appendix explains a high-order convergence phenomenon observed in some numerical benchmark examples for the lowest-order Arnold-Winther method. Adopt the notation from this paper for $k=1$ and let $\sigma$ solve (1.1) and let $\sigma_{h}\in AW_{k}(\mathcal{T})$ solve (1.3).

Theorem 17** (fourth-order convergence).**

Suppose $f=f_{h}\in P_{1}(\mathcal{T};\mathbb{R}^{2})$ and $g=g_{h}\in G(\mathcal{T})$ and suppose that the stress $\sigma\in H^{4}(\Omega;\mathbb{S})$. Then the $L^{2}$ stress error satisfies (with the maximal mesh-size $h$)

[TABLE]

Proof.

Since the stress error $\sigma-\sigma_{h}\in H^{4}(\mathcal{T};\mathbb{S})$ is divergence-free, $\alpha$ vanishes in (2.1) and $\sigma-\sigma_{h}=\operatorname{Curl}^{2}\beta\in H^{4}(\mathcal{T};\mathbb{S})$. Since $\beta\in H^{2}(\mathcal{T})$ is piecewise in $C^{2}$, it follows $\beta\in C^{1}(\overline{\Omega})$. The Arnold-Winther finite elements have nodal degrees of freedom at the vertices and hence $\sigma_{h}$ is continuous at each vertex $z\in\mathcal{N}$. Hence the second derivatives of $\beta\in C^{2}(\mathcal{T})\cap C^{1}(\overline{\Omega})$ are continuous at each vertex $z\in\mathcal{N}$. It follows that the nodal interpolation operator $I_{A}$ associated to the Argyris finite element space $A(\mathcal{T})\subset C^{1}(\Omega)\cap P_{5}(\mathcal{T})$ exists for $\beta$ in the classical sense and is composed of the piecewise local interpolation. This defines $\beta_{h}=I_{A}\beta$ and the divergence-free $\tau_{h}:=\operatorname{Curl}^{2}\beta_{h}\in AW(\mathcal{T})$ test function in (1.1) and in (1.3). Consequently,

[TABLE]

This and standard local interpolation error estimates for the nodal interpolation of the quintic Argyris finite elements [BS08, Bra07, Cia78] show

[TABLE]

(With $\sigma-\sigma_{h}=\operatorname{Curl}^{2}\beta$ and $|\sigma_{h}|_{H^{4}(\mathcal{T})}=0$ for piecewise cubic $\sigma_{h}$ in the last step.) ∎