A workaround would be to adapt this to:
if len(facets) == 0:
indices = np.empty([], dtype=np.int32)
else:
indices = mesh.entities_to_geometry(domain, domain.topology.dim-1, facets, False)
indices = np.unique(indices)
However, the code:
is not really safe in its current status, as you have some heavy implicit assumptions on the mesh geometry and the function space.
I am assuming you are using a first order Lagrange space for your computations?
If so, one should rather use a vertex_to_dofmap construction, as for instance shown in Extract the coordinates of the boundary points and the value of a function on them - #4 by dokken