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))
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)
gradu_x = GridFunction(V2)
gradu_y = GridFunction(V2)
gradu_x.Set(u.Deriv()[0])
gradu_y.Set(u.Deriv()[1])
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)
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)
