Defining a subdomaim from cell data

I would suggest using meshio, as for instance explained in this post: Accessing and marking imported boundaries - #8 by dokken
Here, you can replace: cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
with the read in array (assuming that it is ordered in the same way as the VTK file and cells), and convert it to XDMF (including the cell markers), which in turn can be read into Dolfin as shown for instance in:
Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #3 by dokken