Integration of dolfin.Mesh and ufl.Domain

[Martin Sandve Alnæs]

(note: subclassing as discussed below wont be necessary, I've found another way)

I'm contemplating the side effects of making dolfin.Mesh a subclass of ufl.Domain, which would in many respects make for a smoother ufl/dolfin integration. However I found this comment in dolfins meshes.py:

Replacement constructor for Mesh class. We cannot create a Python subclass of cpp.Mesh since then isinstance(UnitSquare, Mesh) would not be true.

Is there a way around that? If not I'll have to find another way to handle this integration.

[Johan Hake]

I think we could find a way around that. Maybe by subclassing all built-in meshes, like we, to a certain extent do with specialfunctions.py.

Johan

[Martin Sandve Alnæs]

I'm not sure we want that anyway. It gets messy to subclass across ufl/dolfin interfaces and severely restricts the design freedom. In particular, to support a multidomain sublanguage in ufl there will be several domain related types in ufl. Something like (this is not the final design):

class Domain:

class DomainDiscretization(Domain): # Can reference a dolfin.Mesh

class DomainOperator(Domain):
class DomainRegion(DomainOperator):
class DomainIntersection(DomainOperator):
class DomainUnion(DomainOperator):
class DomainInterface(DomainOperator):

From one point of view it is sufficient that DomainDiscretization holds a reference to the mesh. The point is that you can get the dolfin.Mesh through a ufl.Form -> ufl.Measure -> ufl.Domain -> mesh reference.

However I want to be able to write:

vol = 1*dx(mesh)
assemble(vol)

and maybe later (with the help of Andres work on overlapping meshes):

vol = 1*dx(intersection(mesh1, mesh2))
assemble(vol)

We can achieve this through either letting dolfin.Mesh subclass ufl.DomainDiscretization, or through a looser pythonic coupling:

# Add a ufl_domain function to dolfin.Mesh:
def _ufl_domain(self):
    return ufl.DomainDiscretization(
            gdim=self.geometry().dim(), 
            tdim=gdim=self.topology().dim(),
            cell=self.ufl_cell(),
            metadata=self)
dolfin.Mesh.ufl_domain = _ufl_domain

# "Casting" dolfin.Mesh to ufl.DomainDiscretization
# without actually subclassing:
def as_domain(domain):
    if isinstance(domain, ufl.Domain):
        return domain
    elif hasattr(domain, "ufl_domain"):
        return domain.ufl_domain()
    else:
        ufl_error("Invalid domain type.")

Does this look ok to our swig expert?

[Prof Garth Wells]

Having UnitSquareMesh, etc, as a subclass of Mesh is pain - it complicates code on the C++ side in cases too.

The only difference between a Mesh and a UnitSquareMesh is the constructor. A better design would be to have mesh builder/mesh factory functions to build and return a unit foo Mesh.

[Johan Hake]

It is nice to be able to write:

UnitSquareMesh mesh(10,10);

in C++, but the implementation could be done in a factory function, which is called in the constructor of FooMesh. In Python we could then ignore the classes and just call the factory function. Maybe also rename the factory function to FooMesh so user code wont break.

[Martin Sandve Alnæs]

I made a separate issue for that.

[Martin Sandve Alnæs]

After an hour long discussion with Johan, we may want to keep the ufl side a bit simpler than sketched above with all the DomainFoo types, and rather implement easy to use set operations on MeshFunctions in dolfin. Some enhancement of the domain descriptions in ufl is needed though. I'll have to think more about this.

[Martin Sandve Alnæs]

What do people think about this syntax for integration over a particular mesh?

dx = mesh.dx()
ds = mesh.ds(facet_domains)

volume = 1*dx()
surface = 1*ds()
surface23 = 1*ds((2,3))

print assemble(volume)
print assemble(surface)
print assemble(surface23)

[Johan Hake]

I like it. What with subdomains contained in the mesh? Should we automatically generate a MeshFunction and pass that to the returned measure when:

dx = mesh.dx()

if mesh contains marked cells? Then that would be the way to populate the form with mesh-stored domains. One can then over rule that by:

dx = mesh.dx(some_other_domain)

[Martin Sandve Alnæs]

Yes, that makes sense. Loading a mesh with markers from file could then be used like this:

mesh = Mesh("mymesh.xml.gz")
dx = mesh.dx()
ds = mesh.ds()

form = 1*dx(2) + 3*ds((4,5))
print assemble(form)

Later maybe we could get this working with names as well:

form = 1*dx("left") + 3*ds(("wall", "inlet"))

[Prof Garth Wells]

Looks ok. As long as using domain data stored in a mesh file is explicit, otherwise a user can get in a tangle if the mesh file unknowingly contains domain data.

[Martin Sandve Alnæs]

Assembly with integrals on subdomains and no subdomains provided should be an error. The old 0=everywhere behaviour is gone so this should be fine.

[Anders Logg (Chalmers)]

dx = mesh.dx() looks good.

[Martin Sandve Alnæs]

Garth: I don't agree, what tangle? My example above your last post does exactly what you don't want.

Basically you have two situations:

• A form contains only integrals over everywhere. Then subdomains are not used.
• A form contains integrals over subdomains. Then surely the user expects subdomains to be used.

[Martin Sandve Alnæs]

I suspect you're thinking about the old situation where f*dx meant f*dx(0), but that is no longer a problem since f*dx now means integration over the entire domain.

[Prof Garth Wells]

What I'm saying is that a mesh file should not contain domain data that is magically used in the assembly without me as a user knowing that it's being used.

[Martin Sandve Alnæs]

I don't see the problem, but ok. I can do this:

mesh = Mesh("mymesh.xml.gz")
cell_domains = MeshFunction("size_t", mesh, d, mesh.domains())
facet_domains = MeshFunction("size_t", mesh, d-1, mesh.domains())
dx = mesh.dx(cell_domains)
ds = mesh.ds(facet_domains)

form = 1*dx(2) + 3*ds((4,5))
print assemble(form)

It's often useful to have the domains as explicit objects anyway.

[Johan Hake]

I also do not see the problem. With some swig magic we could facilitate:

mesh = Mesh("mymesh.xml.gz")
dx = mesh.domains().dx()
ds = mesh.domains().ds()

form = 1*dx(2) + 3*ds((4,5))
print assemble(form)

Then it is more explicit that the measure includes domain data.

[Prof Garth Wells]

It's not a problem but a prerequisite. Both of the above two code snippets look fine to me.

[Johan Hake]

Or we could hide the generation of the MeshFunction by:

mesh = Mesh("mymesh.xml.gz")
dx = mesh.dx(mesh.domains())
ds = mesh.ds(mesh.domains())

[Martin Sandve Alnæs]

I'm working on something similar, for now assigning myself here to keep track of the discussion in this issue.

[Martin Sandve Alnæs]

Issue #127 was marked as a duplicate of this issue.

[Martin Sandve Alnæs]

Duplicate of #190.

[Martin Sandve Alnæs]

Note: I consider this issue closed, but keep in mind that everything above is a discussion and not exactly how it turned out.