Writing a P2 function to VTK should use VTK's quadratic elements

[Patrick Farrell]

VTK supports P2 fields on triangles and tetrahedra (http://www.vtk.org/doc/nightly/html/classvtkQuadraticTriangle.html).
However, when I dump a P2 field in Dolfin to VTK, I get a P1 representation:

----------------------------------------------------------------
from dolfin import *

mesh = UnitSquare(5, 5)
P2 = FunctionSpace(mesh, "CG", 2)
P1 = FunctionSpace(mesh, "CG", 1)

print "P2.dim(): ", P2.dim()
print "P1.dim(): ", P1.dim()

v = Function(P2)
v.vector()[:] = 1.0
f = File("p2.pvd")
f << v

import os
os.system("grep NumberOfPoints p2000000.vtu")
----------------------------------------------------------------

The output of this script is:
P2.dim(): 121
P1.dim(): 36
<Piece NumberOfPoints="36" NumberOfCells="50">

i.e. the VTK file produced has 36 dof, not 121.

Would it be possible to make VTK output P2 functions without any projections to P1?

[Garth Wells (garth-wells)]

I implemented this a long time ago, but VTK/ParaView did not draw
triangular cells when using high-order triangles, which was a bit
misleading for visualisation.

Garth

On 6 June 2012 13:52, Patrick Farrell <email address hidden> wrote:
> Public bug reported:
>
> VTK supports P2 fields on triangles and tetrahedra (http://www.vtk.org/doc/nightly/html/classvtkQuadraticTriangle.html).
> However, when I dump a P2 field in Dolfin to VTK, I get a P1 representation:
>
> ----------------------------------------------------------------
> from dolfin import *
>
> mesh = UnitSquare(5, 5)
> P2 = FunctionSpace(mesh, "CG", 2)
> P1 = FunctionSpace(mesh, "CG", 1)
>
> print "P2.dim(): ", P2.dim()
> print "P1.dim(): ", P1.dim()
>
> v = Function(P2)
> v.vector()[:] = 1.0
> f = File("p2.pvd")
> f << v
>
> import os
> os.system("grep NumberOfPoints p2000000.vtu")
> ----------------------------------------------------------------
>
> The output of this script is:
> P2.dim():  121
> P1.dim():  36
> <Piece  NumberOfPoints="36" NumberOfCells="50">
>
> i.e. the VTK file produced has 36 dof, not 121.
>
> Would it be possible to make VTK output P2 functions without any
> projections to P1?
>
> ** Affects: dolfin
>     Importance: Undecided
>         Status: New
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/1009441
>
> Title:
>  Writing a P2 function to VTK should use VTK's quadratic elements
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1009441/+subscriptions

[Patrick Farrell (pefarrell)]

How inconvenient. I guess I'll just stick with P1 for now then ..

[Gerard (g-gorman)]

I just bumped into this as well. I think it really should be an option. If you cap from a higher order element to a lower order element the visualisation is also misleading, particularly for narrow channels. Separately we can also chase the VTK developers about improving higher order element visualisation as this is of wide interest.