Solver crashes after upgrade

Dear all,
I have found out recently, that my script for solving magnetostatic problems does not work any more in the current version (1.0) of dolfin on ubuntu (12.04 and 11.10). I am not sure, if it is problem of dolfin since I upgraded other packages as well. The script is as follows

from dolfin import *

set_log_level(DEBUG)
parameters['linear_algebra_backend'] = 'uBLAS'
mesh = Mesh("MAC.xml")
subdomains = MeshFunction("uint", mesh, "MAC_mat.xml")

P0 = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
P1 = FunctionSpace(mesh, "Lagrange", 1)
PN = MixedFunctionSpace((P0, P1))

(v0, q0) = TestFunctions(PN)
(u0, p0) = TrialFunctions(PN)

# Define functions
ZeroVector = Constant((0.0, 0.0, 0.0))
ZeroScalar = Constant(0.0)

# Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

# Boundary condition
bc1 = DirichletBC(PN.sub(0), ZeroVector, DirichletBoundary())
bc2 = DirichletBC(PN.sub(1), ZeroScalar, DirichletBoundary())
bcs = [bc1, bc2]

MAC_chamber_length = 270e-3

class CurrentExpression0(Expression):
    #def __init__(self, mesh):
    # self.mesh = mesh
    def eval_cell(self, values, x, ufc_cell):
        #cell = Cell(self.mesh, ufc_cell.index)
        index = subdomains[ufc_cell.index]
        x, y, z = x[0:3]
        if index == 1 or index == 2:
            SC = 15e-3*30e-3 #coil cross section
            AT = 300*3.5 #number of ampere-turns
            J = AT/SC
            r = sqrt(x**2 + y**2)
            values[0] = -y/r*J
            values[1] = x/r*J
            values[2] = 0.0

        elif index ==3:
            SC = 16e-3*18e-3 #coil cross section
            AT = 200*3.5 #number of ampere-turns
            J = AT/SC
            r = sqrt(x**2 + y**2)
            values[0] = -y/r*J
            values[1] = x/r*J
            values[2] = 0.0
        elif index ==4:
            x0, z0 = 0, 373.5e-3
            x, z = x-x0, z-z0
            SC = 10e-3*18e-3 #coil cross section
            AT = -166*5. #number of ampere-turns
            J = AT/SC*1
            r = sqrt(x**2 + z**2)
            values[0] = -z/r*J
            values[1] = 0.0
            values[2] = x/r*J
        else:
            values[0] = 0.0
            values[1] = 0.0
            values[2] = 0.0

    def value_shape(self):
        return (3,)

J = CurrentExpression0()

class MuExpression0(Expression):
    def eval_cell(self, value, x, ufc_cell):
        mu0 = 4*pi*1e-7
        index = subdomains[ufc_cell.index]
        if index == 5 or index == 6:
            value[0] = mu0*2333
            #permitivity of 1117 steel according to FEMM
        else:
            value[0] = mu0

    def value_shape(self):
        return ()

Mu = MuExpression0()

a = (inner(curl(v0), curl(u0)/Mu) + inner(v0, grad(p0)) + inner(grad(q0), u0))*dx
L = -inner(J, v0)*dx

U = Function(PN)
problem = LinearVariationalProblem(a, L, U, bcs)
solver = LinearVariationalSolver(problem)
solver.solve()
print U
u, p = split(U)

# calculate curl(u)
V = VectorFunctionSpace(mesh, "CG", 1)
g = project(curl(u), V)

outfile = File("MAC1.pvd")
outfile << g

The mesh was generated from the following netgen geometry "MAC.geo" using commands
netgen -geofile=MAC.geo -meshfile=MAC.grid -meshfiletype="DIFFPACK Format" -batchmode
dolfin-convert MAC.grid MAC.xml

algebraic3d
solid box_outer = sphere (0, 0, 320e-3; 1000e-3) -bc=1;

hcoil = 2e-3

