XDMF output is scrambled for function spaces with more than 3 components

[Francesco Ballarin]

Dear FEniCS developers, the following MWE shows an issue in exporting a Function to XDMF for function spaces with more than 3 scalar components.

from dolfin import *

mesh = UnitSquareMesh(3,3)

P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1)
P1_3 = MixedElement(P1, P1, P1)
V_3 = FunctionSpace(mesh, P1_3)
P1_4 = MixedElement(P1, P1, P1, P1)
V_4 = FunctionSpace(mesh, P1_4)

u_3 = Function(V_3)
assign(u_3.sub(0), interpolate(Constant(1), V))
assign(u_3.sub(1), interpolate(Constant(2), V))
assign(u_3.sub(2), interpolate(Constant(3), V))
XDMFFile("out_3.xdmf").write(u_3)

u_4 = Function(V_4)
assign(u_4.sub(0), interpolate(Constant(1), V))
assign(u_4.sub(1), interpolate(Constant(2), V))
assign(u_4.sub(2), interpolate(Constant(3), V))
assign(u_4.sub(3), interpolate(Constant(4), V))
XDMFFile("out_4.xdmf").write(u_4)

I can open correctly out_3.xdmf in Paraview, while apparentely only 3 fields out of 4 are saved in out_4.xdmf, with incorrect values.

Of course an easy workaround is to split before exporting. I am opening this ticket so that you can keep track of this issue.

[Michal Habera]

Hi Francesco. Saving MixedElement and EnrichedElement functions is in general not supported in XDMF (and FEniCS). Split the function

u1, u2, u3, u4 = u_4.split(True)

and save one function per one file.

[Francesco Ballarin]

Hi Michal. Thanks for confirming that the workaround of splitting before exporting is the suggested way to go. However, I still think that the user should at least be warned that the output he/she is going to get in the case of u_4 is wrong. Would it be meaningful to raise an error as in VTKFile::results_write, line 257--275? Thanks!

[Michal Habera]

I understand your point. What should be proper politics for IO of mixed functions, @chris_richardson , @garth-wells ? We are able to save multiple fctions in one file, that should not be a problem. But enriched functions (not nodal in general) might be problem.

[Chris Richardson]

I'd be in favour of an error message for Mixed/Enriched Functions.

[Jan Blechta]

EnrichedElement is only non-nodal element in FEniCS and it is non-nodal everytime unless trivially of one subelement.

[Nathan Sime]

Isn't this just happening because Paraview doesn't know how to interpret the 4D field u_4? Consider the following with the same result:

from dolfin import *

mesh = UnitSquareMesh(3,3)

P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
VV = VectorElement("CG", mesh.ufl_cell(), 1)

P1_3 = MixedElement(P1, VV)
V_3 = FunctionSpace(mesh, P1_3)

P1_4 = MixedElement(P1, P1, VV)
V_4 = FunctionSpace(mesh, P1_4)

u_3 = interpolate(Constant((1, 2, 3)), V_3)
XDMFFile("out_3.xdmf").write(u_3)

u_4 = interpolate(Constant((1, 2, 3, 4)), V_4)
XDMFFile("out_4.xdmf").write(u_4)

There are extension modules for Paraview and VisIt for visualisation of 4D data. This is useful for space-time elements.