Hello,

I am loading my mesh into FEniCS through .xml files, after creating them with dolfin-convert from a GMSH .msh file. As an example, here the mesh is a square and I would like to solve the Poisson equation in its left half. I do this by marking the desired subdomain an then creating a submesh on it. The subdomain is marked with the following code (as the original rectangle was identified by marker 1):

```
subdomains_left = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomains_left.array()[subdomains.array()==1] = 1
```

As I need to apply Dirichlet boundary conditions in some boundaries of this submesh, I would need to preserve the same markers for the boundaries as I had on the original mesh. So I expect that the previous procedure for the subdomains could also be used for “subboundaries”. Say I want to impose boundary conditions on the original boundary marked as 3. Therefore, my aim is to implement something like this (that is not working):

```
boundaries_left = MeshFunction('size_t', mesh_left, mesh_left.topology().dim()-1)
boundaries_left.array()[boundaries.array()==3] = 1
```

I post my actual code, which is working but without aplying the Dirichlet BC that I want.

```
from dolfin import *
# Read mesh, subdomains and boundaries created with GMSH
mesh = Mesh('square.xml')
subdomains = MeshFunction('size_t', mesh, 'square_physical_region.xml')
boundaries = MeshFunction('size_t', mesh, 'square_facet_region.xml')
# Mark left rectangle (with label 1) from original subdomains
subdomains_left = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomains_left.array()[subdomains.array()==1] = 1
# Create submesh for left rectangle
mesh_left = SubMesh(mesh, subdomains_left, 1)
# Define finite element space for submesh
V_left = FunctionSpace(mesh_left, "CG", 1)
# Define variational problem
u = TrialFunction(V_left)
v = TestFunction(V_left)
f = Expression('-2', degree=0)
a = -dot(grad(u), grad(v))*dx
L = - f*v*dx
# Boundary conditions
bc = []
# Compute solution
u = Function(V_left)
solve(a == L, u, bc)
# Save solution to file in XDMF format
with XDMFFile(mesh_left.mpi_comm(), 'u.xdmf') as XFile:
XFile.write_checkpoint(project(u, V_left), 'u', 0)
# I would need to apply Dirichlet boundary condition on specific
# boundaries of the original mesh. For instance, DirichletBC
# should be used to impose sol_teo on boundary marked as 3 in
# the original mesh
# Something like:
#boundaries_left = MeshFunction('size_t', mesh_left, mesh_left.topology().dim()-1)
#boundaries_left.array()[boundaries.array()==3] = 1
#sol_teo = Expression('pow(x[0],2)', degree=2)
#bc = DirichletBC(V_left, sol_teo, boundaries_left, 1)
```

I also attach my square.geo file:

```
//+
SetFactory("OpenCASCADE");
Rectangle(1) = {0, 0, 0, 1, 2, 0};
//+
Point(5) = {2, 0, 0, 1.0};
//+
Point(6) = {2, 2, 0, 1.0};
//+
Line(5) = {3, 6};
//+
Line(6) = {6, 5};
//+
Line(7) = {5, 2};
//+
Curve Loop(3) = {2, 5, 6, 7};
//+
Plane Surface(3) = {3};
//+
Physical Surface("Left_surface") = {1};
//+
Physical Surface("Right_surface") = {3};
//+
Physical Curve("Left") = {4};
//+
Physical Curve("Top_left") = {3};
//+
Physical Curve("Top_right") = {5};
//+
Physical Curve("Right") = {6};
//+
Physical Curve("Bottom_left") = {1};
//+
Physical Curve("Bottom_right") = {7};
//+
Physical Curve("Middle") = {2};
//+
MeshSize {4, 3, 1, 2} = 1;
```

I would appreciate any help in order to solve my problem. Thanks in advance.