A strange error in the latest fenics with import a dolfin converting .xml file

Question

I debug code following:

mesh = Mesh("test.xml")
element = VectorElement('P', 'triangle', 1, dim=4)
ME = FunctionSpace(mesh, element)

# Define test functions
du = TrialFunction(ME)
v_c, v_bs, v_bm, v_bd = TestFunctions(ME)

# Define functions
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

# # Split system functions to access components
c, bs, bm, bd = split(u)
c0, bs0, bm0, bd0 = split(u0)

grad(c) #debug test!

where test.xml is converted from a gmsh test.msh file using dilfin-convert(the msh file is created by PointWise and saved as .msh format, it shows a cube with a very small circle zone at the center of bottom face). when it debug down to the last line, it throw a exception: "cannot determine the geometric dimension from experession". I traced it and seems the Grad(f) function return type has some unkown problems, but c is recognized well. It is strange that 1) within 1.6 version, the code is fine; 2) if I change the first line with mesh = UnitSquareMesh(16, 16), no error.
So I have no idea that whether the error due to the incompatible problem between the new version fenics and dolfin-convert module, or duo to the .msh file itself, or I missed some fenics grammer?
Here are test.xml and test.msh.

Comments

You should double check this--I fixed the error, but I am not sure if this is the right way to solve it.

I just specified a VectorFunctionSpace instead of a FunctionSpace. I didn't see any other use for element in the code, so I commented it out. If you need it for something, you might be able to extract it with ME.element().

from dolfin import *

mesh = Mesh("test.xml")

# element = VectorElement('P', 'triangle', 1, dim=4)
# ME = FunctionSpace(mesh, element)

ME = VectorFunctionSpace(mesh, 'P', 1, dim=4)

# Define test functions
du = TrialFunction(ME)
v_c, v_bs, v_bm, v_bd = TestFunctions(ME)

# Define functions
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

# # Split system functions to access components
c, bs, bm, bd = split(u)
c0, bs0, bm0, bd0 = split(u0)

grad(c) #debug test!

---

Thank you very much! I will try it. The problem may be caused by inertial mechanism, hopefully fenics officially check the problem and fix it.

Answer 1

OK, I found the problem, may be. In the demo "Cahn-Hilliard equation", it shows

mesh = UnitSquareMesh(96, 96)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

So, I thought in my code, it should be

mesh = Mesh("test.xml")
element = VectorElement('P', mesh.ufl_cell(), 1, dim=4)

Also, it's wrong to use "triangle" here because in the mesh cell is 3D.