# Empty Neumann boundary measures from mesh function

On a 3D domain, I want to explicitly assemble Neumann boundary conditions using

• `gmsh` with `Physical Surfaces`
• `dolfin-convert` to convert from `.msh` to the `dolfin` input files
• and a `dolfin.MeshFunction` to assess boundary measures `ds`

What works well on one example, simply fails on another example (see the code below). The returned measure is just empty.

Even more difficult for me to understand, if I use the same `MeshFunction` in `dolfin.DirichletBC` definitions, the resulting values are not empty.

I suspect that there are many pitfalls in treating meshes and boundary conditions. But maybe someone has an idea how to go about this issue.

Here is the example. The mesh `karman3D...` does not work, whereas the `4-segs-...` mesh works (in the sense that the Neumann boundaries do not appear empty)

Btw. I’m using `dolfin 2019`

``````import dolfin
import numpy as np
import scipy.sparse as sps

meshlevel = 1

meshprfx = 'karman3D_lvl{0}'.format(meshlevel)   # not working
meshprfx = '4-segs-3D_lvl{0}'.format(meshlevel)  # working

meshfile = meshprfx + '.xml.gz'
physregs = meshprfx + '_physical_region.xml.gz'
fctsregs = meshprfx + '_facet_region.xml.gz'

mesh = dolfin.Mesh(meshfile)

boundaries = dolfin.MeshFunction('size_t', mesh, fctsregs)
subdomains = dolfin.MeshFunction('size_t', mesh, physregs)

# a Paraview file to check the facet numbering
pvd = dolfin.File('mf.pvd')
pvd << boundaries

ds = dolfin.Measure('ds', subdomain_data=boundaries)

V = dolfin.FunctionSpace(mesh, 'CG', 1)

# checking as Dirichlet type
bcexp = dolfin.Expression("1", degree=1)
for diripe in range(15):
diribcs = dolfin.DirichletBC(V, bcexp, boundaries, diripe)
bcdict = diribcs.get_boundary_values()
bcvals = []
bcvals.extend(bcdict.values())
print(f'id:{diripe} -- Dirichlet valsum: {np.array(bcvals).sum()}')

# checking as Neumann type
v = dolfin.TestFunction(V)
for pe in range(15):
obsop = dolfin.assemble(v*ds(pe)).get_local().reshape((1, V.dim()))
obsopmat = sps.csc_matrix(obsop)
print(f'id:{pe} -- Neumann valsum: {obsopmat.sum()}')

``````

The issues seems to be that I have two `physical volumes` defined in the mesh file (which I used to apply different diffusion parameters in my code).

If I only have one volume – everything works fine. Anyone here who had similar experiences?