Interface between two domains

Question

Hello,

I wonder if there is a more pythonish fenicsy way of marking the edges that make up the interface between two different subdomains, using only topology, because the geometry is unknown a priori.

Proving once again that you can write Fortran in any language, I came up with:

mesh.init()
for f in facets(mesh):
 neighs =  f.entities(d)
 if neighs.size == d:
    if submarker[neighs[0]] != submarker[neighs[1]]:
       interface.array()[f.index()] = 1

Even though this works, I know this can't be right :)
In particular, the loop of all internal facets is pretty horrific. I know about BondaryMesh, but I couldn't figure out how to get the mesh corresponding to each subdomain.

Comments

Rather late comment, but do you have a fully working minimal example to play with by any chance?