Featured

## 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.Assemble()

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

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)

h = specialcf.mesh_size

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)
asig.Assemble()

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)

Attachments
residual ee.py [1.84Kb]
Uploaded Wednesday, 11 January 2017 by Philip Lederer