Periodic Boundary conditions in Poisson 2D

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