[petsc-users] Visualizing structured cell-centered data VTK

[Jed Brown]

"Buesing, Henrik" <hbuesing at eonerc.rwth-aachen.de> writes:

>> >> > I have structured cell-centered data and would like to visualize
>> >> > this with
>> >> Paraview. Up to now I use PetscViewerVTKOpen and VecView to write
>> >> data in *.vts format. I would like to tell PETSc that the fieldtype
>> >> is PETSC_VTK_CELL_FIELD. I have found PetscViewerVTKAddField.
>> >> >
>> >> > Is this the way to go? I was thinking maybe a DMDASetFieldType
>> >> > exists, but
>> >> did not find any. If yes, what is the PetscViewerVTKWriteFunction I
>> >> need to provide?
>> >>
>> >> DMDA does not explicitly support distinguishing between cell and
>> >> point values.  PetscViewerVTKAddField is a developer level routine
>> >> and you would need to implement a function similar to
>> >> DMDAVTKWriteAll_VTS (not at all trivial and you need to read the code
>> >> because it is responsible for almost everything).
>> >
>> > I am looking at src/sys/classes/viewer/impls/vtk/vtkv.c. There is a
>> reference to PETSC_VTK_POINT_FIELD vs. PETSC_VTK_CELL_FIELD. Judging
>> from the output I get, I was assuming fieldtype=PETSC_VTK_POINT_FIELD. I
>> would be totally fine with replacing POINT by CELL everywhere, since all my
>> data is cell-centered.
>> 
>> I think coordinates need to be PointData, not coordinates of cell centroids.
>> DMDA doesn't have that concept.  (Maybe it should, but adding it is no small
>> task and hacking the output is likely to create a lot of edge cases.)
>> 
> I had a look at DMDAVTKWriteAll_VTS. You are totally right that this
> is not generic. www.vtk.org/VTK/img/file-formats.pdf describes the
> StructuredGrid format. "Points" contain only one DataArray specifying
> the coordinates whereas Cells contain three DataArrays with
> connectivity, offsets and types.

Are you mixing up CellData from the StructuredGrid XML spec with Cells
from the UnstructuredGrid?

> What if I use DMDAVTKWriteAll_VTR and write a RectliniearGrid? Then I could just replace PointData with CellData and "hack" only the extent which should leave no edge cases...

I think it's feasible.  You could use DMDAGetInterpolationType() to
determine which variant to use.