1D periodic boundary condition not satisfied

variational formulation
1

I am new to FEniCS and I am trying to solve a type of inviscid Burgers equation with periodic boundary conditions on the unit interval. The code I have written does not throw any errors, but the resulting functions that I plotted did not have periodic boundary conditions.

I have supplied my code below. I have installed FEniCS 2018.1.0 from Anaconda.
I have added a print statement to the map function and noticed that it doesn’t execute.
Am I doing something wrong?

from fenics import *
import matplotlib.pyplot as plt

t0 = 0
num_elements = 300
p_lbd = 2
dt = 0.001

# define periodic boundary
class PeriodicBoundary(SubDomain):

    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS or x[0] > 1 - DOLFIN_EPS and on_boundary)

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


mesh = UnitIntervalMesh(num_elements)
V = FunctionSpace(mesh, 'P', 1, constrained_domain=PeriodicBoundary())

u = Function(V)
v = TestFunction(V)

u0 = Expression("p_lbd + sin(DOLFIN_PI*(2*x[0]))", p_lbd=p_lbd, element=V.ufl_element())

u0 = project(u0, V)

# implicit euler formulation
F = u*v*dx + dt*inner(grad(u),grad(v))*u*dx + dt*0.01*inner(grad(u), grad(v))*dx - u0*v*dx

T = 0.1
t = dt

# steps
while t <= T:
    solve(F == 0, u)
    t += dt
    u0.assign(u)

plot(u)
plt.show()
2

In your inside function you should specify the “target” domain only.
It works for me with

return bool(- DOLFIN_EPS < x[0] < DOLFIN_EPS and on_boundary)
3

Thanks for the hint! I changed that line but I get really strange looking plots now, for example for the initial function I get this

periodicu0

However, evaluating the function at given points seems to give a correct result, so I guess this solved the problem.

4

Hi, thanks for this remark, which saves me a lot of time debugging. Do you know why this happens? Using compute_vertex_values gives the desired results but using .vector.().get_local() gives the strange looking plot.

5

The dofs in the function space is not ordered in the same was as the points in the geometry. To Get the coordinates of the dofs in the same order as u.vector().get_local() returns do x=V.tabulate_dof_coordinates()