#trap coil 1
solid coil1 = cylinder (0, 0, -66.5e-3; 0, 0, -51.5e-3; 77e-3)
 and not cylinder (0, 0, -66.5e-3; 0, 0, -51.5e-3; 47e-3)
 and plane (0, 0, -51.5e-3; 0, 0, 1)
 and plane (0, 0, -66.5e-3; 0, 0, -1) -bc=3 -maxh=10e-3;

#trap coil2
solid coil2 = cylinder (0, 0, -16e-3; 0, 0, -1e-3; 77e-3)
 and not cylinder (0, 0, -16e-3; 0, 0, -1e-3; 47e-3)
 and plane (0, 0, -1e-3; 0, 0, 1)
 and plane (0, 0, -16e-3; 0, 0, -1) -bc=3 -maxh=10e-3;

#detector focusing coil
solid coil3 = cylinder (0, 0, 373.5e-3; 0, 0, 389.5e-3; 98e-3)
 and not cylinder (0, 0, 373.5e-3; 0, 0, 389.5e-3; 80e-3)
 and plane (0, 0, 389.5e-3; 0, 0, 1)
 and plane (0, 0, 373.5e-3; 0, 0, -1) -bc=3 -maxh=10e-3;

#detector deflecting coil
solid coil4 = cylinder (0, 106e-3, 373.5e-3; 0, 118e-3, 373.5e-3; 70e-3)
 and not cylinder (0, 106e-3, 373.5e-3; 0, 118e-3, 373.5e-3; 50e-3)
 and plane (0, 118e-3, 373.5e-3; 0, 1, 0)
 and plane (0, 106e-3, 373.5e-3; 0, -1, 0) -bc=3 -maxh=10e-3;

solid shield1 = cylinder (0, 0, -84e-3; 0, 0, -69e-3; 95e-3)
 and not cylinder (0, 0, -84e-3; 0, 0, -69e-3; 37e-3)
 and plane (0, 0, -69e-3; 0, 0, 1)
 and plane (0, 0, -84e-3; 0, 0, -1) -bc=3 -maxh=10e-3;

solid shield2 = cylinder (0, 0, 2e-3; 0, 0, 17e-3; 95e-3)
 and not cylinder (0, 0, 2e-3; 0, 0, 17e-3; 37e-3)
 and plane (0, 0, 17e-3; 0, 0, 1)
 and plane (0, 0, 2e-3; 0, 0, -1) -bc=3 -maxh=10e-3;

solid detector = cylinder (0, 25e-3, 384.9e-3; 0, 25e-3, 390e-3; 7.25e-3)
 and plane (0, 0, 390e-3; 0, 0, 1)
 and plane (0, 0, 384.9e-3; 0, 0, -1) -bc=3 -maxh=10e-3;

solid axial_area = cylinder (0, 0, -100e-3; 0, 0, 330e-3; 20e-3)
 and plane (0, 0, 330e-3; 0, 0, 1)
 and plane (0, 0, -100e-3; 0, 0, -1) -bc=3 -maxh=10e-3;

solid air = box_outer and not coil1 and not coil2 and not coil3 and not coil4 and not shield1 and not shield2 and not detector and not axial_area -bc=3;

tlo coil1 -col=[0, 1, 0];
tlo coil2 -col=[0, 1, 0];
tlo coil3 -col=[0, 1, 0];
tlo coil4 -col=[0, 1, 0];
tlo shield1 -col=[1, 0, 0];
tlo shield2 -col=[1, 0, 0];
tlo detector -col=[1, 1, 0];
tlo axial_area -col=[0, 0, 1] -transparent;
tlo air -col=[0, 0, 1] -transparent;

Result of running the script:

