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