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