[petsc-users] VTK format?

[Berend van Wachem]

Dear PETSc,

What is the best way to write data from a DMPLEX vector to file, so it 
can be viewed with paraview?
I've found that the standard VTK format works for a serial job, but if 
there is more than 1 processor, the geometry data gets messed up.
I've attached a small working example for a cylinder and the visualised 
geometry with paraview for 1 processors and 4 processors.
Any pointers or "best practice" very much appreciated.
Best regards,

Berend.


#include "petsc.h"

static char help[] = "Simple example to write to .vtk format";

int main(int argc, char **argv)
{
  PetscErrorCode ierr;
  DM dm;
  DM pdm = NULL;
  DMBoundaryType BoundaryType = DM_BOUNDARY_NONE;
  PetscInt RefinementSteps = 3;
  PetscScalar Radius = 0.5, Length = 4.0;
  PetscInt overlap = 1; /* Overlap for neighbour access */
  PetscInt i;
  PetscInt nCoords;
  PetscScalar *coords;
  Vec coordinates;
  Vec xGlobalVector;
  PetscInt dimEmbed;
  const PetscInt dim = 3;
  const PetscInt numFields = 4;
  PetscInt *numComp;
  PetscInt *numDof;
  PetscSection section;
  PetscViewer viewer;


  PetscFunctionBegin;

  ierr = PetscInitialize(&argc, &argv, (char *) 0, help); CHKERRQ(ierr);

  /* Create a Cylinder */
  ierr = DMPlexCreateHexCylinderMesh(PETSC_COMM_WORLD, RefinementSteps, BoundaryType, &dm); CHKERRQ(ierr);

  /* Scale the Cylinder */
  ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr);
  ierr = DMGetCoordinateDim(dm, &dimEmbed); CHKERRQ(ierr);
  ierr = VecGetLocalSize(coordinates, &nCoords); CHKERRQ(ierr);
  ierr = VecGetArray(coordinates, &coords); CHKERRQ(ierr);
  for (i = 0; i < nCoords; i += dimEmbed)
  {
    PetscScalar *coord = &coords[i];
    coord[0] *= Radius;
    coord[1] *= Radius;
    coord[2] *= Length;
  }
  ierr = VecRestoreArray(coordinates, &coords); CHKERRQ(ierr);
  ierr = DMSetCoordinatesLocal(dm, coordinates); CHKERRQ(ierr);


  /* Distribute */
  ierr = DMPlexDistribute(dm, overlap, NULL, &pdm); CHKERRQ(ierr);
  if (pdm)
  {
    ierr = DMDestroy(&dm); CHKERRQ(ierr);
    dm = pdm;
  }

  /* Create a section with 4 variables at cell centers */
  ierr = PetscCalloc2(numFields, &numComp, numFields * (dim + 1), &numDof);
  CHKERRQ(ierr);

  /* Set all primary variables as cell center */
  for (i = 0; i < numFields; i++)
  {
    numComp[i] = 1;
    numDof[i * (dim + 1) + dim] = 1;
  }
  ierr = DMSetNumFields(dm, numFields);
  CHKERRQ(ierr);
  ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, 0, NULL, NULL, NULL, NULL, &section);
  CHKERRQ(ierr);

  /* Push section onto the DM */
  ierr = DMSetLocalSection(dm, section);
  CHKERRQ(ierr);
  ierr = PetscFree2(numComp, numDof);
  CHKERRQ(ierr);
  ierr = PetscSectionDestroy(&section);
  CHKERRQ(ierr);


  /* Create a global vector over the geometry with the section: */
  ierr = DMCreateGlobalVector(dm, &xGlobalVector);CHKERRQ(ierr);
  ierr = VecSet(xGlobalVector, 0.0); CHKERRQ(ierr);


  /* Now write the file to disk in .VTK format */
  ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
  CHKERRQ(ierr);
  ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);
  CHKERRQ(ierr);
  PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
  ierr=PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
  ierr = PetscViewerFileSetName(viewer, "cylindermesh.vtk");
  CHKERRQ(ierr);
  DMView(dm, viewer);
  ierr = VecView(xGlobalVector, viewer);
  CHKERRQ(ierr);



  /* Clean up */
  ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
  ierr = VecDestroy(&xGlobalVector); CHKERRQ(ierr);
  ierr = DMDestroy(&dm); CHKERRQ(ierr);
  ierr = PetscFinalize(); CHKERRQ(ierr);

  PetscFunctionReturn(0);

}