Facet markers do not transfer after SubMesh

[Navid Mirzaei]

When created a SubMesh from a mesh generated by GMSH the facet markers do not work any more. Look at the code below:

from dolfin import *
import numpy as np

Domain=Mesh('eccentric_domain_3D_NC.xml')
Volume = MeshFunction('size_t' , Domain , 'eccentric_domain_3D_NC_physical_region.xml' )  
Boundary = MeshFunction('size_t', Domain , 'eccentric_domain_3D_NC_facet_region.xml')




d_L = Constant(10)
d_H = Constant(10)
D_L = Constant(10000)
D_H = Constant(10000)
mu_L = Constant(100000)
sigma_L = Constant(1000)
alpha_L = Constant(0.1)
gamma_L = Constant(0.1)
nu_H = Constant(10)
kappa = Constant(1)
sigma_H = Constant(1000)

Sub_Mesh=SubMesh(Domain,Bulk,7)
submesh_regions = CellFunction('size_t', Sub_Mesh)
submesh_boundaries = FacetFunction('size_t', Sub_Mesh)
ds = Measure('ds', subdomain_data=submesh_boundaries)
dx = Measure('dx', subdomain_data= submesh_regions)

P1 = FiniteElement('P', tetrahedron, 1)
element = MixedElement([P1, P1])
Mixed_Space = FunctionSpace(Sub_Mesh, element)
v1 , v2=TestFunction(Mixed_Space)
U = Function(Mixed_Space)
L, H= split(U)

F = d_L*L*v1*dx+mu_L*((L*Constant(2))/(Constant(1)+L))*v1*dx-D_L*sigma_L*((Constant(1)+H/alpha_L)/(Constant(1)+H/gamma_L))*v1*ds(3)+D_L*dot(grad(L),grad(v1))*dx(7)+d_H*H*v2*dx+nu_H*(H*Constant(3)/(kappa+H))*v2*dx-D_H*sigma_H*v2*ds(3)+D_H*dot(grad(H),grad(v2))*dx

solve(F == 0, U, solver_parameters={"newton_solver": {"relative_tolerance": 5e-10,   
   "absolute_tolerance": 5e-10,"maximum_iterations": 10}}) 

L_,H_ = U.split()

The .xml files are attached.