How to check that the Neumann boundary condition is correctly applied

variational formulation
4

Consider the following example:

from dolfin import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 2)
u, v = TrialFunction(V), TestFunction(V)

a = dot(grad(u), grad(v))*dx
L = Constant(1.0)*v*dx

LEFT, RIGHT, TOP, BOTTOM = 1, 2, 3, 4
ff = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
AutoSubDomain(lambda x: x[0] < DOLFIN_EPS).mark(ff, LEFT)
AutoSubDomain(lambda x: x[0] > 1.0 - DOLFIN_EPS).mark(ff, RIGHT)
AutoSubDomain(lambda x: x[1] < DOLFIN_EPS).mark(ff, BOTTOM)
AutoSubDomain(lambda x: x[1] > 1.0 - DOLFIN_EPS).mark(ff, TOP)

bcs = [DirichletBC(V, Constant(0.0), ff, BOTTOM)]

uh = Function(V)
solve(a == L, uh, bcs)

n = FacetNormal(mesh)
ds = Measure("ds", subdomain_data=ff)
solution_trace_norm = assemble(dot(grad(uh), n)**2*ds((LEFT, RIGHT, TOP)))**0.5
print(f"Homogeneous Neumann BC trace: {solution_trace_norm}")

which gives me:

Homogeneous Neumann BC trace: 1.0071881730182207e-13