Actual source code: plexhdf5xdmf.c
petsc-3.11.4 2019-09-28
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;
13: PetscStrrchr(path,'/',&tmp);
14: PetscStrcpy(name,tmp);
15: if (tmp != path) {
16: /* '/' found, name is substring of path after last occurence 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: PetscStrcpy(path,"/");
23: }
24: return(0);
25: }
27: PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM dm, PetscViewer viewer)
28: {
29: Vec coordinates;
30: IS cells;
31: PetscInt spatialDim, numCells, numVertices, numCorners;
32: PetscMPIInt rank;
33: MPI_Comm comm;
34: PetscErrorCode ierr;
35: char topo_path[PETSC_MAX_PATH_LEN]="/viz/topology/cells", topo_name[PETSC_MAX_PATH_LEN];
36: char geom_path[PETSC_MAX_PATH_LEN]="/geometry/vertices", geom_name[PETSC_MAX_PATH_LEN];
37: PetscBool seq = PETSC_FALSE;
40: PetscObjectGetComm((PetscObject)dm, &comm);
41: MPI_Comm_rank(comm, &rank);
43: PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"DMPlex HDF5/XDMF Loader Options","PetscViewer");
44: PetscOptionsString("-dm_plex_hdf5_topology_path","HDF5 path of topology dataset",NULL,topo_path,topo_path,sizeof(topo_path),NULL);
45: PetscOptionsString("-dm_plex_hdf5_geometry_path","HDF5 path to geometry dataset",NULL,geom_path,geom_path,sizeof(geom_path),NULL);
46: PetscOptionsBool( "-dm_plex_hdf5_force_sequential","force sequential loading",NULL,seq,&seq,NULL);
47: PetscOptionsEnd();
49: SplitPath_Private(topo_path, topo_name);
50: SplitPath_Private(geom_path, geom_name);
51: PetscInfo2(dm, "Topology group %s, name %s\n", topo_path, topo_name);
52: PetscInfo2(dm, "Geometry group %s, name %s\n", geom_path, geom_name);
54: /* Read topology */
55: PetscViewerHDF5PushGroup(viewer, topo_path);
56: ISCreate(comm, &cells);
57: PetscObjectSetName((PetscObject) cells, topo_name);
58: if (seq) {
59: PetscViewerHDF5ReadSizes(viewer, topo_name, NULL, &numCells);
60: PetscLayoutSetSize(cells->map, numCells);
61: numCells = !rank ? numCells : 0;
62: PetscLayoutSetLocalSize(cells->map, numCells);
63: }
64: ISLoad(cells, viewer);
65: ISGetLocalSize(cells, &numCells);
66: ISGetBlockSize(cells, &numCorners);
67: PetscViewerHDF5PopGroup(viewer);
68: numCells /= numCorners;
70: /* Read geometry */
71: PetscViewerHDF5PushGroup(viewer, geom_path);
72: VecCreate(comm, &coordinates);
73: PetscObjectSetName((PetscObject) coordinates, geom_name);
74: if (seq) {
75: PetscViewerHDF5ReadSizes(viewer, geom_name, NULL, &numVertices);
76: PetscLayoutSetSize(coordinates->map, numVertices);
77: numVertices = !rank ? numVertices : 0;
78: PetscLayoutSetLocalSize(coordinates->map, numVertices);
79: }
80: VecLoad(coordinates, viewer);
81: VecGetLocalSize(coordinates, &numVertices);
82: VecGetBlockSize(coordinates, &spatialDim);
83: PetscViewerHDF5PopGroup(viewer);
84: numVertices /= spatialDim;
86: PetscInfo4(NULL, "Loaded mesh dimensions: numCells %d numCorners %d numVertices %d spatialDim %d\n", numCells, numCorners, numVertices, spatialDim);
88: #if defined(PETSC_USE_DEBUG)
89: /* Check that maximum index referred in cells is in line with global number of vertices */
90: {
91: PetscInt max1, max2;
92: ISGetMinMax(cells, NULL, &max1);
93: VecGetSize(coordinates, &max2);
94: max2 /= spatialDim; max2--;
95: if (max1 > max2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "maximum index in cells = %d > %d = total number of vertices - 1", max1, max2);
96: }
97: #endif
99: {
100: const PetscScalar *coordinates_arr;
101: PetscReal *coordinates_arr_real;
102: const PetscInt *cells_arr;
103: int *cells_arr_int;
104: PetscSF sfVert=NULL;
105: #if defined(PETSC_USE_64BIT_INDICES) || defined(PETSC_USE_COMPLEX)
106: PetscInt i;
107: #endif
109: VecGetArrayRead(coordinates, &coordinates_arr);
110: ISGetIndices(cells, &cells_arr);
112: #if defined(PETSC_USE_64BIT_INDICES)
113: /* convert to 32-bit integers if PetscInt is 64-bit */
114: /*TODO More systematic would be to change all the function arguments to PetscInt */
115: PetscMalloc1(numCells*numCorners, &cells_arr_int);
116: for (i = 0; i < numCells*numCorners; ++i) {
117: PetscMPIIntCast(cells_arr[i], &cells_arr_int[i]);
118: }
119: #else
120: cells_arr_int = (int*)cells_arr;
121: #endif
123: #if defined(PETSC_USE_COMPLEX)
124: /* convert to real numbers if PetscScalar is complex */
125: /*TODO More systematic would be to change all the function arguments to PetscScalar */
126: PetscMalloc1(numVertices*spatialDim, &coordinates_arr_real);
127: for (i = 0; i < numVertices*spatialDim; ++i) {
128: coordinates_arr_real[i] = PetscRealPart(coordinates_arr[i]);
129: #if defined(PETSC_USE_DEBUG)
130: if (PetscImaginaryPart(coordinates_arr[i])) {
131: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector of coordinates contains complex numbers but only real vectors are currently supported.");
132: }
133: #endif /* defined(PETSC_USE_DEBUG) */
134: }
135: #else
136: coordinates_arr_real = (PetscReal*)coordinates_arr;
137: #endif
139: DMSetDimension(dm, spatialDim);
140: DMPlexBuildFromCellList_Parallel_Internal(dm, spatialDim, numCells, numVertices, numCorners, cells_arr_int, PETSC_TRUE, &sfVert);
141: DMPlexBuildCoordinates_Parallel_Internal( dm, spatialDim, numCells, numVertices, sfVert, coordinates_arr_real);
142: VecRestoreArrayRead(coordinates, &coordinates_arr);
143: ISRestoreIndices(cells, &cells_arr);
144: PetscSFDestroy(&sfVert);
145: #if defined(PETSC_USE_64BIT_INDICES)
146: PetscFree(cells_arr_int);
147: #endif
148: #if defined(PETSC_USE_COMPLEX)
149: PetscFree(coordinates_arr_real);
150: #endif
151: }
152: ISDestroy(&cells);
153: VecDestroy(&coordinates);
155: /* scale coordinates - unlike in DMPlexLoad_HDF5_Internal, this can only be done after DM is populated */
156: {
157: PetscReal lengthScale;
159: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
160: DMGetCoordinates(dm, &coordinates);
161: VecScale(coordinates, 1.0/lengthScale);
162: }
164: /* Read Labels */
165: /* TODO: this probably does not work as elements get permuted */
166: /* DMPlexLoadLabels_HDF5_Internal(dm, viewer); */
167: return(0);
168: }
169: #endif