How to define BCs on boundaries of submeshes of a parent mesh?

Hi there! I am experiencing problem defining BCs on submeshes extracted from a parent mesh.
My geometry consists of two concentric circles. I am solving one PDE for the smaller disk and another for the annulus. I successfully recovered the submeshes but having trouble defining the boundaries to my submeshes. Below is the gmsh script followed by the Python script.

gmsh script:

Circle(1) = {0, 0, 0, 5, 0, 2*Pi};
Circle(2) = {0, 0, 0, 3, 0, 2*Pi};
Curve Loop(1) = {2};
Plane Surface(1) = {1};
Curve Loop(2) = {1};
Curve Loop(3) = {2};
Plane Surface(2) = {2, 3};
Physical Curve(1) = {2};
Physical Curve(2) = {1};
Physical Surface(5) = {1};
Physical Surface(4) = {2};

Python script:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
    return out_mesh

import meshio
msh ="example.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)

s5 = MeshView.create(cf, 5) # Extracting subspace 5
s4 = MeshView.create(cf, 4) # Extracting subspace 4

W1 = FunctionSpace(s5,"Lagrange",1) # Defining function space on s5
W2 = FunctionSpace(s4,"Lagrange",1) # Defining function space on s4

def bs5(x):  # Defining boundary for subspace s5, curve 1
  return mf == 1 
def b2s4(x): # Defining boundary for subspace s4, curve 2
  return mf == 2
# Define boundary condition
G , mu = 1, 0.1

bc = DirichletBC(W1, u_D, bs5)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
us = Function(W1)
solve(a == L, us, bc)

# Save solution to file in PVD format for Paraview
file = File("examplelumenvelocity.pvd");
file << us;


D, k, t = 0.7, 0.3, 3
f0, f1, beta = 5, 2, 0.3
gamma1 = (-1/D)*Constant(f0 + f1*(1-exp(-beta*(t**2)))) # Neumann condition on inner curve -dot(grad(c),n)=gamma1
gamma2 = Constant(0.0)    # Dirichlet condition on outer circle 

bc2 = DirichletBC(W2, gamma2, b2s4)


# Define variational problem
w = TestFunction(W2)
c = TrialFunction(W2)
a = inner(grad(c), grad(w))*dx + (((k/D)*c)*w)*dx
a -= - inner(gamma1, w)*ds_custom(bs5)

# Compute solution
cs= Function(W2)
solve(a == 0, cs, bcs)

# Save solution to file in PVD format for Paraview
file = File("example.pvd");
file << cs;

Error message:

*** Warning: Found no facets matching domain for boundary condition.
Invalid subdomain_id <function bs5 at 0x7fc435e3f1f0>.

How to fix this? I am using Dolfin version 2019.2.0.dev0.

I modified the BCs as

bc = DirichletBC(W1, u_D, bs5, method="pointwise") and similarly, the other one as
bc2 = DirichletBC(W2, gamma2, b2s4, method="pointwise")

Now, the error reads

ufl.log.UFLException: Invalid subdomain_id <function bs5 at 0x7f91ac0111f0>.

Does importing the parent mesh from gmsh and then extracting the subdomains mess up the physical tags or is there something wrong with my code? Any suggestion would be extremely helpful.

I think I had a similar issue. To my knowledge there is no direct way unfortunately to convert a MeshFunction defined on the parent mesh to the submesh. At least I could not find an elegant solution.
Below a code, that does it manually: not the fastest but it works.
It basically iterates over every cell/facet to copy the names of the subdomains and facets defined on a MeshFunction on your parent mesh to a new MeshFunction which is defined on your new mesh of the subdomain defined via MeshView.create(…).

def extract_SubMesh_FacetDomain(mesh, mesh_Subdomain, mf_2d, mf_1d, indices):
        Extracts the submesh with the indices 'indices' from mesh 
        similar to
            - mesh: mesh of the system
            - mesh_Subdomain: MeshView object of the subdomain ( created with extract_SubMesh(mesh, mf_2d, indices) )
            - mf_2d: MeshFunction labeling the subdomains of the mesh
            - mf_1d: MeshFunction labeling the facets of the mesh
            - indices: infices of the subdomains stored in mf_2d that should be extracted

            - submesh_boundaries: facets of the submesh extracted from mf_1d
            - submesh_domains:  subdomains of the submesh extracted from mf_2d
    # Step 1: Create a submesh with SubMesh(). This is needed because I need to access the 'parent_vertex_indices' and 'parent_cell_indices', which apparently is not possible via MeshView.create (to my knowledge)
    subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0) 
    for k in indices:
        subdomains.array()[mf_2d.array()==k] = 1

    temp_SubMesh = SubMesh(mesh,  subdomains, 1)
    # Step 2: Iterate over every cell and copy the indices from mf_2d/mf_2d in their new equivalent for the subdomain
    submesh_domains = MeshFunction('size_t', mesh_Subdomain, mesh_Subdomain.topology().dim(), 0)
    submesh_boundaries = MeshFunction('size_t', mesh_Subdomain, mesh_Subdomain.topology().dim()-1, 0)
    vmap ='parent_vertex_indices', 0)
    cmap ='parent_cell_indices', 2)
    for submesh_cell in cells(mesh_Subdomain): # loop for copying the names of the subdomains
        parent_cell = cmap[submesh_cell.index()]
        submesh_domains.array()[submesh_cell.index()] = mf_2d.array()[parent_cell]
    for submesh_facet in facets(mesh_Subdomain): # loop for copying the names of the facets (boundaries)
        submesh_facet_vertices = vmap[submesh_facet.entities(0)]
        for facet in facets(mesh):
            if set(facet.entities(0)) == set(submesh_facet_vertices):
                submesh_boundaries.array()[submesh_facet.index()] = mf_1d.array()[facet.index()]

    return submesh_boundaries, submesh_domains