I have created a mesh in `gmsh`

for the Haut Glacier d’Arolla (www.the-cryosphere.net/2/95/2008/) and converted to `.xml`

format for use in `FEniCS`

using `dolfin-convert`

.

I decided to test my mesh by solving the following Poisson problem.

After solving the discretized problem, checking that the solution is reasonable by visualizing in `ParaView`

I then check the H1 semi-norm of the solution, defined as

in two ways. First by using the `assemble`

command in `dolfin`

, the numerical value returned is `nan`

. There is an alternative way to compute the H1 semi-norm. That is

Where **A** is the finite element stiffness matrix. **u**^T**A u** evaluates to a reasonable nonnegative finite number.

I have included the full code below., but of course you will not be able to test without the mesh.

I am hoping to get insight into why this would be the case, what goes on internally in `FEniCS`

when one assembles a finite element system matrix, such as **A**, with Dirichlet conditions that causes the discrepancy that I am seeing.

import dolfin as dl

import numpy as np

```
# mesh of PDE domain
meshdir = "../"
meshname = "Arolla2"
mesh = dl.Mesh(meshdir+meshname+".xml")
boundary_markers = dl.MeshFunction("size_t", mesh,\
meshdir+meshname+"_facet_region.xml")
# discrete function spaces
P1 = dl.FiniteElement("CG", mesh.ufl_cell(), 1)
Vh = dl.FunctionSpace(mesh, P1)
# Homogeneous Dirichlet conditions
bc = [dl.DirichletBC(Vh, dl.Constant(0.0), boundary_markers, 1),\
dl.DirichletBC(Vh, dl.Constant(0.0), boundary_markers, 2)]
# Right hand side forcing term
f = dl.Constant(1.0)
vh = dl.TestFunction(Vh)
uh = dl.TrialFunction(Vh)
a = dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx
L = f*vh*dl.dx
A, b = dl.assemble_system(a, L, bc)
Asolver = dl.PETScKrylovSolver()
dl.PETScOptions.set("ksp_monitor")
dl.PETScOptions.set("ksp_monitor_true_residual")
dl.PETScOptions.set("ksp_type", "cg")
dl.PETScOptions.set("ksp_rtol", "1.e-12")
dl.PETScOptions.set("ksp_atol", "1.e-12")
Asolver.set_from_options()
Asolver.set_operator(A)
u = dl.Function(Vh)
Asolver.solve(u.vector(), b)
# for the determined solution,
# compute sqrt(int_(Omega) ||grad(u)||^2 dV)
u_H1norm_1 = np.sqrt(dl.assemble(dl.inner(dl.grad(u), dl.grad(u))*dl.dx(mesh)))
print("||u||_H^1(Omega)^2 = %1.3e (directly assembled)" %(u_H1norm_1))
# vector to hold the product of the solution with
# the finite element system matrix A
Au = dl.Vector()
A.init_vector(Au, 0)
# Au = A*u
A.mult(u.vector(), Au)
# sqrt(u^T Au)
u_H1norm_2 = np.sqrt(Au.inner(u.vector()))
print("||u||_H^1(Omega) = %1.3e (from assembled operators)" \
%(u_H1norm_2))
```