1: #if !defined(PETSC_USE_CXX_COMPLEX_FLOAT_WORKAROUND)
  2: #define PETSC_USE_CXX_COMPLEX_FLOAT_WORKAROUND 1
  3: #endif
  4: /*
  5:    Demonstrates using the HDF5 viewer with a DMDA Vec
  6:  - create a global vector containing a gauss profile (exp(-x^2-y^2))
  7:  - write the global vector in a hdf5 file

  9:    The resulting file gauss.h5 can be viewed with Visit (an open source visualization package)
 10:    Or with some versions of MATLAB with data=hdfread('gauss.h5','pressure'); mesh(data);

 12:    The file storage of the vector is independent of the number of processes used.
 13:  */

 15:  #include <petscdm.h>
 16:  #include <petscdmda.h>
 17:  #include <petscsys.h>
 18: #include <petscviewerhdf5.h>

 20: static char help[] = "Test to write HDF5 file from PETSc DMDA Vec.\n\n";

 22: int main(int argc,char **argv)
 23: {
 25:   DM             da2D;
 26:   PetscInt       i,j,ixs, ixm, iys, iym;;
 27:   PetscViewer    H5viewer;
 28:   PetscScalar    xm    = -1.0, xp=1.0;
 29:   PetscScalar    ym    = -1.0, yp=1.0;
 30:   PetscScalar    value = 1.0,dx,dy;
 31:   PetscInt       Nx    = 40, Ny=40;
 32:   Vec            gauss,input;
 33:   PetscScalar    **gauss_ptr;
 34:   PetscReal      norm;
 35:   const char     *vecname;

 37:   dx=(xp-xm)/(Nx-1);
 38:   dy=(yp-ym)/(Ny-1);

 40:   /* Initialize the Petsc context */
 41:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 42:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,Nx,Ny,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da2D);
 43:   DMSetFromOptions(da2D);
 44:   DMSetUp(da2D);

 46:   /* Set the coordinates */
 47:   DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);

 49:   /* Declare gauss as a DMDA component */
 50:   DMCreateGlobalVector(da2D,&gauss);
 51:   PetscObjectSetName((PetscObject) gauss, "pressure");

 53:   /* Initialize vector gauss with a constant value (=1) */
 54:   VecSet(gauss,value);

 56:   /* Get the coordinates of the corners for each process */
 57:   DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0);

 59:   /* Build the gaussian profile (exp(-x^2-y^2)) */
 60:   DMDAVecGetArray(da2D,gauss,&gauss_ptr);
 61:   for (j=iys; j<iys+iym; j++) {
 62:     for (i=ixs; i<ixs+ixm; i++) {
 63:       gauss_ptr[j][i]=PetscExpScalar(-(xm+i*dx)*(xm+i*dx)-(ym+j*dy)*(ym+j*dy));
 64:     }
 65:   }
 66:   DMDAVecRestoreArray(da2D,gauss,&gauss_ptr);

 68:   /* Create the HDF5 viewer */
 69:   PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_WRITE,&H5viewer);
 70:   PetscViewerSetFromOptions(H5viewer);

 72:   /* Write the H5 file */
 73:   VecView(gauss,H5viewer);

 75:   /* Close the viewer */
 76:   PetscViewerDestroy(&H5viewer);

 78:   VecDuplicate(gauss,&input);
 79:   PetscObjectGetName((PetscObject)gauss,&vecname);
 80:   PetscObjectSetName((PetscObject)input,vecname);

 82:   /* Create the HDF5 viewer for reading */
 83:   PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_READ,&H5viewer);
 84:   PetscViewerSetFromOptions(H5viewer);
 85:   VecLoad(input,H5viewer);
 86:   PetscViewerDestroy(&H5viewer);

 88:   VecAXPY(input,-1.0,gauss);
 89:   VecNorm(input,NORM_2,&norm);
 90:   if (norm > 1.e-6) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Vec read in does not match vector written out");

 92:   VecDestroy(&input);
 93:   VecDestroy(&gauss);
 94:   DMDestroy(&da2D);
 95:   PetscFinalize();
 96:   return ierr;
 97: }


100: /*TEST

102:       build:
103:          requires: hdf5 !define(PETSC_USE_CXXCOMPLEX)

105:       test:
106:          nsize: 4

108:       test:
109:          nsize: 4
110:          suffix: 2
111:          args: -viewer_hdf5_base_dimension2
112:          output_file: output/ex10_1.out

114:       test:
115:          nsize: 4
116:          suffix: 3
117:          args: -viewer_hdf5_sp_output
118:          output_file: output/ex10_1.out

120: TEST*/