Strange behaviour with write_checkpoint / read_checkpoint

I/O
7

I cannot reproduce your issue. I separated the mesh generation to a separate file, calling the mesh i load in mesh_mesh.xdmf.

import dolfin as dol
import numpy as np
import meshio

class Delta_Function(dol.UserExpression):
    '''Class for mollified Dirac delta function'''
    def __init__(self, x_0, eps):
        '''parameters
           x_0 : 1D array like. Center
           eps : Float. Standard deviation'''
        self.x_0 = x_0
        self.eps2 = eps**2
        self.dim = len(x_0)
        self.norm = (np.sqrt(2. *np.pi * self.eps2))**self.dim
        dol.UserExpression.__init__(self)
    def eval(self, values, x):
        r2 = np.sum([(x[i]-self.x_0[i])**2 for i in range(self.dim)])
        values[0] = np.exp(-0.5 * r2 /self.eps2)/self.norm
    def value_shape(self):
        return ()


def load_XDMF_files(filename):
    '''load XDMF files'''
    mesh = dol.Mesh()
    with dol.XDMFFile(filename) as infile:
        infile.read(mesh)
    #create subdomain cell function
    mvc = dol.MeshValueCollection('size_t', mesh, 3)
    with dol.XDMFFile(''.join(filename.split(sep ='_')[:-1]) + '_subdomains.xdmf') as infile:
        infile.read(mvc, 'name_to_read')
    subdomain_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    #create boundary subdomain facet function
    mvc = dol.MeshValueCollection('size_t', mesh, 2)
    with dol.XDMFFile(''.join(filename.split(sep ='_')[:-1]) + '_boundary_subdomains.xdmf') as infile:
        infile.read(mvc, 'name_to_read')
    boundary_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    del(mvc)
    return mesh, subdomain_marker, boundary_marker


mesh, subdomain_marker, boundary_marker = load_XDMF_files('mesh_mesh.xdmf')   #define measures
dx = dol.Measure('dx', domain = mesh, subdomain_data = subdomain_marker)
ds = dol.Measure('ds', domain = mesh, subdomain_data = boundary_marker)

element = dol.FiniteElement('P', mesh.ufl_cell(), degree = 1)
W = dol.FunctionSpace(mesh, element)
u = dol.TrialFunction(W)
v = dol.TestFunction(W)
a = dol.inner(dol.grad(u), dol.grad(v)) * dx
L = 4. * np.pi * Delta_Function([0.0, 0.0, 0.52], 0.05) * dol.inner(dol.grad(v), dol.Constant((0,0,1))) * dx
bc = [dol.DirichletBC(W, dol.Constant(0.), boundary_marker, 1),
          dol.DirichletBC(W, dol.Constant(0.), boundary_marker, 2)]

u = dol.Function(W)
problem = dol.LinearVariationalProblem(a, L, u, bc)
solver = dol.LinearVariationalSolver(problem)
solver.solve()

with dol.XDMFFile('potential.xdmf') as outfile:
            outfile.write_checkpoint(u, 'phi')
v = dol.Function(W)
with dol.XDMFFile('potential.xdmf') as infile:
        infile.read_checkpoint(v, 'phi')

index = np.argwhere(u.compute_vertex_values() != v.compute_vertex_values())
print(index)
print(u.compute_vertex_values()[index])
print(v.compute_vertex_values()[index])
print(v.compute_vertex_values().shape)
print(v.compute_vertex_values().shape)

gives me

[]
[]
[]
(50129,)
(50129,)

So are you sure you are running the latest version of dolfin?