Periodic Boundary conditions in Poisson 2D

1

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, [])
2

First of all, note that we now encourage all new users to select DOLFINx instead of DOLFIN, ref: The new DOLFINx solver is now recommended over DOLFIN
DOLFINx has documentation at: https://docs.fenicsproject.org/
Demos: Demos — DOLFINx 0.4.1 documentation
Tutorial: The FEniCSx tutorial — FEniCSx tutorial
Extension for periodic conditons (dolfinx_mpc): dolfinx_mpc/demo_periodic_geometrical.py at master · jorgensd/dolfinx_mpc · GitHub

As for your problem above, running your code throws the following error:

Traceback (most recent call last):
  File "script.py", line 5, in <module>
    mesh = BoxMesh(Point(0.0, 0.0), Point(a, a), a, a)
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfin.cpp.generation.BoxMesh(arg0: dolfin::Point, arg1: dolfin::Point, arg2: int, arg3: int, arg4: int)
    2. dolfin.cpp.generation.BoxMesh(comm: MPICommWrapper, p0: dolfin::Point, p1: dolfin::Point, nx: int, ny: int, nz: int)

Invoked with: <dolfin.cpp.geometry.Point object at 0x7fdff4934458>, <dolfin.cpp.geometry.Point object at 0x7fdff492cbc8>, 10, 10

which makes sense as a box cannot be two dimensional.
I’ve made various changes to your code to make it compile:

from IPython import embed
from dolfin import *

P = 10
# Create mesh and define function space
mesh = RectangleMesh(Point(0.0, 0.0), Point(P, P), 50, 50)

# 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 near(x[0], 0) and on_boundary

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


class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        eps = 1e-5
        return on_boundary and (x[0] > eps and x[0] < P-eps) and (near(x[1], 0) or near(x[1], P))


# 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


bc = DirichletBC(V, u_D, DirichletBoundary())

# Compute solution
u = Function(V)
solve(a == L, u, [bc])
File("u.pvd") << u
1 Like