Dirichlet BCs on an Annulus

mesh
1

I was trying to modify the Poisson demo example for an annular domain in 2D, but I’m struggling with getting it to apply Dirichlet BCs properly. When I run the following code:

import numpy as np
from fenics import *
from mshr import *

r0 = 0.5
r1 = 1.0

mesh = generate_mesh(Circle(Point(0,0),r1)-Circle(Point(0,0),r0), 10)

def inner(x, on_boundary):
    return near(x[0]**2 + x[1]**2, r0**2) and on_boundary
def outer(x, on_boundary):
    return near(x[0]**2 + x[1]**2, r1**2) and on_boundary

u_D = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=1)
V = FunctionSpace(mesh, "CG", 1)

bc = DirichletBC(V, u_D, inner)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

u = Function(V)
solve(a == L, u, bc)

But when I run this, I get the warning:

*** Warning: Found no facets matching domain for boundary condition.

and a zero solution. It seems that I’m not properly tagging the part of the boundary, but I’m not sure what I"ve done wrong.

2

There are two reasons why the boundary facets are not meeting the criteria set by inner:

  1. The default behavior of DirichletBC is to check the midpoints of facets as well as the vertices. If the vertices are on the circular boundary, the midpoint necessarily won’t be. You can turn this off by passing the optional keyword argument check_midpoint=False to the DirichletBC constructor. However, this still won’t work in general, because of the following.

  2. The Circle CSG primitive defined via mshr is actually a many-sided polygon. There is not necessarily one mesh facet per side of this polygon, so, even if you only check mesh vertices, there may still be some mesh vertices on the interiors of the sides making up the “circle”, which will not meet the radius criterion.

The simple way around these issues is to just check radius with inequalities, using a smaller radius on the outside and a larger radius on the inside. If you want to make sure all of the mesh vertices really do lie on the circular boundaries during refinement, you need to use the boundary snapping functionality. This is demonstrated in the following undocumented demo

https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/undocumented/smoothing/demo_smoothing.py

However, this functionality is actually broken in the current stable release. I fixed it with pull request 521 in the development version, but you can also use a quick work-around in the Python interface with 2019.1:

def py_snap_boundary(mesh, sub_domain):
    boundary = BoundaryMesh(mesh,"exterior")
    dim = mesh.geometry().dim()
    x = boundary.coordinates()
    for i in range(0,boundary.num_vertices()):
        sub_domain.snap(x[i,:])
    ALE.move(mesh,boundary)
1 Like
3

I am also working on this problem, and I ran into the same warning regarding no facets matching domain for boundary condition. Could you please share the complete workaround?

Thank You

4

Further to @kamensky’s comment regarding the linear discretisation of the circle geometry, there’s a third argument to near().

def inner(x, on_boundary):
    return near(x[0]**2 + x[1]**2, r0**2, TOL) and on_boundary

where TOL is some acceptable tolerance which will find the facets.

1 Like