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
SplitPath_Private(char path[],char name[])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 */
DMPlexInvertCells_XDMF_Private(DM dm)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
DMPlexLoad_HDF5_Xdmf_Internal(DM dm,PetscViewer viewer)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