Hi,

I used hexahedral elements to mesh a spherical shell (I uploaded the mesh file to github). There was a problem adding MeshTags to boundaries as suggested by Dokken. Some boundary facets are not marked.

The code and result are as follows:

```
boundaries = [(1, lambda x: np.isclose(x[0]**2+x[1]**2+x[2]**2, 1)),
(2, lambda x: np.isclose(x[0]**2+x[1]**2+x[2]**2, 0.35**2))]
facet_indices, facet_markers = [], []
fdim = msh.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(msh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tag)
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
with io.XDMFFile(msh.comm, "facet_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(msh)
xdmf.write_meshtags(facet_tag, msh.geometry)
```

There is no problem with tetrahedral elements. So I want to know how to create identifying the facets for hexahedral elements correctly?

Thanks!