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;
10: size_t len;
12: PetscFunctionBegin;
13: PetscCall(PetscStrrchr(path, '/', &tmp));
14: PetscCall(PetscStrlen(tmp, &len));
15: PetscCall(PetscStrncpy(name, tmp, len + 1)); /* assuming adequate buffer */
16: if (tmp != path) {
17: /* '/' found, name is substring of path after last occurrence of '/'. */
18: /* Trim the '/name' part from path just by inserting null character. */
19: tmp--;
20: *tmp = '\0';
21: } else {
22: /* '/' not found, name = path, path = "/". */
23: PetscCall(PetscStrncpy(path, "/", 2)); /* assuming adequate buffer */
24: }
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: /*
29: - invert (involute) cells of some types according to XDMF/VTK numbering of vertices in a cells
30: - cell type is identified using the number of vertices
31: */
32: static PetscErrorCode DMPlexInvertCells_XDMF_Private(DM dm)
33: {
34: PetscInt dim, *cones, cHeight, cStart, cEnd, p;
35: PetscSection cs;
37: PetscFunctionBegin;
38: PetscCall(DMGetDimension(dm, &dim));
39: if (dim != 3) PetscFunctionReturn(PETSC_SUCCESS);
40: PetscCall(DMPlexGetCones(dm, &cones));
41: PetscCall(DMPlexGetConeSection(dm, &cs));
42: PetscCall(DMPlexGetVTKCellHeight(dm, &cHeight));
43: PetscCall(DMPlexGetHeightStratum(dm, cHeight, &cStart, &cEnd));
44: for (p = cStart; p < cEnd; p++) {
45: PetscInt numCorners, o;
47: PetscCall(PetscSectionGetDof(cs, p, &numCorners));
48: PetscCall(PetscSectionGetOffset(cs, p, &o));
49: switch (numCorners) {
50: case 4:
51: PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cones[o]));
52: break;
53: case 6:
54: PetscCall(DMPlexInvertCell(DM_POLYTOPE_TRI_PRISM, &cones[o]));
55: break;
56: case 8:
57: PetscCall(DMPlexInvertCell(DM_POLYTOPE_HEXAHEDRON, &cones[o]));
58: break;
59: }
60: }
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM dm, PetscViewer viewer)
65: {
66: Vec coordinates;
67: IS cells;
68: PetscInt spatialDim, topoDim = -1, numCells, numVertices, NVertices, numCorners;
69: PetscMPIInt rank;
70: MPI_Comm comm;
71: char topo_path[PETSC_MAX_PATH_LEN] = "/viz/topology/cells", topo_name[PETSC_MAX_PATH_LEN];
72: char geom_path[PETSC_MAX_PATH_LEN] = "/geometry/vertices", geom_name[PETSC_MAX_PATH_LEN];
73: PetscBool seq = PETSC_FALSE;
75: PetscFunctionBegin;
76: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
77: PetscCallMPI(MPI_Comm_rank(comm, &rank));
79: PetscOptionsBegin(PetscObjectComm((PetscObject)dm), ((PetscObject)dm)->prefix, "DMPlex HDF5/XDMF Loader Options", "PetscViewer");
80: PetscCall(PetscOptionsString("-dm_plex_hdf5_topology_path", "HDF5 path of topology dataset", NULL, topo_path, topo_path, sizeof(topo_path), NULL));
81: PetscCall(PetscOptionsString("-dm_plex_hdf5_geometry_path", "HDF5 path to geometry dataset", NULL, geom_path, geom_path, sizeof(geom_path), NULL));
82: PetscCall(PetscOptionsBool("-dm_plex_hdf5_force_sequential", "force sequential loading", NULL, seq, &seq, NULL));
83: PetscOptionsEnd();
85: PetscCall(SplitPath_Private(topo_path, topo_name));
86: PetscCall(SplitPath_Private(geom_path, geom_name));
87: PetscCall(PetscInfo(dm, "Topology group %s, name %s\n", topo_path, topo_name));
88: PetscCall(PetscInfo(dm, "Geometry group %s, name %s\n", geom_path, geom_name));
90: /* Read topology */
91: PetscCall(PetscViewerHDF5PushGroup(viewer, topo_path));
92: PetscCall(ISCreate(comm, &cells));
93: PetscCall(PetscObjectSetName((PetscObject)cells, topo_name));
94: if (seq) {
95: PetscCall(PetscViewerHDF5ReadSizes(viewer, topo_name, NULL, &numCells));
96: PetscCall(PetscLayoutSetSize(cells->map, numCells));
97: numCells = rank == 0 ? numCells : 0;
98: PetscCall(PetscLayoutSetLocalSize(cells->map, numCells));
99: }
100: PetscCall(ISLoad(cells, viewer));
101: PetscCall(ISGetLocalSize(cells, &numCells));
102: PetscCall(ISGetBlockSize(cells, &numCorners));
103: PetscCall(PetscViewerHDF5ReadAttribute(viewer, topo_name, "cell_dim", PETSC_INT, &topoDim, &topoDim));
104: PetscCall(PetscViewerHDF5PopGroup(viewer));
105: numCells /= numCorners;
107: /* Read geometry */
108: PetscCall(PetscViewerHDF5PushGroup(viewer, geom_path));
109: PetscCall(VecCreate(comm, &coordinates));
110: PetscCall(PetscObjectSetName((PetscObject)coordinates, geom_name));
111: if (seq) {
112: PetscCall(PetscViewerHDF5ReadSizes(viewer, geom_name, NULL, &numVertices));
113: PetscCall(PetscLayoutSetSize(coordinates->map, numVertices));
114: numVertices = rank == 0 ? numVertices : 0;
115: PetscCall(PetscLayoutSetLocalSize(coordinates->map, numVertices));
116: }
117: PetscCall(VecLoad(coordinates, viewer));
118: PetscCall(VecGetLocalSize(coordinates, &numVertices));
119: PetscCall(VecGetSize(coordinates, &NVertices));
120: PetscCall(VecGetBlockSize(coordinates, &spatialDim));
121: PetscCall(PetscViewerHDF5PopGroup(viewer));
122: numVertices /= spatialDim;
123: NVertices /= spatialDim;
125: PetscCall(PetscInfo(NULL, "Loaded mesh dimensions: numCells %" PetscInt_FMT " numCorners %" PetscInt_FMT " numVertices %" PetscInt_FMT " spatialDim %" PetscInt_FMT "\n", numCells, numCorners, numVertices, spatialDim));
126: {
127: const PetscScalar *coordinates_arr;
128: PetscReal *coordinates_arr_real;
129: const PetscInt *cells_arr;
130: PetscSF sfVert = NULL;
131: PetscInt i;
133: PetscCall(VecGetArrayRead(coordinates, &coordinates_arr));
134: PetscCall(ISGetIndices(cells, &cells_arr));
136: if (PetscDefined(USE_COMPLEX)) {
137: /* convert to real numbers if PetscScalar is complex */
138: /*TODO More systematic would be to change all the function arguments to PetscScalar */
139: PetscCall(PetscMalloc1(numVertices * spatialDim, &coordinates_arr_real));
140: for (i = 0; i < numVertices * spatialDim; ++i) {
141: coordinates_arr_real[i] = PetscRealPart(coordinates_arr[i]);
142: 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.");
143: }
144: } else coordinates_arr_real = (PetscReal *)coordinates_arr;
146: PetscCall(DMSetDimension(dm, topoDim < 0 ? spatialDim : topoDim));
147: PetscCall(DMPlexBuildFromCellListParallel(dm, numCells, numVertices, NVertices, numCorners, cells_arr, &sfVert, NULL));
148: PetscCall(DMPlexInvertCells_XDMF_Private(dm));
149: PetscCall(DMPlexBuildCoordinatesFromCellListParallel(dm, spatialDim, sfVert, coordinates_arr_real));
150: PetscCall(VecRestoreArrayRead(coordinates, &coordinates_arr));
151: PetscCall(ISRestoreIndices(cells, &cells_arr));
152: PetscCall(PetscSFDestroy(&sfVert));
153: if (PetscDefined(USE_COMPLEX)) PetscCall(PetscFree(coordinates_arr_real));
154: }
155: PetscCall(ISDestroy(&cells));
156: PetscCall(VecDestroy(&coordinates));
158: /* scale coordinates - unlike in DMPlexLoad_HDF5_Internal, this can only be done after DM is populated */
159: {
160: PetscReal lengthScale;
162: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
163: PetscCall(DMGetCoordinates(dm, &coordinates));
164: PetscCall(VecScale(coordinates, 1.0 / lengthScale));
165: }
167: /* Read Labels */
168: /* TODO: this probably does not work as elements get permuted */
169: /* PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer)); */
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
172: #endif