I want to solve the poisson equation (electrostatics problem) with a neumann boundary condition

\left[\frac{\partial u}{\partial n}\right]=800 (constant surface charge density) on some interior facets.

However, the boundary condition seems to be completely ignored and all I get is a constant solution.

```
#-----------Generate Mesh------------
rank = MPI.COMM_WORLD.rank
NEUMANN_TAG = 2137
gmsh.initialize()
gdim = 2
model_rank = 0
mesh_comm = MPI.COMM_WORLD
d = 2 # width of capacitor
h = 1 # height of capacitor plate
r = 0.02 # radius of sphere
sphere = gmsh.model.occ.add_disk(d/2, h/2, 0, r, r)
gmsh.model.occ.synchronize()
box = gmsh.model.occ.add_rectangle(0,0,0, d, h)
gmsh.model.occ.synchronize()
background = gmsh.model.occ.fragment([(2, box)], [(2, sphere)])
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(gdim, [box, sphere], 1, "air")
faces = gmsh.model.getBoundary([(2, sphere)], oriented=False)
gmsh.model.addPhysicalGroup(1, [faces[0][1]], NEUMANN_TAG)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(gdim)
#----Actual solving---------
msh, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
V = fem.functionspace(msh, ("DG", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
neumann_ds = Measure("dS", domain=msh, subdomain_data=ft)
f = x[0]-x[0]
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx - 800* v("+")* neumann_ds(NEUMANN_TAG)
problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
```