Cahn-Hillard Equation Periodic Boundary Condition In 3D?

I want to implement periodic boundary condition in Z planes and different constant fluxes on other boundaries . Can you take a Look at my Code and see what changes i need to do, in order to obtain the results.

python
import random
from dolfin import *


        

# Sub domain for Periodic boundary condition

class PeriodicBoundary(SubDomain):

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

    def map(self, x, y):
        y[0] = x[0]
        y[1] = x[1]
        y[2] = x[2] - 100.0
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 100.0)

class Inside(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Outside(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 100.0)        

left = Left()
inside = Inside()
right = Right()
outside = Outside()
pbc = PeriodicBoundary()


boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)

boundaries.set_all(0)
left.mark(boundaries, 1)
inside.mark(boundaries,5)
right.mark(boundaries, 3)
outside.mark(boundaries, 6)
pbc.mark(boundaries,2)
pbc.mark(boundaries,4)
# Class representing the intial conditions

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world))
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 0.50 + 0.002*(0.5 - random.random())
        values[1] = 0.0
    def value_shape(self):
        return (2,)
# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)
mesh = BoxMesh(Point(0.0,0.0,0.0),Point(100.0,100.0,100.0), 10,10,10)
P1 = FiniteElement("Lagrange",mesh.ufl_cell(), 1, constrained_domain=pbc)
ME = FunctionSpace(mesh, P1*P1,)        

# Model parameters
lmbda  = 5    # surface parameter
dt     = 0.1  # time step
theta  = 0.5  # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson


# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True


# Create mesh and define function space


# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)

# Compute the chemical potential df/dc
c = variable(c)
f    = 50*(c-0.3)**2*(0.7-c)**2
dfdc = diff(f, c)
gL=Constant(0.02)
gR=Constant(-0.03)
gI=Constant(-0.05)
gO=Constant(0.01)
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh,)
# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(5*grad(mu_mid), grad(q))*dx             +gL*q*ds(1)+gR*q*ds(3)+gI*q*ds(5)+gO*q*ds(6)
L1 = mu*v*dx - dfdc*v*dx - dot(lmbda*grad(c), grad(v))*dx
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
file = File("sol6/prpf3.pvd", "compressed")

# Step in time
t = 0.0
T = 100*dt
while (t < T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    file << (u.split()[0], t)