$ python curl.py
Requesting connectivity 3 - 3.
Requesting connectivity 3 - 0.
Requesting connectivity 0 - 3.
Requesting connectivity 3 - 0.
Computing mesh connectivity 0 - 3 from transpose.
Elapsed time: 0.0074029 (compute connectivity 0 - 3)
Computing mesh connectivity 3 - 3 from intersection 3 - 0 - 3.
Elapsed time: 0.267105 (compute connectivity 3 - 3)
Elapsed time: 0.787659 (compute entities dim = 1)
Elapsed time: 0.0127592 (Init dofmap)
Elapsed time: 0.012361 (Init dofmap)
Elapsed time: 0.0136759 (Init dofmap)
Extracted finite element for sub system: FiniteElement('Nedelec 1st kind H(curl)', Cell('tetrahedron', Space(3)), 1, None)
Requesting connectivity 3 - 3.
Elapsed time: 0.899081 (compute entities dim = 2)
Computing sub domain markers for sub domain 0.
Requesting connectivity 2 - 3.
Requesting connectivity 3 - 2.
Computing mesh connectivity 2 - 3 from transpose.
Elapsed time: 0.00682902 (compute connectivity 2 - 3)
Extracted finite element for sub system: FiniteElement('Lagrange', Cell('tetrahedron', Space(3)), 1, None)
Computing sub domain markers for sub domain 0.
Debug: Automatic selection of expression element: <Lagrange vector element of degree ? on a <? cell in ?>: 3 x <CG? on a <? cell in ?>>> [at /usr/lib/python2.7/dist-packages/dolfin/functions/expression.py:688 in _auto_select_element_from_shape()]
Debug: Automatic selection of expression element: <CG? on a <? cell in ?>> [at /usr/lib/python2.7/dist-packages/dolfin/functions/expression.py:688 in _auto_select_element_from_shape()]
Solving linear variational problem.
  Building sparsity pattern over cells [===========================> ] 92.4%
  Building sparsity pattern over cells [==============================] 100.0%
  Matrix of size 95328 x 95328 has 2723628 (0.547471%) nonzero entries.
  Elapsed time: 0.55211 (Build sparsity)
  Elapsed time: 0.193381 (Init tensor)
  Elapsed time: 1.90735e-06 (Delete sparsity)
  Assembling matrix over cells [========> ] 23.1%
  Assembling matrix over cells [=================> ] 46.2%
  Assembling matrix over cells [==============================> ] 80.9%
  Assembling matrix over cells [======================================] 100.0%
  Elapsed time: 3.35065 (Assemble cells)
  Elapsed time: 7.86781e-06 (Build sparsity)
  Elapsed time: 0.00014782 (Init tensor)
  Elapsed time: 9.53674e-07 (Delete sparsity)
  Assembling vector over cells [> ] 1.5%
  Assembling vector over cells [=> ] 2.9%
  Assembling vector over cells [===> ] 8.0%
  Assembling vector over cells [=========> ] 25.3%
  Assembling vector over cells [==================> ] 48.4%
  Assembling vector over cells [===============================> ] 83.0%
  Assembling vector over cells [======================================] 100.0%
  Elapsed time: 5.26788 (Assemble cells)
  Applying boundary conditions to linear system.
  Applying boundary conditions to linear system.
  Symbolic factorization of a matrix of size 95328 x 95328 (UMFPACK).
  LU factorization of a matrix of size 95328 x 95328 (UMFPACK).
Killed

This script used to work with the fenics ppa version from February 2012 (it was 1.0 series, not sure if it was some pre-release...). I also tried now, that it works with python-dolfin 0.9.10 in ubuntu 11.04 (with minor modifications for the VariationalProblem interface etc...). I am sorry for posting the full script with netgen geometry. However, I was not able to reproduce the problem on simple meshes generated by fenics. Before being killed, the script uses less than 30% of memory of my 8GB, so it should not be "out of memory" problem.

So, do you have any idea, what could be the cause of this problem? Maybe some default solver parameters have changed? I tried to change to petsc backend, but it crashes with "out of memory" error in umfpack. I also tried some other solvers and preconditioners without success. I don't have much experience with fem solvers, maybe you could recommend some more suitable solver?

Thank you for any advice
Stepan Roucka