periodic boundary condition for meshes that are not built-in

Question

I'm trying to impose a periodic boundary condition on a cubic domain with a cylinder cut out of it. However, running the code below results in a solution which is not periodic. Interestingly, if I change the mesh to one without a cylinder, namely mesh = BoxMesh(0,0,0,1,1,1,5,5,5), then the solution is periodic.

Here is my code:

from dolfin import *
# Specify periodic boundary conditions
class PeriodicBoundary(SubDomain):

def inside(self, x, on_boundary):
    return bool(near(x[1],0) and on_boundary)

def map(self, x, y):
        y[0] = x[0]
        y[1] = x[1] - 1.
        y[2] = x[2]

# Define 2D geometry
domain = Box(0, 0, 0, 1., 1., 1.)-Cylinder(Point(0.5, 0, 0.5),Point(0.5, 1.0, 0.5), .25, 20)
mesh = Mesh(domain, 40)

#mesh = BoxMesh(0,0,0,1,1,1,5,5,5)

# define function spaces
V = FunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary())
R = FunctionSpace(mesh, "R", 0, constrained_domain=PeriodicBoundary())
W = V * R

f = Expression('sin(x[1])')

# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)

a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = f*v*dx

# Solve variational problem
w = Function(W)
solve(a == L, w)
(u, c) = w.split()

# Plot function
plot(u)
interactive()

Interestingly, even if I use the mesh

domain = Box(0, 0, 0, 1., 1., 1.)
mesh = Mesh(domain, 60)

the solution isn't periodic, so it appears that the issue is related to using a mesh that is not one of the built-in ones, but I can't figure out what it is.

Answer 1

For periodic meshes, the points must match up exactly.
When you use the CGAL mesh generator, this does not happen, whereas with the
built-in BoxMesh, the points are on a regular grid, so they do.

You need to build your mesh with a generator which can constrain the points on the boundaries.

Comments

Thanks a lot!