Hi,

I am having issues solving a Poisson equation with Neumann BCs using external mesh files.

The code doesn’t show any errors but Neumann BCs seem not to be applied. When solving Poisson equations using a FEnics mesh function (e.g. UnitSquareMesh) the simulation works fine by marking boundaries, instantiating DirichletBC() objects as well as generating expressions for NeumannBCs as described in tutorial 14.

(https://fenicsproject.org/pub/tutorial/html/._ftut1014.html)

When using an external mesh files, Neumann BCs don’t seem to work. As an example, I am considering a geometry (generated in gmsh and converted using dolfin-convert tool) that consists of a rectangle with a hole and 5 boundaries defined as follows; Upper part noted *up*, down part noted *dp*, left part noted *lp*, right part noted *rp* and the hole noted *hp*. Constant values are defined for each boundary parts and Neumann BC for the hole.

```
from dolfin import*
mesh = Mesh("mesh.xml")
subdomains = MeshFunction("size_t", mesh, "mesh_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "mesh_facet_region.xml")
V = FunctionSpace(mesh, 'P', 1)
C_up = Constant("20")
C_dp = Constant("200")
C_lp = Constant("60")
C_rp = Constant("120")
C_hp = Constant("150")
kappa = Constant("2")
f = Constant("0")
boundary_conditions = {0: {'Dirichlet': C_up},
1: {'Dirichlet': C_dp},
2: {'Dirichlet': C_lp},
3: {'Dirichlet': C_rp},
4: {'Neumann': C_hp}}
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
boundaries, i)
bcs.append(bc)
u = TrialFunction(V)
v = TestFunction(V)
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
if boundary_conditions[i]['Neumann'] != 0:
g = boundary_conditions[i]['Neumann']
integrals_N.append(g*v*ds(subdomain_id=i, domain=mesh, subdomain_data=boundaries))
a = kappa*inner(grad(u), grad(v))*dx
L = f*v*dx - sum(integrals_N)
u = Function(V)
solve(a == L, u, bcs)
file = File("poisson.pvd")
file << u
```

Solving the problem with only Dirichlet BCs is working.

```
boundary_conditions = {0: {'Dirichlet': C_up},
1: {'Dirichlet': C_dp},
2: {'Dirichlet': C_lp},
3: {'Dirichlet': C_rp},
4: {'Dirichlet': C_hp}}
```

Is there anything missing when applying Neumann BCs?

Many thanks in advance.