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$μ .
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.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.
1. Why the two computed H^2 norms are not the same?
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.
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.