1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/isimpl.h> 3 #include <petsc/private/vecimpl.h> 4 #include <petsclayouthdf5.h> 5 6 static PetscErrorCode SplitPath_Private(char path[], char name[]) 7 { 8 char *tmp; 9 size_t len; 10 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 } 26 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; 35 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; 45 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 } 62 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; 73 74 PetscFunctionBegin; 75 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 76 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 77 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(); 83 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)); 88 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; 105 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; 123 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; 131 132 PetscCall(VecGetArrayRead(coordinates, &coordinates_arr)); 133 PetscCall(ISGetIndices(cells, &cells_arr)); 134 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; 144 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)); 156 157 /* scale coordinates - unlike in DMPlexLoad_HDF5_Internal, this can only be done after DM is populated */ 158 { 159 PetscReal lengthScale; 160 161 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale)); 162 PetscCall(DMGetCoordinates(dm, &coordinates)); 163 PetscCall(VecScale(coordinates, 1.0 / lengthScale)); 164 } 165 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 } 171