Implementation of a residual error estimator for the Poisson problem - a lifting technique

In this example we want to show how one can easily implement an error estimator for the Poisson problem $$ - \Delta u = f \quad \textrm{on } \Omega, $$ with \( f = \sin(\pi x) \sin(\pi y) \) and \(\Omega = (0,1) \times (0,1) \).
# solve primal problem ...
V = H1(mesh, order=order, dirichlet="bottom|right|left|top")
u = V.TrialFunction()
v = V.TestFunction()

n = specialcf.normal(2)

a = BilinearForm(V)
a += SymbolicBFI (grad(u)*grad(v))

f = sin(pi*x)*sin(pi*y)
f_rhs = LinearForm(V)
f_rhs += SymbolicLFI (f*v)

u = GridFunction(V)
u.vec.data = a.mat.Inverse(freedofs=V.FreeDofs()) * f_rhs.vec
The residual error estimator is then given by $$ \eta(u_h)^2 := \sum\limits_{T \in \mathcal{T}} h_T^2 \| f + \Delta u_h\|_{L^2(T)}^2 + \sum\limits_{E \in \mathcal{F}} h_E \| [ \! [ \nabla u_h \cdot n ] \! ] \|_{L^2(E)}^2. $$ For the element contributions of the error estimator, we need the second derivative of our solution on each triangle. For this we simply define a new space V2 and create two gridfunctions gradu_x and gradu_y and use the Set-function of NGSolve to set the new functions to the derivative of the solution. We choose V2 as a \( L^2 \) space as the Set-function then is a simple elementwise interpolation, i.e. exact. Note that, although gradu_x and gradu_y are \( L^2 \) functions, they can still be derived (as this is needed for DG-methods!). Using the integrate function we get the first term.
# calc element terms of error estimator:
# create a L^2 space to set the gradient to the derivative of the solution
V2 = L2(mesh, order = order-1)

gradu_x = GridFunction(V2)
gradu_y = GridFunction(V2)


h = specialcf.mesh_size

laplace_u = gradu_x.Deriv()[0] + gradu_y.Deriv()[1]

eta_element = Integrate(h*h*(f + laplace_u)*(f + laplace_u), mesh)
For the edge contributions we are going to use a lifting techniuqe. The idea is to solve a little problem which reads as a \( L^2 \) minimization problem with a constraint on the edges, thus we seek for \( (u_2,\sigma) \in H^1 \times H(\textrm{div}) \) that fullfills $$ \int_\Omega u_2 v_ 2 + \varepsilon \int_\Omega \sigma \cdot \tau + \sum\limits_{E \in \mathcal{F}} \int_E \sigma \cdot n ~ \tau \cdot n = \sum\limits_{E \in \mathcal{F}} \int_E [ \! [ \nabla u_h \cdot n ] \! ] \tau \cdot n \qquad \forall (v_2,\tau), $$ where \( \varepsilon \) is a small regularisation parameter. For the implementation of this problem we define two bilinearforms. The first one represents the left hand side of our problem.
Sigma = HDiv(mesh, order=order-1)
X = FESpace([V,Sigma])

u2, sigma = X.TrialFunction()
v2, tau = X.TestFunction()

asig = BilinearForm(X)
asig += SymbolicBFI (u2*v2)
asig += SymbolicBFI (1e-4 * sigma*tau)
asig += SymbolicBFI ( 0.5*(sigma*n)*(tau*n), element_boundary=True)
The 0.5 factor appears as the integral on the edge using the element_boundary flag is evaluated on each side of the triangle (but \( \sigma \cdot n \) is continuous!). For the right had side we also define a bilinearform as we are going to use the apply function of NGSolve, thus a matrix vector muliplication.
amixed = BilinearForm(X)
amixed += SymbolicBFI ((tau*n) * (n* (grad(u2)-grad(u2).Other())), VOL, skeleton=True)
The jump on the edge is implemented using the Other() function of the TrialFunction. Note that we do not assemble this matrix! Now we can simply solve this problem:
sigma = GridFunction(X, name="lifted_jumps")
hv1 = sigma.vec.CreateVector()
hv2 = sigma.vec.CreateVector()

hv1[:] = 0
hv1[X.Range(0)] = u.vec
amixed.Apply (hv1, hv2)
sigma.vec.data = asig.mat.Inverse() * hv2
Using the trace inequality we get the desired edge contributions by
sigma2 = sigma.components[1]
eta_edge = Integrate (sigma2*sigma2, mesh)
residual ee.py [1.84Kb]
Uploaded Wednesday, 11 January 2017 by Philip Lederer