SubMesh of a BoundaryMesh

Question

Hi all,

I'm trying to create a SubMesh of a BoundaryMesh. For instance, given a UnitSquareMesh, I would like to get only the vertices on the top boundary of the mesh. I've tried this code:

from dolfin import *

mesh = UnitSquareMesh(10, 10)

class OmegaTop(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] == 1
omega_top = OmegaTop()

subdomains = CellFunction('size_t', mesh)
subdomains.set_all(0)
omega_top.mark(subdomains, 1)

submesh_omega_top = SubMesh(mesh, subdomains, 1)

However, the output is:

 In  [1]: print submesh_omega_top
 Out [1]: <Mesh of topological dimension 2 (triangles) with 0 vertices and 0 cells, ordered>

The point is, that I don't understand, why the SubMesh has no vertices and no cells? I've expected an intervall with 11 vertices and 10 cells. How can I correct the code to get only the vertiecs on the top boundary of the UnitSquareMesh?

By the way, using the folling code leads to the same result:

submesh_omega_top = SubMesh( BoundaryMesh(mesh, "exterior") , subdomains, 1)

Thank you for your support in advance.

Answer 1

Two errors:

First, you need to first extract a BoundaryMesh so that you get a mesh of one dimension lower (in this case 1D).

Second, use near(x[1], 1) rather than x[1] == 1 since round-off errors may otherwise prevent those points from being filtered out.

Comments

Thank you very much. The code

submesh_omega_top = SubMesh(BoundaryMesh(mesh, 'exterior'), omega_top)

provides

<Mesh of topological dimension 1 (intervals) with 11 vertices and 10 cells, ordered>

(even without the use of near(x[1], 1) ) and that's exactly what I wanted. I guess

submesh_omega_top = SubMesh( BoundaryMesh(mesh, "exterior") , subdomains, 1)

was wrong, because the CellFunction is defined on the mesh and not on the BoundaryMesh.