I've been reading through the documentation and I haven't been able to successfully solve the test problem I've been trying to implement using a mesh I generated in GMSH. Rather than renaming the surfaces and volumes in the .geo file, I want to do this post-mesh generation and conversion using dolfin-convert. Is it ok to rename these subdomains and boundaries and define the surface/volume differentials as I have below? I'm trying to determine if this is the issue in my code or if I must look elsewhere for additional problems. If it helps, my domain has two Dirichlet BC's on domain 34 (30 and 31 from the .msh file) and an interior interface between subdomains 33 and 34 with a flux continuity Neumann BC on it. Subdomain 33 is entirely encompassed by subdomain 34.
mesh = Mesh(ofilename)
subdomains = MeshFunction("size_t", mesh, "%s_physical_region.xml"%ofilename.split('.')[0])
boundaries = MeshFunction("size_t", mesh, "%s_facet_region.xml"%ofilename.split('.')[0])
fix_subdomains = CellFunction("size_t", mesh)
fix_subdomains.set_all(0)
fix_subdomains.array()[subdomains.array()==34] = 1
fix_subdomains.array()[subdomains.array()==33] = 2
fix_boundaries = FacetFunction("size_t", mesh)
fix_boundaries.set_all(0)
fix_boundaries.array()[boundaries.array()==30] = 1
fix_boundaries.array()[boundaries.array()==31] = 2
fix_boundaries.array()[boundaries.array()==32] = 3
fix_interface = FacetFunction("size_t", mesh)
fix_interface.set_all(0)
fix_interface.array()[boundaries.array()==32] = 1
dx = Measure("dx", subdomain_data=fix_subdomains)
ds = Measure("ds", subdomain_data=fix_boundaries)
dS = Measure("dS", subdomain_data=fix_interface)