Writing results only from a partial region in the mesh

Hi all,

I am solving a two field problem (one vector field, one scalar field) in 3D in parallel. The geometry has two materials in it as shown in the attached picture (PDMS and glass).
problem%20figure

Now what I want to do at the output stage is, only write the field data from the PDMS region, to be read in paraview so that glass region appear as holes. I am trying to do this in the following way:

mesh=Mesh("mesh.xml") #mesh defined over the complete geometry: PDMS+glass
Y = FunctionSpace(mesh, "CG", 1)   
z  = Function(Y)   #this is the actual field I solve for

mesh_PDMS=Mesh("mesh_0.xml")   #mesh defined only over PDMS region
YR = FunctionSpace(mesh_PDMS, "CG", 1)
zR=Function(YR)

**Solve for z defined in Y**

z.set_allow_extrapolation(True)
assign(zR, project(z,YR, solver_type='cg', preconditioner_type='jacobi'))

file_results = XDMFFile( "result.xdmf" )
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(zR,step)

So in short, I am using the projection function (i have tried interpolation function too) to extract the values of the field z only in the region containing PDMS. But this process is not working properly. I am not getting the correct values. What is it that I might be doing wrong? And is there a better way to do this?

Thanks,
Aditya