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

[Buesing, Henrik]

> >> >> > 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?
Ah, yes. You are right! Ok wrong reasoning, but maybe good conclusion nevertheless. 

> > 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.

I will try it and see where I can go with it.