### Compute H^2 norm

I would like to compute the H^2 norm of a field. Let say $w$

`w`is the numerical solution of a Poisson problem with homogenous Dirichlet boundary condition, a constant source term $f$

`ƒ`and a piecewise constant coefficient $\mu$

`μ`.

$\mu\Delta w=f$

`μ`Δ

`w`=

`ƒ`

This solution is computed with Lagrange finite element of order 2.

```
N = 128
mesh = UnitSquareMesh(N,N)
V2 = FunctionSpace(mesh,'CG',2)
def no_slip_bound(x,on_boundary):
return on_boundary
bc = DirichletBC(V2,Constant(0.),no_slip_bound)
# Parameters
f = Constant(1.)
mu1 = 1.
mu2 = 1.
mu = Expression("mu1 + mu2*(x[0]>=0.5)",mu1=mu1,mu2=mu2,degree=0)
# Weak formulation
u = TrialFunction(V2)
v = TestFunction(V2)
w = Function(V2)
a = mu*inner(grad(u),grad(v))*dx
L = dot(f,v)*dx
solve(a==L,w,bc)
```

I found two possible ways to compute its H^2 norm, but they do not give the same result :

First solution :

```
V1 = FunctionSpace(mesh,"Lagrange",1)
V0 = FunctionSpace(mesh,"DG",0)
# Compute derivatives
dxw = w.dx(0); dxw = project(dxw,V1)
dyw = w.dx(1); dyw = project(dyw,V1)
dxxw = dxw.dx(0); dxxw = project(dxxw,V0)
dxyw = dxw.dx(1); dxyw = project(dxyw,V0)
dyxw = dyw.dx(0); dyxw = project(dyxw,V0)
dyyw = dyw.dx(1); dyyw = project(dyyw,V0)
# Compute L2, H1 and H2 norms
normL2 = sqrt(assemble(dot(w,w)*dx))
normH1 = sqrt(assemble(dot(w,w)*dx) + assemble(dot(dxw,dxw)*dx + dot(dyw,dyw)*dx))
normH2 = sqrt(assemble((dot(w,w) + dot(dxw,dxw) + dot(dyw,dyw) + dot(dxxw,dxxw) + dot(dxyw,dxyw) + dot(dyxw,dyxw) + dot(dyyw,dyyw))*dx))
```

Second solution :

```
normL2 = sqrt(assemble(dot(w,w)*dx))
normH1 = sqrt(assemble(dot(w,w)*dx) + assemble((dot(w.dx(0),w.dx(0)) + dot(w.dx(1),w.dx(1)))*dx))
normH2 = sqrt(assemble((dot(w,w) + dot(w.dx(0),w.dx(0)) + dot(w.dx(1),w.dx(1)) + dot(w.dx(0).dx(0),w.dx(0).dx(0)) + dot(w.dx(0).dx(1),w.dx(0).dx(1)) + dot(w.dx(1).dx(0),w.dx(1).dx(0)) + dot(w.dx(1).dx(1),w.dx(1).dx(1)))*dx))
```

In the first case, the H^2 norm is not bounded, which is what we should expect, since the coefficient in the Poisson equation is not continuous.

However, in the second case, the H^2 norm is bounded.

**Questions :**

1.Why the two computed H^2 norms are not the same?

1.

**2.**Can we really compute the H^2 norm with Lagrange elements, since they are not C^1 ? ( I am a bit confused about that, sorry...)

Thank you for your help.

Cheers,

Fabien

### 1 Answer

`V0 = FonctionSpace(mesh,"DG",0)`

instead of the one you are using, surely the norms will look similar. There is also the chance that you are having numerical artifacts, so maybe compare this methods to calculating the norm with

`normH2 = norm(u, 'H1') + norm(grad(u), 'H10', mesh = mesh)`

where the 0 indicates that you are not getting the \( L^2 \) term (https://fenicsproject.org/olddocs/dolfin/1.5.0/python/programmers-reference/fem/norms/norm.html) and the second one is receiving the mesh because grad(u) is no longer a Function object.... but maybe it won't run anyway. If it doesn't, try projecting first the gradient as you do with the partial derivatives to \( H^1 \). Let me know how it goes.

Best regards

```
from fenics import *
N = 100
mesh = UnitSquareMesh(N,N)
f = Expression("x[0] + x[1]", degree = 1, domain = mesh)
dx = Measure('dx', domain = mesh)
Df = grad(f)
DDf = grad(grad(f))
norm1 = sqrt(assemble( (f*f + dot(Df, Df) + inner(DDf, DDf))*dx ))
V1 = FunctionSpace(mesh, 'CG', 1)
V0 = FunctionSpace(mesh, 'DG', 0)
dxf = f.dx(0); dxf = project(dxf,V1)
dyf = f.dx(1); dyf = project(dyf,V1)
dxxf = dxf.dx(0); dxxf = project(dxxf,V0)
dxyf = dxf.dx(1); dxyf = project(dxyf,V0)
dyxf = dyf.dx(0); dyxf = project(dyxf,V0)
dyyf = dyf.dx(1); dyyf = project(dyyf,V0)
norm2 = sqrt(assemble( (f*f + dxf*dxf + dyf*dyf + dxxf*dxxf + \
dxyf*dxyf + dyxf*dyxf + dyyf*dyyf)*dx))
gradf = as_vector([dxf, dyf])
gradf = project(gradf, VectorFunctionSpace(mesh, 'CG', 1))
norm3 = norm(f, 'H1') + norm(gradf, 'H10')
print norm1, norm2, norm3, sqrt(2./3 + 1./2 + 2.)
```

where the last term is the exact solution, and I get the same value for all three. In any case, I recommend using the first one I did here, of maybe doing

`norm(f, 'H1') + sqrt(assemble(inner(DDf, DDf)*dx))`

because it is cleaner. Also, it avoids projecting, which will change a little bit your numbers (but shouldn't be wrong!), as fenics will (behiind curtains) use the exact values per element, because it just takes derivatives of the local polynomial expression. Try creating functions and comparing the projections vs interpolations to see this effect.

Best regards!

Concerning your first remark, I made a mistake in writing my question... (sorry about that ! I edited my question). In my code I have

Concerning your second remark, calculting the norm with

do not run. However, projecting first the gradient with the partial derivatives to $H^1$

H^{1}does,The bad news is that it gives a different results from the other two ways of computing the $H^2$

H^{2}norm and I do not know which one is correct.Regards,

Fabien