FEniCS: how to specify boundary conditions on a circle inside 2D mesh

Question

I would like to numerically find a mutual capacitance of two stripes of metal on the opposites sides of a cylinder. The problem is obviously a 2D Laplace equation. I would like to find the potential outside the cylinder as well. Therefore I have something like this:

mesh = UnitCircleMesh(50)
V = FunctionSpace(mesh, 'Lagrange', 1)

u_L = Constant(-1)
def left_boundary(x, on_boundary):
    r = math.sqrt(x[0] * x[0] + x[1] * x[1])
    return near(r, 0.7) and x[0] < 0 and between(x[1], (-0.5, 0.5))

u_R = Constant(1)
def right_boundary(x, on_boundary):
    r = math.sqrt(x[0] * x[0] + x[1] * x[1])
    return near(r, 0.7) and x[0] > 0 and between(x[1], (-0.5, 0.5))

leftPlate = DirichletBC(V, u_L, left_boundary)
rightPlate = DirichletBC(V, u_R, right_boundary)
bcs = [leftPlate, rightPlate]

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx

u = Function(V)
solve(a == Constant(0) * v * dx, u, bcs)

When I run it, I receive *** Warning: Found no facets matching domain for boundary condition.. And the found solution is 0. What is wrong?

Answer 1

As the error message suggest the DirichletBC does not hit any mesh entities with corresponding dofs. You need to examine your subdomains. One way to debug these are to apply the DirichletBC to a Function and plot the result:

def left_boundary(x, on_boundary):
    r = math.sqrt(x[0] * x[0] + x[1] * x[1])
    return r < 0.5

leftPlate = DirichletBC(V, u_L, left_boundary)
u = Function(V)
u.vector()[:] = 0
leftPlate.apply(u.vector())
plot(u, interactive=True)

This will plot a nice hat. Now you can start tweak your left_boundary function until you reach the result you want.

Answer 2

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')
```