Actual source code: plexhdf5xdmf.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/isimpl.h>
  3: #include <petsc/private/vecimpl.h>
  4: #include <petsclayouthdf5.h>

  6: #if defined(PETSC_HAVE_HDF5)
  7: static PetscErrorCode SplitPath_Private(char path[], char name[])
  8: {
  9:   char *tmp;

 11:   PetscStrrchr(path, '/', &tmp);
 12:   PetscStrcpy(name, tmp);
 13:   if (tmp != path) {
 14:     /* '/' found, name is substring of path after last occurrence of '/'. */
 15:     /* Trim the '/name' part from path just by inserting null character. */
 16:     tmp--;
 17:     *tmp = '\0';
 18:   } else {
 19:     /* '/' not found, name = path, path = "/". */
 20:     PetscStrcpy(path, "/");
 21:   }
 22:   return 0;
 23: }

 25: /*
 26:   - invert (involute) cells of some types according to XDMF/VTK numbering of vertices in a cells
 27:   - cell type is identified using the number of vertices
 28: */
 29: static PetscErrorCode DMPlexInvertCells_XDMF_Private(DM dm)
 30: {
 31:   PetscInt     dim, *cones, cHeight, cStart, cEnd, p;
 32:   PetscSection cs;

 34:   DMGetDimension(dm, &dim);
 35:   if (dim != 3) return 0;
 36:   DMPlexGetCones(dm, &cones);
 37:   DMPlexGetConeSection(dm, &cs);
 38:   DMPlexGetVTKCellHeight(dm, &cHeight);
 39:   DMPlexGetHeightStratum(dm, cHeight, &cStart, &cEnd);
 40:   for (p = cStart; p < cEnd; p++) {
 41:     PetscInt numCorners, o;

 43:     PetscSectionGetDof(cs, p, &numCorners);
 44:     PetscSectionGetOffset(cs, p, &o);
 45:     switch (numCorners) {
 46:     case 4:
 47:       DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cones[o]);
 48:       break;
 49:     case 6:
 50:       DMPlexInvertCell(DM_POLYTOPE_TRI_PRISM, &cones[o]);
 51:       break;
 52:     case 8:
 53:       DMPlexInvertCell(DM_POLYTOPE_HEXAHEDRON, &cones[o]);
 54:       break;
 55:     }
 56:   }
 57:   return 0;
 58: }

 60: PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM dm, PetscViewer viewer)
 61: {
 62:   Vec         coordinates;
 63:   IS          cells;
 64:   PetscInt    spatialDim, topoDim = -1, numCells, numVertices, NVertices, numCorners;
 65:   PetscMPIInt rank;
 66:   MPI_Comm    comm;
 67:   char        topo_path[PETSC_MAX_PATH_LEN] = "/viz/topology/cells", topo_name[PETSC_MAX_PATH_LEN];
 68:   char        geom_path[PETSC_MAX_PATH_LEN] = "/geometry/vertices", geom_name[PETSC_MAX_PATH_LEN];
 69:   PetscBool   seq                           = PETSC_FALSE;

 71:   PetscObjectGetComm((PetscObject)dm, &comm);
 72:   MPI_Comm_rank(comm, &rank);

 74:   PetscOptionsBegin(PetscObjectComm((PetscObject)dm), ((PetscObject)dm)->prefix, "DMPlex HDF5/XDMF Loader Options", "PetscViewer");
 75:   PetscOptionsString("-dm_plex_hdf5_topology_path", "HDF5 path of topology dataset", NULL, topo_path, topo_path, sizeof(topo_path), NULL);
 76:   PetscOptionsString("-dm_plex_hdf5_geometry_path", "HDF5 path to geometry dataset", NULL, geom_path, geom_path, sizeof(geom_path), NULL);
 77:   PetscOptionsBool("-dm_plex_hdf5_force_sequential", "force sequential loading", NULL, seq, &seq, NULL);
 78:   PetscOptionsEnd();

 80:   SplitPath_Private(topo_path, topo_name);
 81:   SplitPath_Private(geom_path, geom_name);
 82:   PetscInfo(dm, "Topology group %s, name %s\n", topo_path, topo_name);
 83:   PetscInfo(dm, "Geometry group %s, name %s\n", geom_path, geom_name);

 85:   /* Read topology */
 86:   PetscViewerHDF5PushGroup(viewer, topo_path);
 87:   ISCreate(comm, &cells);
 88:   PetscObjectSetName((PetscObject)cells, topo_name);
 89:   if (seq) {
 90:     PetscViewerHDF5ReadSizes(viewer, topo_name, NULL, &numCells);
 91:     PetscLayoutSetSize(cells->map, numCells);
 92:     numCells = rank == 0 ? numCells : 0;
 93:     PetscLayoutSetLocalSize(cells->map, numCells);
 94:   }
 95:   ISLoad(cells, viewer);
 96:   ISGetLocalSize(cells, &numCells);
 97:   ISGetBlockSize(cells, &numCorners);
 98:   PetscViewerHDF5ReadAttribute(viewer, topo_name, "cell_dim", PETSC_INT, &topoDim, &topoDim);
 99:   PetscViewerHDF5PopGroup(viewer);
100:   numCells /= numCorners;

102:   /* Read geometry */
103:   PetscViewerHDF5PushGroup(viewer, geom_path);
104:   VecCreate(comm, &coordinates);
105:   PetscObjectSetName((PetscObject)coordinates, geom_name);
106:   if (seq) {
107:     PetscViewerHDF5ReadSizes(viewer, geom_name, NULL, &numVertices);
108:     PetscLayoutSetSize(coordinates->map, numVertices);
109:     numVertices = rank == 0 ? numVertices : 0;
110:     PetscLayoutSetLocalSize(coordinates->map, numVertices);
111:   }
112:   VecLoad(coordinates, viewer);
113:   VecGetLocalSize(coordinates, &numVertices);
114:   VecGetSize(coordinates, &NVertices);
115:   VecGetBlockSize(coordinates, &spatialDim);
116:   PetscViewerHDF5PopGroup(viewer);
117:   numVertices /= spatialDim;
118:   NVertices /= spatialDim;

120:   PetscInfo(NULL, "Loaded mesh dimensions: numCells %" PetscInt_FMT " numCorners %" PetscInt_FMT " numVertices %" PetscInt_FMT " spatialDim %" PetscInt_FMT "\n", numCells, numCorners, numVertices, spatialDim);
121:   {
122:     const PetscScalar *coordinates_arr;
123:     PetscReal         *coordinates_arr_real;
124:     const PetscInt    *cells_arr;
125:     PetscSF            sfVert = NULL;
126:     PetscInt           i;

128:     VecGetArrayRead(coordinates, &coordinates_arr);
129:     ISGetIndices(cells, &cells_arr);

131:     if (PetscDefined(USE_COMPLEX)) {
132:       /* convert to real numbers if PetscScalar is complex */
133:       /*TODO More systematic would be to change all the function arguments to PetscScalar */
134:       PetscMalloc1(numVertices * spatialDim, &coordinates_arr_real);
135:       for (i = 0; i < numVertices * spatialDim; ++i) {
136:         coordinates_arr_real[i] = PetscRealPart(coordinates_arr[i]);
137:         PetscAssert(PetscImaginaryPart(coordinates_arr[i]) == 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector of coordinates contains complex numbers but only real vectors are currently supported.");
138:       }
139:     } else coordinates_arr_real = (PetscReal *)coordinates_arr;

141:     DMSetDimension(dm, topoDim < 0 ? spatialDim : topoDim);
142:     DMPlexBuildFromCellListParallel(dm, numCells, numVertices, NVertices, numCorners, cells_arr, &sfVert, NULL);
143:     DMPlexInvertCells_XDMF_Private(dm);
144:     DMPlexBuildCoordinatesFromCellListParallel(dm, spatialDim, sfVert, coordinates_arr_real);
145:     VecRestoreArrayRead(coordinates, &coordinates_arr);
146:     ISRestoreIndices(cells, &cells_arr);
147:     PetscSFDestroy(&sfVert);
148:     if (PetscDefined(USE_COMPLEX)) PetscFree(coordinates_arr_real);
149:   }
150:   ISDestroy(&cells);
151:   VecDestroy(&coordinates);

153:   /* scale coordinates - unlike in DMPlexLoad_HDF5_Internal, this can only be done after DM is populated */
154:   {
155:     PetscReal lengthScale;

157:     DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
158:     DMGetCoordinates(dm, &coordinates);
159:     VecScale(coordinates, 1.0 / lengthScale);
160:   }

162:   /* Read Labels */
163:   /* TODO: this probably does not work as elements get permuted */
164:   /* DMPlexLabelsLoad_HDF5_Internal(dm, viewer); */
165:   return 0;
166: }
167: #endif