Write VTK and XDMF formats for quadratic and other non-vertex based elements

[Johann Heller]

Writing a file to vtk using:

file << u

interpolates all data to the vertices even though the actual DOFs might be located somewhere else (e.g. on the edges). I guess it would be nicer to write to vtk without any interpolation. At least there should be an option for no interpolation.

[Prof Garth Wells]

How do you propose doing this? VTK works with vertex data.

[Johann Heller]

I did it for myself by extracting the location of the DOFs (in my case on the edges) by:

dofmap.tabulate_all_coordinates(mesh)

and making a new "Dof-mesh" where the vertices are located at the DOF-locations. It works for me for an extremely small example. In this way the interpolation is avoided.

Is ask this because the forced interpolation might lead to artifacts that imply a wrong solution.

This question is related to issue #488

[Chris Richardson]

Other formats (e.g. XDMF) also support non-vertex data. The main problem has been visualisation - what software are you using to visualise your Edge-based function?

Sorry - after re-reading your reply, I see you are placing extra vertices. I guess this is possible, though may be difficult for higher order elements - it would require calculating a submesh for each possible triangle or tetrahedron (e.g. what do you do with a CG4 tetrahedron?). At least for quadratic elements, it is not too difficult, and might be considered in the future.

[Johann Heller]

In the end i dit it in paraview. I have exported my stuff to Matlab (since i am a little more familiar with that than Python) and wrote a script that writes the edge data to .vtk making a new mesh where the vertices are located at the edges ... pretty under-developed but it works :-)

[Chris Richardson]

It is possible, in fact, to get ParaView to display quadratic elements by supplying the edge values on each cell. It could be done with some effort. Other challenges exist for different element types (e.g. RT) where the vector values are defined at Edges.

[Johann Heller]

Thanks a lot !

It is nice to see that the community is so active and that normal users can make suggestions for the development.

[Martin Sandve Alnæs]

I thought I made an issue for this before but I can't find it so maybe I just talked about it.

It should be fairly simple to extend VTK support for both file and plot to discontinuous functions by duplicating vertices for each cell. All linear non-Lagrange finite elements can be represented exactly in a DG1 basis.

An orthogonal step is quadratic elements by providing edge values as Chris suggested.

With higher order geometries in place later, we can add edge coordinate values as well.

Together these features should make the plotting quite flexible.

[Prof Garth Wells]

@martinal Duplicating vertices can lead to artefacts in ParaView visualisation, especially with transparency.

[Johan Hake]

I am pretty sure Simone @peppu has solved this in his fenicshotools repo. Not sure exactly how to use it as it is sparse on documentation. It is also written to work with meshes of higher order geometries.

[Simone Pezzuto]

@johanhake The package has at least three purposes:

1. Generate high order geometries with Gmsh, saving everything into a .h5 file that can be loaded (using load_geometry) to obtain a Domain and boundary markers;

2. Linearise the Domain (by adding extra vertices), so that you get ordinary dolfin Mesh and Functions from a high-order Lagrangian Domain of any order. This is done by constructing a new mesh using the dofs as vertices (not trivial for tetrahedral meshes);

3. Convert a dolfin Domain to a vtkUnstructuredGrid of the right type, so with quadratic elements if the domain is quadratic (cubic or more are not supported). Then you can add fields and save it so that Paraview is able to read it. It is also possible to export a linear Mesh plus a quadratic FunctionSpace, not just isoparametric.

Of course this is valid only for Lagrangian functions and mappings, so edge-based elements won't work. Also for DG elements is not trivial (except DG0) because when you duplicate the vertices you also lose the capability to correctly perform some post-processing in Paraview, since the mesh is not connected.

I found some comments on VTK mailing-list about this, and one solution is to either define you own VTK cell type (so that it is able to handle geometry queries), or export the dofs as quadrature points (see here).

[Prof Garth Wells]

The issue should be split in two. Quadratic elements should be straightforward since it's supported by ParaView. Edge-based element support is another story.

[Chris Richardson]

Just made some quadratic elements with XDMF. Now waiting on Issue #321 to integrate into dolfin.

[Martin Sandve Alnæs]

Nice. I updated #321 to a be a bit more constructive description of the latest plans.

[Prof Garth Wells]

Nice, but no need to wait for #321 since we can put quadratic solution data on the mesh rather than using nodal (P1) interpolations. Will need a way to robustly determine that an element is P2 Lagrange.

[Martin Sandve Alnæs]

I don't understand what you mean by "put quadratic solution data on the mesh rather than using nodal (P1) interpolations"? Did you read my updated #321 description (CoordinateFunction is gone)?

[Prof Garth Wells]

@martinal Plot the solution data for a P2 element at the six nodes of an element (in 2D) rather than just the vertex values, which is the topic of this Issue. This can be done independently of any cell mapping.

[Martin Sandve Alnæs]

Ok, so just for P2 visualisation, just sampling the edge midpoints on an affine mesh to get nodes to attach P2 dofs to for the xdmf file writing?

[Chris Richardson]

OK, there are multiple issues which could be split from here.

1. Plot P2 data on a regular mesh
2. Plot something on an isoparametric mesh (e.g. quadratic)
3. Plot other things, like RT, which might be edge based, or any other weird FunctionSpace

I've probably confused things by heading for (2) here Maybe I'll split off a few more issues.

[Martin Sandve Alnæs]

I also assumed you went for 2 here, given the quadratic edges in the plot.

But Garth is right, this can be useful before quadratic mesh is in place by just writing the facet midpoints of the affine cells and the facet midpoint values of the function to the vtk file for all elements of degree >= 2.

[Prof Garth Wells]

For clarity, let's split this report into the three points listed by @chris_richardson with some links to each report in the description.

[Prof Garth Wells]

@martinal We might want a new UFC function to get mid-point values efficiently (akin to interpolate_vertex_values), rather than by evaluating the FE basis.

[Martin Sandve Alnæs]

@garth-wells I've previously considered adding evaluation of FE tables in facet and cell midpoints to ffc. Then midpoint values can easily be made available in UFL as well.

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)

[Chris Richardson]

This issue is now being picked up by a project to improve XDMF support for visualising higher order elements.