Unable to get extract the solution at the surface for fine meshes

Errors
1

Hello,

I’m running fenics 2019.1.0 with singularity 3.6.0. I’ve successfully run the physics simulation and am just trying to get the data on the boundary. The script I have to do this works fine with the coarser meshes, but when I attempt it with the finer mesh, I get the following error:

Process 14: *** Warning: Mesh is empty, unable to create entities of dimension 0.
Process 14: *** Warning: Mesh is empty, unable to create entities of dimension 0.
Process 14: Number mesh entities for distributed mesh (for specified vertex ids).
Traceback (most recent call last):
  File "extract_surfaces.py", line 38, in <module>
    Vb = VectorFunctionSpace(surface,'P',1)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/functionspace.py", line 222, in VectorF
unctionSpace
    return FunctionSpace(mesh, element, constrained_domain=constrained_domain)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/functionspace.py", line 31, in __init__
    self._init_from_ufl(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/functionspace.py", line 50, in _init_from_ufl
    dolfin_dofmap = cpp.fem.DofMap(ufc_dofmap, mesh)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to number mesh entities.
*** Reason:  Global vertex indices exist at input. Cannot be renumbered.
*** Where:   This error was encountered inside MeshPartitioning.cpp.
*** Process: 14
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

This is the script that I have:

# Environment setup
import os

# FEniCS imports
from fenics import *
from dolfin import *
import numpy

# Define communicator and rank for this process, in case we are running with MPI
mpi_comm = MPI.comm_world
mpi_rank = mpi_comm.Get_rank()

PROGRESS = 16
TRACE = 13
DBG = 10

set_log_level(PROGRESS)

mesh = Mesh(mpi_comm)
with XDMFFile(mpi_comm, "mesh.xdmf") as infile:
    infile.read(mesh)
# Read boundaries
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
surface = BoundaryMesh(mesh,"exterior")


# Try reading checkpoint
V = VectorFunctionSpace(mesh,'P',1)
Vb = VectorFunctionSpace(surface,'P',1)

u = Function(V)
ub = Function(Vb)

xdmf_vol = XDMFFile(mpi_comm, "displacement.xdmf")
xdmf_sur = XDMFFile(mpi_comm, "displacement_surf.xdmf")

for i in range(0,5):
    xdmf_vol.read_checkpoint(u, "u", i)
    ub = interpolate(u,Vb)
    ub.rename("ub", "displacement")
    #xdmf_sur.write(ub,i)
    if i < 1:
        xdmf_sur.write_checkpoint(ub,"ub",i,XDMFFile.Encoding.HDF5,False)
    else:
        xdmf_sur.write_checkpoint(ub,"ub",i,XDMFFile.Encoding.HDF5,True)

Seeing as the script works for the coarser meshes, I’m not sure what to do. Would there be any workarounds for this? Thanks!