integration of Constant over subdomain fails

[Nico Schlömer]

Another submesh issue (compare with issue #240 and http://fenicsproject.org/qa/2714/integration-subdomain-matching-function-spaces-measures, http://fenicsproject.org/qa/2719/assemlbing-over-subdomain-expecting-completed-domains-point):

Plain integration of 1.0 over a submesh fails with unhelpful error messages.

MWE:

from dolfin import *

mesh = RectangleMesh(0.0, 0.0, 1.0, 1.0, 100, 100)


class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))

obstacle = Obstacle()

subdomains = CellFunction('size_t', mesh)
subdomains.set_all(0)
obstacle.mark(subdomains, 1)

dx = Measure('dx')[subdomains]

submesh = SubMesh(mesh, subdomains, 1)
# Fails:
assemble(Constant(1.0) * dx(1))
# Fails:
assemble(Constant(1.0) * dx(1), mesh=submesh)
# Works:
assemble(Constant(1.0) * dx(1), mesh=mesh)

[Jan Blechta]

This is same problem as with #240. You're integrating on non-matching meshes. If you want to integrate on on submesh you need to provide a measure defined on submesh, not on mesh.

[Martin Sandve Alnæs]

So I have a fix on the way, but there's a bug in your example. Your 'subdomains' relate to the 'mesh' and should not be used with the 'submesh'.

[Martin Sandve Alnæs]

240 is invalid, this one shows an actual bug as well.

[Jan Blechta]

Ok, you mean first assemble. This is not a bug per se but rather a missing feature.

[Martin Sandve Alnæs]

Fix is in next (need to update ufl too to next).

[Nico Schlömer]

I updated to next and still get the same errors with the MWE in the original post.

[Jan Blechta]

Yes, assemble(Constant(1.0) * dx(1)) still does not work as form Constant(1.0) * dx(1) can't be initialized by mesh given within measure.

assemble(Constant(1.0) * dx(1), mesh=submesh) now raises appropriate error message

*** Error:   Unable to extract mesh from form.
*** Reason:  Non-matching meshes for function spaces and/or measures.
*** Where:   This error was encountered inside Form.cpp.

[Nico Schlömer]

I understand now why there is a mismatch, but find that rather unintuitive. After all, dx(1) is the measure belonging to subdomain 1, and so is the constructed submesh. Also, Constant(1.0) * dx(1) looks legit since all information about integrand and domain is given.

[Martin Sandve Alnæs]

It's subdomains and submesh that don't match, the former relating to cells of the mesh not of the submesh. I guess the error message could be clearer on that.

[Martin Sandve Alnæs]

Feel free to contribute with better error checking and error messages but I consider this particular issue closed.

[Prof Garth Wells]

Removing milestone: 1.4 (automated comment)