[firedrake] Problem with paraview and vector functions

[Lawrence Mitchell]

Dear Luca,

> On 12 Mar 2021, at 10:37, luca.mechelli at uni-konstanz.de wrote:
> 
> Dear all,
> 
> I encountered a problem when visualizing a function of a vector function space in paraview, if the order of the finite element is higher than 1. For example, if I run the following code, paraview shows the error message: "Cell type UnknownClass(69) is not a 3D cell. Missing original point id arrays. "
> 
> 
> mesh = UnitSquareMesh(10, 10)
> V = VectorFunctionSpace(mesh, "CG", 2)
> v = Function(V)
> outfile = File("test.pvd")
> outfile.write(v)
> 
> 
> If I replace the instruction "V = VectorFunctionSpace(mesh, "CG", 2)" with "V = VectorFunctionSpace(mesh, "CG", 1)", then I can visualize correctly the function in paraview. Could you please explain me which is my mistake?
> 
> As temporary solution I have defined an additional space V_2 = VectorFunctionSpace(mesh, "CG", 1), on which I interpolate the function v, although I am not really satisfied with the outcome.


I suspect this is because your version of paraview is a bit old. For quadratic (and higher degree) elements, we write VTU files using Paraview's (relatively) recent native support for these elements [this kitware post describes things https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/]. See also https://www.firedrakeproject.org/visualisation.html#visualising-high-order-data

I think support for this was introduced in Paraview 5.5, so if your version is older then things won't work.

One way around this is to interpolate the solution to the low order space. You can have the File object manage this by saying:

outfile = File("test.pvd", target_degree=1)

but this will lose information.

Thanks,

Lawrence