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: /* Logging support */
8: PetscLogEvent DMPLEX_DistributionView, DMPLEX_DistributionLoad;
10: #if defined(PETSC_HAVE_HDF5)
11: static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *);
12: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer, DMPlexStorageVersion);
13: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer, const char[], DMPlexStorageVersion);
14: static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *);
16: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
18: static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer viewer, const char str[], DMPlexStorageVersion *version)
19: {
20: PetscToken t;
21: char *ts;
22: PetscInt i;
23: PetscInt ti[3];
24: DMPlexStorageVersion v;
26: PetscFunctionBegin;
27: PetscCall(PetscTokenCreate(str, '.', &t));
28: for (i = 0; i < 3; i++) {
29: PetscCall(PetscTokenFind(t, &ts));
30: PetscCheck(ts, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Malformed version string %s", str);
31: PetscCall(PetscOptionsStringToInt(ts, &ti[i]));
32: }
33: PetscCall(PetscTokenFind(t, &ts));
34: PetscCheck(!ts, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Malformed version string %s", str);
35: PetscCall(PetscTokenDestroy(&t));
36: PetscCall(PetscNew(&v));
37: v->major = ti[0];
38: v->minor = ti[1];
39: v->subminor = ti[2];
40: PetscCall(PetscViewerCheckVersion_Private(viewer, v));
41: *version = v;
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion v)
46: {
47: PetscContainer cont;
49: PetscFunctionBegin;
50: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)viewer), &cont));
51: PetscCall(PetscContainerSetPointer(cont, v));
52: PetscCall(PetscContainerSetUserDestroy(cont, PetscContainerUserDestroyDefault));
53: PetscCall(PetscObjectCompose((PetscObject)viewer, key, (PetscObject)cont));
54: PetscCall(PetscContainerDestroy(&cont));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion *v)
59: {
60: PetscContainer cont;
62: PetscFunctionBegin;
63: PetscCall(PetscObjectQuery((PetscObject)viewer, key, (PetscObject *)&cont));
64: *v = NULL;
65: if (cont) PetscCall(PetscContainerGetPointer(cont, (void **)v));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: /*
70: Version log:
71: 1.0.0 legacy version (default if no "dmplex_storage_version" attribute found in file)
72: 2.0.0 introduce versioning and multiple topologies
73: 2.1.0 introduce distributions
74: */
75: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer viewer, DMPlexStorageVersion version)
76: {
77: PetscBool valid = PETSC_FALSE;
79: PetscFunctionBegin;
80: switch (version->major) {
81: case 1:
82: switch (version->minor) {
83: case 0:
84: switch (version->subminor) {
85: case 0:
86: valid = PETSC_TRUE;
87: break;
88: }
89: break;
90: }
91: break;
92: case 2:
93: switch (version->minor) {
94: case 0:
95: switch (version->subminor) {
96: case 0:
97: valid = PETSC_TRUE;
98: break;
99: }
100: break;
101: case 1:
102: switch (version->subminor) {
103: case 0:
104: valid = PETSC_TRUE;
105: break;
106: }
107: break;
108: }
109: break;
110: case 3:
111: switch (version->minor) {
112: case 0:
113: switch (version->subminor) {
114: case 0:
115: valid = PETSC_TRUE;
116: break;
117: }
118: break;
119: }
120: break;
121: }
122: PetscCheck(valid, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "DMPlexStorageVersion %d.%d.%d not supported", version->major, version->minor, version->subminor);
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static inline PetscBool DMPlexStorageVersionGE(DMPlexStorageVersion version, int major, int minor, int subminor)
127: {
128: return (PetscBool)((version->major == major && version->minor == minor && version->subminor >= subminor) || (version->major == major && version->minor > minor) || (version->major > major));
129: }
131: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer viewer, DMPlexStorageVersion *version)
132: {
133: const char ATTR_NAME[] = "dmplex_storage_version";
134: PetscBool fileHasVersion;
135: char optVersion[16], fileVersion[16];
137: PetscFunctionBegin;
139: PetscAssertPointer(version, 2);
140: PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, version));
141: if (*version) PetscFunctionReturn(PETSC_SUCCESS);
143: PetscCall(PetscStrncpy(fileVersion, DMPLEX_STORAGE_VERSION_STABLE, sizeof(fileVersion)));
144: PetscCall(PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion));
145: if (fileHasVersion) {
146: char *tmp;
148: PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &tmp));
149: PetscCall(PetscStrncpy(fileVersion, tmp, sizeof(fileVersion)));
150: PetscCall(PetscFree(tmp));
151: }
152: PetscCall(PetscStrncpy(optVersion, fileVersion, sizeof(optVersion)));
153: PetscObjectOptionsBegin((PetscObject)viewer);
154: PetscCall(PetscOptionsString("-dm_plex_view_hdf5_storage_version", "DMPlex HDF5 viewer storage version", NULL, optVersion, optVersion, sizeof(optVersion), NULL));
155: PetscOptionsEnd();
156: if (!fileHasVersion) {
157: PetscCall(PetscViewerHDF5WriteAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, optVersion));
158: } else {
159: PetscBool flg;
161: PetscCall(PetscStrcmp(fileVersion, optVersion, &flg));
162: PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but file already has version %s - cannot mix versions", optVersion, fileVersion);
163: }
164: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT));
165: PetscCall(PetscViewerParseVersion_Private(viewer, optVersion, version));
166: PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, *version));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer viewer, DMPlexStorageVersion *version)
171: {
172: const char ATTR_NAME[] = "dmplex_storage_version";
173: char *defaultVersion;
174: char *versionString;
176: PetscFunctionBegin;
178: PetscAssertPointer(version, 2);
179: PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, version));
180: if (*version) PetscFunctionReturn(PETSC_SUCCESS);
182: //TODO string HDF5 attribute handling is terrible and should be redesigned
183: PetscCall(PetscStrallocpy(DMPLEX_STORAGE_VERSION_FIRST, &defaultVersion));
184: PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, &defaultVersion, &versionString));
185: PetscCall(PetscViewerParseVersion_Private(viewer, versionString, version));
186: PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, *version));
187: PetscCall(PetscFree(versionString));
188: PetscCall(PetscFree(defaultVersion));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode DMPlexGetHDF5Name_Private(DM dm, const char *name[])
193: {
194: PetscFunctionBegin;
195: if (((PetscObject)dm)->name) {
196: PetscCall(PetscObjectGetName((PetscObject)dm, name));
197: } else {
198: *name = "plex";
199: }
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
204: {
205: Vec stamp;
206: PetscMPIInt rank;
208: PetscFunctionBegin;
209: if (seqnum < 0) PetscFunctionReturn(PETSC_SUCCESS);
210: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
211: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp));
212: PetscCall(VecSetBlockSize(stamp, 1));
213: PetscCall(PetscObjectSetName((PetscObject)stamp, seqname));
214: if (rank == 0) {
215: PetscReal timeScale;
216: PetscBool istime;
218: PetscCall(PetscStrncmp(seqname, "time", 5, &istime));
219: if (istime) {
220: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale));
221: value *= timeScale;
222: }
223: PetscCall(VecSetValue(stamp, 0, value, INSERT_VALUES));
224: }
225: PetscCall(VecAssemblyBegin(stamp));
226: PetscCall(VecAssemblyEnd(stamp));
227: PetscCall(PetscViewerHDF5PushGroup(viewer, "/"));
228: PetscCall(PetscViewerHDF5PushTimestepping(viewer));
229: PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); /* seqnum < 0 jumps out above */
230: PetscCall(VecView(stamp, viewer));
231: PetscCall(PetscViewerHDF5PopTimestepping(viewer));
232: PetscCall(PetscViewerHDF5PopGroup(viewer));
233: PetscCall(VecDestroy(&stamp));
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
238: {
239: Vec stamp;
240: PetscMPIInt rank;
242: PetscFunctionBegin;
243: if (seqnum < 0) PetscFunctionReturn(PETSC_SUCCESS);
244: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
245: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp));
246: PetscCall(VecSetBlockSize(stamp, 1));
247: PetscCall(PetscObjectSetName((PetscObject)stamp, seqname));
248: PetscCall(PetscViewerHDF5PushGroup(viewer, "/"));
249: PetscCall(PetscViewerHDF5PushTimestepping(viewer));
250: PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); /* seqnum < 0 jumps out above */
251: PetscCall(VecLoad(stamp, viewer));
252: PetscCall(PetscViewerHDF5PopTimestepping(viewer));
253: PetscCall(PetscViewerHDF5PopGroup(viewer));
254: if (rank == 0) {
255: const PetscScalar *a;
256: PetscReal timeScale;
257: PetscBool istime;
259: PetscCall(VecGetArrayRead(stamp, &a));
260: *value = a[0];
261: PetscCall(VecRestoreArrayRead(stamp, &a));
262: PetscCall(PetscStrncmp(seqname, "time", 5, &istime));
263: if (istime) {
264: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale));
265: *value /= timeScale;
266: }
267: }
268: PetscCall(VecDestroy(&stamp));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
273: {
274: IS cutcells = NULL;
275: const PetscInt *cutc;
276: PetscInt cellHeight, vStart, vEnd, cStart, cEnd, c;
278: PetscFunctionBegin;
279: if (!cutLabel) PetscFunctionReturn(PETSC_SUCCESS);
280: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
281: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
282: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
283: /* Label vertices that should be duplicated */
284: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel));
285: PetscCall(DMLabelGetStratumIS(cutLabel, 2, &cutcells));
286: if (cutcells) {
287: PetscInt n;
289: PetscCall(ISGetIndices(cutcells, &cutc));
290: PetscCall(ISGetLocalSize(cutcells, &n));
291: for (c = 0; c < n; ++c) {
292: if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
293: PetscInt *closure = NULL;
294: PetscInt closureSize, cl, value;
296: PetscCall(DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure));
297: for (cl = 0; cl < closureSize * 2; cl += 2) {
298: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
299: PetscCall(DMLabelGetValue(cutLabel, closure[cl], &value));
300: if (value == 1) PetscCall(DMLabelSetValue(*cutVertexLabel, closure[cl], 1));
301: }
302: }
303: PetscCall(DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure));
304: }
305: }
306: PetscCall(ISRestoreIndices(cutcells, &cutc));
307: }
308: PetscCall(ISDestroy(&cutcells));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
313: {
314: DM dm;
315: DM dmBC;
316: PetscSection section, sectionGlobal;
317: Vec gv;
318: const char *name;
319: PetscViewerFormat format;
320: PetscInt seqnum;
321: PetscReal seqval;
322: PetscBool isseq;
324: PetscFunctionBegin;
325: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
326: PetscCall(VecGetDM(v, &dm));
327: PetscCall(DMGetLocalSection(dm, §ion));
328: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
329: PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar)seqval, viewer));
330: if (seqnum >= 0) {
331: PetscCall(PetscViewerHDF5PushTimestepping(viewer));
332: PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
333: }
334: PetscCall(PetscViewerGetFormat(viewer, &format));
335: PetscCall(DMGetOutputDM(dm, &dmBC));
336: PetscCall(DMGetGlobalSection(dmBC, §ionGlobal));
337: PetscCall(DMGetGlobalVector(dmBC, &gv));
338: PetscCall(PetscObjectGetName((PetscObject)v, &name));
339: PetscCall(PetscObjectSetName((PetscObject)gv, name));
340: PetscCall(DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv));
341: PetscCall(DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv));
342: PetscCall(PetscObjectTypeCompare((PetscObject)gv, VECSEQ, &isseq));
343: if (isseq) PetscCall(VecView_Seq(gv, viewer));
344: else PetscCall(VecView_MPI(gv, viewer));
345: if (format == PETSC_VIEWER_HDF5_VIZ) {
346: /* Output visualization representation */
347: PetscInt numFields, f;
348: DMLabel cutLabel, cutVertexLabel = NULL;
350: PetscCall(PetscSectionGetNumFields(section, &numFields));
351: PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
352: for (f = 0; f < numFields; ++f) {
353: Vec subv;
354: IS is;
355: const char *fname, *fgroup, *componentName;
356: char subname[PETSC_MAX_PATH_LEN];
357: PetscInt Nc, c, Nt, t;
358: PetscInt *pStart, *pEnd;
359: PetscViewerVTKFieldType *ft;
361: PetscCall(DMPlexGetFieldTypes_Internal(dm, section, f, &Nt, &pStart, &pEnd, &ft));
362: for (t = 0; t < Nt; ++t) {
363: fgroup = (ft[t] == PETSC_VTK_POINT_VECTOR_FIELD) || (ft[t] == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
364: PetscCall(PetscSectionGetFieldName(section, f, &fname));
365: if (!fname || ft[t] == PETSC_VTK_INVALID) continue;
367: if (!t) {
368: PetscCall(PetscViewerHDF5PushGroup(viewer, fgroup));
369: } else {
370: char group[PETSC_MAX_PATH_LEN];
372: PetscCall(PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "%s_%" PetscInt_FMT, fgroup, t));
373: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
374: }
376: if (cutLabel) {
377: const PetscScalar *ga;
378: PetscScalar *suba;
379: PetscInt gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;
381: PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
382: PetscCall(PetscSectionGetFieldComponents(section, f, &Nc));
383: for (p = pStart[t]; p < pEnd[t]; ++p) {
384: PetscInt gdof, fdof = 0, val;
386: PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
387: if (gdof > 0) PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
388: subSize += fdof;
389: PetscCall(DMLabelGetValue(cutVertexLabel, p, &val));
390: if (val == 1) extSize += fdof;
391: }
392: PetscCall(VecCreate(PetscObjectComm((PetscObject)gv), &subv));
393: PetscCall(VecSetSizes(subv, subSize + extSize, PETSC_DETERMINE));
394: PetscCall(VecSetBlockSize(subv, Nc));
395: PetscCall(VecSetType(subv, VECSTANDARD));
396: PetscCall(VecGetOwnershipRange(gv, &gstart, NULL));
397: PetscCall(VecGetArrayRead(gv, &ga));
398: PetscCall(VecGetArray(subv, &suba));
399: for (p = pStart[t]; p < pEnd[t]; ++p) {
400: PetscInt gdof, goff, val;
402: PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
403: if (gdof > 0) {
404: PetscInt fdof, fc, f2, poff = 0;
406: PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
407: /* Can get rid of this loop by storing field information in the global section */
408: for (f2 = 0; f2 < f; ++f2) {
409: PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof));
410: poff += fdof;
411: }
412: PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
413: for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff + poff + fc - gstart];
414: PetscCall(DMLabelGetValue(cutVertexLabel, p, &val));
415: if (val == 1) {
416: for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize + newOff] = ga[goff + poff + fc - gstart];
417: }
418: }
419: }
420: PetscCall(VecRestoreArrayRead(gv, &ga));
421: PetscCall(VecRestoreArray(subv, &suba));
422: PetscCall(DMLabelDestroy(&cutVertexLabel));
423: } else {
424: PetscCall(PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart[t], pEnd[t], &is, &subv));
425: }
426: PetscCall(PetscStrncpy(subname, name, sizeof(subname)));
427: PetscCall(PetscStrlcat(subname, "_", sizeof(subname)));
428: PetscCall(PetscStrlcat(subname, fname, sizeof(subname)));
429: PetscCall(PetscObjectSetName((PetscObject)subv, subname));
430: if (isseq) PetscCall(VecView_Seq(subv, viewer));
431: else PetscCall(VecView_MPI(subv, viewer));
432: if ((ft[t] == PETSC_VTK_POINT_VECTOR_FIELD) || (ft[t] == PETSC_VTK_CELL_VECTOR_FIELD)) {
433: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "vector"));
434: } else {
435: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "scalar"));
436: }
438: /* Output the component names in the field if available */
439: PetscCall(PetscSectionGetFieldComponents(section, f, &Nc));
440: for (c = 0; c < Nc; ++c) {
441: char componentNameLabel[PETSC_MAX_PATH_LEN];
442: PetscCall(PetscSectionGetComponentName(section, f, c, &componentName));
443: PetscCall(PetscSNPrintf(componentNameLabel, sizeof(componentNameLabel), "componentName%" PetscInt_FMT, c));
444: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, componentNameLabel, PETSC_STRING, componentName));
445: }
447: if (cutLabel) PetscCall(VecDestroy(&subv));
448: else PetscCall(PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart[t], pEnd[t], &is, &subv));
449: PetscCall(PetscViewerHDF5PopGroup(viewer));
450: }
451: PetscCall(DMPlexRestoreFieldTypes_Internal(dm, section, f, &Nt, &pStart, &pEnd, &ft));
452: }
453: }
454: if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
455: PetscCall(DMRestoreGlobalVector(dmBC, &gv));
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
460: {
461: DM dm;
462: Vec locv;
463: PetscObject isZero;
464: const char *name;
465: PetscReal time;
467: PetscFunctionBegin;
468: PetscCall(VecGetDM(v, &dm));
469: PetscCall(DMGetLocalVector(dm, &locv));
470: PetscCall(PetscObjectGetName((PetscObject)v, &name));
471: PetscCall(PetscObjectSetName((PetscObject)locv, name));
472: PetscCall(PetscObjectQuery((PetscObject)v, "__Vec_bc_zero__", &isZero));
473: PetscCall(PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", isZero));
474: PetscCall(DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv));
475: PetscCall(DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv));
476: PetscCall(DMGetOutputSequenceNumber(dm, NULL, &time));
477: PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL));
478: PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
479: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
480: PetscCall(VecView_Plex_Local_HDF5_Internal(locv, viewer));
481: PetscCall(PetscViewerPopFormat(viewer));
482: PetscCall(PetscViewerHDF5PopGroup(viewer));
483: PetscCall(PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", NULL));
484: PetscCall(DMRestoreLocalVector(dm, &locv));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
489: {
490: PetscBool isseq;
492: PetscFunctionBegin;
493: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
494: PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
495: if (isseq) PetscCall(VecView_Seq(v, viewer));
496: else PetscCall(VecView_MPI(v, viewer));
497: PetscCall(PetscViewerHDF5PopGroup(viewer));
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
502: {
503: DM dm;
504: Vec locv;
505: const char *name;
506: PetscInt seqnum;
508: PetscFunctionBegin;
509: PetscCall(VecGetDM(v, &dm));
510: PetscCall(DMGetLocalVector(dm, &locv));
511: PetscCall(PetscObjectGetName((PetscObject)v, &name));
512: PetscCall(PetscObjectSetName((PetscObject)locv, name));
513: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, NULL));
514: PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
515: if (seqnum >= 0) {
516: PetscCall(PetscViewerHDF5PushTimestepping(viewer));
517: PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
518: }
519: PetscCall(VecLoad_Plex_Local(locv, viewer));
520: if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
521: PetscCall(PetscViewerHDF5PopGroup(viewer));
522: PetscCall(DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v));
523: PetscCall(DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v));
524: PetscCall(DMRestoreLocalVector(dm, &locv));
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
529: {
530: DM dm;
531: PetscInt seqnum;
533: PetscFunctionBegin;
534: PetscCall(VecGetDM(v, &dm));
535: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, NULL));
536: PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
537: if (seqnum >= 0) {
538: PetscCall(PetscViewerHDF5PushTimestepping(viewer));
539: PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
540: }
541: PetscCall(VecLoad_Default(v, viewer));
542: if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
543: PetscCall(PetscViewerHDF5PopGroup(viewer));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: static PetscErrorCode DMPlexDistributionView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
548: {
549: MPI_Comm comm;
550: PetscMPIInt size, rank;
551: PetscInt size_petsc_int;
552: const char *topologydm_name, *distribution_name;
553: const PetscInt *gpoint;
554: PetscInt pStart, pEnd, p;
555: PetscSF pointSF;
556: PetscInt nroots, nleaves;
557: const PetscInt *ilocal;
558: const PetscSFNode *iremote;
559: IS chartSizesIS, ownersIS, gpointsIS;
560: PetscInt *chartSize, *owners, *gpoints;
562: PetscFunctionBegin;
563: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
564: PetscCallMPI(MPI_Comm_size(comm, &size));
565: PetscCallMPI(MPI_Comm_rank(comm, &rank));
566: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
567: PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
568: if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS);
569: PetscCall(PetscLogEventBegin(DMPLEX_DistributionView, viewer, 0, 0, 0));
570: size_petsc_int = (PetscInt)size;
571: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "comm_size", PETSC_INT, (void *)&size_petsc_int));
572: PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
573: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
574: PetscCall(PetscMalloc1(1, &chartSize));
575: *chartSize = pEnd - pStart;
576: PetscCall(PetscMalloc1(*chartSize, &owners));
577: PetscCall(PetscMalloc1(*chartSize, &gpoints));
578: PetscCall(DMGetPointSF(dm, &pointSF));
579: PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote));
580: for (p = pStart; p < pEnd; ++p) {
581: PetscInt gp = gpoint[p - pStart];
583: owners[p - pStart] = rank;
584: gpoints[p - pStart] = (gp < 0 ? -(gp + 1) : gp);
585: }
586: for (p = 0; p < nleaves; ++p) {
587: PetscInt ilocalp = (ilocal ? ilocal[p] : p);
589: owners[ilocalp] = iremote[p].rank;
590: }
591: PetscCall(ISCreateGeneral(comm, 1, chartSize, PETSC_OWN_POINTER, &chartSizesIS));
592: PetscCall(ISCreateGeneral(comm, *chartSize, owners, PETSC_OWN_POINTER, &ownersIS));
593: PetscCall(ISCreateGeneral(comm, *chartSize, gpoints, PETSC_OWN_POINTER, &gpointsIS));
594: PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes"));
595: PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners"));
596: PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers"));
597: PetscCall(ISView(chartSizesIS, viewer));
598: PetscCall(ISView(ownersIS, viewer));
599: PetscCall(ISView(gpointsIS, viewer));
600: PetscCall(ISDestroy(&chartSizesIS));
601: PetscCall(ISDestroy(&ownersIS));
602: PetscCall(ISDestroy(&gpointsIS));
603: PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
604: PetscCall(PetscLogEventEnd(DMPLEX_DistributionView, viewer, 0, 0, 0));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: static PetscErrorCode DMPlexTopologyView_HDF5_Inner_Private(DM dm, IS globalPointNumbers, PetscViewer viewer, PetscInt pStart, PetscInt pEnd, const char pointsName[], const char coneSizesName[], const char conesName[], const char orientationsName[])
609: {
610: IS coneSizesIS, conesIS, orientationsIS;
611: PetscInt *coneSizes, *cones, *orientations;
612: const PetscInt *gpoint;
613: PetscInt nPoints = 0, conesSize = 0;
614: PetscInt p, c, s;
615: MPI_Comm comm;
617: PetscFunctionBegin;
618: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
619: PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
620: for (p = pStart; p < pEnd; ++p) {
621: if (gpoint[p] >= 0) {
622: PetscInt coneSize;
624: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
625: nPoints += 1;
626: conesSize += coneSize;
627: }
628: }
629: PetscCall(PetscMalloc1(nPoints, &coneSizes));
630: PetscCall(PetscMalloc1(conesSize, &cones));
631: PetscCall(PetscMalloc1(conesSize, &orientations));
632: for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
633: if (gpoint[p] >= 0) {
634: const PetscInt *cone, *ornt;
635: PetscInt coneSize, cp;
637: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
638: PetscCall(DMPlexGetCone(dm, p, &cone));
639: PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
640: coneSizes[s] = coneSize;
641: for (cp = 0; cp < coneSize; ++cp, ++c) {
642: cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]] + 1) : gpoint[cone[cp]];
643: orientations[c] = ornt[cp];
644: }
645: ++s;
646: }
647: }
648: PetscCheck(s == nPoints, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %" PetscInt_FMT " != %" PetscInt_FMT, s, nPoints);
649: PetscCheck(c == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %" PetscInt_FMT " != %" PetscInt_FMT, c, conesSize);
650: PetscCall(ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS));
651: PetscCall(ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS));
652: PetscCall(ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS));
653: PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
654: PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
655: PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
656: PetscCall(ISView(coneSizesIS, viewer));
657: PetscCall(ISView(conesIS, viewer));
658: PetscCall(ISView(orientationsIS, viewer));
659: PetscCall(ISDestroy(&coneSizesIS));
660: PetscCall(ISDestroy(&conesIS));
661: PetscCall(ISDestroy(&orientationsIS));
662: if (pointsName) {
663: IS pointsIS;
664: PetscInt *points;
666: PetscCall(PetscMalloc1(nPoints, &points));
667: for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
668: if (gpoint[p] >= 0) {
669: points[s] = gpoint[p];
670: ++s;
671: }
672: }
673: PetscCall(ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS));
674: PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName));
675: PetscCall(ISView(pointsIS, viewer));
676: PetscCall(ISDestroy(&pointsIS));
677: }
678: PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode DMPlexTopologyView_HDF5_Legacy_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
683: {
684: const char *pointsName, *coneSizesName, *conesName, *orientationsName;
685: PetscInt pStart, pEnd, dim;
687: PetscFunctionBegin;
688: pointsName = "order";
689: coneSizesName = "cones";
690: conesName = "cells";
691: orientationsName = "orientation";
692: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
693: PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers, viewer, pStart, pEnd, pointsName, coneSizesName, conesName, orientationsName));
694: PetscCall(DMGetDimension(dm, &dim));
695: PetscCall(PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *)&dim));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: //TODO get this numbering right away without needing this function
700: /* Renumber global point numbers so that they are 0-based per stratum */
701: static PetscErrorCode RenumberGlobalPointNumbersPerStratum_Private(DM dm, IS globalPointNumbers, IS *newGlobalPointNumbers, IS *strataPermutation)
702: {
703: PetscInt d, depth, p, n;
704: PetscInt *offsets;
705: const PetscInt *gpn;
706: PetscInt *ngpn;
707: MPI_Comm comm;
709: PetscFunctionBegin;
710: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
711: PetscCall(ISGetLocalSize(globalPointNumbers, &n));
712: PetscCall(ISGetIndices(globalPointNumbers, &gpn));
713: PetscCall(PetscMalloc1(n, &ngpn));
714: PetscCall(DMPlexGetDepth(dm, &depth));
715: PetscCall(PetscMalloc1(depth + 1, &offsets));
716: for (d = 0; d <= depth; d++) {
717: PetscInt pStart, pEnd;
719: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
720: offsets[d] = PETSC_MAX_INT;
721: for (p = pStart; p < pEnd; p++) {
722: if (gpn[p] >= 0 && gpn[p] < offsets[d]) offsets[d] = gpn[p];
723: }
724: }
725: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, offsets, depth + 1, MPIU_INT, MPI_MIN, comm));
726: for (d = 0; d <= depth; d++) {
727: PetscInt pStart, pEnd;
729: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
730: for (p = pStart; p < pEnd; p++) ngpn[p] = gpn[p] - PetscSign(gpn[p]) * offsets[d];
731: }
732: PetscCall(ISRestoreIndices(globalPointNumbers, &gpn));
733: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)globalPointNumbers), n, ngpn, PETSC_OWN_POINTER, newGlobalPointNumbers));
734: {
735: PetscInt *perm;
737: PetscCall(PetscMalloc1(depth + 1, &perm));
738: for (d = 0; d <= depth; d++) perm[d] = d;
739: PetscCall(PetscSortIntWithPermutation(depth + 1, offsets, perm));
740: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, depth + 1, perm, PETSC_OWN_POINTER, strataPermutation));
741: }
742: PetscCall(PetscFree(offsets));
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static PetscErrorCode DMPlexTopologyView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
747: {
748: IS globalPointNumbers0, strataPermutation;
749: const char *coneSizesName, *conesName, *orientationsName;
750: PetscInt depth, d;
751: MPI_Comm comm;
753: PetscFunctionBegin;
754: coneSizesName = "cone_sizes";
755: conesName = "cones";
756: orientationsName = "orientations";
757: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
758: PetscCall(DMPlexGetDepth(dm, &depth));
759: {
760: PetscInt dim;
762: PetscCall(DMGetDimension(dm, &dim));
763: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "cell_dim", PETSC_INT, &dim));
764: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "depth", PETSC_INT, &depth));
765: }
767: PetscCall(RenumberGlobalPointNumbersPerStratum_Private(dm, globalPointNumbers, &globalPointNumbers0, &strataPermutation));
768: /* TODO dirty trick to save serial IS using the same parallel viewer */
769: {
770: IS spOnComm;
771: PetscInt n = 0, N;
772: const PetscInt *idx = NULL;
773: const PetscInt *old;
774: PetscMPIInt rank;
776: PetscCallMPI(MPI_Comm_rank(comm, &rank));
777: PetscCall(ISGetLocalSize(strataPermutation, &N));
778: PetscCall(ISGetIndices(strataPermutation, &old));
779: if (!rank) {
780: n = N;
781: idx = old;
782: }
783: PetscCall(ISCreateGeneral(comm, n, idx, PETSC_COPY_VALUES, &spOnComm));
784: PetscCall(ISRestoreIndices(strataPermutation, &old));
785: PetscCall(ISDestroy(&strataPermutation));
786: strataPermutation = spOnComm;
787: }
788: PetscCall(PetscObjectSetName((PetscObject)strataPermutation, "permutation"));
789: PetscCall(ISView(strataPermutation, viewer));
790: PetscCall(PetscViewerHDF5PushGroup(viewer, "strata"));
791: for (d = 0; d <= depth; d++) {
792: PetscInt pStart, pEnd;
793: char group[128];
795: PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, d));
796: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
797: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
798: PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers0, viewer, pStart, pEnd, NULL, coneSizesName, conesName, orientationsName));
799: PetscCall(PetscViewerHDF5PopGroup(viewer));
800: }
801: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */
802: PetscCall(ISDestroy(&globalPointNumbers0));
803: PetscCall(ISDestroy(&strataPermutation));
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
808: {
809: DMPlexStorageVersion version;
810: const char *topologydm_name;
811: char group[PETSC_MAX_PATH_LEN];
813: PetscFunctionBegin;
814: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
815: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
816: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
817: PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name));
818: } else {
819: PetscCall(PetscStrncpy(group, "/", sizeof(group)));
820: }
821: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
823: PetscCall(PetscViewerHDF5PushGroup(viewer, "topology"));
824: if (version->major < 3) {
825: PetscCall(DMPlexTopologyView_HDF5_Legacy_Private(dm, globalPointNumbers, viewer));
826: } else {
827: /* since DMPlexStorageVersion 3.0.0 */
828: PetscCall(DMPlexTopologyView_HDF5_Private(dm, globalPointNumbers, viewer));
829: }
830: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */
832: if (DMPlexStorageVersionGE(version, 2, 1, 0)) {
833: const char *distribution_name;
835: PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
836: PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions"));
837: PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
838: PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name));
839: PetscCall(DMPlexDistributionView_HDF5_Private(dm, globalPointNumbers, viewer));
840: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */
841: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */
842: }
844: PetscCall(PetscViewerHDF5PopGroup(viewer));
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
849: {
850: PetscSF sfPoint;
851: DMLabel cutLabel, cutVertexLabel = NULL;
852: IS globalVertexNumbers, cutvertices = NULL;
853: const PetscInt *gcell, *gvertex, *cutverts = NULL;
854: PetscInt *vertices;
855: PetscInt conesSize = 0;
856: PetscInt dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
858: PetscFunctionBegin;
859: *numCorners = 0;
860: PetscCall(DMGetDimension(dm, &dim));
861: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
862: PetscCall(ISGetIndices(globalCellNumbers, &gcell));
864: for (cell = cStart; cell < cEnd; ++cell) {
865: PetscInt *closure = NULL;
866: PetscInt closureSize, v, Nc = 0;
868: if (gcell[cell] < 0) continue;
869: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
870: for (v = 0; v < closureSize * 2; v += 2) {
871: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
872: }
873: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
874: conesSize += Nc;
875: if (!numCornersLocal) numCornersLocal = Nc;
876: else if (numCornersLocal != Nc) numCornersLocal = 1;
877: }
878: PetscCall(MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
879: PetscCheck(!numCornersLocal || !(numCornersLocal != *numCorners || *numCorners == 0), PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
880: /* Handle periodic cuts by identifying vertices which should be duplicated */
881: PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
882: PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
883: if (cutVertexLabel) PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices));
884: if (cutvertices) {
885: PetscCall(ISGetIndices(cutvertices, &cutverts));
886: PetscCall(ISGetLocalSize(cutvertices, &vExtra));
887: }
888: PetscCall(DMGetPointSF(dm, &sfPoint));
889: if (cutLabel) {
890: const PetscInt *ilocal;
891: const PetscSFNode *iremote;
892: PetscInt nroots, nleaves;
894: PetscCall(PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote));
895: if (nleaves < 0) {
896: PetscCall(PetscObjectReference((PetscObject)sfPoint));
897: } else {
898: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), &sfPoint));
899: PetscCall(PetscSFSetGraph(sfPoint, nroots + vExtra, nleaves, (PetscInt *)ilocal, PETSC_COPY_VALUES, (PetscSFNode *)iremote, PETSC_COPY_VALUES));
900: }
901: } else {
902: PetscCall(PetscObjectReference((PetscObject)sfPoint));
903: }
904: /* Number all vertices */
905: PetscCall(DMPlexCreateNumbering_Plex(dm, vStart, vEnd + vExtra, 0, NULL, sfPoint, &globalVertexNumbers));
906: PetscCall(PetscSFDestroy(&sfPoint));
907: /* Create cones */
908: PetscCall(ISGetIndices(globalVertexNumbers, &gvertex));
909: PetscCall(PetscMalloc1(conesSize, &vertices));
910: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
911: PetscInt *closure = NULL;
912: PetscInt closureSize, Nc = 0, p, value = -1;
913: PetscBool replace;
915: if (gcell[cell] < 0) continue;
916: if (cutLabel) PetscCall(DMLabelGetValue(cutLabel, cell, &value));
917: replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
918: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
919: for (p = 0; p < closureSize * 2; p += 2) {
920: if ((closure[p] >= vStart) && (closure[p] < vEnd)) closure[Nc++] = closure[p];
921: }
922: PetscCall(DMPlexReorderCell(dm, cell, closure));
923: for (p = 0; p < Nc; ++p) {
924: PetscInt nv, gv = gvertex[closure[p] - vStart];
926: if (replace) {
927: PetscCall(PetscFindInt(closure[p], vExtra, cutverts, &nv));
928: if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
929: }
930: vertices[v++] = gv < 0 ? -(gv + 1) : gv;
931: }
932: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
933: }
934: PetscCall(ISRestoreIndices(globalVertexNumbers, &gvertex));
935: PetscCall(ISDestroy(&globalVertexNumbers));
936: PetscCall(ISRestoreIndices(globalCellNumbers, &gcell));
937: if (cutvertices) PetscCall(ISRestoreIndices(cutvertices, &cutverts));
938: PetscCall(ISDestroy(&cutvertices));
939: PetscCall(DMLabelDestroy(&cutVertexLabel));
940: PetscCheck(v == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %" PetscInt_FMT " != %" PetscInt_FMT, v, conesSize);
941: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS));
942: PetscCall(PetscLayoutSetBlockSize((*cellIS)->map, *numCorners));
943: PetscCall(PetscObjectSetName((PetscObject)*cellIS, "cells"));
944: PetscFunctionReturn(PETSC_SUCCESS);
945: }
947: static PetscErrorCode DMPlexTopologyView_HDF5_XDMF_Private(DM dm, IS globalCellNumbers, PetscViewer viewer)
948: {
949: DM cdm;
950: DMLabel depthLabel, ctLabel;
951: IS cellIS;
952: PetscInt dim, depth, cellHeight, c, n = 0;
954: PetscFunctionBegin;
955: PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz"));
956: PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
958: PetscCall(PetscViewerHDF5PopGroup(viewer));
959: PetscCall(DMGetDimension(dm, &dim));
960: PetscCall(DMPlexGetDepth(dm, &depth));
961: PetscCall(DMGetCoordinateDM(dm, &cdm));
962: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
963: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
964: PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel));
965: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
966: const DMPolytopeType ict = (DMPolytopeType)c;
967: PetscInt pStart, pEnd, dep, numCorners;
968: PetscBool output = PETSC_FALSE, doOutput;
970: if (ict == DM_POLYTOPE_FV_GHOST) continue;
971: PetscCall(DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd));
972: if (pStart >= 0) {
973: PetscCall(DMLabelGetValue(depthLabel, pStart, &dep));
974: if (dep == depth - cellHeight) output = PETSC_TRUE;
975: }
976: PetscCall(MPIU_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
977: if (!doOutput) continue;
978: PetscCall(CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS));
979: if (!n) {
980: PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/topology"));
981: } else {
982: char group[PETSC_MAX_PATH_LEN];
984: PetscCall(PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%" PetscInt_FMT, n));
985: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
986: }
987: PetscCall(ISView(cellIS, viewer));
988: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_corners", PETSC_INT, (void *)&numCorners));
989: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_dim", PETSC_INT, (void *)&dim));
990: PetscCall(ISDestroy(&cellIS));
991: PetscCall(PetscViewerHDF5PopGroup(viewer));
992: ++n;
993: }
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
998: {
999: DM cdm;
1000: Vec coordinates, newcoords;
1001: PetscReal lengthScale;
1002: PetscInt m, M, bs;
1004: PetscFunctionBegin;
1005: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
1006: PetscCall(DMGetCoordinateDM(dm, &cdm));
1007: PetscCall(DMGetCoordinates(dm, &coordinates));
1008: PetscCall(VecCreate(PetscObjectComm((PetscObject)coordinates), &newcoords));
1009: PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices"));
1010: PetscCall(VecGetSize(coordinates, &M));
1011: PetscCall(VecGetLocalSize(coordinates, &m));
1012: PetscCall(VecSetSizes(newcoords, m, M));
1013: PetscCall(VecGetBlockSize(coordinates, &bs));
1014: PetscCall(VecSetBlockSize(newcoords, bs));
1015: PetscCall(VecSetType(newcoords, VECSTANDARD));
1016: PetscCall(VecCopy(coordinates, newcoords));
1017: PetscCall(VecScale(newcoords, lengthScale));
1018: /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
1019: PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry"));
1020: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
1021: PetscCall(VecView(newcoords, viewer));
1022: PetscCall(PetscViewerPopFormat(viewer));
1023: PetscCall(PetscViewerHDF5PopGroup(viewer));
1024: PetscCall(VecDestroy(&newcoords));
1025: PetscFunctionReturn(PETSC_SUCCESS);
1026: }
1028: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
1029: {
1030: DM cdm;
1031: Vec coords, newcoords;
1032: PetscInt m, M, bs;
1033: PetscReal lengthScale;
1034: const char *topologydm_name, *coordinatedm_name, *coordinates_name;
1036: PetscFunctionBegin;
1037: {
1038: PetscViewerFormat format;
1039: DMPlexStorageVersion version;
1041: PetscCall(PetscViewerGetFormat(viewer, &format));
1042: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
1043: if (!DMPlexStorageVersionGE(version, 2, 0, 0) || format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ) {
1044: PetscCall(DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer));
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1047: }
1048: /* since 2.0.0 */
1049: PetscCall(DMGetCoordinateDM(dm, &cdm));
1050: PetscCall(DMGetCoordinates(dm, &coords));
1051: PetscCall(PetscObjectGetName((PetscObject)cdm, &coordinatedm_name));
1052: PetscCall(PetscObjectGetName((PetscObject)coords, &coordinates_name));
1053: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1054: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1055: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1056: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name));
1057: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name));
1058: PetscCall(PetscViewerHDF5PopGroup(viewer));
1059: PetscCall(PetscViewerHDF5PopGroup(viewer));
1060: PetscCall(DMPlexSectionView(dm, viewer, cdm));
1061: PetscCall(VecCreate(PetscObjectComm((PetscObject)coords), &newcoords));
1062: PetscCall(PetscObjectSetName((PetscObject)newcoords, coordinates_name));
1063: PetscCall(VecGetSize(coords, &M));
1064: PetscCall(VecGetLocalSize(coords, &m));
1065: PetscCall(VecSetSizes(newcoords, m, M));
1066: PetscCall(VecGetBlockSize(coords, &bs));
1067: PetscCall(VecSetBlockSize(newcoords, bs));
1068: PetscCall(VecSetType(newcoords, VECSTANDARD));
1069: PetscCall(VecCopy(coords, newcoords));
1070: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
1071: PetscCall(VecScale(newcoords, lengthScale));
1072: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
1073: PetscCall(DMPlexGlobalVectorView(dm, viewer, cdm, newcoords));
1074: PetscCall(PetscViewerPopFormat(viewer));
1075: PetscCall(VecDestroy(&newcoords));
1076: PetscFunctionReturn(PETSC_SUCCESS);
1077: }
1079: static PetscErrorCode DMPlexCoordinatesView_HDF5_XDMF_Private(DM dm, PetscViewer viewer)
1080: {
1081: DM cdm;
1082: Vec coordinatesLocal, newcoords;
1083: PetscSection cSection, cGlobalSection;
1084: PetscScalar *coords, *ncoords;
1085: DMLabel cutLabel, cutVertexLabel = NULL;
1086: const PetscReal *L;
1087: PetscReal lengthScale;
1088: PetscInt vStart, vEnd, v, bs, N, coordSize, dof, off, d;
1089: PetscBool localized, embedded;
1091: PetscFunctionBegin;
1092: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1093: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
1094: PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal));
1095: PetscCall(VecGetBlockSize(coordinatesLocal, &bs));
1096: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
1097: if (localized == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1098: PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L));
1099: PetscCall(DMGetCoordinateDM(dm, &cdm));
1100: PetscCall(DMGetLocalSection(cdm, &cSection));
1101: PetscCall(DMGetGlobalSection(cdm, &cGlobalSection));
1102: PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
1103: N = 0;
1105: PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
1106: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &newcoords));
1107: PetscCall(PetscSectionGetDof(cSection, vStart, &dof));
1108: PetscCall(PetscPrintf(PETSC_COMM_SELF, "DOF: %" PetscInt_FMT "\n", dof));
1109: embedded = (PetscBool)(L && dof == 2 && !cutLabel);
1110: if (cutVertexLabel) {
1111: PetscCall(DMLabelGetStratumSize(cutVertexLabel, 1, &v));
1112: N += dof * v;
1113: }
1114: for (v = vStart; v < vEnd; ++v) {
1115: PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof));
1116: if (dof < 0) continue;
1117: if (embedded) N += dof + 1;
1118: else N += dof;
1119: }
1120: if (embedded) PetscCall(VecSetBlockSize(newcoords, bs + 1));
1121: else PetscCall(VecSetBlockSize(newcoords, bs));
1122: PetscCall(VecSetSizes(newcoords, N, PETSC_DETERMINE));
1123: PetscCall(VecSetType(newcoords, VECSTANDARD));
1124: PetscCall(VecGetArray(coordinatesLocal, &coords));
1125: PetscCall(VecGetArray(newcoords, &ncoords));
1126: coordSize = 0;
1127: for (v = vStart; v < vEnd; ++v) {
1128: PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof));
1129: PetscCall(PetscSectionGetOffset(cSection, v, &off));
1130: if (dof < 0) continue;
1131: if (embedded) {
1132: if (L && (L[0] > 0.0) && (L[1] > 0.0)) {
1133: PetscReal theta, phi, r, R;
1134: /* XY-periodic */
1135: /* Suppose its an y-z circle, then
1136: \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
1137: and the circle in that plane is
1138: \hat r cos(phi) + \hat x sin(phi) */
1139: theta = 2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1];
1140: phi = 2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0];
1141: r = L[0] / (2.0 * PETSC_PI * 2.0 * L[1]);
1142: R = L[1] / (2.0 * PETSC_PI);
1143: ncoords[coordSize++] = PetscSinReal(phi) * r;
1144: ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
1145: ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
1146: } else if (L && (L[0] > 0.0)) {
1147: /* X-periodic */
1148: ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
1149: ncoords[coordSize++] = coords[off + 1];
1150: ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
1151: } else if (L && (L[1] > 0.0)) {
1152: /* Y-periodic */
1153: ncoords[coordSize++] = coords[off + 0];
1154: ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
1155: ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
1156: #if 0
1157: } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
1158: PetscReal phi, r, R;
1159: /* Mobius strip */
1160: /* Suppose its an x-z circle, then
1161: \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
1162: and in that plane we rotate by pi as we go around the circle
1163: \hat r cos(phi/2) + \hat y sin(phi/2) */
1164: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
1165: R = L[0];
1166: r = PetscRealPart(coords[off+1]) - L[1]/2.0;
1167: ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
1168: ncoords[coordSize++] = PetscSinReal(phi/2.0) * r;
1169: ncoords[coordSize++] = PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
1170: #endif
1171: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
1172: } else {
1173: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off + d];
1174: }
1175: }
1176: if (cutVertexLabel) {
1177: IS vertices;
1178: const PetscInt *verts;
1179: PetscInt n;
1181: PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &vertices));
1182: if (vertices) {
1183: PetscCall(ISGetIndices(vertices, &verts));
1184: PetscCall(ISGetLocalSize(vertices, &n));
1185: for (v = 0; v < n; ++v) {
1186: PetscCall(PetscSectionGetDof(cSection, verts[v], &dof));
1187: PetscCall(PetscSectionGetOffset(cSection, verts[v], &off));
1188: for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off + d] + ((L[d] > 0.) ? L[d] : 0.0);
1189: }
1190: PetscCall(ISRestoreIndices(vertices, &verts));
1191: PetscCall(ISDestroy(&vertices));
1192: }
1193: }
1194: PetscCheck(coordSize == N, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %" PetscInt_FMT " != %" PetscInt_FMT, coordSize, N);
1195: PetscCall(DMLabelDestroy(&cutVertexLabel));
1196: PetscCall(VecRestoreArray(coordinatesLocal, &coords));
1197: PetscCall(VecRestoreArray(newcoords, &ncoords));
1198: PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices"));
1199: PetscCall(VecScale(newcoords, lengthScale));
1200: PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz"));
1201: PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
1202: PetscCall(PetscViewerHDF5PopGroup(viewer));
1203: PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/geometry"));
1204: PetscCall(VecView(newcoords, viewer));
1205: PetscCall(PetscViewerHDF5PopGroup(viewer));
1206: PetscCall(VecDestroy(&newcoords));
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
1211: {
1212: const char *topologydm_name;
1213: const PetscInt *gpoint;
1214: PetscInt numLabels, l;
1215: DMPlexStorageVersion version;
1216: char group[PETSC_MAX_PATH_LEN];
1218: PetscFunctionBegin;
1219: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
1220: PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
1221: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1222: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1223: PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name));
1224: } else {
1225: PetscCall(PetscStrncpy(group, "/labels", sizeof(group)));
1226: }
1227: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1228: PetscCall(DMGetNumLabels(dm, &numLabels));
1229: for (l = 0; l < numLabels; ++l) {
1230: DMLabel label;
1231: const char *name;
1232: IS valueIS, pvalueIS, globalValueIS;
1233: const PetscInt *values;
1234: PetscInt numValues, v;
1235: PetscBool isDepth, output;
1237: PetscCall(DMGetLabelByNum(dm, l, &label));
1238: PetscCall(PetscObjectGetName((PetscObject)label, &name));
1239: PetscCall(DMGetLabelOutput(dm, name, &output));
1240: PetscCall(PetscStrncmp(name, "depth", 10, &isDepth));
1241: if (isDepth || !output) continue;
1242: PetscCall(PetscViewerHDF5PushGroup(viewer, name));
1243: PetscCall(DMLabelGetValueIS(label, &valueIS));
1244: /* Must copy to a new IS on the global comm */
1245: PetscCall(ISGetLocalSize(valueIS, &numValues));
1246: PetscCall(ISGetIndices(valueIS, &values));
1247: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS));
1248: PetscCall(ISRestoreIndices(valueIS, &values));
1249: PetscCall(ISAllGather(pvalueIS, &globalValueIS));
1250: PetscCall(ISDestroy(&pvalueIS));
1251: PetscCall(ISSortRemoveDups(globalValueIS));
1252: PetscCall(ISGetLocalSize(globalValueIS, &numValues));
1253: PetscCall(ISGetIndices(globalValueIS, &values));
1254: for (v = 0; v < numValues; ++v) {
1255: IS stratumIS, globalStratumIS;
1256: const PetscInt *spoints = NULL;
1257: PetscInt *gspoints, n = 0, gn, p;
1258: const char *iname = "indices";
1259: char group[PETSC_MAX_PATH_LEN];
1261: PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, values[v]));
1262: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1263: PetscCall(DMLabelGetStratumIS(label, values[v], &stratumIS));
1265: if (stratumIS) PetscCall(ISGetLocalSize(stratumIS, &n));
1266: if (stratumIS) PetscCall(ISGetIndices(stratumIS, &spoints));
1267: for (gn = 0, p = 0; p < n; ++p)
1268: if (gpoint[spoints[p]] >= 0) ++gn;
1269: PetscCall(PetscMalloc1(gn, &gspoints));
1270: for (gn = 0, p = 0; p < n; ++p)
1271: if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
1272: if (stratumIS) PetscCall(ISRestoreIndices(stratumIS, &spoints));
1273: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS));
1274: PetscCall(PetscObjectSetName((PetscObject)globalStratumIS, iname));
1276: PetscCall(ISView(globalStratumIS, viewer));
1277: PetscCall(ISDestroy(&globalStratumIS));
1278: PetscCall(ISDestroy(&stratumIS));
1279: PetscCall(PetscViewerHDF5PopGroup(viewer));
1280: }
1281: PetscCall(ISRestoreIndices(globalValueIS, &values));
1282: PetscCall(ISDestroy(&globalValueIS));
1283: PetscCall(ISDestroy(&valueIS));
1284: PetscCall(PetscViewerHDF5PopGroup(viewer));
1285: }
1286: PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
1287: PetscCall(PetscViewerHDF5PopGroup(viewer));
1288: PetscFunctionReturn(PETSC_SUCCESS);
1289: }
1291: /* We only write cells and vertices. Does this screw up parallel reading? */
1292: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
1293: {
1294: IS globalPointNumbers;
1295: PetscViewerFormat format;
1296: PetscBool viz_geom = PETSC_FALSE, xdmf_topo = PETSC_FALSE, petsc_topo = PETSC_FALSE;
1298: PetscFunctionBegin;
1299: PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers));
1300: PetscCall(DMPlexCoordinatesView_HDF5_Internal(dm, viewer));
1302: PetscCall(PetscViewerGetFormat(viewer, &format));
1303: switch (format) {
1304: case PETSC_VIEWER_HDF5_VIZ:
1305: viz_geom = PETSC_TRUE;
1306: xdmf_topo = PETSC_TRUE;
1307: break;
1308: case PETSC_VIEWER_HDF5_XDMF:
1309: xdmf_topo = PETSC_TRUE;
1310: break;
1311: case PETSC_VIEWER_HDF5_PETSC:
1312: petsc_topo = PETSC_TRUE;
1313: break;
1314: case PETSC_VIEWER_DEFAULT:
1315: case PETSC_VIEWER_NATIVE:
1316: viz_geom = PETSC_TRUE;
1317: xdmf_topo = PETSC_TRUE;
1318: petsc_topo = PETSC_TRUE;
1319: break;
1320: default:
1321: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
1322: }
1324: if (viz_geom) PetscCall(DMPlexCoordinatesView_HDF5_XDMF_Private(dm, viewer));
1325: if (xdmf_topo) PetscCall(DMPlexTopologyView_HDF5_XDMF_Private(dm, globalPointNumbers, viewer));
1326: if (petsc_topo) {
1327: PetscCall(DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer));
1328: PetscCall(DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer));
1329: }
1331: PetscCall(ISDestroy(&globalPointNumbers));
1332: PetscFunctionReturn(PETSC_SUCCESS);
1333: }
1335: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
1336: {
1337: MPI_Comm comm;
1338: const char *topologydm_name;
1339: const char *sectiondm_name;
1340: PetscSection gsection;
1342: PetscFunctionBegin;
1343: PetscCall(PetscObjectGetComm((PetscObject)sectiondm, &comm));
1344: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1345: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
1346: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1347: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1348: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1349: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1350: PetscCall(DMGetGlobalSection(sectiondm, &gsection));
1351: /* Save raw section */
1352: PetscCall(PetscSectionView(gsection, viewer));
1353: /* Save plex wrapper */
1354: {
1355: PetscInt pStart, pEnd, p, n;
1356: IS globalPointNumbers;
1357: const PetscInt *gpoints;
1358: IS orderIS;
1359: PetscInt *order;
1361: PetscCall(PetscSectionGetChart(gsection, &pStart, &pEnd));
1362: PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers));
1363: PetscCall(ISGetIndices(globalPointNumbers, &gpoints));
1364: for (p = pStart, n = 0; p < pEnd; ++p)
1365: if (gpoints[p] >= 0) n++;
1366: /* "order" is an array of global point numbers.
1367: When loading, it is used with topology/order array
1368: to match section points with plex topology points. */
1369: PetscCall(PetscMalloc1(n, &order));
1370: for (p = pStart, n = 0; p < pEnd; ++p)
1371: if (gpoints[p] >= 0) order[n++] = gpoints[p];
1372: PetscCall(ISRestoreIndices(globalPointNumbers, &gpoints));
1373: PetscCall(ISDestroy(&globalPointNumbers));
1374: PetscCall(ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS));
1375: PetscCall(PetscObjectSetName((PetscObject)orderIS, "order"));
1376: PetscCall(ISView(orderIS, viewer));
1377: PetscCall(ISDestroy(&orderIS));
1378: }
1379: PetscCall(PetscViewerHDF5PopGroup(viewer));
1380: PetscCall(PetscViewerHDF5PopGroup(viewer));
1381: PetscCall(PetscViewerHDF5PopGroup(viewer));
1382: PetscCall(PetscViewerHDF5PopGroup(viewer));
1383: PetscFunctionReturn(PETSC_SUCCESS);
1384: }
1386: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1387: {
1388: const char *topologydm_name;
1389: const char *sectiondm_name;
1390: const char *vec_name;
1391: PetscInt bs;
1393: PetscFunctionBegin;
1394: /* Check consistency */
1395: {
1396: PetscSF pointsf, pointsf1;
1398: PetscCall(DMGetPointSF(dm, &pointsf));
1399: PetscCall(DMGetPointSF(sectiondm, &pointsf1));
1400: PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
1401: }
1402: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1403: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
1404: PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
1405: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1406: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1407: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1408: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1409: PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
1410: PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
1411: PetscCall(VecGetBlockSize(vec, &bs));
1412: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs));
1413: PetscCall(VecSetBlockSize(vec, 1));
1414: /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but, */
1415: /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view */
1416: /* is set to VecView_Plex, which would save vec in a predefined location. */
1417: /* To save vec in where we want, we create a new Vec (temp) with */
1418: /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
1419: {
1420: Vec temp;
1421: const PetscScalar *array;
1422: PetscLayout map;
1424: PetscCall(VecCreate(PetscObjectComm((PetscObject)vec), &temp));
1425: PetscCall(PetscObjectSetName((PetscObject)temp, vec_name));
1426: PetscCall(VecGetLayout(vec, &map));
1427: PetscCall(VecSetLayout(temp, map));
1428: PetscCall(VecSetUp(temp));
1429: PetscCall(VecGetArrayRead(vec, &array));
1430: PetscCall(VecPlaceArray(temp, array));
1431: PetscCall(VecView(temp, viewer));
1432: PetscCall(VecResetArray(temp));
1433: PetscCall(VecRestoreArrayRead(vec, &array));
1434: PetscCall(VecDestroy(&temp));
1435: }
1436: PetscCall(VecSetBlockSize(vec, bs));
1437: PetscCall(PetscViewerHDF5PopGroup(viewer));
1438: PetscCall(PetscViewerHDF5PopGroup(viewer));
1439: PetscCall(PetscViewerHDF5PopGroup(viewer));
1440: PetscCall(PetscViewerHDF5PopGroup(viewer));
1441: PetscCall(PetscViewerHDF5PopGroup(viewer));
1442: PetscCall(PetscViewerHDF5PopGroup(viewer));
1443: PetscFunctionReturn(PETSC_SUCCESS);
1444: }
1446: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1447: {
1448: MPI_Comm comm;
1449: const char *topologydm_name;
1450: const char *sectiondm_name;
1451: const char *vec_name;
1452: PetscSection section;
1453: PetscBool includesConstraints;
1454: Vec gvec;
1455: PetscInt m, bs;
1457: PetscFunctionBegin;
1458: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1459: /* Check consistency */
1460: {
1461: PetscSF pointsf, pointsf1;
1463: PetscCall(DMGetPointSF(dm, &pointsf));
1464: PetscCall(DMGetPointSF(sectiondm, &pointsf1));
1465: PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
1466: }
1467: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1468: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
1469: PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
1470: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1471: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1472: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1473: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1474: PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
1475: PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
1476: PetscCall(VecGetBlockSize(vec, &bs));
1477: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs));
1478: PetscCall(VecCreate(comm, &gvec));
1479: PetscCall(PetscObjectSetName((PetscObject)gvec, vec_name));
1480: PetscCall(DMGetGlobalSection(sectiondm, §ion));
1481: PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints));
1482: if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m));
1483: else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m));
1484: PetscCall(VecSetSizes(gvec, m, PETSC_DECIDE));
1485: PetscCall(VecSetUp(gvec));
1486: PetscCall(DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec));
1487: PetscCall(DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec));
1488: PetscCall(VecView(gvec, viewer));
1489: PetscCall(VecDestroy(&gvec));
1490: PetscCall(PetscViewerHDF5PopGroup(viewer));
1491: PetscCall(PetscViewerHDF5PopGroup(viewer));
1492: PetscCall(PetscViewerHDF5PopGroup(viewer));
1493: PetscCall(PetscViewerHDF5PopGroup(viewer));
1494: PetscCall(PetscViewerHDF5PopGroup(viewer));
1495: PetscCall(PetscViewerHDF5PopGroup(viewer));
1496: PetscFunctionReturn(PETSC_SUCCESS);
1497: }
1499: struct _n_LoadLabelsCtx {
1500: MPI_Comm comm;
1501: PetscMPIInt rank;
1502: DM dm;
1503: PetscViewer viewer;
1504: DMLabel label;
1505: PetscSF sfXC;
1506: PetscLayout layoutX;
1507: };
1508: typedef struct _n_LoadLabelsCtx *LoadLabelsCtx;
1510: static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx)
1511: {
1512: PetscFunctionBegin;
1513: PetscCall(PetscNew(ctx));
1514: PetscCall(PetscObjectReference((PetscObject)((*ctx)->dm = dm)));
1515: PetscCall(PetscObjectReference((PetscObject)((*ctx)->viewer = viewer)));
1516: PetscCall(PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm));
1517: PetscCallMPI(MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank));
1518: (*ctx)->sfXC = sfXC;
1519: if (sfXC) {
1520: PetscInt nX;
1522: PetscCall(PetscObjectReference((PetscObject)sfXC));
1523: PetscCall(PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL));
1524: PetscCall(PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX));
1525: }
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx)
1530: {
1531: PetscFunctionBegin;
1532: if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
1533: PetscCall(DMDestroy(&(*ctx)->dm));
1534: PetscCall(PetscViewerDestroy(&(*ctx)->viewer));
1535: PetscCall(PetscSFDestroy(&(*ctx)->sfXC));
1536: PetscCall(PetscLayoutDestroy(&(*ctx)->layoutX));
1537: PetscCall(PetscFree(*ctx));
1538: PetscFunctionReturn(PETSC_SUCCESS);
1539: }
1541: /*
1542: A: on-disk points
1543: X: global points [0, NX)
1544: C: distributed plex points
1545: */
1546: static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS)
1547: {
1548: MPI_Comm comm = ctx->comm;
1549: PetscSF sfXC = ctx->sfXC;
1550: PetscLayout layoutX = ctx->layoutX;
1551: PetscSF sfXA;
1552: const PetscInt *A_points;
1553: PetscInt nX, nC;
1554: PetscInt n;
1556: PetscFunctionBegin;
1557: PetscCall(PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL));
1558: PetscCall(ISGetLocalSize(stratumIS, &n));
1559: PetscCall(ISGetIndices(stratumIS, &A_points));
1560: PetscCall(PetscSFCreate(comm, &sfXA));
1561: PetscCall(PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points));
1562: PetscCall(ISCreate(comm, newStratumIS));
1563: PetscCall(ISSetType(*newStratumIS, ISGENERAL));
1564: {
1565: PetscInt i;
1566: PetscBool *A_mask, *X_mask, *C_mask;
1568: PetscCall(PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask));
1569: for (i = 0; i < n; i++) A_mask[i] = PETSC_TRUE;
1570: PetscCall(PetscSFReduceBegin(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE));
1571: PetscCall(PetscSFReduceEnd(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE));
1572: PetscCall(PetscSFBcastBegin(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR));
1573: PetscCall(PetscSFBcastEnd(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR));
1574: PetscCall(ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask));
1575: PetscCall(PetscFree3(A_mask, X_mask, C_mask));
1576: }
1577: PetscCall(PetscSFDestroy(&sfXA));
1578: PetscCall(ISRestoreIndices(stratumIS, &A_points));
1579: PetscFunctionReturn(PETSC_SUCCESS);
1580: }
1582: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data)
1583: {
1584: LoadLabelsCtx ctx = (LoadLabelsCtx)op_data;
1585: PetscViewer viewer = ctx->viewer;
1586: DMLabel label = ctx->label;
1587: MPI_Comm comm = ctx->comm;
1588: IS stratumIS;
1589: const PetscInt *ind;
1590: PetscInt value, N, i;
1592: PetscCall(PetscOptionsStringToInt(vname, &value));
1593: PetscCall(ISCreate(comm, &stratumIS));
1594: PetscCall(PetscObjectSetName((PetscObject)stratumIS, "indices"));
1595: PetscCall(PetscViewerHDF5PushGroup(viewer, vname)); /* labels/<lname>/<vname> */
1597: if (!ctx->sfXC) {
1598: /* Force serial load */
1599: PetscCall(PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N));
1600: PetscCall(PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0));
1601: PetscCall(PetscLayoutSetSize(stratumIS->map, N));
1602: }
1603: PetscCall(ISLoad(stratumIS, viewer));
1605: if (ctx->sfXC) {
1606: IS newStratumIS;
1608: PetscCallHDF5(ReadLabelStratumHDF5_Distribute_Private, (stratumIS, ctx, &newStratumIS));
1609: PetscCall(ISDestroy(&stratumIS));
1610: stratumIS = newStratumIS;
1611: }
1613: PetscCall(PetscViewerHDF5PopGroup(viewer));
1614: PetscCall(ISGetLocalSize(stratumIS, &N));
1615: PetscCall(ISGetIndices(stratumIS, &ind));
1616: for (i = 0; i < N; ++i) PetscCall(DMLabelSetValue(label, ind[i], value));
1617: PetscCall(ISRestoreIndices(stratumIS, &ind));
1618: PetscCall(ISDestroy(&stratumIS));
1619: return 0;
1620: }
1622: /* TODO: Fix this code, it is returning PETSc error codes when it should be translating them to herr_t codes */
1623: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data)
1624: {
1625: LoadLabelsCtx ctx = (LoadLabelsCtx)op_data;
1626: DM dm = ctx->dm;
1627: hsize_t idx = 0;
1628: PetscErrorCode ierr;
1629: PetscBool flg;
1630: herr_t err;
1632: PetscCall(DMHasLabel(dm, lname, &flg));
1633: if (flg) PetscCall(DMRemoveLabel(dm, lname, NULL));
1634: ierr = DMCreateLabel(dm, lname);
1635: if (ierr) return (herr_t)ierr;
1636: ierr = DMGetLabel(dm, lname, &ctx->label);
1637: if (ierr) return (herr_t)ierr;
1638: ierr = PetscViewerHDF5PushGroup(ctx->viewer, lname);
1639: if (ierr) return (herr_t)ierr;
1640: /* Iterate over the label's strata */
1641: PetscCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1642: ierr = PetscViewerHDF5PopGroup(ctx->viewer);
1643: if (ierr) return (herr_t)ierr;
1644: return err;
1645: }
1647: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1648: {
1649: const char *topologydm_name;
1650: LoadLabelsCtx ctx;
1651: hsize_t idx = 0;
1652: char group[PETSC_MAX_PATH_LEN];
1653: DMPlexStorageVersion version;
1654: PetscBool distributed, hasGroup;
1656: PetscFunctionBegin;
1657: PetscCall(DMPlexIsDistributed(dm, &distributed));
1658: if (distributed) PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load");
1659: PetscCall(LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx));
1660: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1661: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
1662: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1663: PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name));
1664: } else {
1665: PetscCall(PetscStrncpy(group, "labels", sizeof(group)));
1666: }
1667: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1668: PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup));
1669: if (hasGroup) {
1670: hid_t fileId, groupId;
1672: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &fileId, &groupId));
1673: /* Iterate over labels */
1674: PetscCallHDF5(H5Literate, (groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx));
1675: PetscCallHDF5(H5Gclose, (groupId));
1676: }
1677: PetscCall(PetscViewerHDF5PopGroup(viewer));
1678: PetscCall(LoadLabelsCtxDestroy(&ctx));
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: static PetscErrorCode DMPlexDistributionLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF sf, PetscSF *distsf, DM *distdm)
1683: {
1684: MPI_Comm comm;
1685: PetscMPIInt size, rank;
1686: PetscInt dist_size;
1687: const char *distribution_name;
1688: PetscInt p, lsize;
1689: IS chartSizesIS, ownersIS, gpointsIS;
1690: const PetscInt *chartSize, *owners, *gpoints;
1691: PetscLayout layout;
1692: PetscBool has;
1694: PetscFunctionBegin;
1695: *distsf = NULL;
1696: *distdm = NULL;
1697: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1698: PetscCallMPI(MPI_Comm_size(comm, &size));
1699: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1700: PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
1701: if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS);
1702: PetscCall(PetscLogEventBegin(DMPLEX_DistributionLoad, viewer, 0, 0, 0));
1703: PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &has));
1704: if (!has) {
1705: char *full_group;
1707: PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &full_group));
1708: PetscCheck(has, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Distribution %s cannot be found: HDF5 group %s not found in file", distribution_name, full_group);
1709: }
1710: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "comm_size", PETSC_INT, NULL, (void *)&dist_size));
1711: PetscCheck(dist_size == (PetscInt)size, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mismatching comm sizes: comm size of this session (%d) != comm size used for %s (%" PetscInt_FMT ")", size, distribution_name, dist_size);
1712: PetscCall(ISCreate(comm, &chartSizesIS));
1713: PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes"));
1714: PetscCall(ISCreate(comm, &ownersIS));
1715: PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners"));
1716: PetscCall(ISCreate(comm, &gpointsIS));
1717: PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers"));
1718: PetscCall(PetscLayoutSetLocalSize(chartSizesIS->map, 1));
1719: PetscCall(ISLoad(chartSizesIS, viewer));
1720: PetscCall(ISGetIndices(chartSizesIS, &chartSize));
1721: PetscCall(PetscLayoutSetLocalSize(ownersIS->map, *chartSize));
1722: PetscCall(PetscLayoutSetLocalSize(gpointsIS->map, *chartSize));
1723: PetscCall(ISLoad(ownersIS, viewer));
1724: PetscCall(ISLoad(gpointsIS, viewer));
1725: PetscCall(ISGetIndices(ownersIS, &owners));
1726: PetscCall(ISGetIndices(gpointsIS, &gpoints));
1727: PetscCall(PetscSFCreate(comm, distsf));
1728: PetscCall(PetscSFSetFromOptions(*distsf));
1729: PetscCall(PetscLayoutCreate(comm, &layout));
1730: PetscCall(PetscSFGetGraph(sf, &lsize, NULL, NULL, NULL));
1731: PetscCall(PetscLayoutSetLocalSize(layout, lsize));
1732: PetscCall(PetscLayoutSetBlockSize(layout, 1));
1733: PetscCall(PetscLayoutSetUp(layout));
1734: PetscCall(PetscSFSetGraphLayout(*distsf, layout, *chartSize, NULL, PETSC_OWN_POINTER, gpoints));
1735: PetscCall(PetscLayoutDestroy(&layout));
1736: /* Migrate DM */
1737: {
1738: PetscInt pStart, pEnd;
1739: PetscSFNode *buffer0, *buffer1, *buffer2;
1741: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1742: PetscCall(PetscMalloc2(pEnd - pStart, &buffer0, lsize, &buffer1));
1743: PetscCall(PetscMalloc1(*chartSize, &buffer2));
1744: {
1745: PetscSF workPointSF;
1746: PetscInt workNroots, workNleaves;
1747: const PetscInt *workIlocal;
1748: const PetscSFNode *workIremote;
1750: for (p = pStart; p < pEnd; ++p) {
1751: buffer0[p - pStart].rank = rank;
1752: buffer0[p - pStart].index = p - pStart;
1753: }
1754: PetscCall(DMGetPointSF(dm, &workPointSF));
1755: PetscCall(PetscSFGetGraph(workPointSF, &workNroots, &workNleaves, &workIlocal, &workIremote));
1756: for (p = 0; p < workNleaves; ++p) {
1757: PetscInt workIlocalp = (workIlocal ? workIlocal[p] : p);
1759: buffer0[workIlocalp].rank = -1;
1760: }
1761: }
1762: for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
1763: for (p = 0; p < *chartSize; ++p) buffer2[p].rank = -1;
1764: PetscCall(PetscSFReduceBegin(sf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC));
1765: PetscCall(PetscSFReduceEnd(sf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC));
1766: PetscCall(PetscSFBcastBegin(*distsf, MPIU_2INT, buffer1, buffer2, MPI_REPLACE));
1767: PetscCall(PetscSFBcastEnd(*distsf, MPIU_2INT, buffer1, buffer2, MPI_REPLACE));
1768: if (PetscDefined(USE_DEBUG)) {
1769: for (p = 0; p < *chartSize; ++p) {
1770: PetscCheck(buffer2[p].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Found negative root rank %" PetscInt_FMT " at local point %" PetscInt_FMT " on rank %d when making migrationSF", buffer2[p].rank, p, rank);
1771: }
1772: }
1773: PetscCall(PetscFree2(buffer0, buffer1));
1774: PetscCall(DMCreate(comm, distdm));
1775: PetscCall(DMSetType(*distdm, DMPLEX));
1776: PetscCall(PetscObjectSetName((PetscObject)*distdm, ((PetscObject)dm)->name));
1777: PetscCall(DMPlexDistributionSetName(*distdm, distribution_name));
1778: {
1779: PetscSF migrationSF;
1781: PetscCall(PetscSFCreate(comm, &migrationSF));
1782: PetscCall(PetscSFSetFromOptions(migrationSF));
1783: PetscCall(PetscSFSetGraph(migrationSF, pEnd - pStart, *chartSize, NULL, PETSC_OWN_POINTER, buffer2, PETSC_OWN_POINTER));
1784: PetscCall(PetscSFSetUp(migrationSF));
1785: PetscCall(DMPlexMigrate(dm, migrationSF, *distdm));
1786: PetscCall(PetscSFDestroy(&migrationSF));
1787: }
1788: }
1789: /* Set pointSF */
1790: {
1791: PetscSF pointSF;
1792: PetscInt *ilocal, nleaves, q;
1793: PetscSFNode *iremote, *buffer0, *buffer1;
1795: PetscCall(PetscMalloc2(*chartSize, &buffer0, lsize, &buffer1));
1796: for (p = 0, nleaves = 0; p < *chartSize; ++p) {
1797: if (owners[p] == rank) {
1798: buffer0[p].rank = rank;
1799: } else {
1800: buffer0[p].rank = -1;
1801: nleaves++;
1802: }
1803: buffer0[p].index = p;
1804: }
1805: for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
1806: PetscCall(PetscSFReduceBegin(*distsf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC));
1807: PetscCall(PetscSFReduceEnd(*distsf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC));
1808: for (p = 0; p < *chartSize; ++p) buffer0[p].rank = -1;
1809: PetscCall(PetscSFBcastBegin(*distsf, MPIU_2INT, buffer1, buffer0, MPI_REPLACE));
1810: PetscCall(PetscSFBcastEnd(*distsf, MPIU_2INT, buffer1, buffer0, MPI_REPLACE));
1811: if (PetscDefined(USE_DEBUG)) {
1812: for (p = 0; p < *chartSize; ++p) PetscCheck(buffer0[p].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Found negative root rank %" PetscInt_FMT " at local point %" PetscInt_FMT " on rank %d when making pointSF", buffer0[p].rank, p, rank);
1813: }
1814: PetscCall(PetscMalloc1(nleaves, &ilocal));
1815: PetscCall(PetscMalloc1(nleaves, &iremote));
1816: for (p = 0, q = 0; p < *chartSize; ++p) {
1817: if (buffer0[p].rank != rank) {
1818: ilocal[q] = p;
1819: iremote[q].rank = buffer0[p].rank;
1820: iremote[q].index = buffer0[p].index;
1821: q++;
1822: }
1823: }
1824: PetscCheck(q == nleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching leaf sizes: %" PetscInt_FMT " != %" PetscInt_FMT, q, nleaves);
1825: PetscCall(PetscFree2(buffer0, buffer1));
1826: PetscCall(PetscSFCreate(comm, &pointSF));
1827: PetscCall(PetscSFSetFromOptions(pointSF));
1828: PetscCall(PetscSFSetGraph(pointSF, *chartSize, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1829: PetscCall(DMSetPointSF(*distdm, pointSF));
1830: {
1831: DM cdm;
1833: PetscCall(DMGetCoordinateDM(*distdm, &cdm));
1834: PetscCall(DMSetPointSF(cdm, pointSF));
1835: }
1836: PetscCall(PetscSFDestroy(&pointSF));
1837: }
1838: PetscCall(ISRestoreIndices(chartSizesIS, &chartSize));
1839: PetscCall(ISRestoreIndices(ownersIS, &owners));
1840: PetscCall(ISRestoreIndices(gpointsIS, &gpoints));
1841: PetscCall(ISDestroy(&chartSizesIS));
1842: PetscCall(ISDestroy(&ownersIS));
1843: PetscCall(ISDestroy(&gpointsIS));
1844: /* Record that overlap has been manually created. */
1845: /* This is to pass `DMPlexCheckPointSF()`, which checks that */
1846: /* pointSF does not contain cells in the leaves if overlap = 0. */
1847: PetscCall(DMPlexSetOverlap_Plex(*distdm, NULL, DMPLEX_OVERLAP_MANUAL));
1848: PetscCall(DMPlexDistributeSetDefault(*distdm, PETSC_FALSE));
1849: PetscCall(DMPlexReorderSetDefault(*distdm, DM_REORDER_DEFAULT_FALSE));
1850: PetscCall(PetscLogEventEnd(DMPLEX_DistributionLoad, viewer, 0, 0, 0));
1851: PetscFunctionReturn(PETSC_SUCCESS);
1852: }
1854: static PetscErrorCode DMPlexTopologyLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer, PetscSF *sf)
1855: {
1856: MPI_Comm comm;
1857: const char *pointsName, *coneSizesName, *conesName, *orientationsName;
1858: IS pointsIS, coneSizesIS, conesIS, orientationsIS;
1859: const PetscInt *points, *coneSizes, *cones, *orientations;
1860: PetscInt *cone, *ornt;
1861: PetscInt dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
1862: PetscMPIInt size, rank;
1864: PetscFunctionBegin;
1865: pointsName = "order";
1866: coneSizesName = "cones";
1867: conesName = "cells";
1868: orientationsName = "orientation";
1869: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1870: PetscCallMPI(MPI_Comm_size(comm, &size));
1871: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1872: PetscCall(ISCreate(comm, &pointsIS));
1873: PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName));
1874: PetscCall(ISCreate(comm, &coneSizesIS));
1875: PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
1876: PetscCall(ISCreate(comm, &conesIS));
1877: PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
1878: PetscCall(ISCreate(comm, &orientationsIS));
1879: PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
1880: PetscCall(PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject)conesIS, "cell_dim", PETSC_INT, NULL, &dim));
1881: PetscCall(DMSetDimension(dm, dim));
1882: {
1883: /* Force serial load */
1884: PetscCall(PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np));
1885: PetscCall(PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0));
1886: PetscCall(PetscLayoutSetSize(pointsIS->map, Np));
1887: pEnd = rank == 0 ? Np : 0;
1888: PetscCall(PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np));
1889: PetscCall(PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0));
1890: PetscCall(PetscLayoutSetSize(coneSizesIS->map, Np));
1891: PetscCall(PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N));
1892: PetscCall(PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0));
1893: PetscCall(PetscLayoutSetSize(conesIS->map, N));
1894: PetscCall(PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N));
1895: PetscCall(PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0));
1896: PetscCall(PetscLayoutSetSize(orientationsIS->map, N));
1897: }
1898: PetscCall(ISLoad(pointsIS, viewer));
1899: PetscCall(ISLoad(coneSizesIS, viewer));
1900: PetscCall(ISLoad(conesIS, viewer));
1901: PetscCall(ISLoad(orientationsIS, viewer));
1902: /* Create Plex */
1903: PetscCall(DMPlexSetChart(dm, 0, pEnd));
1904: PetscCall(ISGetIndices(pointsIS, &points));
1905: PetscCall(ISGetIndices(coneSizesIS, &coneSizes));
1906: for (p = 0; p < pEnd; ++p) {
1907: PetscCall(DMPlexSetConeSize(dm, points[p], coneSizes[p]));
1908: maxConeSize = PetscMax(maxConeSize, coneSizes[p]);
1909: }
1910: PetscCall(DMSetUp(dm));
1911: PetscCall(ISGetIndices(conesIS, &cones));
1912: PetscCall(ISGetIndices(orientationsIS, &orientations));
1913: PetscCall(PetscMalloc2(maxConeSize, &cone, maxConeSize, &ornt));
1914: for (p = 0, q = 0; p < pEnd; ++p) {
1915: for (c = 0; c < coneSizes[p]; ++c, ++q) {
1916: cone[c] = cones[q];
1917: ornt[c] = orientations[q];
1918: }
1919: PetscCall(DMPlexSetCone(dm, points[p], cone));
1920: PetscCall(DMPlexSetConeOrientation(dm, points[p], ornt));
1921: }
1922: PetscCall(PetscFree2(cone, ornt));
1923: /* Create global section migration SF */
1924: if (sf) {
1925: PetscLayout layout;
1926: PetscInt *globalIndices;
1928: PetscCall(PetscMalloc1(pEnd, &globalIndices));
1929: /* plex point == globalPointNumber in this case */
1930: for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
1931: PetscCall(PetscLayoutCreate(comm, &layout));
1932: PetscCall(PetscLayoutSetSize(layout, Np));
1933: PetscCall(PetscLayoutSetBlockSize(layout, 1));
1934: PetscCall(PetscLayoutSetUp(layout));
1935: PetscCall(PetscSFCreate(comm, sf));
1936: PetscCall(PetscSFSetFromOptions(*sf));
1937: PetscCall(PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices));
1938: PetscCall(PetscLayoutDestroy(&layout));
1939: PetscCall(PetscFree(globalIndices));
1940: }
1941: /* Clean-up */
1942: PetscCall(ISRestoreIndices(pointsIS, &points));
1943: PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes));
1944: PetscCall(ISRestoreIndices(conesIS, &cones));
1945: PetscCall(ISRestoreIndices(orientationsIS, &orientations));
1946: PetscCall(ISDestroy(&pointsIS));
1947: PetscCall(ISDestroy(&coneSizesIS));
1948: PetscCall(ISDestroy(&conesIS));
1949: PetscCall(ISDestroy(&orientationsIS));
1950: /* Fill in the rest of the topology structure */
1951: PetscCall(DMPlexSymmetrize(dm));
1952: PetscCall(DMPlexStratify(dm));
1953: PetscFunctionReturn(PETSC_SUCCESS);
1954: }
1956: /* Representation of two DMPlex strata in 0-based global numbering */
1957: struct _n_PlexLayer {
1958: PetscInt d;
1959: IS conesIS, orientationsIS;
1960: PetscSection coneSizesSection;
1961: PetscLayout vertexLayout;
1962: PetscSF overlapSF, l2gSF; //TODO maybe confusing names (in DMPlex in general)
1963: PetscInt offset, conesOffset, leafOffset;
1964: };
1965: typedef struct _n_PlexLayer *PlexLayer;
1967: static PetscErrorCode PlexLayerDestroy(PlexLayer *layer)
1968: {
1969: PetscFunctionBegin;
1970: if (!*layer) PetscFunctionReturn(PETSC_SUCCESS);
1971: PetscCall(PetscSectionDestroy(&(*layer)->coneSizesSection));
1972: PetscCall(ISDestroy(&(*layer)->conesIS));
1973: PetscCall(ISDestroy(&(*layer)->orientationsIS));
1974: PetscCall(PetscSFDestroy(&(*layer)->overlapSF));
1975: PetscCall(PetscSFDestroy(&(*layer)->l2gSF));
1976: PetscCall(PetscLayoutDestroy(&(*layer)->vertexLayout));
1977: PetscCall(PetscFree(*layer));
1978: PetscFunctionReturn(PETSC_SUCCESS);
1979: }
1981: static PetscErrorCode PlexLayerCreate_Private(PlexLayer *layer)
1982: {
1983: PetscFunctionBegin;
1984: PetscCall(PetscNew(layer));
1985: (*layer)->d = -1;
1986: (*layer)->offset = -1;
1987: (*layer)->conesOffset = -1;
1988: (*layer)->leafOffset = -1;
1989: PetscFunctionReturn(PETSC_SUCCESS);
1990: }
1992: static PetscErrorCode PlexLayerLoad_Private(PlexLayer layer, PetscViewer viewer, PetscInt d, PetscLayout pointsLayout)
1993: {
1994: char path[128];
1995: MPI_Comm comm;
1996: const char *coneSizesName, *conesName, *orientationsName;
1997: IS coneSizesIS, conesIS, orientationsIS;
1998: PetscSection coneSizesSection;
1999: PetscLayout vertexLayout = NULL;
2000: PetscInt s;
2002: PetscFunctionBegin;
2003: coneSizesName = "cone_sizes";
2004: conesName = "cones";
2005: orientationsName = "orientations";
2006: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
2008: /* query size of next lower depth stratum (next lower dimension) */
2009: if (d > 0) {
2010: PetscInt NVertices;
2012: PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT "/%s", d - 1, coneSizesName));
2013: PetscCall(PetscViewerHDF5ReadSizes(viewer, path, NULL, &NVertices));
2014: PetscCall(PetscLayoutCreate(comm, &vertexLayout));
2015: PetscCall(PetscLayoutSetSize(vertexLayout, NVertices));
2016: PetscCall(PetscLayoutSetUp(vertexLayout));
2017: }
2019: PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT, d));
2020: PetscCall(PetscViewerHDF5PushGroup(viewer, path));
2022: /* create coneSizesSection from stored IS coneSizes */
2023: {
2024: const PetscInt *coneSizes;
2026: PetscCall(ISCreate(comm, &coneSizesIS));
2027: PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
2028: if (pointsLayout) PetscCall(ISSetLayout(coneSizesIS, pointsLayout));
2029: PetscCall(ISLoad(coneSizesIS, viewer));
2030: if (!pointsLayout) PetscCall(ISGetLayout(coneSizesIS, &pointsLayout));
2031: PetscCall(ISGetIndices(coneSizesIS, &coneSizes));
2032: PetscCall(PetscSectionCreate(comm, &coneSizesSection));
2033: //TODO different start ?
2034: PetscCall(PetscSectionSetChart(coneSizesSection, 0, pointsLayout->n));
2035: for (s = 0; s < pointsLayout->n; ++s) PetscCall(PetscSectionSetDof(coneSizesSection, s, coneSizes[s]));
2036: PetscCall(PetscSectionSetUp(coneSizesSection));
2037: PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes));
2038: {
2039: PetscLayout tmp = NULL;
2041: /* We need to keep the layout until the end of function */
2042: PetscCall(PetscLayoutReference((PetscLayout)pointsLayout, &tmp));
2043: }
2044: PetscCall(ISDestroy(&coneSizesIS));
2045: }
2047: /* use value layout of coneSizesSection as layout of cones and orientations */
2048: {
2049: PetscLayout conesLayout;
2051: PetscCall(PetscSectionGetValueLayout(comm, coneSizesSection, &conesLayout));
2052: PetscCall(ISCreate(comm, &conesIS));
2053: PetscCall(ISCreate(comm, &orientationsIS));
2054: PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
2055: PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
2056: PetscCall(PetscLayoutDuplicate(conesLayout, &conesIS->map));
2057: PetscCall(PetscLayoutDuplicate(conesLayout, &orientationsIS->map));
2058: PetscCall(ISLoad(conesIS, viewer));
2059: PetscCall(ISLoad(orientationsIS, viewer));
2060: PetscCall(PetscLayoutDestroy(&conesLayout));
2061: }
2063: /* check assertion that layout of points is the same as point layout of coneSizesSection */
2064: {
2065: PetscLayout pointsLayout0;
2066: PetscBool flg;
2068: PetscCall(PetscSectionGetPointLayout(comm, coneSizesSection, &pointsLayout0));
2069: PetscCall(PetscLayoutCompare(pointsLayout, pointsLayout0, &flg));
2070: PetscCheck(flg, comm, PETSC_ERR_PLIB, "points layout != coneSizesSection point layout");
2071: PetscCall(PetscLayoutDestroy(&pointsLayout0));
2072: }
2073: PetscCall(PetscViewerHDF5PopGroup(viewer));
2074: PetscCall(PetscLayoutDestroy(&pointsLayout));
2076: layer->d = d;
2077: layer->conesIS = conesIS;
2078: layer->coneSizesSection = coneSizesSection;
2079: layer->orientationsIS = orientationsIS;
2080: layer->vertexLayout = vertexLayout;
2081: PetscFunctionReturn(PETSC_SUCCESS);
2082: }
2084: static PetscErrorCode PlexLayerDistribute_Private(PlexLayer layer, PetscSF cellLocalToGlobalSF)
2085: {
2086: IS newConesIS, newOrientationsIS;
2087: PetscSection newConeSizesSection;
2088: MPI_Comm comm;
2090: PetscFunctionBegin;
2091: PetscCall(PetscObjectGetComm((PetscObject)cellLocalToGlobalSF, &comm));
2092: PetscCall(PetscSectionCreate(comm, &newConeSizesSection));
2093: //TODO rename to something like ISDistribute() with optional PetscSection argument
2094: PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->conesIS, newConeSizesSection, &newConesIS));
2095: PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->orientationsIS, newConeSizesSection, &newOrientationsIS));
2097: PetscCall(PetscObjectSetName((PetscObject)newConeSizesSection, ((PetscObject)layer->coneSizesSection)->name));
2098: PetscCall(PetscObjectSetName((PetscObject)newConesIS, ((PetscObject)layer->conesIS)->name));
2099: PetscCall(PetscObjectSetName((PetscObject)newOrientationsIS, ((PetscObject)layer->orientationsIS)->name));
2100: PetscCall(PetscSectionDestroy(&layer->coneSizesSection));
2101: PetscCall(ISDestroy(&layer->conesIS));
2102: PetscCall(ISDestroy(&layer->orientationsIS));
2103: layer->coneSizesSection = newConeSizesSection;
2104: layer->conesIS = newConesIS;
2105: layer->orientationsIS = newOrientationsIS;
2106: PetscFunctionReturn(PETSC_SUCCESS);
2107: }
2109: //TODO share code with DMPlexBuildFromCellListParallel()
2110: #include <petsc/private/hashseti.h>
2111: static PetscErrorCode PlexLayerCreateSFs_Private(PlexLayer layer, PetscSF *vertexOverlapSF, PetscSF *sfXC)
2112: {
2113: PetscLayout vertexLayout = layer->vertexLayout;
2114: PetscSection coneSection = layer->coneSizesSection;
2115: IS cellVertexData = layer->conesIS;
2116: IS coneOrientations = layer->orientationsIS;
2117: PetscSF vl2gSF, vOverlapSF;
2118: PetscInt *verticesAdj;
2119: PetscInt i, n, numVerticesAdj;
2120: const PetscInt *cvd, *co = NULL;
2121: MPI_Comm comm;
2123: PetscFunctionBegin;
2124: PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm));
2125: PetscCall(PetscSectionGetStorageSize(coneSection, &n));
2126: {
2127: PetscInt n0;
2129: PetscCall(ISGetLocalSize(cellVertexData, &n0));
2130: PetscCheck(n == n0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local size of IS cellVertexData = %" PetscInt_FMT " != %" PetscInt_FMT " = storage size of PetscSection coneSection", n0, n);
2131: PetscCall(ISGetIndices(cellVertexData, &cvd));
2132: }
2133: if (coneOrientations) {
2134: PetscInt n0;
2136: PetscCall(ISGetLocalSize(coneOrientations, &n0));
2137: PetscCheck(n == n0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local size of IS coneOrientations = %" PetscInt_FMT " != %" PetscInt_FMT " = storage size of PetscSection coneSection", n0, n);
2138: PetscCall(ISGetIndices(coneOrientations, &co));
2139: }
2140: /* Get/check global number of vertices */
2141: {
2142: PetscInt NVerticesInCells = PETSC_MIN_INT;
2144: /* NVerticesInCells = max(cellVertexData) + 1 */
2145: for (i = 0; i < n; i++)
2146: if (cvd[i] > NVerticesInCells) NVerticesInCells = cvd[i];
2147: ++NVerticesInCells;
2148: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, comm));
2150: if (vertexLayout->n == PETSC_DECIDE && vertexLayout->N == PETSC_DECIDE) vertexLayout->N = NVerticesInCells;
2151: else
2152: PetscCheck(vertexLayout->N == PETSC_DECIDE || vertexLayout->N >= NVerticesInCells, comm, PETSC_ERR_ARG_SIZ, "Specified global number of vertices %" PetscInt_FMT " must be greater than or equal to the number of unique vertices in the cell-vertex dataset %" PetscInt_FMT,
2153: vertexLayout->N, NVerticesInCells);
2154: PetscCall(PetscLayoutSetUp(vertexLayout));
2155: }
2156: /* Find locally unique vertices in cellVertexData */
2157: {
2158: PetscHSetI vhash;
2159: PetscInt off = 0;
2161: PetscCall(PetscHSetICreate(&vhash));
2162: for (i = 0; i < n; ++i) PetscCall(PetscHSetIAdd(vhash, cvd[i]));
2163: PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj));
2164: PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj));
2165: PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj));
2166: PetscAssert(off == numVerticesAdj, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assertion failed: off == numVerticesAdj (%" PetscInt_FMT " != %" PetscInt_FMT ")", off, numVerticesAdj);
2167: PetscCall(PetscHSetIDestroy(&vhash));
2168: }
2169: /* We must sort vertices to preserve numbering */
2170: PetscCall(PetscSortInt(numVerticesAdj, verticesAdj));
2171: /* Connect locally unique vertices with star forest */
2172: PetscCall(PetscSFCreateByMatchingIndices(vertexLayout, numVerticesAdj, verticesAdj, NULL, 0, numVerticesAdj, verticesAdj, NULL, 0, &vl2gSF, &vOverlapSF));
2173: PetscCall(PetscObjectSetName((PetscObject)vOverlapSF, "overlapSF"));
2174: PetscCall(PetscObjectSetName((PetscObject)vl2gSF, "localToGlobalSF"));
2176: PetscCall(PetscFree(verticesAdj));
2177: *vertexOverlapSF = vOverlapSF;
2178: *sfXC = vl2gSF;
2179: PetscFunctionReturn(PETSC_SUCCESS);
2180: }
2182: static PetscErrorCode PlexLayerCreateCellSFs_Private(PlexLayer layer, PetscSF *cellOverlapSF, PetscSF *cellLocalToGlobalSF)
2183: {
2184: PetscSection coneSection = layer->coneSizesSection;
2185: PetscInt nCells;
2186: MPI_Comm comm;
2188: PetscFunctionBegin;
2189: PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm));
2190: {
2191: PetscInt cStart;
2193: PetscCall(PetscSectionGetChart(coneSection, &cStart, &nCells));
2194: PetscCheck(cStart == 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "coneSection must start at 0");
2195: }
2196: /* Create overlapSF as empty SF with the right number of roots */
2197: PetscCall(PetscSFCreate(comm, cellOverlapSF));
2198: PetscCall(PetscSFSetGraph(*cellOverlapSF, nCells, 0, NULL, PETSC_USE_POINTER, NULL, PETSC_USE_POINTER));
2199: PetscCall(PetscSFSetUp(*cellOverlapSF));
2200: /* Create localToGlobalSF as identity mapping */
2201: {
2202: PetscLayout map;
2204: PetscCall(PetscLayoutCreateFromSizes(comm, nCells, PETSC_DECIDE, 1, &map));
2205: PetscCall(PetscSFCreateFromLayouts(map, map, cellLocalToGlobalSF));
2206: PetscCall(PetscSFSetUp(*cellLocalToGlobalSF));
2207: PetscCall(PetscLayoutDestroy(&map));
2208: }
2209: PetscFunctionReturn(PETSC_SUCCESS);
2210: }
2212: static PetscErrorCode DMPlexTopologyBuildFromLayers_Private(DM dm, PetscInt depth, PlexLayer *layers, IS strataPermutation)
2213: {
2214: const PetscInt *permArr;
2215: PetscInt d, nPoints;
2216: MPI_Comm comm;
2218: PetscFunctionBegin;
2219: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2220: PetscCall(ISGetIndices(strataPermutation, &permArr));
2222: /* Count points, strata offsets and cones offsets (taking strataPermutation into account) */
2223: {
2224: PetscInt stratumOffset = 0;
2225: PetscInt conesOffset = 0;
2227: for (d = 0; d <= depth; d++) {
2228: const PetscInt e = permArr[d];
2229: const PlexLayer l = layers[e];
2230: PetscInt lo, n, size;
2232: PetscCall(PetscSectionGetChart(l->coneSizesSection, &lo, &n));
2233: PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &size));
2234: PetscCheck(lo == 0, comm, PETSC_ERR_PLIB, "starting point should be 0 in coneSizesSection %" PetscInt_FMT, d);
2235: l->offset = stratumOffset;
2236: l->conesOffset = conesOffset;
2237: stratumOffset += n;
2238: conesOffset += size;
2239: }
2240: nPoints = stratumOffset;
2241: }
2243: /* Set interval for all plex points */
2244: //TODO we should store starting point of plex
2245: PetscCall(DMPlexSetChart(dm, 0, nPoints));
2247: /* Set up plex coneSection from layer coneSections */
2248: {
2249: PetscSection coneSection;
2251: PetscCall(DMPlexGetConeSection(dm, &coneSection));
2252: for (d = 0; d <= depth; d++) {
2253: const PlexLayer l = layers[d];
2254: PetscInt n, q;
2256: PetscCall(PetscSectionGetChart(l->coneSizesSection, NULL, &n));
2257: for (q = 0; q < n; q++) {
2258: const PetscInt p = l->offset + q;
2259: PetscInt coneSize;
2261: PetscCall(PetscSectionGetDof(l->coneSizesSection, q, &coneSize));
2262: PetscCall(PetscSectionSetDof(coneSection, p, coneSize));
2263: }
2264: }
2265: }
2266: //TODO this is terrible, DMSetUp_Plex() should be DMPlexSetUpSections() or so
2267: PetscCall(DMSetUp(dm));
2269: /* Renumber cones points from layer-global numbering to plex-local numbering */
2270: {
2271: PetscInt *cones, *ornts;
2273: PetscCall(DMPlexGetCones(dm, &cones));
2274: PetscCall(DMPlexGetConeOrientations(dm, &ornts));
2275: for (d = 1; d <= depth; d++) {
2276: const PlexLayer l = layers[d];
2277: PetscInt i, lConesSize;
2278: PetscInt *lCones;
2279: const PetscInt *lOrnts;
2280: PetscInt *pCones = &cones[l->conesOffset];
2281: PetscInt *pOrnts = &ornts[l->conesOffset];
2283: PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &lConesSize));
2284: /* Get cones in local plex numbering */
2285: {
2286: ISLocalToGlobalMapping l2g;
2287: PetscLayout vertexLayout = l->vertexLayout;
2288: PetscSF vertexSF = layers[d - 1]->l2gSF; /* vertices of this layer are cells of previous layer */
2289: const PetscInt *gCones;
2290: PetscInt lConesSize0;
2292: PetscCall(ISGetLocalSize(l->conesIS, &lConesSize0));
2293: PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(conesIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize);
2294: PetscCall(ISGetLocalSize(l->orientationsIS, &lConesSize0));
2295: PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(orientationsIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize);
2297: PetscCall(PetscMalloc1(lConesSize, &lCones));
2298: PetscCall(ISGetIndices(l->conesIS, &gCones));
2299: PetscCall(ISLocalToGlobalMappingCreateSF(vertexSF, vertexLayout->rstart, &l2g));
2300: PetscCall(ISGlobalToLocalMappingApply(l2g, IS_GTOLM_MASK, lConesSize, gCones, &lConesSize0, lCones));
2301: PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "global to local does not cover all indices (%" PetscInt_FMT " of %" PetscInt_FMT ")", lConesSize0, lConesSize);
2302: PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
2303: PetscCall(ISRestoreIndices(l->conesIS, &gCones));
2304: }
2305: PetscCall(ISGetIndices(l->orientationsIS, &lOrnts));
2306: /* Set cones, need to add stratum offset */
2307: for (i = 0; i < lConesSize; i++) {
2308: pCones[i] = lCones[i] + layers[d - 1]->offset; /* cone points of current layer are points of previous layer */
2309: pOrnts[i] = lOrnts[i];
2310: }
2311: PetscCall(PetscFree(lCones));
2312: PetscCall(ISRestoreIndices(l->orientationsIS, &lOrnts));
2313: }
2314: }
2315: PetscCall(DMPlexSymmetrize(dm));
2316: PetscCall(DMPlexStratify(dm));
2317: PetscCall(ISRestoreIndices(strataPermutation, &permArr));
2318: PetscFunctionReturn(PETSC_SUCCESS);
2319: }
2321: static PetscErrorCode PlexLayerConcatenateSFs_Private(MPI_Comm comm, PetscInt depth, PlexLayer layers[], IS strataPermutation, PetscSF *overlapSF, PetscSF *l2gSF)
2322: {
2323: PetscInt d;
2324: PetscSF *osfs, *lsfs;
2325: PetscInt *leafOffsets;
2326: const PetscInt *permArr;
2328: PetscFunctionBegin;
2329: PetscCall(ISGetIndices(strataPermutation, &permArr));
2330: PetscCall(PetscCalloc3(depth + 1, &osfs, depth + 1, &lsfs, depth + 1, &leafOffsets));
2331: for (d = 0; d <= depth; d++) {
2332: const PetscInt e = permArr[d];
2334: PetscAssert(e == layers[e]->d, PETSC_COMM_SELF, PETSC_ERR_PLIB, "assertion: e == layers[e]->d");
2335: osfs[d] = layers[e]->overlapSF;
2336: lsfs[d] = layers[e]->l2gSF;
2337: leafOffsets[d] = layers[e]->offset;
2338: }
2339: PetscCall(PetscSFConcatenate(comm, depth + 1, osfs, PETSCSF_CONCATENATE_ROOTMODE_LOCAL, leafOffsets, overlapSF));
2340: PetscCall(PetscSFConcatenate(comm, depth + 1, lsfs, PETSCSF_CONCATENATE_ROOTMODE_GLOBAL, leafOffsets, l2gSF));
2341: PetscCall(PetscFree3(osfs, lsfs, leafOffsets));
2342: PetscCall(ISRestoreIndices(strataPermutation, &permArr));
2343: PetscFunctionReturn(PETSC_SUCCESS);
2344: }
2346: static PetscErrorCode DMPlexTopologyLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF *sfXC)
2347: {
2348: PlexLayer *layers;
2349: IS strataPermutation;
2350: PetscLayout pointsLayout;
2351: PetscInt depth;
2352: PetscInt d;
2353: MPI_Comm comm;
2355: PetscFunctionBegin;
2356: {
2357: PetscInt dim;
2359: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "depth", PETSC_INT, NULL, &depth));
2360: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "cell_dim", PETSC_INT, NULL, &dim));
2361: PetscCall(DMSetDimension(dm, dim));
2362: }
2363: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2365: {
2366: IS spOnComm;
2368: PetscCall(ISCreate(comm, &spOnComm));
2369: PetscCall(PetscObjectSetName((PetscObject)spOnComm, "permutation"));
2370: PetscCall(ISLoad(spOnComm, viewer));
2371: /* have the same serial IS on every rank */
2372: PetscCall(ISAllGather(spOnComm, &strataPermutation));
2373: PetscCall(PetscObjectSetName((PetscObject)strataPermutation, ((PetscObject)spOnComm)->name));
2374: PetscCall(ISDestroy(&spOnComm));
2375: }
2377: /* Create layers, load raw data for each layer */
2378: PetscCall(PetscViewerHDF5PushGroup(viewer, "strata"));
2379: PetscCall(PetscMalloc1(depth + 1, &layers));
2380: for (d = depth, pointsLayout = NULL; d >= 0; pointsLayout = layers[d]->vertexLayout, d--) {
2381: PetscCall(PlexLayerCreate_Private(&layers[d]));
2382: PetscCall(PlexLayerLoad_Private(layers[d], viewer, d, pointsLayout));
2383: }
2384: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */
2386: for (d = depth; d >= 0; d--) {
2387: /* Redistribute cells and vertices for each applicable layer */
2388: if (d < depth) PetscCall(PlexLayerDistribute_Private(layers[d], layers[d]->l2gSF));
2389: /* Create vertex overlap SF and vertex localToGlobal SF for each applicable layer */
2390: if (d > 0) PetscCall(PlexLayerCreateSFs_Private(layers[d], &layers[d - 1]->overlapSF, &layers[d - 1]->l2gSF));
2391: }
2392: /* Build trivial SFs for the cell layer as well */
2393: PetscCall(PlexLayerCreateCellSFs_Private(layers[depth], &layers[depth]->overlapSF, &layers[depth]->l2gSF));
2395: /* Build DMPlex topology from the layers */
2396: PetscCall(DMPlexTopologyBuildFromLayers_Private(dm, depth, layers, strataPermutation));
2398: /* Build overall point SF alias overlap SF */
2399: {
2400: PetscSF overlapSF;
2402: PetscCall(PlexLayerConcatenateSFs_Private(comm, depth, layers, strataPermutation, &overlapSF, sfXC));
2403: PetscCall(DMSetPointSF(dm, overlapSF));
2404: PetscCall(PetscSFDestroy(&overlapSF));
2405: }
2407: for (d = depth; d >= 0; d--) PetscCall(PlexLayerDestroy(&layers[d]));
2408: PetscCall(PetscFree(layers));
2409: PetscCall(ISDestroy(&strataPermutation));
2410: PetscFunctionReturn(PETSC_SUCCESS);
2411: }
2413: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sfXC)
2414: {
2415: DMPlexStorageVersion version;
2416: const char *topologydm_name;
2417: char group[PETSC_MAX_PATH_LEN];
2418: PetscSF sfwork = NULL;
2420: PetscFunctionBegin;
2421: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2422: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2423: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
2424: PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name));
2425: } else {
2426: PetscCall(PetscStrncpy(group, "/", sizeof(group)));
2427: }
2428: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
2430: PetscCall(PetscViewerHDF5PushGroup(viewer, "topology"));
2431: if (version->major < 3) {
2432: PetscCall(DMPlexTopologyLoad_HDF5_Legacy_Private(dm, viewer, &sfwork));
2433: } else {
2434: /* since DMPlexStorageVersion 3.0.0 */
2435: PetscCall(DMPlexTopologyLoad_HDF5_Private(dm, viewer, &sfwork));
2436: }
2437: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */
2439: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
2440: DM distdm;
2441: PetscSF distsf;
2442: const char *distribution_name;
2443: PetscBool exists;
2445: PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
2446: /* this group is guaranteed to be present since DMPlexStorageVersion 2.1.0 */
2447: PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions"));
2448: PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &exists));
2449: if (exists) {
2450: PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name));
2451: PetscCall(DMPlexDistributionLoad_HDF5_Private(dm, viewer, sfwork, &distsf, &distdm));
2452: if (distdm) {
2453: PetscCall(DMPlexReplace_Internal(dm, &distdm));
2454: PetscCall(PetscSFDestroy(&sfwork));
2455: sfwork = distsf;
2456: }
2457: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */
2458: }
2459: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */
2460: }
2461: if (sfXC) {
2462: *sfXC = sfwork;
2463: } else {
2464: PetscCall(PetscSFDestroy(&sfwork));
2465: }
2467: PetscCall(PetscViewerHDF5PopGroup(viewer));
2468: PetscFunctionReturn(PETSC_SUCCESS);
2469: }
2471: /* If the file is old, it not only has different path to the coordinates, but */
2472: /* does not contain coordinateDMs, so must fall back to the old implementation. */
2473: static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
2474: {
2475: PetscSection coordSection;
2476: Vec coordinates;
2477: PetscReal lengthScale;
2478: PetscInt spatialDim, N, numVertices, vStart, vEnd, v;
2479: PetscMPIInt rank;
2481: PetscFunctionBegin;
2482: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
2483: /* Read geometry */
2484: PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry"));
2485: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates));
2486: PetscCall(PetscObjectSetName((PetscObject)coordinates, "vertices"));
2487: {
2488: /* Force serial load */
2489: PetscCall(PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N));
2490: PetscCall(VecSetSizes(coordinates, rank == 0 ? N : 0, N));
2491: PetscCall(VecSetBlockSize(coordinates, spatialDim));
2492: }
2493: PetscCall(VecLoad(coordinates, viewer));
2494: PetscCall(PetscViewerHDF5PopGroup(viewer));
2495: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
2496: PetscCall(VecScale(coordinates, 1.0 / lengthScale));
2497: PetscCall(VecGetLocalSize(coordinates, &numVertices));
2498: PetscCall(VecGetBlockSize(coordinates, &spatialDim));
2499: numVertices /= spatialDim;
2500: /* Create coordinates */
2501: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
2502: PetscCheck(numVertices == vEnd - vStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of coordinates loaded %" PetscInt_FMT " does not match number of vertices %" PetscInt_FMT, numVertices, vEnd - vStart);
2503: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2504: PetscCall(PetscSectionSetNumFields(coordSection, 1));
2505: PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spatialDim));
2506: PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd));
2507: for (v = vStart; v < vEnd; ++v) {
2508: PetscCall(PetscSectionSetDof(coordSection, v, spatialDim));
2509: PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spatialDim));
2510: }
2511: PetscCall(PetscSectionSetUp(coordSection));
2512: PetscCall(DMSetCoordinates(dm, coordinates));
2513: PetscCall(VecDestroy(&coordinates));
2514: PetscFunctionReturn(PETSC_SUCCESS);
2515: }
2517: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
2518: {
2519: DMPlexStorageVersion version;
2520: DM cdm;
2521: Vec coords;
2522: PetscInt blockSize;
2523: PetscReal lengthScale;
2524: PetscSF lsf;
2525: const char *topologydm_name;
2526: char *coordinatedm_name, *coordinates_name;
2528: PetscFunctionBegin;
2529: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2530: if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
2531: PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer));
2532: PetscFunctionReturn(PETSC_SUCCESS);
2533: }
2534: /* else: since DMPlexStorageVersion 2.0.0 */
2535: PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load");
2536: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2537: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
2538: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
2539: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, NULL, &coordinatedm_name));
2540: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, NULL, &coordinates_name));
2541: PetscCall(PetscViewerHDF5PopGroup(viewer));
2542: PetscCall(PetscViewerHDF5PopGroup(viewer));
2543: PetscCall(DMGetCoordinateDM(dm, &cdm));
2544: PetscCall(PetscObjectSetName((PetscObject)cdm, coordinatedm_name));
2545: PetscCall(PetscFree(coordinatedm_name));
2546: /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
2547: PetscCall(DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf));
2548: PetscCall(DMCreateLocalVector(cdm, &coords));
2549: PetscCall(PetscObjectSetName((PetscObject)coords, coordinates_name));
2550: PetscCall(PetscFree(coordinates_name));
2551: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
2552: PetscCall(DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords));
2553: PetscCall(PetscViewerPopFormat(viewer));
2554: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
2555: PetscCall(VecScale(coords, 1.0 / lengthScale));
2556: PetscCall(DMSetCoordinatesLocal(dm, coords));
2557: PetscCall(VecGetBlockSize(coords, &blockSize));
2558: PetscCall(DMSetCoordinateDim(dm, blockSize));
2559: PetscCall(VecDestroy(&coords));
2560: PetscCall(PetscSFDestroy(&lsf));
2561: PetscFunctionReturn(PETSC_SUCCESS);
2562: }
2564: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
2565: {
2566: DMPlexStorageVersion version;
2568: PetscFunctionBegin;
2569: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2570: if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
2571: PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL));
2572: PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL));
2573: PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer));
2574: } else {
2575: PetscSF sfXC;
2577: /* since DMPlexStorageVersion 2.0.0 */
2578: PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC));
2579: PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC));
2580: PetscCall(DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC));
2581: PetscCall(PetscSFDestroy(&sfXC));
2582: }
2583: PetscFunctionReturn(PETSC_SUCCESS);
2584: }
2586: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2587: {
2588: MPI_Comm comm;
2589: PetscInt pStart, pEnd, p, m;
2590: PetscInt *goffs, *ilocal;
2591: PetscBool rootIncludeConstraints, leafIncludeConstraints;
2593: PetscFunctionBegin;
2594: PetscCall(PetscObjectGetComm((PetscObject)leafSection, &comm));
2595: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2596: PetscCall(PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints));
2597: PetscCall(PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints));
2598: if (rootIncludeConstraints && leafIncludeConstraints) PetscCall(PetscSectionGetStorageSize(leafSection, &m));
2599: else PetscCall(PetscSectionGetConstrainedStorageSize(leafSection, &m));
2600: PetscCall(PetscMalloc1(m, &ilocal));
2601: PetscCall(PetscMalloc1(m, &goffs));
2602: /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
2603: /* for the top-level section (not for each field), so one must have */
2604: /* rootSection->pointMajor == PETSC_TRUE. */
2605: PetscCheck(rootSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
2606: /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
2607: PetscCheck(leafSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
2608: for (p = pStart, m = 0; p < pEnd; ++p) {
2609: PetscInt dof, cdof, i, j, off, goff;
2610: const PetscInt *cinds;
2612: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2613: if (dof < 0) continue;
2614: goff = globalOffsets[p - pStart];
2615: PetscCall(PetscSectionGetOffset(leafSection, p, &off));
2616: PetscCall(PetscSectionGetConstraintDof(leafSection, p, &cdof));
2617: PetscCall(PetscSectionGetConstraintIndices(leafSection, p, &cinds));
2618: for (i = 0, j = 0; i < dof; ++i) {
2619: PetscBool constrained = (PetscBool)(j < cdof && i == cinds[j]);
2621: if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {
2622: ilocal[m] = off++;
2623: goffs[m++] = goff++;
2624: } else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
2625: else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff;
2626: if (constrained) ++j;
2627: }
2628: }
2629: PetscCall(PetscSFCreate(comm, sectionSF));
2630: PetscCall(PetscSFSetFromOptions(*sectionSF));
2631: PetscCall(PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs));
2632: PetscCall(PetscFree(goffs));
2633: PetscFunctionReturn(PETSC_SUCCESS);
2634: }
2636: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
2637: {
2638: MPI_Comm comm;
2639: PetscMPIInt size, rank;
2640: const char *topologydm_name;
2641: const char *sectiondm_name;
2642: PetscSection sectionA, sectionB;
2643: PetscInt nX, n, i;
2644: PetscSF sfAB;
2646: PetscFunctionBegin;
2647: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2648: PetscCallMPI(MPI_Comm_size(comm, &size));
2649: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2650: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2651: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
2652: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
2653: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
2654: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
2655: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
2656: /* A: on-disk points */
2657: /* X: list of global point numbers, [0, NX) */
2658: /* B: plex points */
2659: /* Load raw section (sectionA) */
2660: PetscCall(PetscSectionCreate(comm, §ionA));
2661: PetscCall(PetscSectionLoad(sectionA, viewer));
2662: PetscCall(PetscSectionGetChart(sectionA, NULL, &n));
2663: /* Create sfAB: A -> B */
2664: #if defined(PETSC_USE_DEBUG)
2665: {
2666: PetscInt N, N1;
2668: PetscCall(PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1));
2669: PetscCall(MPIU_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm));
2670: PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Mismatching sizes: on-disk order array size (%" PetscInt_FMT ") != number of loaded section points (%" PetscInt_FMT ")", N1, N);
2671: }
2672: #endif
2673: {
2674: IS orderIS;
2675: const PetscInt *gpoints;
2676: PetscSF sfXA, sfAX;
2677: PetscLayout layout;
2678: PetscSFNode *owners, *buffer;
2679: PetscInt nleaves;
2680: PetscInt *ilocal;
2681: PetscSFNode *iremote;
2683: /* Create sfAX: A -> X */
2684: PetscCall(ISCreate(comm, &orderIS));
2685: PetscCall(PetscObjectSetName((PetscObject)orderIS, "order"));
2686: PetscCall(PetscLayoutSetLocalSize(orderIS->map, n));
2687: PetscCall(ISLoad(orderIS, viewer));
2688: PetscCall(PetscLayoutCreate(comm, &layout));
2689: PetscCall(PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL));
2690: PetscCall(PetscLayoutSetLocalSize(layout, nX));
2691: PetscCall(PetscLayoutSetBlockSize(layout, 1));
2692: PetscCall(PetscLayoutSetUp(layout));
2693: PetscCall(PetscSFCreate(comm, &sfXA));
2694: PetscCall(ISGetIndices(orderIS, &gpoints));
2695: PetscCall(PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints));
2696: PetscCall(ISRestoreIndices(orderIS, &gpoints));
2697: PetscCall(ISDestroy(&orderIS));
2698: PetscCall(PetscLayoutDestroy(&layout));
2699: PetscCall(PetscMalloc1(n, &owners));
2700: PetscCall(PetscMalloc1(nX, &buffer));
2701: for (i = 0; i < n; ++i) {
2702: owners[i].rank = rank;
2703: owners[i].index = i;
2704: }
2705: for (i = 0; i < nX; ++i) {
2706: buffer[i].rank = -1;
2707: buffer[i].index = -1;
2708: }
2709: PetscCall(PetscSFReduceBegin(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC));
2710: PetscCall(PetscSFReduceEnd(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC));
2711: PetscCall(PetscSFDestroy(&sfXA));
2712: PetscCall(PetscFree(owners));
2713: for (i = 0, nleaves = 0; i < nX; ++i)
2714: if (buffer[i].rank >= 0) nleaves++;
2715: PetscCall(PetscMalloc1(nleaves, &ilocal));
2716: PetscCall(PetscMalloc1(nleaves, &iremote));
2717: for (i = 0, nleaves = 0; i < nX; ++i) {
2718: if (buffer[i].rank >= 0) {
2719: ilocal[nleaves] = i;
2720: iremote[nleaves].rank = buffer[i].rank;
2721: iremote[nleaves].index = buffer[i].index;
2722: nleaves++;
2723: }
2724: }
2725: PetscCall(PetscFree(buffer));
2726: PetscCall(PetscSFCreate(comm, &sfAX));
2727: PetscCall(PetscSFSetFromOptions(sfAX));
2728: PetscCall(PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2729: PetscCall(PetscSFCompose(sfAX, sfXB, &sfAB));
2730: PetscCall(PetscSFDestroy(&sfAX));
2731: }
2732: PetscCall(PetscViewerHDF5PopGroup(viewer));
2733: PetscCall(PetscViewerHDF5PopGroup(viewer));
2734: PetscCall(PetscViewerHDF5PopGroup(viewer));
2735: PetscCall(PetscViewerHDF5PopGroup(viewer));
2736: /* Create plex section (sectionB) */
2737: PetscCall(DMGetLocalSection(sectiondm, §ionB));
2738: if (lsf || gsf) {
2739: PetscLayout layout;
2740: PetscInt M, m;
2741: PetscInt *offsetsA;
2742: PetscBool includesConstraintsA;
2744: PetscCall(PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB));
2745: PetscCall(PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA));
2746: if (includesConstraintsA) PetscCall(PetscSectionGetStorageSize(sectionA, &m));
2747: else PetscCall(PetscSectionGetConstrainedStorageSize(sectionA, &m));
2748: PetscCall(MPIU_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm));
2749: PetscCall(PetscLayoutCreate(comm, &layout));
2750: PetscCall(PetscLayoutSetSize(layout, M));
2751: PetscCall(PetscLayoutSetUp(layout));
2752: if (lsf) {
2753: PetscSF lsfABdata;
2755: PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata));
2756: *lsf = lsfABdata;
2757: }
2758: if (gsf) {
2759: PetscSection gsectionB, gsectionB1;
2760: PetscBool includesConstraintsB;
2761: PetscSF gsfABdata, pointsf;
2763: PetscCall(DMGetGlobalSection(sectiondm, &gsectionB1));
2764: PetscCall(PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB));
2765: PetscCall(DMGetPointSF(sectiondm, &pointsf));
2766: PetscCall(PetscSectionCreateGlobalSection(sectionB, pointsf, PETSC_TRUE, includesConstraintsB, PETSC_TRUE, &gsectionB));
2767: PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata));
2768: PetscCall(PetscSectionDestroy(&gsectionB));
2769: *gsf = gsfABdata;
2770: }
2771: PetscCall(PetscLayoutDestroy(&layout));
2772: PetscCall(PetscFree(offsetsA));
2773: } else {
2774: PetscCall(PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB));
2775: }
2776: PetscCall(PetscSFDestroy(&sfAB));
2777: PetscCall(PetscSectionDestroy(§ionA));
2778: PetscFunctionReturn(PETSC_SUCCESS);
2779: }
2781: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
2782: {
2783: MPI_Comm comm;
2784: const char *topologydm_name;
2785: const char *sectiondm_name;
2786: const char *vec_name;
2787: Vec vecA;
2788: PetscInt mA, m, bs;
2789: const PetscInt *ilocal;
2790: const PetscScalar *src;
2791: PetscScalar *dest;
2793: PetscFunctionBegin;
2794: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2795: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2796: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
2797: PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
2798: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
2799: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
2800: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
2801: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
2802: PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
2803: PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
2804: PetscCall(VecCreate(comm, &vecA));
2805: PetscCall(PetscObjectSetName((PetscObject)vecA, vec_name));
2806: PetscCall(PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL));
2807: /* Check consistency */
2808: {
2809: PetscSF pointsf, pointsf1;
2810: PetscInt m1, i, j;
2812: PetscCall(DMGetPointSF(dm, &pointsf));
2813: PetscCall(DMGetPointSF(sectiondm, &pointsf1));
2814: PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
2815: #if defined(PETSC_USE_DEBUG)
2816: {
2817: PetscInt MA, MA1;
2819: PetscCall(MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm));
2820: PetscCall(PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1));
2821: PetscCheck(MA1 == MA, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total SF root size (%" PetscInt_FMT ") != On-disk vector data size (%" PetscInt_FMT ")", MA, MA1);
2822: }
2823: #endif
2824: PetscCall(VecGetLocalSize(vec, &m1));
2825: PetscCheck(m1 >= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Target vector size (%" PetscInt_FMT ") < SF leaf size (%" PetscInt_FMT ")", m1, m);
2826: for (i = 0; i < m; ++i) {
2827: j = ilocal ? ilocal[i] : i;
2828: PetscCheck(j >= 0 && j < m1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Leaf's %" PetscInt_FMT "-th index, %" PetscInt_FMT ", not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", i, j, (PetscInt)0, m1);
2829: }
2830: }
2831: PetscCall(VecSetSizes(vecA, mA, PETSC_DECIDE));
2832: PetscCall(VecLoad(vecA, viewer));
2833: PetscCall(VecGetArrayRead(vecA, &src));
2834: PetscCall(VecGetArray(vec, &dest));
2835: PetscCall(PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE));
2836: PetscCall(PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE));
2837: PetscCall(VecRestoreArray(vec, &dest));
2838: PetscCall(VecRestoreArrayRead(vecA, &src));
2839: PetscCall(VecDestroy(&vecA));
2840: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *)&bs));
2841: PetscCall(VecSetBlockSize(vec, bs));
2842: PetscCall(PetscViewerHDF5PopGroup(viewer));
2843: PetscCall(PetscViewerHDF5PopGroup(viewer));
2844: PetscCall(PetscViewerHDF5PopGroup(viewer));
2845: PetscCall(PetscViewerHDF5PopGroup(viewer));
2846: PetscCall(PetscViewerHDF5PopGroup(viewer));
2847: PetscCall(PetscViewerHDF5PopGroup(viewer));
2848: PetscFunctionReturn(PETSC_SUCCESS);
2849: }
2850: #endif