How do I check if a Facet is part of the boundary?

Question

I would like to define a function that's zero everywhere except on the boundary of my mesh. How do I do this? What I have so far is:

mesh = Mesh('mesh.xml')
bmesh = BoundaryMesh(mesh, 'exterior')
my_function = FacetFunction('int', mesh)
my_function.set_all(0)
for facet in facets(mesh):
    if facet in bmesh: # What's the proper check here?
        my_function[facet] = 42

Not surprisingly, the "facet in bmesh" expression doesn't work. But what do I have to use instead?

Answer 1

Found the answer in question 927: the function that I need is facet.exterior().

Still, if someone can tell me the route I should have followed to get this from the documentation, please go ahead. I certainly feel a bit lost..

Answer 2

The following might be faster than the solution from the first answer:

mapping = bmesh.entity_map(2)
boundary_facets = [ Facet(mesh, mapping[cell.index()])
                    for cell in cells(bmesh) ]
for facet in boundary_facets:
    my_function[facet] = 1