Domain handling for forms with multiple domains

[Lawrence Mitchell]

We're developing a coastal flow solver and want an ALE formulation. It turns out the natural way to do this in the timestepping is to have a single mesh topology and two coordinate fields (effectively two UFL domains). Since the coordinate field is obtained from the UFL domain object most things fall through. There are a few minor glitches which I'd like to iron out here.

First the proposed usage:

mesh1 = Mesh(...) # Get this mesh from anywhere.
mesh2_coords = Function(VectorFunctionSpace(mesh1, "CG", 1))
mesh2_coords.interpolate(SpatialCoordinate(mesh1)) # start out the same values

V = FunctionSpace(mesh1, ...)
V2 = FunctionSpace(mesh2, ...)

u = TestFunction(V)
u2 = TestFunction(V2)

F = u*dx(domain=mesh1) + u2*dx(domain=mesh2)

The intention here is that the local kernel is the same, the domain just indicates which coordinate field should be used to pull back to reference space.

This currently doesn't work, because the form has two test functions defined on different domains (and UFL doesn't like this), particularly the problem arises from analysis.extract_coefficients_and_arguments.

An alternative way of doing this, which nearly works in all cases is this:

F = u*dx(domain=mesh1) + u*dx(domain=mesh2)

That is, only ever use test and trial functions from one of the domains, the domain argument to the measure indicates where to get the coordinate from.

This works, except if I just write:

F = u*dx(domain=mesh2)

Because the domains of the form are only determined from the measure, and used to calculate a domain renumbering. But u has a different domain. I note the following TODO in form.py:

    def _analyze_domains(self):
        from ufl.domain import join_domains, sort_domains

        # Collect unique integration domains
        integration_domains = join_domains([itg.ufl_domain() for itg in self._integrals])

        # Make canonically ordered list of the domains
        self._integration_domains = sort_domains(integration_domains)

        # TODO: Not including domains from coefficients and arguments here, may need that later
        self._domain_numbering = dict((d, i) for i, d in enumerate(self._integration_domains))

If I just incorporate the domains from coefficients and arguments into the domain_numbering then everything appears to work.

Thoughts?

[Martin Sandve Alnæs]

I'm not sure, I've burned my fingers on this before and I'm worried it might break something. In particular dolfin does funky stuff with functions from other meshes that is not grounded in the intentions in ufl.

See this issue:

https://bitbucket.org/fenics-project/uflacs-deprecated/issues/49/representation-fails-with-coefficient-from

[Lawrence Mitchell]

Hmmm. That's a tricky one. The ALE case is a slightly weird one since you have the same mesh topology you just want to use a different coordinate field (so you're guaranteed to have matching meshes). Another thing we'd consider is being able to provide to the measure the coordinate field to use for integration (but this just looks like a special case of separate domains). Another option might be to introduce a terminal modifier that pushes an expression forward from one domain to another. Then you could write:

stuff_on_mesh1*dx(domain=mesh1) + PushForward(other_stuff_on_mesh1, mesh1, mesh2)*dx(domain=mesh2)

Which is just explicitly annotating all the terminals that their domain has changed. I think this approach shouldn't break existing code, because there are no existing uses. And should allow the domain type-checking to work.