Computing functions on boundaries gives incorrect result in parallel

[Garth Wells]

The MeshFunction that stores exterior facets and which is used in parallel, is computed when a DofMap is created. This is a problem when assembling functionals on exterior facets and passing the mesh because the mesh never computes the "exterior facets" MeshFunction. The below code reproduced the error.

from dolfin import *
mesh = UnitSquare(24, 24)
f = Constant(1.0)
M = f*ds

# This is a hack to get around the bug
#if MPI.num_processes() > 1:
# cpp.MeshPartitioning.number_entities(mesh, mesh.topology().dim() - 1);

# This should print 4.0
print assemble(M, mesh=mesh)

[Marie Rognes (meg-simula)]

This is probably the same problem:

from dolfin import *

mesh = UnitSquare(4, 4)

boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
boundaries[0] = 1

V = FunctionSpace(mesh, "CG", 2)
f = Expression("x[0] + 2")
f = interpolate(f, V)

M = f*ds(1)

# This should print 0 but prints 0.75
print assemble(M, exterior_facet_domains=boundaries)

[Anders Logg (logg)]

Why don't we just call MeshPartitioning.number_entities?