Do the periodic boundary condition solvable? (pure periodic boundary condition)

13 months ago by
Hi, all.

My new question comes now.

There is a demo here for a 2D square where the left and right boundaries are periodic and the top and bottom boundaries are Dirichlet boundary conditions. I want to solve a pure periodic problem which the left and right boundaries are periodic and the top and bottom boundaries are periodic boundaries as well. So, I change the demo, as follows.

from dolfin import *

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

# Left and Bottom boundaries are "target domain" G
def inside(self, x, on_boundary):
return bool((near(x[0],0) or near(x[1],0)) and on_boundary)

def map(self, x, y):
# Map right boundary (H) to left boundary (G)
if near(x[0],0):
y[0] = x[0] - 1.0
y[1] = x[1]
# Map Top boundary to Bottom boundary
elif near(x[1],0):
y[0] = x[0]
y[1] = x[1] - 1.0

# Create mesh and finite element
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary())
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

f = Expression(["x[0]*sin(5.0*pi*x[1]) + 1.0*exp(-(pow((x[0]-0.5),2) + pow((x[1]-0.5),2))/0.02)"])
a = inner(grad(u), grad(v))*dx
L = inner(f,v)*dx

# Compute solution
u = Function(V)
solve(a == L, u)

# Plot solution
plot(u, interactive=True)

*** Error: Unable to periodic boundary mapping.
*** Reason: Need to set coordinate 0 in
*** Where: This error was encountered inside PeriodicBoundaryComputation.cpp.
*** Process: 0

According to the ERROR, it seems the definition of periodic subdomain is not correct. Could you tell me more about the logic of class PeriodicBoundary(SubDomain) ? and the problem left and right are periodic boundaries, top and bottom are periodic boundaries without Dirichlet boundary condition are solvable?
One more question. Does FEniCS handle the matrice of periodic problems as does in interiors?  

Thank !

Best, Hamilton
Community: FEniCS Project

1 Answer

13 months ago by
You want to map right boundary to left. So you must check


and similarly

well, it works.
written 13 months ago by Hamilton  
Could you give more explanations of the principle of periodic boundary conditions. I am confused the definition about it.
For example:

1, If the boundary condition of my problem is pure periodic problem, without other conditions, is it solvable? If so, we need to add pure Neumann boundary condition?

2, class PeriodicBoundary(SubDomain): give two def. The first one is bool type. In my sense, if the left one is a target, the right one supposed to be: x_{right} = x_{left} + 1. However, the demo shows that :

# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(near(x[0],0) and on_boundary)

# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y[0] = x[0] - 1.0
y[1] = x[1]

If i change y[0] = x[0]+1.0; it works, however the solution is a bit differnt from the demo result.

What is wrong?
written 12 months ago by Hamilton  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »