[petsc-users] Pulling submeshes out of a DMPlex object

[Lawrence Mitchell]

Dear all,

for various reasons I'd like, given an existing DMPlex mesh, to be
able to pull submeshes out.  These might take various forms:

1. The locally visible part of a distributed mesh (so that I can
discretise and solve independent problems on each local part)

2. Some subset of entities of a given codimension (for example, the
boundary of a domain, or the skeleton of the mesh).

So I think conceptually what I want is:

Create a DMLabel that marks entities of a given codimension.

Say

DMPlexCreateSubmesh(dm, label, value, height, comm, &subdm)

This should walk all entities of the specified height and gather them
if they are labelled with the provided value, then the subdm is
created by treating those entities as the cells of a new mesh and
adding their transitive closure.

There appears to be some submesh functionality in plex right now, but
it doesn't appear to do what I want:

  DMPlexCreateSubmesh - Extract a hypersurface from the mesh using
vertices defined by a label
  DMPlexFilter - Extract a subset of mesh cells defined by a label as
a separate mesh

The latter looks close to what I want for 1, except the change of
communicator.

The former seems somewhat odd.  I would have thought the way to
extract a hypersurface would be to mark the faces of the hypersurface
and build a mesh as a I described above.

It looks like it should be possible to reuse some of the existing
submesh routines to do what I want to do, but I'm a little lost.

What I tried was just adding an option to
DMPlexCreateSubmeshGeneric_Interpolated to mark a submesh by
traversing the closure of points marked with a given value.

But that swiftly ran into problems because the resulting submesh
didn't have depth strata as expected, and looking at the existing code
I'm quite confused.  The subpointMap that a submesh DM has labels the
points of the parents with their depth (in the original mesh?).  But
this seems to be used to determine the depth of the points in the
submesh.  Which seems odd to me.

I realise this is all rather woolly right now, perhaps in attempting
to clarify any queries I will explain better what I want to do.

Cheers,

Lawrence