Actual source code: plexhdf5.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/isimpl.h>
3: #include <petsc/private/vecimpl.h>
4: #include <petsc/private/viewerhdf5impl.h>
5: #include <petsclayouthdf5.h>
7: #if defined(PETSC_HAVE_HDF5)
9: typedef struct DMPlexStorageVersion {
10: PetscInt major, minor, subminor;
11: } DMPlexStorageVersion;
13: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
15: static PetscErrorCode DMPlexStorageVersionParseString_Private(DM dm, const char str[], DMPlexStorageVersion *v)
16: {
17: PetscToken t;
18: char *ts;
19: PetscInt i;
20: PetscInt ti[3];
22: PetscTokenCreate(str, '.', &t);
23: for (i=0; i<3; i++) {
24: PetscTokenFind(t, &ts);
26: PetscOptionsStringToInt(ts, &ti[i]);
27: }
28: PetscTokenFind(t, &ts);
30: PetscTokenDestroy(&t);
31: v->major = ti[0];
32: v->minor = ti[1];
33: v->subminor = ti[2];
34: return 0;
35: }
37: static PetscErrorCode DMPlexStorageVersionSetUpWriting_Private(DM dm, PetscViewer viewer, DMPlexStorageVersion *version)
38: {
39: const char ATTR_NAME[] = "dmplex_storage_version";
40: PetscBool fileHasVersion;
41: char optVersion[16], fileVersion[16];
42: PetscErrorCode ierr;
44: PetscStrcpy(fileVersion, DMPLEX_STORAGE_VERSION_STABLE);
45: PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion);
46: if (fileHasVersion) {
47: char *tmp;
49: PetscViewerHDF5ReadAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, NULL, &tmp);
50: PetscStrcpy(fileVersion, tmp);
51: PetscFree(tmp);
52: }
53: PetscStrcpy(optVersion, fileVersion);
54: PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"DMPlex HDF5 Viewer Options","PetscViewer");
55: PetscOptionsString("-dm_plex_view_hdf5_storage_version","DMPlex HDF5 viewer storage version",NULL,optVersion,optVersion,sizeof(optVersion),NULL);
56: PetscOptionsEnd();
57: if (!fileHasVersion) {
58: PetscViewerHDF5WriteAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, optVersion);
59: } else {
60: PetscBool flg;
62: PetscStrcmp(fileVersion, optVersion, &flg);
64: }
65: PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT);
66: DMPlexStorageVersionParseString_Private(dm, optVersion, version);
67: return 0;
68: }
70: static PetscErrorCode DMPlexStorageVersionGet_Private(DM dm, PetscViewer viewer, DMPlexStorageVersion *version)
71: {
72: const char ATTR_NAME[] = "dmplex_storage_version";
73: char *defaultVersion;
74: char *versionString;
76: //TODO string HDF5 attribute handling is terrible and should be redesigned
77: PetscStrallocpy("1.0.0", &defaultVersion);
78: PetscViewerHDF5ReadAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, &defaultVersion, &versionString);
79: DMPlexStorageVersionParseString_Private(dm, versionString, version);
80: PetscFree(versionString);
81: PetscFree(defaultVersion);
82: return 0;
83: }
85: static PetscErrorCode DMPlexGetHDF5Name_Private(DM dm, const char *name[])
86: {
87: if (((PetscObject)dm)->name) {
88: PetscObjectGetName((PetscObject)dm, name);
89: } else {
90: *name = "plex";
91: }
92: return 0;
93: }
95: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
96: {
97: Vec stamp;
98: PetscMPIInt rank;
100: if (seqnum < 0) return 0;
101: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
102: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
103: VecSetBlockSize(stamp, 1);
104: PetscObjectSetName((PetscObject) stamp, seqname);
105: if (rank == 0) {
106: PetscReal timeScale;
107: PetscBool istime;
109: PetscStrncmp(seqname, "time", 5, &istime);
110: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); value *= timeScale;}
111: VecSetValue(stamp, 0, value, INSERT_VALUES);
112: }
113: VecAssemblyBegin(stamp);
114: VecAssemblyEnd(stamp);
115: PetscViewerHDF5PushGroup(viewer, "/");
116: PetscViewerHDF5PushTimestepping(viewer);
117: PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
118: VecView(stamp, viewer);
119: PetscViewerHDF5PopTimestepping(viewer);
120: PetscViewerHDF5PopGroup(viewer);
121: VecDestroy(&stamp);
122: return 0;
123: }
125: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
126: {
127: Vec stamp;
128: PetscMPIInt rank;
130: if (seqnum < 0) return 0;
131: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
132: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
133: VecSetBlockSize(stamp, 1);
134: PetscObjectSetName((PetscObject) stamp, seqname);
135: PetscViewerHDF5PushGroup(viewer, "/");
136: PetscViewerHDF5PushTimestepping(viewer);
137: PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
138: VecLoad(stamp, viewer);
139: PetscViewerHDF5PopTimestepping(viewer);
140: PetscViewerHDF5PopGroup(viewer);
141: if (rank == 0) {
142: const PetscScalar *a;
143: PetscReal timeScale;
144: PetscBool istime;
146: VecGetArrayRead(stamp, &a);
147: *value = a[0];
148: VecRestoreArrayRead(stamp, &a);
149: PetscStrncmp(seqname, "time", 5, &istime);
150: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); *value /= timeScale;}
151: }
152: VecDestroy(&stamp);
153: return 0;
154: }
156: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
157: {
158: IS cutcells = NULL;
159: const PetscInt *cutc;
160: PetscInt cellHeight, vStart, vEnd, cStart, cEnd, c;
162: if (!cutLabel) return 0;
163: DMPlexGetVTKCellHeight(dm, &cellHeight);
164: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
165: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
166: /* Label vertices that should be duplicated */
167: DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel);
168: DMLabelGetStratumIS(cutLabel, 2, &cutcells);
169: if (cutcells) {
170: PetscInt n;
172: ISGetIndices(cutcells, &cutc);
173: ISGetLocalSize(cutcells, &n);
174: for (c = 0; c < n; ++c) {
175: if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
176: PetscInt *closure = NULL;
177: PetscInt closureSize, cl, value;
179: DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
180: for (cl = 0; cl < closureSize*2; cl += 2) {
181: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
182: DMLabelGetValue(cutLabel, closure[cl], &value);
183: if (value == 1) {
184: DMLabelSetValue(*cutVertexLabel, closure[cl], 1);
185: }
186: }
187: }
188: DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
189: }
190: }
191: ISRestoreIndices(cutcells, &cutc);
192: }
193: ISDestroy(&cutcells);
194: return 0;
195: }
197: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
198: {
199: DM dm;
200: DM dmBC;
201: PetscSection section, sectionGlobal;
202: Vec gv;
203: const char *name;
204: PetscViewerVTKFieldType ft;
205: PetscViewerFormat format;
206: PetscInt seqnum;
207: PetscReal seqval;
208: PetscBool isseq;
210: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
211: VecGetDM(v, &dm);
212: DMGetLocalSection(dm, §ion);
213: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
214: DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
215: if (seqnum >= 0) {
216: PetscViewerHDF5PushTimestepping(viewer);
217: PetscViewerHDF5SetTimestep(viewer, seqnum);
218: }
219: PetscViewerGetFormat(viewer, &format);
220: DMGetOutputDM(dm, &dmBC);
221: DMGetGlobalSection(dmBC, §ionGlobal);
222: DMGetGlobalVector(dmBC, &gv);
223: PetscObjectGetName((PetscObject) v, &name);
224: PetscObjectSetName((PetscObject) gv, name);
225: DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
226: DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
227: PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
228: if (isseq) VecView_Seq(gv, viewer);
229: else VecView_MPI(gv, viewer);
230: if (format == PETSC_VIEWER_HDF5_VIZ) {
231: /* Output visualization representation */
232: PetscInt numFields, f;
233: DMLabel cutLabel, cutVertexLabel = NULL;
235: PetscSectionGetNumFields(section, &numFields);
236: DMGetLabel(dm, "periodic_cut", &cutLabel);
237: for (f = 0; f < numFields; ++f) {
238: Vec subv;
239: IS is;
240: const char *fname, *fgroup, *componentName;
241: char subname[PETSC_MAX_PATH_LEN];
242: PetscInt pStart, pEnd, Nc, c;
244: DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
245: fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
246: PetscSectionGetFieldName(section, f, &fname);
247: if (!fname) continue;
248: PetscViewerHDF5PushGroup(viewer, fgroup);
249: if (cutLabel) {
250: const PetscScalar *ga;
251: PetscScalar *suba;
252: PetscInt gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;
254: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
255: PetscSectionGetFieldComponents(section, f, &Nc);
256: for (p = pStart; p < pEnd; ++p) {
257: PetscInt gdof, fdof = 0, val;
259: PetscSectionGetDof(sectionGlobal, p, &gdof);
260: if (gdof > 0) PetscSectionGetFieldDof(section, p, f, &fdof);
261: subSize += fdof;
262: DMLabelGetValue(cutVertexLabel, p, &val);
263: if (val == 1) extSize += fdof;
264: }
265: VecCreate(PetscObjectComm((PetscObject) gv), &subv);
266: VecSetSizes(subv, subSize+extSize, PETSC_DETERMINE);
267: VecSetBlockSize(subv, Nc);
268: VecSetType(subv, VECSTANDARD);
269: VecGetOwnershipRange(gv, &gstart, NULL);
270: VecGetArrayRead(gv, &ga);
271: VecGetArray(subv, &suba);
272: for (p = pStart; p < pEnd; ++p) {
273: PetscInt gdof, goff, val;
275: PetscSectionGetDof(sectionGlobal, p, &gdof);
276: if (gdof > 0) {
277: PetscInt fdof, fc, f2, poff = 0;
279: PetscSectionGetOffset(sectionGlobal, p, &goff);
280: /* Can get rid of this loop by storing field information in the global section */
281: for (f2 = 0; f2 < f; ++f2) {
282: PetscSectionGetFieldDof(section, p, f2, &fdof);
283: poff += fdof;
284: }
285: PetscSectionGetFieldDof(section, p, f, &fdof);
286: for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff+poff+fc - gstart];
287: DMLabelGetValue(cutVertexLabel, p, &val);
288: if (val == 1) {
289: for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize+newOff] = ga[goff+poff+fc - gstart];
290: }
291: }
292: }
293: VecRestoreArrayRead(gv, &ga);
294: VecRestoreArray(subv, &suba);
295: DMLabelDestroy(&cutVertexLabel);
296: } else {
297: PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
298: }
299: PetscStrncpy(subname, name,sizeof(subname));
300: PetscStrlcat(subname, "_",sizeof(subname));
301: PetscStrlcat(subname, fname,sizeof(subname));
302: PetscObjectSetName((PetscObject) subv, subname);
303: if (isseq) VecView_Seq(subv, viewer);
304: else VecView_MPI(subv, viewer);
305: if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
306: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "vector");
307: } else {
308: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "scalar");
309: }
311: /* Output the component names in the field if available */
312: PetscSectionGetFieldComponents(section, f, &Nc);
313: for (c = 0; c < Nc; ++c){
314: char componentNameLabel[PETSC_MAX_PATH_LEN];
315: PetscSectionGetComponentName(section, f, c, &componentName);
316: PetscSNPrintf(componentNameLabel, sizeof(componentNameLabel), "componentName%D", c);
317: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, componentNameLabel, PETSC_STRING, componentName);
318: }
320: if (cutLabel) VecDestroy(&subv);
321: else PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
322: PetscViewerHDF5PopGroup(viewer);
323: }
324: }
325: if (seqnum >= 0) {
326: PetscViewerHDF5PopTimestepping(viewer);
327: }
328: DMRestoreGlobalVector(dmBC, &gv);
329: return 0;
330: }
332: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
333: {
334: DM dm;
335: Vec locv;
336: PetscObject isZero;
337: const char *name;
338: PetscReal time;
340: VecGetDM(v, &dm);
341: DMGetLocalVector(dm, &locv);
342: PetscObjectGetName((PetscObject) v, &name);
343: PetscObjectSetName((PetscObject) locv, name);
344: PetscObjectQuery((PetscObject) v, "__Vec_bc_zero__", &isZero);
345: PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", isZero);
346: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
347: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
348: DMGetOutputSequenceNumber(dm, NULL, &time);
349: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
350: PetscViewerHDF5PushGroup(viewer, "/fields");
351: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
352: VecView_Plex_Local_HDF5_Internal(locv, viewer);
353: PetscViewerPopFormat(viewer);
354: PetscViewerHDF5PopGroup(viewer);
355: PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", NULL);
356: DMRestoreLocalVector(dm, &locv);
357: return 0;
358: }
360: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
361: {
362: PetscBool isseq;
364: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
365: PetscViewerHDF5PushGroup(viewer, "/fields");
366: if (isseq) VecView_Seq(v, viewer);
367: else VecView_MPI(v, viewer);
368: PetscViewerHDF5PopGroup(viewer);
369: return 0;
370: }
372: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
373: {
374: DM dm;
375: Vec locv;
376: const char *name;
377: PetscInt seqnum;
379: VecGetDM(v, &dm);
380: DMGetLocalVector(dm, &locv);
381: PetscObjectGetName((PetscObject) v, &name);
382: PetscObjectSetName((PetscObject) locv, name);
383: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
384: PetscViewerHDF5PushGroup(viewer, "/fields");
385: if (seqnum >= 0) {
386: PetscViewerHDF5PushTimestepping(viewer);
387: PetscViewerHDF5SetTimestep(viewer, seqnum);
388: }
389: VecLoad_Plex_Local(locv, viewer);
390: if (seqnum >= 0) {
391: PetscViewerHDF5PopTimestepping(viewer);
392: }
393: PetscViewerHDF5PopGroup(viewer);
394: DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
395: DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
396: DMRestoreLocalVector(dm, &locv);
397: return 0;
398: }
400: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
401: {
402: DM dm;
403: PetscInt seqnum;
405: VecGetDM(v, &dm);
406: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
407: PetscViewerHDF5PushGroup(viewer, "/fields");
408: if (seqnum >= 0) {
409: PetscViewerHDF5PushTimestepping(viewer);
410: PetscViewerHDF5SetTimestep(viewer, seqnum);
411: }
412: VecLoad_Default(v, viewer);
413: if (seqnum >= 0) {
414: PetscViewerHDF5PopTimestepping(viewer);
415: }
416: PetscViewerHDF5PopGroup(viewer);
417: return 0;
418: }
420: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
421: {
422: const char *topologydm_name;
423: const char *pointsName, *coneSizesName, *conesName, *orientationsName;
424: IS pointsIS, coneSizesIS, conesIS, orientationsIS;
425: PetscInt *points, *coneSizes, *cones, *orientations;
426: const PetscInt *gpoint;
427: PetscInt pStart, pEnd, nPoints = 0, conesSize = 0;
428: PetscInt p, c, s;
429: DMPlexStorageVersion version;
430: char group[PETSC_MAX_PATH_LEN];
431: MPI_Comm comm;
433: PetscObjectGetComm((PetscObject) dm, &comm);
434: pointsName = "order";
435: coneSizesName = "cones";
436: conesName = "cells";
437: orientationsName = "orientation";
438: DMPlexStorageVersionSetUpWriting_Private(dm, viewer, &version);
439: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
440: if (version.major <= 1) {
441: PetscStrcpy(group, "/topology");
442: } else {
443: PetscSNPrintf(group, sizeof(group), "topologies/%s/topology", topologydm_name);
444: }
445: PetscViewerHDF5PushGroup(viewer, group);
446: DMPlexGetChart(dm, &pStart, &pEnd);
447: ISGetIndices(globalPointNumbers, &gpoint);
448: for (p = pStart; p < pEnd; ++p) {
449: if (gpoint[p] >= 0) {
450: PetscInt coneSize;
452: DMPlexGetConeSize(dm, p, &coneSize);
453: nPoints += 1;
454: conesSize += coneSize;
455: }
456: }
457: PetscMalloc1(nPoints, &points);
458: PetscMalloc1(nPoints, &coneSizes);
459: PetscMalloc1(conesSize, &cones);
460: PetscMalloc1(conesSize, &orientations);
461: for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
462: if (gpoint[p] >= 0) {
463: const PetscInt *cone, *ornt;
464: PetscInt coneSize, cp;
466: DMPlexGetConeSize(dm, p, &coneSize);
467: DMPlexGetCone(dm, p, &cone);
468: DMPlexGetConeOrientation(dm, p, &ornt);
469: points[s] = gpoint[p];
470: coneSizes[s] = coneSize;
471: for (cp = 0; cp < coneSize; ++cp, ++c) {
472: cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]];
473: orientations[c] = ornt[cp];
474: }
475: ++s;
476: }
477: }
480: ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS);
481: ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS);
482: ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS);
483: ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS);
484: PetscObjectSetName((PetscObject) pointsIS, pointsName);
485: PetscObjectSetName((PetscObject) coneSizesIS, coneSizesName);
486: PetscObjectSetName((PetscObject) conesIS, conesName);
487: PetscObjectSetName((PetscObject) orientationsIS, orientationsName);
488: ISView(pointsIS, viewer);
489: ISView(coneSizesIS, viewer);
490: ISView(conesIS, viewer);
491: ISView(orientationsIS, viewer);
492: ISDestroy(&pointsIS);
493: ISDestroy(&coneSizesIS);
494: ISDestroy(&conesIS);
495: ISDestroy(&orientationsIS);
496: ISRestoreIndices(globalPointNumbers, &gpoint);
497: {
498: PetscInt dim;
500: DMGetDimension(dm, &dim);
501: PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *) &dim);
502: }
503: PetscViewerHDF5PopGroup(viewer);
504: return 0;
505: }
507: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
508: {
509: PetscSF sfPoint;
510: DMLabel cutLabel, cutVertexLabel = NULL;
511: IS globalVertexNumbers, cutvertices = NULL;
512: const PetscInt *gcell, *gvertex, *cutverts = NULL;
513: PetscInt *vertices;
514: PetscInt conesSize = 0;
515: PetscInt dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
517: *numCorners = 0;
518: DMGetDimension(dm, &dim);
519: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
520: ISGetIndices(globalCellNumbers, &gcell);
522: for (cell = cStart; cell < cEnd; ++cell) {
523: PetscInt *closure = NULL;
524: PetscInt closureSize, v, Nc = 0;
526: if (gcell[cell] < 0) continue;
527: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
528: for (v = 0; v < closureSize*2; v += 2) {
529: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
530: }
531: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
532: conesSize += Nc;
533: if (!numCornersLocal) numCornersLocal = Nc;
534: else if (numCornersLocal != Nc) numCornersLocal = 1;
535: }
536: MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
538: /* Handle periodic cuts by identifying vertices which should be duplicated */
539: DMGetLabel(dm, "periodic_cut", &cutLabel);
540: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
541: if (cutVertexLabel) DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);
542: if (cutvertices) {
543: ISGetIndices(cutvertices, &cutverts);
544: ISGetLocalSize(cutvertices, &vExtra);
545: }
546: DMGetPointSF(dm, &sfPoint);
547: if (cutLabel) {
548: const PetscInt *ilocal;
549: const PetscSFNode *iremote;
550: PetscInt nroots, nleaves;
552: PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
553: if (nleaves < 0) {
554: PetscObjectReference((PetscObject) sfPoint);
555: } else {
556: PetscSFCreate(PetscObjectComm((PetscObject) sfPoint), &sfPoint);
557: PetscSFSetGraph(sfPoint, nroots+vExtra, nleaves, (PetscInt*)ilocal, PETSC_COPY_VALUES, (PetscSFNode*)iremote, PETSC_COPY_VALUES);
558: }
559: } else {
560: PetscObjectReference((PetscObject) sfPoint);
561: }
562: /* Number all vertices */
563: DMPlexCreateNumbering_Plex(dm, vStart, vEnd+vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
564: PetscSFDestroy(&sfPoint);
565: /* Create cones */
566: ISGetIndices(globalVertexNumbers, &gvertex);
567: PetscMalloc1(conesSize, &vertices);
568: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
569: PetscInt *closure = NULL;
570: PetscInt closureSize, Nc = 0, p, value = -1;
571: PetscBool replace;
573: if (gcell[cell] < 0) continue;
574: if (cutLabel) DMLabelGetValue(cutLabel, cell, &value);
575: replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
576: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
577: for (p = 0; p < closureSize*2; p += 2) {
578: if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
579: closure[Nc++] = closure[p];
580: }
581: }
582: DMPlexReorderCell(dm, cell, closure);
583: for (p = 0; p < Nc; ++p) {
584: PetscInt nv, gv = gvertex[closure[p] - vStart];
586: if (replace) {
587: PetscFindInt(closure[p], vExtra, cutverts, &nv);
588: if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
589: }
590: vertices[v++] = gv < 0 ? -(gv+1) : gv;
591: }
592: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
593: }
594: ISRestoreIndices(globalVertexNumbers, &gvertex);
595: ISDestroy(&globalVertexNumbers);
596: ISRestoreIndices(globalCellNumbers, &gcell);
597: if (cutvertices) ISRestoreIndices(cutvertices, &cutverts);
598: ISDestroy(&cutvertices);
599: DMLabelDestroy(&cutVertexLabel);
601: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS);
602: PetscLayoutSetBlockSize((*cellIS)->map, *numCorners);
603: PetscObjectSetName((PetscObject) *cellIS, "cells");
604: return 0;
605: }
607: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, IS globalCellNumbers, PetscViewer viewer)
608: {
609: DM cdm;
610: DMLabel depthLabel, ctLabel;
611: IS cellIS;
612: PetscInt dim, depth, cellHeight, c;
613: hid_t fileId, groupId;
615: PetscViewerHDF5PushGroup(viewer, "/viz");
616: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
617: PetscStackCallHDF5(H5Gclose,(groupId));
619: PetscViewerHDF5PopGroup(viewer);
620: DMGetDimension(dm, &dim);
621: DMPlexGetDepth(dm, &depth);
622: DMGetCoordinateDM(dm, &cdm);
623: DMPlexGetVTKCellHeight(dm, &cellHeight);
624: DMPlexGetDepthLabel(dm, &depthLabel);
625: DMPlexGetCellTypeLabel(dm, &ctLabel);
626: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
627: const DMPolytopeType ict = (DMPolytopeType) c;
628: PetscInt pStart, pEnd, dep, numCorners, n = 0;
629: PetscBool output = PETSC_FALSE, doOutput;
631: if (ict == DM_POLYTOPE_FV_GHOST) continue;
632: DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
633: if (pStart >= 0) {
634: DMLabelGetValue(depthLabel, pStart, &dep);
635: if (dep == depth - cellHeight) output = PETSC_TRUE;
636: }
637: MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
638: if (!doOutput) continue;
639: CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS);
640: if (!n) {
641: PetscViewerHDF5PushGroup(viewer, "/viz/topology");
642: } else {
643: char group[PETSC_MAX_PATH_LEN];
645: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%D", n);
646: PetscViewerHDF5PushGroup(viewer, group);
647: }
648: ISView(cellIS, viewer);
649: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_corners", PETSC_INT, (void *) &numCorners);
650: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_dim", PETSC_INT, (void *) &dim);
651: ISDestroy(&cellIS);
652: PetscViewerHDF5PopGroup(viewer);
653: ++n;
654: }
655: return 0;
656: }
658: static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
659: {
660: DM cdm;
661: Vec coordinates, newcoords;
662: PetscReal lengthScale;
663: PetscInt m, M, bs;
665: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
666: DMGetCoordinateDM(dm, &cdm);
667: DMGetCoordinates(dm, &coordinates);
668: VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
669: PetscObjectSetName((PetscObject) newcoords, "vertices");
670: VecGetSize(coordinates, &M);
671: VecGetLocalSize(coordinates, &m);
672: VecSetSizes(newcoords, m, M);
673: VecGetBlockSize(coordinates, &bs);
674: VecSetBlockSize(newcoords, bs);
675: VecSetType(newcoords,VECSTANDARD);
676: VecCopy(coordinates, newcoords);
677: VecScale(newcoords, lengthScale);
678: /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
679: PetscViewerHDF5PushGroup(viewer, "/geometry");
680: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
681: VecView(newcoords, viewer);
682: PetscViewerPopFormat(viewer);
683: PetscViewerHDF5PopGroup(viewer);
684: VecDestroy(&newcoords);
685: return 0;
686: }
688: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
689: {
690: DM cdm;
691: Vec coords, newcoords;
692: PetscInt m, M, bs;
693: PetscReal lengthScale;
694: const char *topologydm_name, *coordinatedm_name, *coordinates_name;
696: {
697: PetscViewerFormat format;
698: DMPlexStorageVersion version;
700: PetscViewerGetFormat(viewer, &format);
701: DMPlexStorageVersionSetUpWriting_Private(dm, viewer, &version);
702: if (format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ || version.major <= 1) {
703: DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer);
704: return 0;
705: }
706: }
707: DMGetCoordinateDM(dm, &cdm);
708: DMGetCoordinates(dm, &coords);
709: PetscObjectGetName((PetscObject)cdm, &coordinatedm_name);
710: PetscObjectGetName((PetscObject)coords, &coordinates_name);
711: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
712: PetscViewerHDF5PushGroup(viewer, "topologies");
713: PetscViewerHDF5PushGroup(viewer, topologydm_name);
714: PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name);
715: PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name);
716: PetscViewerHDF5PopGroup(viewer);
717: PetscViewerHDF5PopGroup(viewer);
718: DMPlexSectionView(dm, viewer, cdm);
719: VecCreate(PetscObjectComm((PetscObject)coords), &newcoords);
720: PetscObjectSetName((PetscObject)newcoords, coordinates_name);
721: VecGetSize(coords, &M);
722: VecGetLocalSize(coords, &m);
723: VecSetSizes(newcoords, m, M);
724: VecGetBlockSize(coords, &bs);
725: VecSetBlockSize(newcoords, bs);
726: VecSetType(newcoords,VECSTANDARD);
727: VecCopy(coords, newcoords);
728: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
729: VecScale(newcoords, lengthScale);
730: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
731: DMPlexGlobalVectorView(dm, viewer, cdm, newcoords);
732: PetscViewerPopFormat(viewer);
733: VecDestroy(&newcoords);
734: return 0;
735: }
737: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
738: {
739: DM cdm;
740: Vec coordinatesLocal, newcoords;
741: PetscSection cSection, cGlobalSection;
742: PetscScalar *coords, *ncoords;
743: DMLabel cutLabel, cutVertexLabel = NULL;
744: const PetscReal *L;
745: const DMBoundaryType *bd;
746: PetscReal lengthScale;
747: PetscInt vStart, vEnd, v, bs, N, coordSize, dof, off, d;
748: PetscBool localized, embedded;
749: hid_t fileId, groupId;
751: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
752: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
753: DMGetCoordinatesLocal(dm, &coordinatesLocal);
754: VecGetBlockSize(coordinatesLocal, &bs);
755: DMGetCoordinatesLocalized(dm, &localized);
756: if (localized == PETSC_FALSE) return 0;
757: DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
758: DMGetCoordinateDM(dm, &cdm);
759: DMGetLocalSection(cdm, &cSection);
760: DMGetGlobalSection(cdm, &cGlobalSection);
761: DMGetLabel(dm, "periodic_cut", &cutLabel);
762: N = 0;
764: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
765: VecCreate(PetscObjectComm((PetscObject) dm), &newcoords);
766: PetscSectionGetDof(cSection, vStart, &dof);
767: PetscPrintf(PETSC_COMM_SELF, "DOF: %D\n", dof);
768: embedded = (PetscBool) (L && dof == 2 && !cutLabel);
769: if (cutVertexLabel) {
770: DMLabelGetStratumSize(cutVertexLabel, 1, &v);
771: N += dof*v;
772: }
773: for (v = vStart; v < vEnd; ++v) {
774: PetscSectionGetDof(cGlobalSection, v, &dof);
775: if (dof < 0) continue;
776: if (embedded) N += dof+1;
777: else N += dof;
778: }
779: if (embedded) VecSetBlockSize(newcoords, bs+1);
780: else VecSetBlockSize(newcoords, bs);
781: VecSetSizes(newcoords, N, PETSC_DETERMINE);
782: VecSetType(newcoords, VECSTANDARD);
783: VecGetArray(coordinatesLocal, &coords);
784: VecGetArray(newcoords, &ncoords);
785: coordSize = 0;
786: for (v = vStart; v < vEnd; ++v) {
787: PetscSectionGetDof(cGlobalSection, v, &dof);
788: PetscSectionGetOffset(cSection, v, &off);
789: if (dof < 0) continue;
790: if (embedded) {
791: if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
792: PetscReal theta, phi, r, R;
793: /* XY-periodic */
794: /* Suppose its an y-z circle, then
795: \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
796: and the circle in that plane is
797: \hat r cos(phi) + \hat x sin(phi) */
798: theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
799: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
800: r = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
801: R = L[1]/(2.0*PETSC_PI);
802: ncoords[coordSize++] = PetscSinReal(phi) * r;
803: ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
804: ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
805: } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
806: /* X-periodic */
807: ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
808: ncoords[coordSize++] = coords[off+1];
809: ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
810: } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
811: /* Y-periodic */
812: ncoords[coordSize++] = coords[off+0];
813: ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
814: ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
815: } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
816: PetscReal phi, r, R;
817: /* Mobius strip */
818: /* Suppose its an x-z circle, then
819: \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
820: and in that plane we rotate by pi as we go around the circle
821: \hat r cos(phi/2) + \hat y sin(phi/2) */
822: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
823: R = L[0];
824: r = PetscRealPart(coords[off+1]) - L[1]/2.0;
825: ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
826: ncoords[coordSize++] = PetscSinReal(phi/2.0) * r;
827: ncoords[coordSize++] = PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
828: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
829: } else {
830: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
831: }
832: }
833: if (cutVertexLabel) {
834: IS vertices;
835: const PetscInt *verts;
836: PetscInt n;
838: DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
839: if (vertices) {
840: ISGetIndices(vertices, &verts);
841: ISGetLocalSize(vertices, &n);
842: for (v = 0; v < n; ++v) {
843: PetscSectionGetDof(cSection, verts[v], &dof);
844: PetscSectionGetOffset(cSection, verts[v], &off);
845: for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
846: }
847: ISRestoreIndices(vertices, &verts);
848: ISDestroy(&vertices);
849: }
850: }
852: DMLabelDestroy(&cutVertexLabel);
853: VecRestoreArray(coordinatesLocal, &coords);
854: VecRestoreArray(newcoords, &ncoords);
855: PetscObjectSetName((PetscObject) newcoords, "vertices");
856: VecScale(newcoords, lengthScale);
857: PetscViewerHDF5PushGroup(viewer, "/viz");
858: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
859: PetscStackCallHDF5(H5Gclose,(groupId));
860: PetscViewerHDF5PopGroup(viewer);
861: PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
862: VecView(newcoords, viewer);
863: PetscViewerHDF5PopGroup(viewer);
864: VecDestroy(&newcoords);
865: return 0;
866: }
868: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
869: {
870: const char *topologydm_name;
871: const PetscInt *gpoint;
872: PetscInt numLabels, l;
873: DMPlexStorageVersion version;
874: char group[PETSC_MAX_PATH_LEN];
876: DMPlexStorageVersionSetUpWriting_Private(dm, viewer, &version);
877: ISGetIndices(globalPointNumbers, &gpoint);
878: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
879: if (version.major <= 1) {
880: PetscStrcpy(group, "/labels");
881: } else {
882: PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name);
883: }
884: PetscViewerHDF5PushGroup(viewer, group);
885: DMGetNumLabels(dm, &numLabels);
886: for (l = 0; l < numLabels; ++l) {
887: DMLabel label;
888: const char *name;
889: IS valueIS, pvalueIS, globalValueIS;
890: const PetscInt *values;
891: PetscInt numValues, v;
892: PetscBool isDepth, output;
894: DMGetLabelByNum(dm, l, &label);
895: PetscObjectGetName((PetscObject)label, &name);
896: DMGetLabelOutput(dm, name, &output);
897: PetscStrncmp(name, "depth", 10, &isDepth);
898: if (isDepth || !output) continue;
899: PetscViewerHDF5PushGroup(viewer, name);
900: DMLabelGetValueIS(label, &valueIS);
901: /* Must copy to a new IS on the global comm */
902: ISGetLocalSize(valueIS, &numValues);
903: ISGetIndices(valueIS, &values);
904: ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
905: ISRestoreIndices(valueIS, &values);
906: ISAllGather(pvalueIS, &globalValueIS);
907: ISDestroy(&pvalueIS);
908: ISSortRemoveDups(globalValueIS);
909: ISGetLocalSize(globalValueIS, &numValues);
910: ISGetIndices(globalValueIS, &values);
911: for (v = 0; v < numValues; ++v) {
912: IS stratumIS, globalStratumIS;
913: const PetscInt *spoints = NULL;
914: PetscInt *gspoints, n = 0, gn, p;
915: const char *iname = "indices";
916: char group[PETSC_MAX_PATH_LEN];
918: PetscSNPrintf(group, sizeof(group), "%D", values[v]);
919: PetscViewerHDF5PushGroup(viewer, group);
920: DMLabelGetStratumIS(label, values[v], &stratumIS);
922: if (stratumIS) ISGetLocalSize(stratumIS, &n);
923: if (stratumIS) ISGetIndices(stratumIS, &spoints);
924: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
925: PetscMalloc1(gn,&gspoints);
926: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
927: if (stratumIS) ISRestoreIndices(stratumIS, &spoints);
928: ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
929: PetscObjectSetName((PetscObject) globalStratumIS, iname);
931: ISView(globalStratumIS, viewer);
932: ISDestroy(&globalStratumIS);
933: ISDestroy(&stratumIS);
934: PetscViewerHDF5PopGroup(viewer);
935: }
936: ISRestoreIndices(globalValueIS, &values);
937: ISDestroy(&globalValueIS);
938: ISDestroy(&valueIS);
939: PetscViewerHDF5PopGroup(viewer);
940: }
941: ISRestoreIndices(globalPointNumbers, &gpoint);
942: PetscViewerHDF5PopGroup(viewer);
943: return 0;
944: }
946: /* We only write cells and vertices. Does this screw up parallel reading? */
947: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
948: {
949: IS globalPointNumbers;
950: PetscViewerFormat format;
951: PetscBool viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
953: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
954: DMPlexCoordinatesView_HDF5_Internal(dm, viewer);
956: PetscViewerGetFormat(viewer, &format);
957: switch (format) {
958: case PETSC_VIEWER_HDF5_VIZ:
959: viz_geom = PETSC_TRUE;
960: xdmf_topo = PETSC_TRUE;
961: break;
962: case PETSC_VIEWER_HDF5_XDMF:
963: xdmf_topo = PETSC_TRUE;
964: break;
965: case PETSC_VIEWER_HDF5_PETSC:
966: petsc_topo = PETSC_TRUE;
967: break;
968: case PETSC_VIEWER_DEFAULT:
969: case PETSC_VIEWER_NATIVE:
970: viz_geom = PETSC_TRUE;
971: xdmf_topo = PETSC_TRUE;
972: petsc_topo = PETSC_TRUE;
973: break;
974: default:
975: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
976: }
978: if (viz_geom) DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);
979: if (xdmf_topo) DMPlexWriteTopology_Vertices_HDF5_Static(dm, globalPointNumbers, viewer);
980: if (petsc_topo) {
981: DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer);
982: DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer);
983: }
985: ISDestroy(&globalPointNumbers);
986: return 0;
987: }
989: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
990: {
991: MPI_Comm comm;
992: const char *topologydm_name;
993: const char *sectiondm_name;
994: PetscSection gsection;
996: PetscObjectGetComm((PetscObject)sectiondm, &comm);
997: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
998: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
999: PetscViewerHDF5PushGroup(viewer, "topologies");
1000: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1001: PetscViewerHDF5PushGroup(viewer, "dms");
1002: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1003: DMGetGlobalSection(sectiondm, &gsection);
1004: /* Save raw section */
1005: PetscSectionView(gsection, viewer);
1006: /* Save plex wrapper */
1007: {
1008: PetscInt pStart, pEnd, p, n;
1009: IS globalPointNumbers;
1010: const PetscInt *gpoints;
1011: IS orderIS;
1012: PetscInt *order;
1014: PetscSectionGetChart(gsection, &pStart, &pEnd);
1015: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
1016: ISGetIndices(globalPointNumbers, &gpoints);
1017: for (p = pStart, n = 0; p < pEnd; ++p) if (gpoints[p] >= 0) n++;
1018: /* "order" is an array of global point numbers.
1019: When loading, it is used with topology/order array
1020: to match section points with plex topology points. */
1021: PetscMalloc1(n, &order);
1022: for (p = pStart, n = 0; p < pEnd; ++p) if (gpoints[p] >= 0) order[n++] = gpoints[p];
1023: ISRestoreIndices(globalPointNumbers, &gpoints);
1024: ISDestroy(&globalPointNumbers);
1025: ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS);
1026: PetscObjectSetName((PetscObject)orderIS, "order");
1027: ISView(orderIS, viewer);
1028: ISDestroy(&orderIS);
1029: }
1030: PetscViewerHDF5PopGroup(viewer);
1031: PetscViewerHDF5PopGroup(viewer);
1032: PetscViewerHDF5PopGroup(viewer);
1033: PetscViewerHDF5PopGroup(viewer);
1034: return 0;
1035: }
1037: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1038: {
1039: const char *topologydm_name;
1040: const char *sectiondm_name;
1041: const char *vec_name;
1042: PetscInt bs;
1044: /* Check consistency */
1045: {
1046: PetscSF pointsf, pointsf1;
1048: DMGetPointSF(dm, &pointsf);
1049: DMGetPointSF(sectiondm, &pointsf1);
1051: }
1052: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1053: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1054: PetscObjectGetName((PetscObject)vec, &vec_name);
1055: PetscViewerHDF5PushGroup(viewer, "topologies");
1056: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1057: PetscViewerHDF5PushGroup(viewer, "dms");
1058: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1059: PetscViewerHDF5PushGroup(viewer, "vecs");
1060: PetscViewerHDF5PushGroup(viewer, vec_name);
1061: VecGetBlockSize(vec, &bs);
1062: PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *) &bs);
1063: VecSetBlockSize(vec, 1);
1064: /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but, */
1065: /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view */
1066: /* is set to VecView_Plex, which would save vec in a predefined location. */
1067: /* To save vec in where we want, we create a new Vec (temp) with */
1068: /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
1069: {
1070: Vec temp;
1071: const PetscScalar *array;
1072: PetscLayout map;
1074: VecCreate(PetscObjectComm((PetscObject)vec), &temp);
1075: PetscObjectSetName((PetscObject)temp, vec_name);
1076: VecGetLayout(vec, &map);
1077: VecSetLayout(temp, map);
1078: VecSetUp(temp);
1079: VecGetArrayRead(vec, &array);
1080: VecPlaceArray(temp, array);
1081: VecView(temp, viewer);
1082: VecResetArray(temp);
1083: VecRestoreArrayRead(vec, &array);
1084: VecDestroy(&temp);
1085: }
1086: VecSetBlockSize(vec, bs);
1087: PetscViewerHDF5PopGroup(viewer);
1088: PetscViewerHDF5PopGroup(viewer);
1089: PetscViewerHDF5PopGroup(viewer);
1090: PetscViewerHDF5PopGroup(viewer);
1091: PetscViewerHDF5PopGroup(viewer);
1092: PetscViewerHDF5PopGroup(viewer);
1093: return 0;
1094: }
1096: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1097: {
1098: MPI_Comm comm;
1099: const char *topologydm_name;
1100: const char *sectiondm_name;
1101: const char *vec_name;
1102: PetscSection section;
1103: PetscBool includesConstraints;
1104: Vec gvec;
1105: PetscInt m, bs;
1107: PetscObjectGetComm((PetscObject)dm, &comm);
1108: /* Check consistency */
1109: {
1110: PetscSF pointsf, pointsf1;
1112: DMGetPointSF(dm, &pointsf);
1113: DMGetPointSF(sectiondm, &pointsf1);
1115: }
1116: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1117: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1118: PetscObjectGetName((PetscObject)vec, &vec_name);
1119: PetscViewerHDF5PushGroup(viewer, "topologies");
1120: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1121: PetscViewerHDF5PushGroup(viewer, "dms");
1122: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1123: PetscViewerHDF5PushGroup(viewer, "vecs");
1124: PetscViewerHDF5PushGroup(viewer, vec_name);
1125: VecGetBlockSize(vec, &bs);
1126: PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *) &bs);
1127: VecCreate(comm, &gvec);
1128: PetscObjectSetName((PetscObject)gvec, vec_name);
1129: DMGetGlobalSection(sectiondm, §ion);
1130: PetscSectionGetIncludesConstraints(section, &includesConstraints);
1131: if (includesConstraints) PetscSectionGetStorageSize(section, &m);
1132: else PetscSectionGetConstrainedStorageSize(section, &m);
1133: VecSetSizes(gvec, m, PETSC_DECIDE);
1134: VecSetUp(gvec);
1135: DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec);
1136: DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec);
1137: VecView(gvec, viewer);
1138: VecDestroy(&gvec);
1139: PetscViewerHDF5PopGroup(viewer);
1140: PetscViewerHDF5PopGroup(viewer);
1141: PetscViewerHDF5PopGroup(viewer);
1142: PetscViewerHDF5PopGroup(viewer);
1143: PetscViewerHDF5PopGroup(viewer);
1144: PetscViewerHDF5PopGroup(viewer);
1145: return 0;
1146: }
1148: struct _n_LoadLabelsCtx {
1149: MPI_Comm comm;
1150: PetscMPIInt rank;
1151: DM dm;
1152: PetscViewer viewer;
1153: DMLabel label;
1154: PetscSF sfXC;
1155: PetscLayout layoutX;
1156: };
1157: typedef struct _n_LoadLabelsCtx *LoadLabelsCtx;
1159: static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx)
1160: {
1161: PetscNew(ctx);
1162: PetscObjectReference((PetscObject) ((*ctx)->dm = dm));
1163: PetscObjectReference((PetscObject) ((*ctx)->viewer = viewer));
1164: PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm);
1165: MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank);
1166: (*ctx)->sfXC = sfXC;
1167: if (sfXC) {
1168: PetscInt nX;
1170: PetscObjectReference((PetscObject) sfXC);
1171: PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL);
1172: PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX);
1173: }
1174: return 0;
1175: }
1177: static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx)
1178: {
1179: if (!*ctx) return 0;
1180: DMDestroy(&(*ctx)->dm);
1181: PetscViewerDestroy(&(*ctx)->viewer);
1182: PetscSFDestroy(&(*ctx)->sfXC);
1183: PetscLayoutDestroy(&(*ctx)->layoutX);
1184: PetscFree(*ctx);
1185: return 0;
1186: }
1188: /*
1189: A: on-disk points
1190: X: global points [0, NX)
1191: C: distributed plex points
1192: */
1193: static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS)
1194: {
1195: MPI_Comm comm = ctx->comm;
1196: PetscSF sfXC = ctx->sfXC;
1197: PetscLayout layoutX = ctx->layoutX;
1198: PetscSF sfXA;
1199: const PetscInt *A_points;
1200: PetscInt nX, nC;
1201: PetscInt n;
1203: PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL);
1204: ISGetLocalSize(stratumIS, &n);
1205: ISGetIndices(stratumIS, &A_points);
1206: PetscSFCreate(comm, &sfXA);
1207: PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points);
1208: ISCreate(comm, newStratumIS);
1209: ISSetType(*newStratumIS,ISGENERAL);
1210: {
1211: PetscInt i;
1212: PetscBool *A_mask, *X_mask, *C_mask;
1214: PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask);
1215: for (i=0; i<n; i++) A_mask[i] = PETSC_TRUE;
1216: PetscSFReduceBegin(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE);
1217: PetscSFReduceEnd( sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE);
1218: PetscSFBcastBegin( sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR);
1219: PetscSFBcastEnd( sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR);
1220: ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask);
1221: PetscFree3(A_mask, X_mask, C_mask);
1222: }
1223: PetscSFDestroy(&sfXA);
1224: ISRestoreIndices(stratumIS, &A_points);
1225: return 0;
1226: }
1228: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data)
1229: {
1230: LoadLabelsCtx ctx = (LoadLabelsCtx) op_data;
1231: PetscViewer viewer = ctx->viewer;
1232: DMLabel label = ctx->label;
1233: MPI_Comm comm = ctx->comm;
1234: IS stratumIS;
1235: const PetscInt *ind;
1236: PetscInt value, N, i;
1238: PetscOptionsStringToInt(vname, &value);
1239: ISCreate(comm, &stratumIS);
1240: PetscObjectSetName((PetscObject) stratumIS, "indices");
1241: PetscViewerHDF5PushGroup(viewer, vname); /* labels/<lname>/<vname> */
1243: if (!ctx->sfXC) {
1244: /* Force serial load */
1245: PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
1246: PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0);
1247: PetscLayoutSetSize(stratumIS->map, N);
1248: }
1249: ISLoad(stratumIS, viewer);
1251: if (ctx->sfXC) {
1252: IS newStratumIS;
1254: ReadLabelStratumHDF5_Distribute_Private(stratumIS, ctx, &newStratumIS);
1255: ISDestroy(&stratumIS);
1256: stratumIS = newStratumIS;
1257: }
1259: PetscViewerHDF5PopGroup(viewer);
1260: ISGetLocalSize(stratumIS, &N);
1261: ISGetIndices(stratumIS, &ind);
1262: for (i = 0; i < N; ++i) DMLabelSetValue(label, ind[i], value);
1263: ISRestoreIndices(stratumIS, &ind);
1264: ISDestroy(&stratumIS);
1265: return 0;
1266: }
1268: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data)
1269: {
1270: LoadLabelsCtx ctx = (LoadLabelsCtx) op_data;
1271: DM dm = ctx->dm;
1272: hsize_t idx = 0;
1274: PetscBool flg;
1275: herr_t err;
1277: DMHasLabel(dm, lname, &flg);
1278: if (flg) DMRemoveLabel(dm, lname, NULL);
1279: DMCreateLabel(dm, lname); if (ierr) return (herr_t) ierr;
1280: DMGetLabel(dm, lname, &ctx->label); if (ierr) return (herr_t) ierr;
1281: PetscViewerHDF5PushGroup(ctx->viewer, lname); /* labels/<lname> */
1282: /* Iterate over the label's strata */
1283: PetscStackCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1284: PetscViewerHDF5PopGroup(ctx->viewer);
1285: return err;
1286: }
1288: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1289: {
1290: const char *topologydm_name;
1291: LoadLabelsCtx ctx;
1292: hsize_t idx = 0;
1293: char group[PETSC_MAX_PATH_LEN];
1294: DMPlexStorageVersion version;
1295: PetscBool distributed, hasGroup;
1297: DMPlexIsDistributed(dm, &distributed);
1298: if (distributed) {
1300: }
1301: LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx);
1302: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1303: DMPlexStorageVersionGet_Private(dm, viewer, &version);
1304: if (version.major <= 1) {
1305: PetscStrcpy(group, "labels");
1306: } else {
1307: PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name);
1308: }
1309: PetscViewerHDF5PushGroup(viewer, group);
1310: PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup);
1311: if (hasGroup) {
1312: hid_t fileId, groupId;
1314: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
1315: /* Iterate over labels */
1316: PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx));
1317: PetscStackCallHDF5(H5Gclose,(groupId));
1318: }
1319: PetscViewerHDF5PopGroup(viewer);
1320: LoadLabelsCtxDestroy(&ctx);
1321: return 0;
1322: }
1324: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sf)
1325: {
1326: MPI_Comm comm;
1327: const char *topologydm_name;
1328: const char *pointsName, *coneSizesName, *conesName, *orientationsName;
1329: IS pointsIS, coneSizesIS, conesIS, orientationsIS;
1330: const PetscInt *points, *coneSizes, *cones, *orientations;
1331: PetscInt *cone, *ornt;
1332: PetscInt dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
1333: PetscMPIInt size, rank;
1334: char group[PETSC_MAX_PATH_LEN];
1335: DMPlexStorageVersion version;
1337: pointsName = "order";
1338: coneSizesName = "cones";
1339: conesName = "cells";
1340: orientationsName = "orientation";
1341: PetscObjectGetComm((PetscObject)dm, &comm);
1342: MPI_Comm_size(comm, &size);
1343: MPI_Comm_rank(comm, &rank);
1344: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1345: DMPlexStorageVersionGet_Private(dm, viewer, &version);
1346: if (version.major <= 1) {
1347: PetscStrcpy(group, "/topology");
1348: } else {
1349: PetscSNPrintf(group, sizeof(group), "topologies/%s/topology", topologydm_name);
1350: }
1351: PetscViewerHDF5PushGroup(viewer, group);
1352: ISCreate(comm, &pointsIS);
1353: PetscObjectSetName((PetscObject) pointsIS, pointsName);
1354: ISCreate(comm, &coneSizesIS);
1355: PetscObjectSetName((PetscObject) coneSizesIS, coneSizesName);
1356: ISCreate(comm, &conesIS);
1357: PetscObjectSetName((PetscObject) conesIS, conesName);
1358: ISCreate(comm, &orientationsIS);
1359: PetscObjectSetName((PetscObject) orientationsIS, orientationsName);
1360: PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject) conesIS, "cell_dim", PETSC_INT, NULL, &dim);
1361: DMSetDimension(dm, dim);
1362: {
1363: /* Force serial load */
1364: PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np);
1365: PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0);
1366: PetscLayoutSetSize(pointsIS->map, Np);
1367: pEnd = rank == 0 ? Np : 0;
1368: PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np);
1369: PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0);
1370: PetscLayoutSetSize(coneSizesIS->map, Np);
1371: PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N);
1372: PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0);
1373: PetscLayoutSetSize(conesIS->map, N);
1374: PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N);
1375: PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0);
1376: PetscLayoutSetSize(orientationsIS->map, N);
1377: }
1378: ISLoad(pointsIS, viewer);
1379: ISLoad(coneSizesIS, viewer);
1380: ISLoad(conesIS, viewer);
1381: ISLoad(orientationsIS, viewer);
1382: PetscViewerHDF5PopGroup(viewer);
1383: /* Create Plex */
1384: DMPlexSetChart(dm, 0, pEnd);
1385: ISGetIndices(pointsIS, &points);
1386: ISGetIndices(coneSizesIS, &coneSizes);
1387: for (p = 0; p < pEnd; ++p) {
1388: DMPlexSetConeSize(dm, points[p], coneSizes[p]);
1389: maxConeSize = PetscMax(maxConeSize, coneSizes[p]);
1390: }
1391: DMSetUp(dm);
1392: ISGetIndices(conesIS, &cones);
1393: ISGetIndices(orientationsIS, &orientations);
1394: PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
1395: for (p = 0, q = 0; p < pEnd; ++p) {
1396: for (c = 0; c < coneSizes[p]; ++c, ++q) {
1397: cone[c] = cones[q];
1398: ornt[c] = orientations[q];
1399: }
1400: DMPlexSetCone(dm, points[p], cone);
1401: DMPlexSetConeOrientation(dm, points[p], ornt);
1402: }
1403: PetscFree2(cone,ornt);
1404: /* Create global section migration SF */
1405: if (sf) {
1406: PetscLayout layout;
1407: PetscInt *globalIndices;
1409: PetscMalloc1(pEnd, &globalIndices);
1410: /* plex point == globalPointNumber in this case */
1411: for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
1412: PetscLayoutCreate(comm, &layout);
1413: PetscLayoutSetSize(layout, Np);
1414: PetscLayoutSetBlockSize(layout, 1);
1415: PetscLayoutSetUp(layout);
1416: PetscSFCreate(comm, sf);
1417: PetscSFSetFromOptions(*sf);
1418: PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices);
1419: PetscLayoutDestroy(&layout);
1420: PetscFree(globalIndices);
1421: }
1422: /* Clean-up */
1423: ISRestoreIndices(pointsIS, &points);
1424: ISRestoreIndices(coneSizesIS, &coneSizes);
1425: ISRestoreIndices(conesIS, &cones);
1426: ISRestoreIndices(orientationsIS, &orientations);
1427: ISDestroy(&pointsIS);
1428: ISDestroy(&coneSizesIS);
1429: ISDestroy(&conesIS);
1430: ISDestroy(&orientationsIS);
1431: /* Fill in the rest of the topology structure */
1432: DMPlexSymmetrize(dm);
1433: DMPlexStratify(dm);
1434: return 0;
1435: }
1437: /* If the file is old, it not only has different path to the coordinates, but */
1438: /* does not contain coordinateDMs, so must fall back to the old implementation. */
1439: static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
1440: {
1441: PetscSection coordSection;
1442: Vec coordinates;
1443: PetscReal lengthScale;
1444: PetscInt spatialDim, N, numVertices, vStart, vEnd, v;
1445: PetscMPIInt rank;
1447: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1448: /* Read geometry */
1449: PetscViewerHDF5PushGroup(viewer, "/geometry");
1450: VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
1451: PetscObjectSetName((PetscObject) coordinates, "vertices");
1452: {
1453: /* Force serial load */
1454: PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
1455: VecSetSizes(coordinates, !rank ? N : 0, N);
1456: VecSetBlockSize(coordinates, spatialDim);
1457: }
1458: VecLoad(coordinates, viewer);
1459: PetscViewerHDF5PopGroup(viewer);
1460: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1461: VecScale(coordinates, 1.0/lengthScale);
1462: VecGetLocalSize(coordinates, &numVertices);
1463: VecGetBlockSize(coordinates, &spatialDim);
1464: numVertices /= spatialDim;
1465: /* Create coordinates */
1466: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1468: DMGetCoordinateSection(dm, &coordSection);
1469: PetscSectionSetNumFields(coordSection, 1);
1470: PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
1471: PetscSectionSetChart(coordSection, vStart, vEnd);
1472: for (v = vStart; v < vEnd; ++v) {
1473: PetscSectionSetDof(coordSection, v, spatialDim);
1474: PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
1475: }
1476: PetscSectionSetUp(coordSection);
1477: DMSetCoordinates(dm, coordinates);
1478: VecDestroy(&coordinates);
1479: return 0;
1480: }
1482: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1483: {
1484: DM cdm;
1485: Vec coords;
1486: PetscInt blockSize;
1487: PetscReal lengthScale;
1488: PetscSF lsf;
1489: const char *topologydm_name;
1490: char *coordinatedm_name, *coordinates_name;
1492: {
1493: DMPlexStorageVersion version;
1495: DMPlexStorageVersionGet_Private(dm, viewer, &version);
1496: if (version.major <= 1) {
1497: DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer);
1498: return 0;
1499: }
1500: }
1502: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1503: PetscViewerHDF5PushGroup(viewer, "topologies");
1504: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1505: PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING , NULL, &coordinatedm_name);
1506: PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING , NULL, &coordinates_name);
1507: PetscViewerHDF5PopGroup(viewer);
1508: PetscViewerHDF5PopGroup(viewer);
1509: DMGetCoordinateDM(dm, &cdm);
1510: PetscObjectSetName((PetscObject)cdm, coordinatedm_name);
1511: PetscFree(coordinatedm_name);
1512: /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
1513: DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf);
1514: DMCreateLocalVector(cdm, &coords);
1515: PetscObjectSetName((PetscObject)coords, coordinates_name);
1516: PetscFree(coordinates_name);
1517: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
1518: DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords);
1519: PetscViewerPopFormat(viewer);
1520: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1521: VecScale(coords, 1.0/lengthScale);
1522: DMSetCoordinatesLocal(dm, coords);
1523: VecGetBlockSize(coords, &blockSize);
1524: DMSetCoordinateDim(dm, blockSize);
1525: VecDestroy(&coords);
1526: PetscSFDestroy(&lsf);
1527: return 0;
1528: }
1530: static PetscErrorCode DMPlexLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
1531: {
1532: DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL);
1533: DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL);
1534: DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer);
1535: return 0;
1536: }
1538: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1539: {
1540: PetscSF sfXC;
1542: {
1543: DMPlexStorageVersion version;
1545: DMPlexStorageVersionGet_Private(dm, viewer, &version);
1546: if (version.major <= 1) {
1547: DMPlexLoad_HDF5_Legacy_Private(dm, viewer);
1548: return 0;
1549: }
1550: }
1551: DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC);
1552: DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC);
1553: DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC);
1554: PetscSFDestroy(&sfXC);
1555: return 0;
1556: }
1558: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1559: {
1560: MPI_Comm comm;
1561: PetscInt pStart, pEnd, p, m;
1562: PetscInt *goffs, *ilocal;
1563: PetscBool rootIncludeConstraints, leafIncludeConstraints;
1565: PetscObjectGetComm((PetscObject)leafSection, &comm);
1566: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1567: PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints);
1568: PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints);
1569: if (rootIncludeConstraints && leafIncludeConstraints) PetscSectionGetStorageSize(leafSection, &m);
1570: else PetscSectionGetConstrainedStorageSize(leafSection, &m);
1571: PetscMalloc1(m, &ilocal);
1572: PetscMalloc1(m, &goffs);
1573: /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
1574: /* for the top-level section (not for each field), so one must have */
1575: /* rootSection->pointMajor == PETSC_TRUE. */
1577: /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
1579: for (p = pStart, m = 0; p < pEnd; ++p) {
1580: PetscInt dof, cdof, i, j, off, goff;
1581: const PetscInt *cinds;
1583: PetscSectionGetDof(leafSection, p, &dof);
1584: if (dof < 0) continue;
1585: goff = globalOffsets[p-pStart];
1586: PetscSectionGetOffset(leafSection, p, &off);
1587: PetscSectionGetConstraintDof(leafSection, p, &cdof);
1588: PetscSectionGetConstraintIndices(leafSection, p, &cinds);
1589: for (i = 0, j = 0; i < dof; ++i) {
1590: PetscBool constrained = (PetscBool) (j < cdof && i == cinds[j]);
1592: if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {ilocal[m] = off++; goffs[m++] = goff++;}
1593: else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
1594: else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff;
1595: if (constrained) ++j;
1596: }
1597: }
1598: PetscSFCreate(comm, sectionSF);
1599: PetscSFSetFromOptions(*sectionSF);
1600: PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs);
1601: PetscFree(goffs);
1602: return 0;
1603: }
1605: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
1606: {
1607: MPI_Comm comm;
1608: PetscMPIInt size, rank;
1609: const char *topologydm_name;
1610: const char *sectiondm_name;
1611: PetscSection sectionA, sectionB;
1612: PetscInt nX, n, i;
1613: PetscSF sfAB;
1615: PetscObjectGetComm((PetscObject)dm, &comm);
1616: MPI_Comm_size(comm, &size);
1617: MPI_Comm_rank(comm, &rank);
1618: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1619: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1620: PetscViewerHDF5PushGroup(viewer, "topologies");
1621: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1622: PetscViewerHDF5PushGroup(viewer, "dms");
1623: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1624: /* A: on-disk points */
1625: /* X: list of global point numbers, [0, NX) */
1626: /* B: plex points */
1627: /* Load raw section (sectionA) */
1628: PetscSectionCreate(comm, §ionA);
1629: PetscSectionLoad(sectionA, viewer);
1630: PetscSectionGetChart(sectionA, NULL, &n);
1631: /* Create sfAB: A -> B */
1632: #if defined(PETSC_USE_DEBUG)
1633: {
1634: PetscInt N, N1;
1636: PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1);
1637: MPI_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm);
1639: }
1640: #endif
1641: {
1642: IS orderIS;
1643: const PetscInt *gpoints;
1644: PetscSF sfXA, sfAX;
1645: PetscLayout layout;
1646: PetscSFNode *owners, *buffer;
1647: PetscInt nleaves;
1648: PetscInt *ilocal;
1649: PetscSFNode *iremote;
1651: /* Create sfAX: A -> X */
1652: ISCreate(comm, &orderIS);
1653: PetscObjectSetName((PetscObject)orderIS, "order");
1654: PetscLayoutSetLocalSize(orderIS->map, n);
1655: ISLoad(orderIS, viewer);
1656: PetscLayoutCreate(comm, &layout);
1657: PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL);
1658: PetscLayoutSetLocalSize(layout, nX);
1659: PetscLayoutSetBlockSize(layout, 1);
1660: PetscLayoutSetUp(layout);
1661: PetscSFCreate(comm, &sfXA);
1662: ISGetIndices(orderIS, &gpoints);
1663: PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints);
1664: ISRestoreIndices(orderIS, &gpoints);
1665: ISDestroy(&orderIS);
1666: PetscLayoutDestroy(&layout);
1667: PetscMalloc1(n, &owners);
1668: PetscMalloc1(nX, &buffer);
1669: for (i = 0; i < n; ++i) {owners[i].rank = rank; owners[i].index = i;}
1670: for (i = 0; i < nX; ++i) {buffer[i].rank = -1; buffer[i].index = -1;}
1671: PetscSFReduceBegin(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
1672: PetscSFReduceEnd(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
1673: PetscSFDestroy(&sfXA);
1674: PetscFree(owners);
1675: for (i = 0, nleaves = 0; i < nX; ++i) if (buffer[i].rank >= 0) nleaves++;
1676: PetscMalloc1(nleaves, &ilocal);
1677: PetscMalloc1(nleaves, &iremote);
1678: for (i = 0, nleaves = 0; i < nX; ++i) {
1679: if (buffer[i].rank >= 0) {
1680: ilocal[nleaves] = i;
1681: iremote[nleaves].rank = buffer[i].rank;
1682: iremote[nleaves].index = buffer[i].index;
1683: nleaves++;
1684: }
1685: }
1686: PetscSFCreate(comm, &sfAX);
1687: PetscSFSetFromOptions(sfAX);
1688: PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1689: /* Fix PetscSFCompose() and replace the code-block below with: */
1690: /* PetscSFCompose(sfAX, sfXB, &sfAB); */
1691: /* which currently causes segmentation fault due to sparse map. */
1692: {
1693: PetscInt npoints;
1694: PetscInt mleaves;
1695: PetscInt *jlocal;
1696: PetscSFNode *jremote;
1698: PetscSFGetGraph(sfXB, NULL, &npoints, NULL, NULL);
1699: PetscMalloc1(npoints, &owners);
1700: for (i = 0; i < npoints; ++i) {owners[i].rank = -1; owners[i].index = -1;}
1701: PetscSFBcastBegin(sfXB, MPIU_2INT, buffer, owners, MPI_REPLACE);
1702: PetscSFBcastEnd(sfXB, MPIU_2INT, buffer, owners, MPI_REPLACE);
1703: for (i = 0, mleaves = 0; i < npoints; ++i) if (owners[i].rank >= 0) mleaves++;
1704: PetscMalloc1(mleaves, &jlocal);
1705: PetscMalloc1(mleaves, &jremote);
1706: for (i = 0, mleaves = 0; i < npoints; ++i) {
1707: if (owners[i].rank >= 0) {
1708: jlocal[mleaves] = i;
1709: jremote[mleaves].rank = owners[i].rank;
1710: jremote[mleaves].index = owners[i].index;
1711: mleaves++;
1712: }
1713: }
1714: PetscSFCreate(comm, &sfAB);
1715: PetscSFSetFromOptions(sfAB);
1716: PetscSFSetGraph(sfAB, n, mleaves, jlocal, PETSC_OWN_POINTER, jremote, PETSC_OWN_POINTER);
1717: PetscFree(owners);
1718: }
1719: PetscFree(buffer);
1720: PetscSFDestroy(&sfAX);
1721: }
1722: PetscViewerHDF5PopGroup(viewer);
1723: PetscViewerHDF5PopGroup(viewer);
1724: PetscViewerHDF5PopGroup(viewer);
1725: PetscViewerHDF5PopGroup(viewer);
1726: /* Create plex section (sectionB) */
1727: DMGetLocalSection(sectiondm, §ionB);
1728: if (lsf || gsf) {
1729: PetscLayout layout;
1730: PetscInt M, m;
1731: PetscInt *offsetsA;
1732: PetscBool includesConstraintsA;
1734: PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB);
1735: PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA);
1736: if (includesConstraintsA) PetscSectionGetStorageSize(sectionA, &m);
1737: else PetscSectionGetConstrainedStorageSize(sectionA, &m);
1738: MPI_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm);
1739: PetscLayoutCreate(comm, &layout);
1740: PetscLayoutSetSize(layout, M);
1741: PetscLayoutSetUp(layout);
1742: if (lsf) {
1743: PetscSF lsfABdata;
1745: DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata);
1746: *lsf = lsfABdata;
1747: }
1748: if (gsf) {
1749: PetscSection gsectionB, gsectionB1;
1750: PetscBool includesConstraintsB;
1751: PetscSF gsfABdata, pointsf;
1753: DMGetGlobalSection(sectiondm, &gsectionB1);
1754: PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB);
1755: DMGetPointSF(sectiondm, &pointsf);
1756: PetscSectionCreateGlobalSection(sectionB, pointsf, includesConstraintsB, PETSC_TRUE, &gsectionB);
1757: DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata);
1758: PetscSectionDestroy(&gsectionB);
1759: *gsf = gsfABdata;
1760: }
1761: PetscLayoutDestroy(&layout);
1762: PetscFree(offsetsA);
1763: } else {
1764: PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB);
1765: }
1766: PetscSFDestroy(&sfAB);
1767: PetscSectionDestroy(§ionA);
1768: return 0;
1769: }
1771: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
1772: {
1773: MPI_Comm comm;
1774: const char *topologydm_name;
1775: const char *sectiondm_name;
1776: const char *vec_name;
1777: Vec vecA;
1778: PetscInt mA, m, bs;
1779: const PetscInt *ilocal;
1780: const PetscScalar *src;
1781: PetscScalar *dest;
1783: PetscObjectGetComm((PetscObject)dm, &comm);
1784: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1785: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1786: PetscObjectGetName((PetscObject)vec, &vec_name);
1787: PetscViewerHDF5PushGroup(viewer, "topologies");
1788: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1789: PetscViewerHDF5PushGroup(viewer, "dms");
1790: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1791: PetscViewerHDF5PushGroup(viewer, "vecs");
1792: PetscViewerHDF5PushGroup(viewer, vec_name);
1793: VecCreate(comm, &vecA);
1794: PetscObjectSetName((PetscObject)vecA, vec_name);
1795: PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL);
1796: /* Check consistency */
1797: {
1798: PetscSF pointsf, pointsf1;
1799: PetscInt m1, i, j;
1801: DMGetPointSF(dm, &pointsf);
1802: DMGetPointSF(sectiondm, &pointsf1);
1804: #if defined(PETSC_USE_DEBUG)
1805: {
1806: PetscInt MA, MA1;
1808: MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm);
1809: PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1);
1811: }
1812: #endif
1813: VecGetLocalSize(vec, &m1);
1815: for (i = 0; i < m; ++i) {
1816: j = ilocal ? ilocal[i] : i;
1818: }
1819: }
1820: VecSetSizes(vecA, mA, PETSC_DECIDE);
1821: VecLoad(vecA, viewer);
1822: VecGetArrayRead(vecA, &src);
1823: VecGetArray(vec, &dest);
1824: PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
1825: PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
1826: VecRestoreArray(vec, &dest);
1827: VecRestoreArrayRead(vecA, &src);
1828: VecDestroy(&vecA);
1829: PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *) &bs);
1830: VecSetBlockSize(vec, bs);
1831: PetscViewerHDF5PopGroup(viewer);
1832: PetscViewerHDF5PopGroup(viewer);
1833: PetscViewerHDF5PopGroup(viewer);
1834: PetscViewerHDF5PopGroup(viewer);
1835: PetscViewerHDF5PopGroup(viewer);
1836: PetscViewerHDF5PopGroup(viewer);
1837: return 0;
1838: }
1839: #endif