The issue is that when refining the mesh, new facets are introduced into the mesh (that had no previous parent facet, as they are internal to a cell, these facets get large numbers as values.
I would like to note that your example is a very good one (good way of presenting a minimal working example), but I have simplified it slightly for easier visual interpretation:
from dolfin import *
parameters["refinement_algorithm"] = "plaza_with_parent_facets"
class Top(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and x[1] > 1.0-1e-10
mesh = UnitSquareMesh(2,2)
d = mesh.topology().dim()
boundaries = MeshFunction("size_t", mesh, d-1)
boundaries.set_all(0)
top = Top()
top.mark(boundaries,1)
mesh1 = adapt(mesh)
boundaries1 = adapt(boundaries, mesh1)
File("cube_boundaries.pvd") << boundaries
File("cube_boundaries_refined.pvd") << boundaries1
In the following picture, the left figure is the adapted mesh (with the value==1 marked as a bolder line), and the RHS the original mesh (with similar bold marking).
As you can see, facets that original belonged to a facet has kept its marker, but new ones have gotten an arbitrary large value.
