No error raised when assemble() over Mesh with a MeshFunction from another Mesh is used as subdomains

Here's a bug that was really hard to find and that bit for a long time.

When assemble()ing an Expression over a subdomain, Dolfin needs to have the mesh supplied along with the Expression. When supplying the entire mesh and specifying the subdomain by dx(k), everything works as expected. When supplying the corresponding SubMesh to assemble(), however, Dolfin computes something else.

The following code highlights the behavior.

from dolfin import *

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5 + DOLFIN_EPS
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return 0.5 - DOLFIN_EPS < x[0]

left = Left()
right = Right()

mesh = UnitSquareMesh(20, 20)

domains = CellFunction('size_t', mesh)
left.mark(domains, 1)
right.mark(domains, 2)

dx = Measure('dx')[domains]

f = Expression('x[0]', cell=triangle)
print(assemble(f*dx(2), mesh=mesh))
print(assemble(f*dx(2), mesh=SubMesh(mesh, domains, 2)))

Comments (6)