Project/interpolate constants to different mesh subdomains

I have a Dolfin mesh that defines two or more material/sub-domains like:

mesh = Mesh("test.xml")
subd = MeshFunction("size_t", mesh, mesh.topology().dim(), mesh.domains() )
dx = Measure("dx", domain=mesh, subdomain_data=subd)

E0 = FiniteElement("P", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, E0)
u = Function(V)
v = TestFunction(V)

I can use this (dx measuse) to define PDE forms in different domains like:

c1*u*v*dx(0) + c2*u*v*dx(1)

However, I would like to interpolate/project different constants to
the different subdomains for the initial solution vector, something like

u0 = [ interpolate(Constant(0.0), V.dx(0)), \
       interpolate(Constant(1.0), V.dx(1)) ]

Is this possible somehow (as interpolate(Constant(1.0),V) affects everything)?

See How to define different materials / importet 3D geometry (gmsh)