[firedrake] Problem with paraview and vector functions

[luca.mechelli at uni-konstanz.de]

Dear Lawrence,
 
Il giorno Venerdi, Marzo 12, 2021 11:47 CET, Lawrence Mitchell <wence at gmx.li> ha scritto: 
 
> 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.
> 

You are right, my Paraview version was too old (5.4). Updating it has fixed the problem (sometimes I forget the simplest solution :D sorry!)

> 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 
 
Thank you for the help,

Luca