Re: [petsc-users] VTK format?

Berend van Wachem Mon, 07 Sep 2020 06:59:20 -0700

Dear Jed,

Thanks for the suggestion. I've tried writing to .vtu (see attachedcode), but that doesn't seem to work for a section containing multiplefields? Upon running the attached code, I get the error:



$ mpirun -np 1 visualizemesh

[0]PETSC ERROR: --------------------- Error Message--------------------------------------------------------------

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Total number of field components 1 != block size 4

[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.htmlfor trouble shooting.

[0]PETSC ERROR: Petsc Release Version 3.13.3, Jul 01, 2020

[0]PETSC ERROR: visualizemesh on a linux-gcc-opt namedivt14.mb.uni-magdeburg.de by berend Mon Sep 7 15:55:10 2020[0]PETSC ERROR: Configure options --with-debugging=0 COPTFLAGS="-O3-march=native -mtune=native" CXXOPTFLAGS="-O3 -march=native-mtune=native" FOPTFLAGS="-O3 -march=native -mtune=native" --with-clean--download-metis=yes --download-parmetis=yes --download-hdf5--download-p4est --download-triangle --download-tetgen--with-zlib-lib=/usr/lib64/libz.a --with-zlib-include=/usr/include[0]PETSC ERROR: #1 DMPlexVTKWriteAll_VTU() line 334 in/usr/local/petsc-3.13.3/src/dm/impls/plex/plexvtu.c[0]PETSC ERROR: #2 DMPlexVTKWriteAll() line 688 in/usr/local/petsc-3.13.3/src/dm/impls/plex/plexvtk.c[0]PETSC ERROR: #3 PetscViewerFlush_VTK() line 100 in/usr/local/petsc-3.13.3/src/sys/classes/viewer/impls/vtk/vtkv.c[0]PETSC ERROR: #4 PetscViewerFlush() line 26 in/usr/local/petsc-3.13.3/src/sys/classes/viewer/interface/flush.c[0]PETSC ERROR: #5 PetscViewerDestroy() line 113 in/usr/local/petsc-3.13.3/src/sys/classes/viewer/interface/view.c

[0]PETSC ERROR: #6 main() line 102 in visualizemesh.c
[0]PETSC ERROR: No PETSc Option Table entries

[0]PETSC ERROR: ----------------End of Error Message -------send entireerror message to [email protected]

application called MPI_Abort(MPI_COMM_SELF, 102077) - process 0


Or am I making a mistake?

Thanks, best wishes, Berend.




On 2020-09-03 16:53, Jed Brown wrote:

Use the xml format (not the legacy format) by naming your file.vtuinstead of file.vtk


On Thu, Sep 3, 2020, at 8:17 AM, Berend van Wachem wrote:

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.




*Attachments:*

  * visualisemesh-1proc.png
  * visualisemesh-4proc.png
  * visualizemesh.c

#include "petsc.h"

static char help[] = "Simple example to write to .vtu 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);
  ierr = PetscViewerFileSetName(viewer, "cylindermesh.vtu");
  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);

}

Re: [petsc-users] VTK format?

Reply via email to