Defining boundaries on gmsh subdomains

mesh
1

I have a problem where I am considering multiple domains. I am attempting to formulate the problem in one of the subdomains by conveniently pulling mesh, subdomain and boundary data from .xml files derived from my .msh. I define the submesh and subdomain without a problem, but when trying to define the boundaries of my submesh, I get the following error:

IndexError: boolean index did not match indexed array along dimension 0; dimension is 27115 but corresponding boolean dimension is 29582

I believe that this comes from the inconsistency in the codimensions between facets in the overall mesh, and what would be edges in the subdomain. I can simply define the boundaries geometrically, although this would only be a temporary solution that would need to be fixed later on. The easiest approach may be to define separate meshes, seeing as I’ll need to interpolate regardless. Although, if there is an easy solution to this then I’d rather keep the single mesh.

Any thoughts would be appreciated!

Update: I realize that the two boundary arrays are of different length. Is there any straightforward approach to extracting the relevant information (the boundary data from the larger array with only considering one subdomain)? Thank you.

MWE:

from fenics import *
from mshr import *
import numpy as np
import meshio
from msh2xdmf import import_mesh
#from fenicstools import interpolate_nonmatching_mesh

#Create mesh, domains and boundaries
mesh, bound, subd, association_table= import_mesh(
prefix=‘Air_Domain12’,
dim=2,
subdomains=True,
)

File(“Boundaries.xml”) << bound

fix_subdomains = MeshFunction(“size_t”, mesh, mesh.topology().dim())
fix_subdomains.set_all(0)
fix_subdomains.array()[subd.array()==11] = 1

submesh = SubMesh(mesh, fix_subdomains, 1)

boundaries = MeshFunction(“size_t”, mesh, “Boundaries.xml”)
submesh_boundaries = MeshFunction(‘size_t’, submesh, mesh.topology().dim()-1)
submesh_boundaries.set_all(0)
submesh_boundaries.array()[boundaries.array()==6] = 1

2

It’s far from ideal, but I went ahead and coded up a rough fix:

boundaries = MeshFunction(“size_t”, mesh, “Boundaries.xml”)
sub_boundaries = MeshFunction(“size_t”, submesh, mesh.topology().dim()-1,0)

tol = mesh.hmax()/5
a = 10 # set to your max subdomain index number
for i in range(a):
counter = -1
for facet in SubsetIterator(boundaries,i):
counter = counter + 1
if counter > -1:
x = facet.midpoint().x()
y = facet.midpoint().y()
for sub_facet in SubsetIterator(sub_boundaries,0):
x_sub = sub_facet.midpoint().x()
y_sub = sub_facet.midpoint().y()
if abs(x - x_sub) <= tol * abs(x) and abs(y - y_sub) <= tol * abs(y):
sub_boundaries.array()[sub_facet.index()] = boundaries.array()[facet.index()]

This process sends boundary data from the FacetFunction of the overall mesh to the FacetFunction of a specific subdomain (will need to be slightly modified for more than two subdomains). It’s computationally expensive, but not entirely impractical. If anyone has a better fix, please share!

3

We can directly get the boundary by using these commands in gmsh:

#This is for domian
Physical Surface("1")={9};
#This is for boundary
Physical Line("1") ={1};

and in FEniCS we can directly use for example:

 facets = MeshFunction("size_t", mesh, "filename_facet_region.xml")
 BC_Bot = DirichletBC(V.sub(1), 0, facets,2)