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: static PetscErrorCode SplitPath_Private(char path[], char name[])
      7: {
      8:   char  *tmp;
      9:   size_t len;

     11:   PetscFunctionBegin;
     12:   PetscCall(PetscStrrchr(path, '/', &tmp));
     13:   PetscCall(PetscStrlen(tmp, &len));
     14:   PetscCall(PetscStrncpy(name, tmp, len + 1)); /* assuming adequate buffer */
     15:   if (tmp != path) {
     16:     /* '/' found, name is substring of path after last occurrence of '/'. */
     17:     /* Trim the '/name' part from path just by inserting null character. */
     18:     tmp--;
     19:     *tmp = '\0';
     20:   } else {
     21:     /* '/' not found, name = path, path = "/". */
     22:     PetscCall(PetscStrncpy(path, "/", 2)); /* assuming adequate buffer */
     23:   }
     24:   PetscFunctionReturn(PETSC_SUCCESS);
     25: }

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

     36:   PetscFunctionBegin;
     37:   PetscCall(DMGetDimension(dm, &dim));
     38:   if (dim != 3) PetscFunctionReturn(PETSC_SUCCESS);
     39:   PetscCall(DMPlexGetCones(dm, &cones));
     40:   PetscCall(DMPlexGetConeSection(dm, &cs));
     41:   PetscCall(DMPlexGetVTKCellHeight(dm, &cHeight));
     42:   PetscCall(DMPlexGetHeightStratum(dm, cHeight, &cStart, &cEnd));
     43:   for (p = cStart; p < cEnd; p++) {
     44:     PetscInt numCorners, o;

     46:     PetscCall(PetscSectionGetDof(cs, p, &numCorners));
     47:     PetscCall(PetscSectionGetOffset(cs, p, &o));
     48:     switch (numCorners) {
     49:     case 4:
     50:       PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cones[o]));
     51:       break;
     52:     case 6:
     53:       PetscCall(DMPlexInvertCell(DM_POLYTOPE_TRI_PRISM, &cones[o]));
     54:       break;
     55:     case 8:
     56:       PetscCall(DMPlexInvertCell(DM_POLYTOPE_HEXAHEDRON, &cones[o]));
     57:       break;
     58:     }
     59:   }
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

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

     74:   PetscFunctionBegin;
     75:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
     76:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

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

     84:   PetscCall(SplitPath_Private(topo_path, topo_name));
     85:   PetscCall(SplitPath_Private(geom_path, geom_name));
     86:   PetscCall(PetscInfo(dm, "Topology group %s, name %s\n", topo_path, topo_name));
     87:   PetscCall(PetscInfo(dm, "Geometry group %s, name %s\n", geom_path, geom_name));

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

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

    124:   PetscCall(PetscInfo(NULL, "Loaded mesh dimensions: numCells %" PetscInt_FMT " numCorners %" PetscInt_FMT " numVertices %" PetscInt_FMT " spatialDim %" PetscInt_FMT "\n", numCells, numCorners, numVertices, spatialDim));
    125:   {
    126:     const PetscScalar *coordinates_arr;
    127:     PetscReal         *coordinates_arr_real;
    128:     const PetscInt    *cells_arr;
    129:     PetscSF            sfVert = NULL;
    130:     PetscInt           i;

    132:     PetscCall(VecGetArrayRead(coordinates, &coordinates_arr));
    133:     PetscCall(ISGetIndices(cells, &cells_arr));

    135:     if (PetscDefined(USE_COMPLEX)) {
    136:       /* convert to real numbers if PetscScalar is complex */
    137:       /*TODO More systematic would be to change all the function arguments to PetscScalar */
    138:       PetscCall(PetscMalloc1(numVertices * spatialDim, &coordinates_arr_real));
    139:       for (i = 0; i < numVertices * spatialDim; ++i) {
    140:         coordinates_arr_real[i] = PetscRealPart(coordinates_arr[i]);
    141:         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.");
    142:       }
    143:     } else coordinates_arr_real = (PetscReal *)coordinates_arr;

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

    157:   /* scale coordinates - unlike in DMPlexLoad_HDF5_Internal, this can only be done after DM is populated */
    158:   {
    159:     PetscReal lengthScale;

    161:     PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
    162:     PetscCall(DMGetCoordinates(dm, &coordinates));
    163:     PetscCall(VecScale(coordinates, 1.0 / lengthScale));
    164:   }

    166:   /* Read Labels */
    167:   /* TODO: this probably does not work as elements get permuted */
    168:   /* PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer)); */
    169:   PetscFunctionReturn(PETSC_SUCCESS);
    170: }