I use something like that in my script. I start off with a scalar function that labels the areas that I want to have a given constant vector. The mesh is created externally, but basically just labels material 1 with a 1, material 2 with a 2, material 3 with a 3. In the snippet below, I am trying to create the constant vector field for material 3:
D = FunctionSpace(mesh, "DG", 0)
cell_labels = MeshFunction("size_t", mesh, "material_type.xml.gz")
helper = numpy.asarray(cell_labels.array(), dtype=numpy.int32) - 1
dm = D.dofmap()
for cell in cells(mesh):
helper[dm.cell_dofs(cell.index())] = cell_labels[cell] - 1
material_3_location = Function(D)
material_3_location.vector()[:] = numpy.choose(helper, [0., 0., 1.])
You can then make your piecewise constant function by projecting or interpolating.
V = FunctionSpace(mesh, "N1curl", 1)
V3 = project(Constant((0., 1., 0.))*material_3_location, V)