1: /*
  2:    Demonstrates using the HDF5 viewer with a DMDA Vec
  3:  - create a global vector containing a gauss profile (exp(-x^2-y^2))
  4:  - write the global vector in a hdf5 file 

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

  9:    The file storage of the vector is independent of the number of processes used.
 10:  */

 12: #include <math.h>
 13: #include <petscdmda.h>
 14: #include <petscsys.h>

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

 18: int main(int argc,char **argv)
 19: {
 21:   DM             da2D;
 22:   PetscInt       i,j,ixs, ixm, iys, iym;;
 23:   PetscViewer    H5viewer;
 24:   PetscScalar    xm=-1.0, xp=1.0;
 25:   PetscScalar    ym=-1.0, yp=1.0;
 26:   PetscScalar    value=1.0,dx,dy;
 27:   PetscInt       Nx=40, Ny=40;
 28:   Vec            gauss;
 29:   PetscScalar    **gauss_ptr;
 30: 
 31:   dx=(xp-xm)/(Nx-1);
 32:   dy=(yp-ym)/(Ny-1);
 33: 
 34:   // Initialize the Petsc context
 35:   PetscInitialize(&argc,&argv,(char*)0,help);
 36: 
 37:   // Build of the DMDA
 38:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,Nx,Ny,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da2D);
 39: 
 40:   // Set the coordinates
 41:   DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
 42: 
 43:   // Declare gauss as a DMDA component
 44:   DMCreateGlobalVector(da2D,&gauss);
 45:   PetscObjectSetName((PetscObject) gauss, "pressure");
 46: 
 47:   // Initialize vector gauss with a constant value (=1)
 48:   VecSet(gauss,value);
 49: 
 50:   // Get the coordinates of the corners for each process
 51:   DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0);
 52: 
 53:   /* Build the gaussian profile (exp(-x^2-y^2)) */
 54:   DMDAVecGetArray(da2D,gauss,&gauss_ptr);
 55:   for (j=iys; j<iys+iym; j++){
 56:     for (i=ixs; i<ixs+ixm; i++){
 57:       gauss_ptr[j][i]=exp(-(xm+i*dx)*(xm+i*dx)-(ym+j*dy)*(ym+j*dy));
 58:     }
 59:   }
 60:   DMDAVecRestoreArray(da2D,gauss,&gauss_ptr);
 61: 
 62:   // Create the HDF5 viewer
 63:   PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_WRITE,&H5viewer);
 64: 
 65:   // Write the H5 file
 66:   VecView(gauss,H5viewer);
 67: 
 68:   // Cleaning stage
 69:   PetscViewerDestroy(&H5viewer);
 70:   VecDestroy(&gauss);
 71:   DMDestroy(&da2D);
 72:   PetscFinalize();
 73:     return 0;
 74: }