Auto-adaptivity fails in 3D after first mesh refinement

[Marie Rognes]

from dolfin import *

# Load mesh and define function space
ou = 1
mesh = Box(0,-6,-6,48,6,6,8,2,2)
mesh.order()
mesh.init(2)
V = VectorFunctionSpace(mesh, "Lagrange", ou)

# Sub domain for clamp at left end
def left(x, on_boundary):
 tol = 1.0e-14
 return abs(x[0]) < tol and on_boundary

# Sub domain for shear at right end
class RightBoundary(SubDomain):
 def inside(self, x, on_boundary):
  tol = 1.0e-14
  return abs(x[0]-48) < tol and on_boundary

# Mark boundary parts
domain_parts = MeshFunction("uint",mesh,mesh.topology().dim())
boundary_parts = MeshFunction("uint",mesh,mesh.topology().dim()-1)
dOmega_t = RightBoundary()
dOmega_t.mark(boundary_parts,0)

# Define variational problem
v = TestFunction(V)
u = TrialFunction(V)
f = Constant((0.0, 0.0, 0.0))
tg = Expression(('0.0','-40.0/288.0*1.0/12*(x[1]+6.0)*(x[1]-6.0)','0.0'))
E = 30000.0
nu = 0.25
rho = 1

mu = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1+nu)*(1.0 - 2*nu))

def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d)

a = inner(grad(v), sigma(u))*dx
b = inner(rho*v,u)*dx
L = inner(v,f)*dx+inner(v,tg)*ds(0)
L.exterior_facet_domains = boundary_parts

# Define goal
u = Function(V)
g = u[1]*ds(0)
g.exterior_facet_domains = boundary_parts

u.name = "Displacements"
# Set up boundary condition at left end
c = Constant((0.0, 0.0, 0.0))
bcl = DirichletBC(V, c, left)

# Set up PDE and solve
eps = 1.e-2
problem = VariationalProblem(a, L, bcl)
problem.solve(u,eps,g)