Empty Neumann boundary measures from mesh function

General
1

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)

Please see here for downloading the meshes and this code snippet.

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()}')
2

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?