Implementation of subdomains

[wence-]

Implementing Subdomains

To support solutions of problems where we solve for some set of variables in one region of the domain and a different set in another, for example mesh-aligned FSI, we need to be able to mark regions of
the mesh and define FunctionSpaces (and hence Functions) just on those regions.

Why not mesh coupling?

Issue #714 has more details on surface coupling of two meshes. The difficulty, it turns out, is that the interpolation operator is not the identity for all but the most trivial cases. Even when the meshes are topologically conforming. Since our current needs are for topologically conforming meshes, subdomains are an easier solution. When we do want to revisit surface coupling of different meshes the topologically conforming case will fall out as a special case of the general non-conforming interface.

Data structures

We start with a single Mesh and need a way to create subdomains in it. External mesh generators probably have a way to mark the cells with some ID, we should ensure that the mesh reader can handle these and provide them to use appropriately.

A SubDomain

A subdomain collects a set of cells in the mesh. When constructing function space data structures by interfacing with the DM, it needs to ensure iteration only occurs over the selected cells. The correct
underlying data structure is therefore probably a DMLabel.

DOLFIN already has an API for defining subdomains, we should check what it looks like, I think they start with a MeshFunction which is used to mark entities.

Subdomains should be allowed to overlap.

Function spaces

It must be possible to build a FunctionSpace on a SubDomain. The resulting space should have no degrees of freedom outside the domain, and entity->node maps should similarly only map from the entities in the subdomain.

This immediately means that we can't represent subdomains using the current pyop2 Subset functionality (which just allow iteration of a subset of a full iteration set by explicitly enumerating the selected indices).

Subdomains therefore need to behave basically like Mesh objects.

Functions

If the FunctionSpace is implemented correctly, these just fall out.

Integration

Integration of functions defined on a single subdomain should be straightforward. A few issues arise. What should the dS and ds integrals be? Should they refer to the global mesh's interior and
exterior facets, or should they be subdomain-specific. In other words, should ds be the portion of the subdomain facets that are exterior in the global mesh, or should it be the boundary of the subdomain?

Integration on overlaps

We will need a way of performing integrals over the overlaps between two subdomains. We therefore probably want a way of expressing this. Maybe something like:

overlap = Intersection(a, b)

assemble(f*dS(overlap))

Integrates over all interior facets that are in the intersection of subdomains a and b. This suggests that the answer to the previous question of what dS and ds mean is that they should be the same as in the global mesh.

Integration over the global domain

This should presumably become illegal for coefficients and arguments defined on subdomains. Does this become a problem?

Support in PyOP2

I think we don't need any new functionality here, since the subdomains are going to be responsible for creating the correct PyOP2 Sets and Maps.

[wence-]

This is intended as a discussion to iterate on to get the design and API of subdomains correct before implementing things.

[pefarrell]

Does this plan support subdomains of different topological dimension (e.g. defining a Function on a part of the 2D boundary of the 3D volume mesh)?

A subdomain collects a set of cells in the mesh

Shouldn't a subdomain be a set of geometric entities (cells, facets, lines, etc)?

[wence-]

I think if we get it right, yes.

[krober10nd]

Hello, Is there any progress on this? I'm having difficulty prescribing surfaces in a msh file not created from gmsh. I'm attempting to do this with custom scripts in Python, but so far no success. In general for geophysical applications, gmsh isn't always the best route for mesh generation since you need to have a fantastic boundary description to mesh. Is there an alternative way to prescribe boundary conditions inside FireDrake given that I know which points are boundaries and their locations in space?