I have also had this problem and spent a lot of time on various forums and I have finally come up with a good solution. This solution will allow you to specify a boundary condition on any whole subdomain within your geometry, and each subdomain can have its own separate boundary condition. I am posting this solution here as the old official Fenics QA forum has now closed, and this solution addresses several unanswered questions on those forums as well. The new forum is available here.
In my particular problem I have two circles inside of a rectangle and Lagrangian function space. In this example I set the left wall to be 1, the right wall to be 0, the left circle to be 0, and the right circle to be 3.
My entire code is included, but this is the important block:
marked_cells = fenics.SubsetIterator(CELL_MESHFUNCTION_RESULT, SUBDOMAIN_INDEX_NUMBER)
for cells in marked_cells:
for f in fenics.facets(cells):
FACET_MESHFUNCTION_RESULT[f] = CHOSEN_CONSTANT
# later in code
BC = fenics.DirichletBC(FUNCTION_SPACE,
fenics.Constant(0), FACET_MESHFUNCTION_RESULT, CHOSEN_CONSTANT)
The result yields the figure:

My understanding of what is going on is that the SubsetIterator function provides the subset of cells within the mesh which match a particular subdomain index number. Next, in a loop, we reassign any facet element which is inside any of these cells with a new index number which we can reference later. Then, we use this index to apply a boundary condition to each of these facets.
This method is more efficient than some others I've found since it does not loop over every facet in every cell in the mesh, only those in the desired subdomain. My intuition is that facets are numbered sequentially, and therefore using low index numbers to mark facets may be undesirable as there may be a random facet in the model which is unintentionally marked with a boundary condition. Using indices much larger than the number of facets may alleviate this issue, but this is still unknown. I was unable to find any mismatched facets in my testing.
Here is the full code (Windows 10, Ubuntu 20 through WSL2, Python 3.8, Fenics 2019.1.0)
# Ryan Budde CID 2021
# imports
import fenics as fn
import mshr
from math import sin, cos, pi
import matplotlib.pyplot as plt
import pprint as pprint
# constants
squareW = 3
squareL = 2
circRad = 0.25
# create background geometry
domain = mshr.Rectangle(fn.Point(0,0), fn.Point(squareW, squareL))
# create source and sink circles
circPos = mshr.Circle(fn.Point(1, 1), circRad)
circNeg = mshr.Circle(fn.Point(2, 1), circRad)
# assign the circles to the domain
domain.set_subdomain(1, circPos)
domain.set_subdomain(2, circNeg)
# generate the mesh
mesh = mshr.generate_mesh(domain, 28)
# define subdomain markers and facets
markers = fn.MeshFunction('size_t',mesh, mesh.topology().dim(), mesh.domains())
boundaries = fn.MeshFunction('size_t',mesh, mesh.topology().dim()-1, mesh.domains())
# define function space
V = fn.FunctionSpace(mesh, 'Lagrange', 1) # a first order lagrangian function
# establish BC as C++ command
leftBC = 'near(x[0], 0)' # left wall
rightBC = 'near(x[0], 3)' # right wall
# establish BCs for subdomains
marked_cells = fn.SubsetIterator(markers,1) #left circle
for cells in marked_cells:
for f in fn.facets(cells):
boundaries[f] = 2
marked_cells = fn.SubsetIterator(markers,2) # right circle
for cells in marked_cells:
for f in fn.facets(cells):
boundaries[f] = 3
# assign BC
bc1 = fn.DirichletBC(V, fn.Constant(2), leftBC)
bc2 = fn.DirichletBC(V, fn.Constant(0), rightBC)
bc3 = fn.DirichletBC(V, fn.Constant(0), boundaries, 2)
bc4 = fn.DirichletBC(V, fn.Constant(3), boundaries, 3)
bcs = [bc1, bc2, bc3, bc4]
# assign dx
dx = fn.Measure('dx', domain=mesh, subdomain_data=markers)
# change permittivity of materials
class Permittivity(fn.UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 0:
values[0] = 1 # vacuum
elif self.markers[cell.index] == 1:
values[0] = 100 # small increase
elif self.markers[cell.index] == 2:
values[0] = 1e5 # large increase
mat_perm = Permittivity(markers, degree=1)
# variational problem
u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = mat_perm * fn.dot(fn.grad(u), fn.grad(v)) * dx
L = 1*v*dx(2) + 1*v*dx(3)
# solve problem
u = fn.Function(V)
fn.solve(a == L, u, bcs)
# plot results
plt.subplot(1,2,1)
fn.plot(markers, title='Domains')
plt.subplot(1,2,2)
c = fn.plot(u, title='Fields')
plt.colorbar(c)
plt.show()
print('done')
```
near(r, 0.7, tol)
as target facets would be rough approximation to circle. $\endgroup$mesh = UnitCircleMesh(n)
where $n = \frac{2}{0.7} m$ both $n$, $m$ being natural. $\endgroup$near
function very much and tweakn
so that mesh hits better radius 0.7 as suggested above. $\endgroup$