7
$\begingroup$

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?

$\endgroup$
6
  • $\begingroup$ Are you trying to hit boundary facets or facets on radius 0.7 (in unit circle)? $\endgroup$ – Jan Blechta May 24 '13 at 21:35
  • $\begingroup$ Are you trying to hit exterior or interior (on radius 0.7) boundary? If the latter is your intent consider wheter there are facets with their vertices and midpoint on radius cca 0.7 forming connected interior boundary. If yes, you probably need to increase tolerance by near(r, 0.7, tol) as target facets would be rough approximation to circle. $\endgroup$ – Jan Blechta May 24 '13 at 21:43
  • $\begingroup$ Also fit your mesh = UnitCircleMesh(n) where $n = \frac{2}{0.7} m$ both $n$, $m$ being natural. $\endgroup$ – Jan Blechta May 24 '13 at 21:54
  • $\begingroup$ I want to calculate the potential both in the exterior and the interior of the cylinder. $\endgroup$ – facetus May 24 '13 at 22:38
  • $\begingroup$ Then increase tolerance of near function very much and tweak n so that mesh hits better radius 0.7 as suggested above. $\endgroup$ – Jan Blechta May 24 '13 at 22:44
5
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Thanks! I'll try this. What is a subdomain (excuse me for a silly question)? $\endgroup$ – facetus May 24 '13 at 22:41
  • 1
    $\begingroup$ @noxmetus: SubDomain is DOLFIN class for definition of some subdomain of your spatial domain, see fenicsproject.org/documentation/dolfin/1.2.0/python/…. Function left_boundary(x, on_boundary) is convenience way of defining subdomains without actually subclassing SubDomain. Notice that left_boundary has same interface as SubDomain.inside(). $\endgroup$ – Jan Blechta May 25 '13 at 0:05
  • 1
    $\begingroup$ Note that you can plot boundary conditions directly without creating a Function first: plot(leftPlate) $\endgroup$ – Anders Logg May 26 '13 at 21:55
1
$\begingroup$

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: simulation results

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')
```
$\endgroup$
3
  • $\begingroup$ Is the Fenics Q&A closed now? What about this? $\endgroup$ – nicoguaro Feb 3 at 23:26
  • $\begingroup$ Good to know about this! When old forum QA responses show up on a search engine, like this link, they link to a "new" QA forum and that link is broken. I will update the answer. $\endgroup$ – Ryan Feb 4 at 15:46
  • $\begingroup$ In that case I suggest that you post a question in the new Q&A mentioning this problem. $\endgroup$ – nicoguaro Feb 4 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.