Hello, I am starting with FEniCS and I’m finding very hard to learn this as intead of documentation which exlains each function etc I can only find demos.
I want to learn how to implement Periodic Boundary Conditions on a 2D system where 2 faces (left and right) have PBC and the other two walls have Dirichlet BC.
All the online examples are too old, the only working example found is in 3D and have no idea why it does not work when I “translate” it to 2D.
Anybody has a working demo? The code I tried to modify to make it into a 2D system looks like the following
from fenics
import *
a = 10
## Create mesh and define function space
mesh = BoxMesh(Point(0.0, 0.0), Point(a, a), a, a)
## u_D(x,y) = 1 + x^2 + 2y^2
## Expression object takes mathematical expression in the C++ syntax
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the two slave edges
return bool ((near(x[0], 0) or near(x[1], 0)) and
(not ((near(x[0], a)) or
(near(x[0], a) and near(x[1], a))
)) and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
#### define mapping for a single point in the box, such that 3 mappings are required
if near(x[0], a) and near(x[1], a) and near(x[2],a):
y[0] = x[0] - a
y[1] = x[1] - a
##### define mapping for edges in the box, such that mapping in 2 Cartesian coordinates are required
if near(x[0], a):
y[0] = x[0] - a
y[1] = x[1]
elif near(x[1], a):
y[0] = x[0]
y[1] = x[1] - a
elif near(x[0], a) and near(x[1], a):
y[0] = x[0] - a
y[1] = x[1] - a
#### right maps to left: left/right is defined as the x-direction
elif near(x[0], a):
y[0] = x[0] - a
y[1] = x[1]
### bottom maps to top: is defined as the y-direction
elif near(x[1], a):
y[0] = x[0]
y[1] = x[1] - a
## Define boundary condition
pbc = PeriodicBoundary()
V = FunctionSpace(mesh, 'Lagrange', 1, constrained_domain=pbc)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
f = Constant(-6.0)
L = f*v*dx
## Compute solution
u = Function (V)
solve(a == L, u, [])