Submesh from boundarymesh

Question

How would one properly define a submesh of a boundarymesh, when the boundaries are given by a meshfunction?

from dolfin import *
mesh = Mesh("mesh.xml")
boundaries = MeshFunction('size_t', mesh, "mesh_facet_region.xml")
boundarymesh = BoundaryMesh(mesh, 'exterior')

submesh = SubMesh(boundarymesh, boundaries, 0) # does not work

When using as subdomain to define the submesh in the last line, there is no problem (see also http://fenicsproject.org/qa/7039/submesh-of-a-boundarymesh).

Answer 1

If you are not worried about solving this by using the tags in your mesh/geo file and just want to get done with your problem, then something like this should work:

mesh = Mesh("2DSquare.xml")
o = BoundaryMesh(mesh, 'exterior')
print o
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
top = Top()
submesh = SubMesh(o, top)
print submesh

However, if your mesh is complex or you specifically want to use the tags, then these links might be useful for you:

1) https://answers.dogfood.paddev.net/dolfin/+question/186727

2) https://bitbucket.org/fenics-project/dolfin/issues/156/segmentation-fault-in-boundingboxtree2d (he even has the .xml and physical.xml attached)

3) https://answers.launchpad.net/dolfin/+question/214113

Let us know if your code works.

Comments

Thank you for your quick reply. However, it does not answer my question (or I fail to see how).

I have been using the subdomain approach until now when defining my submesh, and that indeed works. But now I need to generalize to complex meshes. I created a meshfunction that contains and numbers the boundaryfacets, and the goal now is to derive a mesh of one dimension lower corresponding to such a boundary. There are no functionspaces defined yet at this point.

In the last comment from the link I posted, there might be a clue:

submesh = SubMesh(boundarymesh, boundaries, 0)

might not work because boundaries is defined on the original mesh, i.e. not on the boundarymesh. But that leaves me with the question how to make it work then?

---

Ok. I thought you just want to understand how to work with the tags, but you specifically want to do it through 'boundaries'. I tried some stuff. Didn't work.

Awaiting an expert.

Answer 2

The following seems to work. I am sure the kind FEniCS developers provided us with a shorter way to do the following, but I am not aware of it yet.

from dolfin import *
mesh = Mesh("mesh.xml")
boundaries = MeshFunction('size_t', mesh, "mesh_facet_region.xml")
boundarymesh = BoundaryMesh(mesh, 'exterior')

# copy the meshfunction, but defined on the boundarymesh this time
bdim = boundarymesh.topology().dim()
boundary_boundaries = MeshFunction('size_t', boundarymesh, bdim)
boundary_boundaries.set_all(0)
for i, facet in enumerate(entities(boundarymesh, bdim)):
    parent_meshentity = boundarymesh.entity_map(bdim)[i]
    parent_boundarynumber = boundaries.array()[parent_meshentity]
    boundary_boundaries.array()[i] = parent_boundarynumber

submesh = SubMesh(boundarymesh, boundary_boundaries, 0) # works

I only tested this briefly, but it seems to work. If it turns out to be wrong I will change this answer accordingly. If anyone comes along with a shorter solution, I'd be happy to mark it as an answer.