1: static char help[] = "An example of writing a global Vec from a DMPlex with HDF5 format.\n";

  3: #include <petscdmplex.h>
  4: #include <petscviewerhdf5.h>

  8: int main(int argc, char **argv)
  9: {
 10:   MPI_Comm       comm;
 11:   DM             dm;
 12:   Vec            v, nv, rv, coord;
 13:   PetscBool      test_read = PETSC_FALSE, verbose = PETSC_FALSE, flg;
 14:   PetscViewer    hdf5Viewer;
 15:   PetscInt       dim         = 2;
 16:   PetscInt       numFields   = 1;
 17:   PetscInt       numBC       = 0;
 18:   PetscInt       numComp[1]  = {2};
 19:   PetscInt       numDof[3]   = {2, 0, 0};
 20:   PetscInt       bcFields[1] = {0};
 21:   IS             bcPoints[1] = {NULL};
 22:   PetscSection   section;
 23:   PetscReal      norm;

 26:   PetscInitialize(&argc, &argv, (char *) 0, help);
 27:   comm = PETSC_COMM_WORLD;

 29:   PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");
 30:   PetscOptionsBool("-test_read","Test reading from the HDF5 file","",PETSC_FALSE,&test_read,NULL);
 31:   PetscOptionsBool("-verbose","print the Vecs","",PETSC_FALSE,&verbose,NULL);
 32:   PetscOptionsInt("-dim","the dimension of the problem","",2,&dim,NULL);
 33:   PetscOptionsEnd();

 35:   DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, &dm);
 36:   DMGetDimension(dm, &dim);
 37:   numDof[0] = dim;
 38:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, &section);
 39:   DMSetDefaultSection(dm, section);
 40:   PetscSectionDestroy(&section);
 41:   DMSetUseNatural(dm, PETSC_TRUE);
 42:   {
 43:     DM dmDist;

 45:     DMPlexDistribute(dm, 0, NULL, &dmDist);
 46:     if (dmDist) {
 47:       DMDestroy(&dm);
 48:       dm   = dmDist;
 49:     }
 50:   }

 52:   DMCreateGlobalVector(dm, &v);
 53:   PetscObjectSetName((PetscObject) v, "V");
 54:   DMGetCoordinates(dm, &coord);
 55:   VecCopy(coord, v);

 57:   if (verbose) {
 58:     PetscInt size, bs;

 60:     VecGetSize(v, &size);
 61:     VecGetBlockSize(v, &bs);
 62:     PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%d\tblock size=%d\n", size, bs);
 63:     VecView(v, PETSC_VIEWER_STDOUT_WORLD);
 64:   }

 66:   DMCreateGlobalVector(dm, &nv);
 67:   PetscObjectSetName((PetscObject) nv, "NV");
 68:   DMPlexGlobalToNaturalBegin(dm, v, nv);
 69:   DMPlexGlobalToNaturalEnd(dm, v, nv);

 71:   if (verbose) {
 72:     PetscInt size, bs;

 74:     VecGetSize(nv, &size);
 75:     VecGetBlockSize(nv, &bs);
 76:     PetscPrintf(PETSC_COMM_WORLD, "====  V in natural ordering. size==%d\tblock size=%d\n", size, bs);
 77:     VecView(nv, PETSC_VIEWER_STDOUT_WORLD);
 78:   }

 80:   VecViewFromOptions(v, NULL, "-global_vec_view");

 82:   if (test_read) {
 83:     DMCreateGlobalVector(dm, &rv);
 84:     PetscObjectSetName((PetscObject) rv, "V");
 85:     /* Test native read */
 86:     PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
 87:     PetscViewerPushFormat(hdf5Viewer, PETSC_VIEWER_NATIVE);
 88:     VecLoad(rv, hdf5Viewer);
 89:     PetscViewerPopFormat(hdf5Viewer);
 90:     PetscViewerDestroy(&hdf5Viewer);
 91:     if (verbose) {
 92:       PetscInt size, bs;

 94:       VecGetSize(rv, &size);
 95:       VecGetBlockSize(rv, &bs);
 96:       PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
 97:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
 98:     }
 99:     VecEqual(rv, v, &flg);
100:     if (flg) {
101:       PetscPrintf(PETSC_COMM_WORLD, "V and RV are equal\n");
102:     } else {
103:       PetscPrintf(PETSC_COMM_WORLD, "V and RV are not equal\n\n");
104:       VecAXPY(rv, -1.0, v);
105:       VecNorm(rv, NORM_INFINITY, &norm);
106:       PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm);
107:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
108:     }
109:     /* Test raw read */
110:     PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
111:     VecLoad(rv, hdf5Viewer);
112:     PetscViewerDestroy(&hdf5Viewer);
113:     if (verbose) {
114:       PetscInt size, bs;

116:       VecGetSize(rv, &size);
117:       VecGetBlockSize(rv, &bs);
118:       PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
119:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
120:     }
121:     VecEqual(rv, nv, &flg);
122:     if (flg) {
123:       PetscPrintf(PETSC_COMM_WORLD, "NV and RV are equal\n");
124:     } else {
125:       PetscPrintf(PETSC_COMM_WORLD, "NV and RV are not equal\n\n");
126:       VecAXPY(rv, -1.0, v);
127:       VecNorm(rv, NORM_INFINITY, &norm);
128:       PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm);
129:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
130:     }
131:     VecDestroy(&rv);
132:   }
133:   VecDestroy(&nv);
134:   VecDestroy(&v);
135:   DMDestroy(&dm);
136:   PetscFinalize();
137:   return 0;
138: }