How can I implement a combination of Dirichlet und Perodic Boundary Conditions while using Line Tags created with GMSH?

Hello Community,
I am still relative new to FEniCS and this is my first forum post.

I was already able to create meshes with GMSH and import them into FEniCS using the msh to xdmf conversion method. Furthermore I was able to use the boundarys specified in GMSH to define DirechletBC in FEniCS. But I can’t figure out how to apply periodic BCs with the Line Tags created with GMSH.

My simple geo file:

// Mesh resulution
lc = 1e-1;

// Corner points of a rectangle
Point(1) = {0, 0, 0, lc};
Point(2) = {1, 0, 0, lc};
Point(3) = {1, 1, 0, lc};
Point(4) = {0, 1, 0, lc};

// Lines in between the corner points
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};

// Curve Loop and Surface of the rectangle
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};

// Line Tags for BCs. Numeric names for FEniCS
Physical Curve(1111) = {4};
Physical Curve(4444) = {1};
Physical Curve(5555) = {3};
Physical Curve(9999) = {2};

Physical Surface(555) = {1};

Mesh conversion script:

import meshio
import numpy as np

# --------------------------------------------------
# ----------MESH CONVERSION----------
# Directory and name of the mesh which will be transformed
Meshdirectory = 'OrgMeshes/'
MeshFile = Meshdirectory + 'rect_periodic.msh'


def create_mesh(mesh, cell_type, prune_z=False):
    # Read all relevant cells and get their physical data
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    # Read all points of the mesh
    outpoints = mesh.points

    if prune_z:
        # In this case the z axis is deleted from the mesh (gmsh and salome always create 3d meshes)
        outpoints = np.delete(arr=mesh.points, obj=2, axis=1)

    # The relevant points and cells are combined to a meshio mesh and returned
    return meshio.Mesh(points=outpoints, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]})


# Read in the mesh file
msh = meshio.read(MeshFile)

# Extract the mesh and the boundaries from the msh file
mesh_2d = create_mesh(msh, "triangle", True)
boundaries = create_mesh(msh, "line", True)

# Write the extracted mesh and boundaries to xdmf files (the h5 files are created automatically)
meshio.write("mesh.xdmf", mesh_2d)
meshio.write("boundaries.xdmf", boundaries)
# --------------------------------------------------

Mesh loading script:

import os
import csv
import matplotlib.pyplot as plt
import numpy as np
from dolfin import cpp
from fenics import *

# --------------------------------------------------
# ----------MESH LOADING----------
mesh = Mesh()

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1)
with XDMFFile("boundaries.xdmf") as infile:
    infile.read(mvc, "name_to_read")
loadedBoundaryList = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=loadedBoundaryList, subdomain_id=12)

plot(mesh)
# --------------------------------------------------

# --------------------------------------------------
# ----------SIMULATION----------
D = 1                   # Diffusion coefficient
T = 1.0                 # Final time
num_steps = 100         # Number of time steps
num_steps_initial = 10  # Number of initial time steps
dt = T / num_steps      # Time step size
dt_initial = dt / 10    # Initial time step size
dir_name = 'SimTest'   # Folder for the results

V = FunctionSpace(mesh, 'P', 2)

INLET = DirichletBC(V, Constant(1), loadedBoundaryList, 1111)
OUTLET = DirichletBC(V, Constant(0), loadedBoundaryList, 9999)

bcs = [INLET, OUTLET]

# Define initial value
C_0 = interpolate(Constant(0), V)

# Define variational problem
C = TrialFunction(V)
v = TestFunction(V)

F_BE = (C - C_0) * v * dx + dt_initial * D * dot(grad(C), grad(v)) * dx
F_C = (C - C_0) * v * dx + 0.5 * dt * D * dot(grad(C), grad(v)) * dx + 0.5 * dt * D * dot(grad(C_0), grad(v)) * dx

F_BE_l, F_BE_r = lhs(F_BE), rhs(F_BE)
F_C_l, F_C_r = lhs(F_C), rhs(F_C)

How can i change the last script so that i not only apply the DirichletBCs for 1111 and 9999 but also a periodicBC between 4444 and 5555?
I hope I gave enough informations and my question is understandable otherwise please don’t hesitate to correct me or ask further questions.
Thank you all in advance,
Alex