[petsc-dev] HDF5 VecView bug for DMPlex global vectors.

[Ataollah Mesgarnejad]

Dear developers,

I came across this bug when playing with HDF5 VecView. HDF5 VecView crashes
if the default section for DMPlex only has dofs on the edges. Attached is a
small code to reproduce the error.

Best,
Ata

/**
* @file
*
* Created by Ataollah Mesgarnejad on 3/26/15.
*
* This is a test for DMPlex global to natural and natural to global ordering.
*
*/


static char help[] = "An example of the usage of PetscSF to scatter data back and forth between a \
a dist DM and its uniformly distributed parallel counterpart.\n";

#include <petscdmplex.h>
#include <petscsnes.h>
#include <petscsf.h>
#include <exodusII.h>
#include <petsc-private/dmimpl.h>     /*I      "petscdm.h"     I*/
#include <petsc-private/vecimpl.h>
#include <petscsf.h>
#include <petscviewerhdf5.h>
#define VERBOSE

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
  PetscErrorCode      ierr;
  DM                  dm,distDM;
  PetscViewer         hdf5Viewer;
  PetscMPIInt         numproc,rank;
  PetscViewer         stdoutViewer;
  PetscInt            dim;
  PetscInt            numFields = 1;
  PetscInt            numComp[1] = {1};
  PetscInt            numDof[6] = {0,1,0}; /*CRASHES */
  //PetscInt            numDof[6] = {1,1,0}; /*DOES NOT CRASH */
  PetscInt            bcFields[1] = {0},numBC=0;
  IS                  bcPoints[1] = {NULL};
  PetscSection        seqSection,distSection;
  Vec                 V,distV;
  PetscSF             pointSF;
  MPI_Comm            comm;

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
  ierr = PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&stdoutViewer);CHKERRQ(ierr);
  comm = PETSC_COMM_WORLD;

  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","PetscSF Test Options","none");CHKERRQ(ierr);
  ierr = PetscOptionsEnd();

  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&numproc);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
  ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, &dm);
  ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);

  ierr = DMPlexCreateSection(dm,dim,numFields,numComp,numDof,numBC,bcFields,bcPoints,NULL,&seqSection);CHKERRQ(ierr);
  ierr = DMSetDefaultSection(dm,seqSection);CHKERRQ(ierr);

  ierr = DMPlexDistribute(dm,0,&pointSF,&distDM);CHKERRQ(ierr);

  ierr = DMPlexCreateSection(distDM,dim,numFields,numComp,numDof,numBC,bcFields,bcPoints,NULL,&distSection);CHKERRQ(ierr);    
  ierr = DMSetDefaultSection(distDM,distSection);CHKERRQ(ierr);

  ierr = DMGetGlobalVector(dm,&V); CHKERRQ(ierr);
  ierr = DMGetGlobalVector(distDM,&distV);CHKERRQ(ierr);

  ierr = PetscViewerHDF5Open(comm,"test.h5", FILE_MODE_WRITE, &hdf5Viewer);CHKERRQ(ierr);
  ierr = VecView(distV, hdf5Viewer);CHKERRQ(ierr);
  PetscViewerDestroy(&hdf5Viewer);

  ierr = DMRestoreGlobalVector(dm,&V);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(distDM,&distV);CHKERRQ(ierr);

  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  if (distDM) ierr = DMDestroy(&distDM);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return 0;
}