Preserve original mesh boundaries on submesh

mesh
1

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.

2

Hi! If you are creating and tagging submeshes in gmsh itself, you can try the following commands in your script after transporting the mesh. These commands preserve all the submesh (these are 2D submeshes) boundaries.

Sub_Mesh=SubMesh(mesh,cf,5)  # Extracting a submesh tagged 5 in gmsh
boundary_marker = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim()-1, 0)
ncells = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim())
surface_marker = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim(), 0)
vmap = Sub_Mesh.data().array("parent_vertex_indices", 0)
cmap = Sub_Mesh.data().array("parent_cell_indices", Sub_Mesh.topology().dim())

for c in cells(Sub_Mesh):
  parent_cell = Cell(mesh, cmap[c.index()])
  surface_marker.array()[c.index()] = cf.array()[parent_cell.index()]
  for f in facets(parent_cell):
    for g in facets(c):
      g_vertices = vmap[g.entities(0)]
      if set(f.entities(0)) == set(g_vertices):
        boundary_marker.array()[g.index()] = mf.array()[f.index()]   
ds = Measure('ds', subdomain_data=boundary_marker)
dx = Measure('dx', subdomain_data=surface_marker)

submesh=SubMesh(mesh,cf,4) # Extracting submesh tagged 4 in gmsh
boundary_marker = MeshFunction("size_t", submesh, submesh.topology().dim()-1, 0)
ncells = MeshFunction("size_t", submesh, submesh.topology().dim())
surface_marker = MeshFunction("size_t", submesh, submesh.topology().dim(), 0)
vmap = submesh.data().array("parent_vertex_indices", 0)
cmap = submesh.data().array("parent_cell_indices", Sub_Mesh.topology().dim())
        
for c in cells(submesh):
  parent_cell = Cell(mesh, cmap[c.index()])
  surface_marker.array()[c.index()] = cf.array()[parent_cell.index()]
  for f in facets(parent_cell):
    for g in facets(c):
      g_vertices = vmap[g.entities(0)]
      if set(f.entities(0)) == set(g_vertices):
        boundary_marker.array()[g.index()] = mf.array()[f.index()] 
ds1 = Measure('ds', subdomain_data=boundary_marker)
dx1 = Measure('dx', subdomain_data=surface_marker)