29.1. General framework#

Equilibrated residual error estimators provide upper bounds for the discretization error in energy norm without any generic constant. We consider the standard problem: find $u\in V:=H_0^1(\mathrm{\Omega })$ such that

The left hand side defines the bilinear-form $A(\cdot ,\cdot )$, the right hand side the linear-form $f(\cdot )$. We define a finite element sub-space $V_h\subset V$ of order $k$, and the finite element solution

We assume that $f$ is element-wise polynomial of order $k-1$, and $\lambda$ is element-wise constant and positive.

The residual $r(\cdot )\in V^{\ast }$ is

Since

we aim in estimating $\Vert r\Vert$ in the norm dual to $\Vert \cdot {\Vert }_A$, which is essentially the $H^{-1}$-norm. In general, the direct evaluation of this norm is not feasible. Using the structure of the problem, we can represent the residual as

where $r_T$ and $r_E$ are given given as

The element-residual $r_T$ is a polynomial of order $k-1$ on the element $T$, and the edge residual (the normal jump) is a polynomial of order $k-1$ on the edge $E$.

The residual error estimator estimates the residual in terms of weighted $L_2$-norms:

Here, ${\lambda }_E$ is some averaging of the coefficients on the two elements containing the edge $E$. The equivalence holds with constants depending on the shape of elements, the relative jump of the coefficient, and the polynomial order $k$.

The equilibrated residual error estimator ${\eta }^{er}$ is defined in terms of the same data $r_T$ and $r_E$. It satisfies

The lower bound depends on the shape of elements and the coefficient $\lambda$, but is robust with respect to the polynomial order $k$.

The main idea is the following: Instead of calculating the $H^{-1}$-norm of $r$, we compute a lifting ${\sigma }^{\mathrm{\Delta }}$ such that $\mathrm{div}{\sigma }^{\mathrm{\Delta }}=r$, and calculate the $L_2$-norm of ${\sigma }^{\mathrm{\Delta }}$. Since $r$ is not a regular function, the equation must be posed in distributional form:

Then, the residual can be estimated without envolving any generic constant:

The norm $\Vert {\sigma }^{\mathrm{\Delta }}\Vert :=\int {\lambda }^{-1}|{\sigma }^{\mathrm{\Delta }}{|}^2$ can be evaluated easily.

Remark: The flux-postprocessing $\sigma :=\lambda \mathrm{\nabla }{u}_h+{\sigma }^{\mathrm{\Delta }}$ provides a flux $\sigma \in H(\mathrm{div})$ such that $\mathrm{div}\sigma =f$, i.e. the flux is in exact equilibrium with the source $f$. Thus the name.

29.2. Computation of the lifting $\Vert {\sigma }^{\mathrm{\Delta }}\Vert$#

The residual is a functional of the form

where $r_T$ and $r_E$ are polynomials of order $k-1$. We search for ${\sigma }^{\mathrm{\Delta }}$ which is element-wise a vector-valued polynomial of order $k$, and not continuous across edges. Element-wise integration by parts gives

Thus $\mathrm{div}\sigma =r$ in distributional sense reads as

for all elements $T$ and edges $E$. We could now pose the problem

We minimize the weighted-$L_2$ norm since we want to find the smallest possible upper bound for the error. This is already a computable approach. But, the problem is global, and its solution is of comparable cost as the solution of the original finite element system. The existence of a $\sigma$ such that $\mathrm{div}\sigma =r$ also needs a proof.

We want to localize the construction of the flux. Local problems are associated with vertex-patches ${\omega }_V={\cup }_{T:V\in T}T$. We proceed in two steps:

• localization of the residual: $r=\sum _V{r}^V$

• local liftings: find ${\sigma }^V$ such that $\mathrm{div}{\sigma }^V=r^V$ on the vertex patch. Then, for $\sigma :=\sum {\sigma }^V$ there holds $\mathrm{div}\sigma =r$

The localization is given by multiplication of the $P^1$ vertex basis functions (hat-functions) ${\varphi }_V$:

Since $\sum _V{\varphi }_V=1$, there holds $\sum r^V(\cdot )=r(\cdot )$. The localized residual has the same structure of element and edge terms:

with

The local residual vanishes on constants on the patch:

The last equality follows from the Galerkin-orthogonality.

We give an explicit construction of the lifting ${\sigma }^V$ in terms of the Brezzi-Douglas-Marini (BDM) element. The $k^{th}$ order BDM element on a triangle is given by $V_T=[P^k{]}^2$ and the degrees of freedom:

• ${\int }_E\sigma \cdot n{q}_i$ with $q_i$ a basis for $P^k(E)$

• ${\int }_T\mathrm{div}\sigma q_i$ with $q_i$ a basis for $P^{k-1}(T)\cap L_2^0(T)$

• ${\int }_T\sigma \cdot \mathrm{curl}{q}_i$ with $q_i$ a basis for $P_0^{k+1}(T)$

Exercise: Show that these dofs are unisolvent. Count dimensions, and prove that $[\mathrm{\forall }i:{\psi }_i(\sigma )=0]\Rightarrow \sigma =0$.

Now, we give an explicit construction of equilibrated fluxes on a vertex patch. Label elements $T_1,T_2,\dots T_n$ in a counter-clock-wise order. Edge $E_i$ is the common edge between triangle $T_{i-1}$ and $T_i$ (with identifying $T_0=T_n$). We define $\sigma$ by specifying the dofs of the BDM element:

