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