How to solve Poisson with periodic boundary conditions

Errors
1

Hi guys, I am a new Fenics user. I am trying to solve a Poisson equation with PBC. However, I got the same value of my solution. May I ask what kind of mistakes I made?

a = 1e-7
L = 1e-7
num_points = 30
x = np.linspace(-L/2, L/2, num_points, endpoint=False)
y = np.linspace(0, L, num_points, endpoint=False)
z = np.linspace(0, L, num_points, endpoint=False)
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')


# Gaussian width
sigma = L / (1*num_points)  # Adjust as needed

rho_array = np.zeros((num_points, num_points, num_points))

# Coordinates for the 10 point charges on the x=0 plane
charge_positions = np.random.rand(10, 3) * 1e-7

charge_positions[:,0] = 0

#charge_positions = np.array([[0, 0, 0]])

# For each charge, add a Gaussian distribution centered at its position
for charge_pos in charge_positions:
    x_pos, y_pos, z_pos = charge_pos
    r2 = (X-x_pos)**2 + (Y - y_pos)**2 + (Z - z_pos)**2
    delta_approx = (1 / (np.pi**(3/2) * sigma**3)) * np.exp(-r2 / sigma**2)
    rho_array += delta_approx
    
rho_array = np.roll(rho_array, shift=-num_points//2, axis=0)
rho_array = rho_array * 2.26e-10

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) or near(x[2], 0)) and 
            (not ((near(x[0], a) and near(x[2], a)) or 
                  (near(x[0], a) and near(x[1], a)) or
                  (near(x[1], a) and near(x[2], 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
            y[2] = x[2] - a
        ##### define mapping for edges in the box, such that mapping in 2 Cartesian coordinates are required
        if near(x[0], a) and near(x[2], a):
            y[0] = x[0] - a
            y[1] = x[1] 
            y[2] = x[2] - a      
        elif near(x[1], a) and near(x[2], a):
            y[0] = x[0] 
            y[1] = x[1] - a
            y[2] = x[2] - a
        elif near(x[0], a) and near(x[1], a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2]         
        #### 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]
            y[2] = x[2]
        ### back maps to front: front/back is defined as the y-direction    
        elif near(x[1], a):
            y[0] = x[0]
            y[1] = x[1] - a
            y[2] = x[2] 
        #### top maps to bottom: top/bottom is defined as the z-direction        
        elif near(x[2], a):
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - a

mesh = BoxMesh(Point(0, 0, 0), Point(L, L, L), num_points, num_points, num_points)
PBC = PeriodicBoundary()
V = FunctionSpace(mesh, 'Lagrange', 1, constrained_domain=PBC)
phi = TrialFunction(V)
v = TestFunction(V)

rho = Function(V)
rho.vector().set_local(rho_array.flatten())
rho.vector().apply('insert')
#rho = Constant(1)

a = dot(grad(phi), grad(v))*dx
l = rho*v*dx

phi_solution = Function(V)
solve(a == l, phi_solution)

phi_values = phi_solution.vector().get_local()