Parallel run doesn't work with restrictions

We've put together a little saddle-point problem for Poisson. Everything seems to work beautifully in serial, but when I run mpirun -n 2 python poissonSimple.py I get the following output

Process 0: Number of global vertices: 25
Process 0: Number of global cells: 32
Process 0: Solving linear variational problem.
Process 1: Solving linear variational problem.
Traceback (most recent call last):
  File "poissonSimple.py", line 36, in <module>
    solve(a==L, uSol, bcs)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 268, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 297, in _solve_varproblem
    solver.solve()
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
***
*** https://answers.launchpad.net/dolfin
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to create mesh entity.
*** Reason: Mesh entity index 170476920 out of range [0, 16] for entity of dimension 2.
*** Where: This error was encountered inside MeshEntity.cpp.
*** Process: 1
*** -------------------------------------------------------------------------

I am able to run some of the sample codes in parallel. The main difference between our code and some of the samples is that we are using the newer restriction="facet" functionality. Our Poisson solver is below. Thanks for any help.

# homogeneous poisson with -Del u = 1 on unit square
from dolfin import *

# Create mesh and define function space
degree = 1
delta_p = 2
N = 4
mesh = UnitSquareMesh(N,N)

u0 = Constant(0.0)

# define spaces
U = FunctionSpace(mesh, "CG", degree)    # trial functions
V = FunctionSpace(mesh, "DG", degree+delta_p) # discontinuous test functions
F = FunctionSpace(mesh, "RT", degree, restriction="facet") # fluxes
E = MixedFunctionSpace([U,V,F])
(u,e,f) = TrialFunctions(E)
(du,v,df) = TestFunctions(E)
n = FacetNormal(mesh)

def ip(e,v):
 return inner(grad(e),grad(v)) + inner(e,v)

def b(u,v,f):
 return inner(grad(u),grad(v))*dx - inner(dot(f('+'),n('+')),jump(v)) * dS - inner(dot(f,n),v)*ds

a = b(u,v,f) + b(du,e,df) + ip(e,v)*dx

L = inner(Constant(1.0),v)*dx

bcs = []
bcs.append(DirichletBC(E.sub(0), u0, DomainBoundary())) # boundary conditions on u

uSol = Function(E)
solve(a==L, uSol, bcs)
(u,e,f) = uSol.split()

plot(u, title="u")
plot(e, title="e")
plot(mesh, title="mesh")
interactive()

Comments (8)