• Start on $T_1$. We set ${\sigma }_n=-r_{E_1}^V$ on edge $E_1$. On the edge on the patch-boundary we set ${\sigma }_n=0$, and on $E_2$ we set ${\sigma }_n=const$ such that ${\int }_{\partial T_1}{\sigma }_n={\int }_{T_1}{r}_T^V$. We use the dofs of type (ii) to specify ${\int }_T\mathrm{div}\sigma q={\int }_T{r}_T^{V}q\mathrm{\forall }q\in P^{k-1}\cap L_2^0(T)$. Together with get $\mathrm{div}\sigma =r_T$. Dofs of type (iii) are not needed, and set 0. There holds $${\int }_{E_2}{\sigma }_n={\int }_{T_1}{r}_T^V-{\int }_{E_1}{\sigma }_n={\int }_{T_1}{r}_{T_1}^V+{\int }_{E_1}{r}_{E_1}^V$$

• Continue with element $T_2$. On edge $E_2$ common with $T_1$ set ${\sigma }_n$ such that $[\sigma \cdot n{]}_{E_2}=r_{E_2}$. Otherwise, proceed as on $T_1$. Thus $${\int }_{E_3}{\sigma }_n={\int }_{T_1}{r}_{T_1}^V+{\int }_{E_1}{r}_{E_1}^V+{\int }_{T_2}{r}_{T_2}^V+{\int }_{E_2}{r}_{E_2}^V$$

• Continue to element $T_n$. Observe that on $T_n$: $${\int }_{E_1}{\sigma }_n=\sum _{i=1}^n{\int }_{T_i}{r}_{T_i}^V+\sum _{i=1}^n{\int }_{E_i}{r}_{E_i}^V=0,$$whichfollowsfrom$r^V(1) = 0$.Thus,also$[\sigma \cdot n]{E_1} = r{E_1}^V$ is satisfied.

This explicit construction proves the existence of an equilibrated flux. Instead of this explicit construction, one may solve a local constrained optimization problem

This applies also for 3D. Furthoer notes

• mixed boundary conditions are possible

• the efficiency for the h-FEM is shown by scaling arguments, and equivalence to the residual error estimator

• efficiency is also proven to be robust with respect to polynomial order $k$, examples show overestimation less than 1.5

Literature:

• D. Braess and J. Schöberl: Equilibrated Residual Error Estimator for Maxwell’s Equations. Mathematics of Computation, Vol 77(262), 651-672, 2008

• D. Braess, V. Pillwein and J. Schöberl: Equilibrated Residual Error Estimates are p-Robust. Computer Methods in Applied Mechanics and Engineering. Vol 198, 1189-1197, 2009

29.3. Equilibration in NGSolve#

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import *

r1 = Rectangle(1,1).Face()
r1.edges.Max(X).name='right'
r1.edges.Min(X).name='left'
r1.edges.Max(Y).name='top'
r1.edges.Min(Y).name='bottom'
r1.faces.name='air'
r2 = MoveTo(0.6,0.6).Rectangle(0.2,0.2).Face()
r2.faces.name='source'

r1 -= r2
shape = Glue( [r1,r2] )

mesh = Mesh(OCCGeometry(shape,dim=2).GenerateMesh(maxh=0.1))
Draw (mesh);

order = 3
order_equ = 3
dirichlet="left|bottom"
fes = H1(mesh, order=order, dirichlet=dirichlet)
u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
source = mesh.MaterialCF( { "source" : x*y }, default=0 )
f = LinearForm(source*v*dx).Assemble()

gfu = GridFunction(fes, name="solution")
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw (gfu);

construct $\sigma$ such that $\mathrm{div}\sigma =f+\mathrm{\Delta }{u}_h$ from local contributions $\sigma =\sum {\sigma }^V$:

where $\mathrm{div}{\sigma }^V={\varphi }_V(f+\mathrm{\Delta }{u}_h)$ is understood in distributional sense, what is

# variational formulation for patch correction problems
Xsigma = Discontinuous (HDiv(mesh, order=order_equ))
Xw     = L2(mesh,order=order_equ-1)
Xwf    = FacetFESpace(mesh, order=order_equ)
Xlam   = NumberSpace(mesh)
X = Xsigma*Xw*Xwf*Xlam 

sigma, w, wf, lam = X.TrialFunction()
tau,   v, vf, mu  = X.TestFunction()

n = specialcf.normal(mesh.dim)
bfequ = (tau*sigma + w*div(tau) + v*div(sigma) + w*mu + lam*v) * dx 
bfequ += (-sigma*n*vf-tau*n*wf)*dx(element_boundary=True)  
bfequ += (-1e10*wf.Trace()*vf.Trace())*ds(definedon=mesh.Boundaries(dirichlet)) # penalty for Dirichlet
bfequ += (-1e10*mu*lam)*ds(definedon=mesh.Boundaries(dirichlet))  # no constraint at Dir bnd

# the residual: element-term + edge term:
lfequ = (source+Trace(gfu.Operator("hesse")))*v*dx 
lfequ += (-grad(gfu)*n*vf)*dx(element_boundary=True)

gfequ = GridFunction(X)

with TaskManager():
    PatchwiseSolve (bf=bfequ, lf=lfequ, gf=gfequ)

fesflux = HDiv(mesh, order=order_equ)
equflux = GridFunction(fesflux)
equflux.Set (grad(gfu)-gfequ.components[0])

Draw (equflux, mesh);

Check equilibrium:

Draw (div(equflux)+source, mesh, "residuum");