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_SINGLE_LIBRARY_INTERN 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;
778: PetscBool ocompress;
780: PetscFunctionBegin;
781: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
782: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
783: PetscCallMPI(MPI_Comm_size(comm, &size));
784: PetscCallMPI(MPI_Comm_rank(comm, &rank));
785: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
786: PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
787: if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS);
788: PetscCall(PetscLogEventBegin(DMPLEX_DistributionView, viewer, 0, 0, 0));
789: size_petsc_int = (PetscInt)size;
790: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "comm_size", PETSC_INT, (void *)&size_petsc_int));
791: PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
792: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
793: PetscCall(PetscMalloc1(1, &chartSize));
794: *chartSize = pEnd - pStart;
795: PetscCall(PetscMalloc1(*chartSize, &owners));
796: PetscCall(PetscMalloc1(*chartSize, &gpoints));
797: PetscCall(DMGetPointSF(dm, &pointSF));
798: PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote));
799: for (p = pStart; p < pEnd; ++p) {
800: PetscInt gp = gpoint[p - pStart];
802: owners[p - pStart] = rank;
803: gpoints[p - pStart] = (gp < 0 ? -(gp + 1) : gp);
804: }
805: for (p = 0; p < nleaves; ++p) {
806: PetscInt ilocalp = (ilocal ? ilocal[p] : p);
808: owners[ilocalp] = iremote[p].rank;
809: }
810: PetscCall(ISCreateGeneral(comm, 1, chartSize, PETSC_OWN_POINTER, &chartSizesIS));
811: PetscCall(ISCreateGeneral(comm, *chartSize, owners, PETSC_OWN_POINTER, &ownersIS));
812: PetscCall(ISCreateGeneral(comm, *chartSize, gpoints, PETSC_OWN_POINTER, &gpointsIS));
813: if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
814: PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
815: PetscCall(PetscViewerHDF5SetCompress(viewer, 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: if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
828: PetscCall(PetscLogEventEnd(DMPLEX_DistributionView, viewer, 0, 0, 0));
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: 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[])
833: {
834: DMPlexStorageVersion version;
835: IS coneSizesIS, conesIS, orientationsIS;
836: PetscInt *coneSizes, *cones, *orientations;
837: const PetscInt *gpoint;
838: PetscInt nPoints = 0, conesSize = 0;
839: PetscInt p, c, s;
840: PetscBool ocompress;
841: MPI_Comm comm;
843: PetscFunctionBegin;
844: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
845: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
846: PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
847: for (p = pStart; p < pEnd; ++p) {
848: if (gpoint[p] >= 0) {
849: PetscInt coneSize;
851: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
852: nPoints += 1;
853: conesSize += coneSize;
854: }
855: }
856: PetscCall(PetscMalloc1(nPoints, &coneSizes));
857: PetscCall(PetscMalloc1(conesSize, &cones));
858: PetscCall(PetscMalloc1(conesSize, &orientations));
859: for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
860: if (gpoint[p] >= 0) {
861: const PetscInt *cone, *ornt;
862: PetscInt coneSize, cp;
864: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
865: PetscCall(DMPlexGetCone(dm, p, &cone));
866: PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
867: coneSizes[s] = coneSize;
868: for (cp = 0; cp < coneSize; ++cp, ++c) {
869: cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]] + 1) : gpoint[cone[cp]];
870: orientations[c] = ornt[cp];
871: }
872: ++s;
873: }
874: }
875: PetscCheck(s == nPoints, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %" PetscInt_FMT " != %" PetscInt_FMT, s, nPoints);
876: PetscCheck(c == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %" PetscInt_FMT " != %" PetscInt_FMT, c, conesSize);
877: PetscCall(ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS));
878: PetscCall(ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS));
879: PetscCall(ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS));
880: PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
881: PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
882: PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
883: if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
884: PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
885: PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
886: }
887: PetscCall(ISView(coneSizesIS, viewer));
888: PetscCall(ISView(conesIS, viewer));
889: PetscCall(ISView(orientationsIS, viewer));
890: PetscCall(ISDestroy(&coneSizesIS));
891: PetscCall(ISDestroy(&conesIS));
892: PetscCall(ISDestroy(&orientationsIS));
893: if (pointsName) {
894: IS pointsIS;
895: PetscInt *points;
897: PetscCall(PetscMalloc1(nPoints, &points));
898: for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
899: if (gpoint[p] >= 0) {
900: points[s] = gpoint[p];
901: ++s;
902: }
903: }
904: PetscCall(ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS));
905: PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName));
906: PetscCall(ISView(pointsIS, viewer));
907: PetscCall(ISDestroy(&pointsIS));
908: }
909: PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
910: if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: static PetscErrorCode DMPlexTopologyView_HDF5_Legacy_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
915: {
916: const char *pointsName, *coneSizesName, *conesName, *orientationsName;
917: PetscInt pStart, pEnd, dim;
919: PetscFunctionBegin;
920: pointsName = "order";
921: coneSizesName = "cones";
922: conesName = "cells";
923: orientationsName = "orientation";
924: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
925: PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers, viewer, pStart, pEnd, pointsName, coneSizesName, conesName, orientationsName));
926: PetscCall(DMGetDimension(dm, &dim));
927: PetscCall(PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *)&dim));
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: //TODO get this numbering right away without needing this function
932: /* Renumber global point numbers so that they are 0-based per stratum */
933: static PetscErrorCode RenumberGlobalPointNumbersPerStratum_Private(DM dm, IS globalPointNumbers, IS *newGlobalPointNumbers, IS *strataPermutation)
934: {
935: PetscInt d, depth, p, n;
936: PetscInt *offsets;
937: const PetscInt *gpn;
938: PetscInt *ngpn;
939: MPI_Comm comm;
941: PetscFunctionBegin;
942: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
943: PetscCall(ISGetLocalSize(globalPointNumbers, &n));
944: PetscCall(ISGetIndices(globalPointNumbers, &gpn));
945: PetscCall(PetscMalloc1(n, &ngpn));
946: PetscCall(DMPlexGetDepth(dm, &depth));
947: PetscCall(PetscMalloc1(depth + 1, &offsets));
948: for (d = 0; d <= depth; d++) {
949: PetscInt pStart, pEnd;
951: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
952: offsets[d] = PETSC_INT_MAX;
953: for (p = pStart; p < pEnd; p++) {
954: if (gpn[p] >= 0 && gpn[p] < offsets[d]) offsets[d] = gpn[p];
955: }
956: }
957: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, offsets, depth + 1, MPIU_INT, MPI_MIN, comm));
958: for (d = 0; d <= depth; d++) {
959: PetscInt pStart, pEnd;
961: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
962: for (p = pStart; p < pEnd; p++) ngpn[p] = gpn[p] - PetscSign(gpn[p]) * offsets[d];
963: }
964: PetscCall(ISRestoreIndices(globalPointNumbers, &gpn));
965: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)globalPointNumbers), n, ngpn, PETSC_OWN_POINTER, newGlobalPointNumbers));
966: {
967: PetscInt *perm;
969: PetscCall(PetscMalloc1(depth + 1, &perm));
970: for (d = 0; d <= depth; d++) perm[d] = d;
971: PetscCall(PetscSortIntWithPermutation(depth + 1, offsets, perm));
972: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, depth + 1, perm, PETSC_OWN_POINTER, strataPermutation));
973: }
974: PetscCall(PetscFree(offsets));
975: PetscFunctionReturn(PETSC_SUCCESS);
976: }
978: static PetscErrorCode DMPlexTopologyView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
979: {
980: IS globalPointNumbers0, strataPermutation;
981: const char *coneSizesName, *conesName, *orientationsName;
982: PetscInt depth, d;
983: MPI_Comm comm;
985: PetscFunctionBegin;
986: coneSizesName = "cone_sizes";
987: conesName = "cones";
988: orientationsName = "orientations";
989: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
990: PetscCall(DMPlexGetDepth(dm, &depth));
991: {
992: PetscInt dim;
994: PetscCall(DMGetDimension(dm, &dim));
995: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "cell_dim", PETSC_INT, &dim));
996: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "depth", PETSC_INT, &depth));
997: }
999: PetscCall(RenumberGlobalPointNumbersPerStratum_Private(dm, globalPointNumbers, &globalPointNumbers0, &strataPermutation));
1000: /* TODO dirty trick to save serial IS using the same parallel viewer */
1001: {
1002: IS spOnComm;
1003: PetscInt n = 0, N;
1004: const PetscInt *idx = NULL;
1005: const PetscInt *old;
1006: PetscMPIInt rank;
1008: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1009: PetscCall(ISGetLocalSize(strataPermutation, &N));
1010: PetscCall(ISGetIndices(strataPermutation, &old));
1011: if (!rank) {
1012: n = N;
1013: idx = old;
1014: }
1015: PetscCall(ISCreateGeneral(comm, n, idx, PETSC_COPY_VALUES, &spOnComm));
1016: PetscCall(ISRestoreIndices(strataPermutation, &old));
1017: PetscCall(ISDestroy(&strataPermutation));
1018: strataPermutation = spOnComm;
1019: }
1020: PetscCall(PetscObjectSetName((PetscObject)strataPermutation, "permutation"));
1021: PetscCall(ISView(strataPermutation, viewer));
1022: PetscCall(PetscViewerHDF5PushGroup(viewer, "strata"));
1023: for (d = 0; d <= depth; d++) {
1024: PetscInt pStart, pEnd;
1025: char group[128];
1027: PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, d));
1028: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1029: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
1030: PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers0, viewer, pStart, pEnd, NULL, coneSizesName, conesName, orientationsName));
1031: PetscCall(PetscViewerHDF5PopGroup(viewer));
1032: }
1033: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */
1034: PetscCall(ISDestroy(&globalPointNumbers0));
1035: PetscCall(ISDestroy(&strataPermutation));
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
1040: {
1041: DMPlexStorageVersion version;
1042: const char *topologydm_name;
1043: char group[PETSC_MAX_PATH_LEN];
1045: PetscFunctionBegin;
1046: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
1047: PetscCall(PetscInfo(dm, "Writing DM %s storage version %d.%d.%d\n", dm->hdr.name, version->major, version->minor, version->subminor));
1048: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1049: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1050: PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name));
1051: } else {
1052: PetscCall(PetscStrncpy(group, "/", sizeof(group)));
1053: }
1054: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1056: PetscCall(PetscViewerHDF5PushGroup(viewer, "topology"));
1057: if (version->major < 3) {
1058: PetscCall(DMPlexTopologyView_HDF5_Legacy_Private(dm, globalPointNumbers, viewer));
1059: } else {
1060: /* since DMPlexStorageVersion 3.0.0 */
1061: PetscCall(DMPlexTopologyView_HDF5_Private(dm, globalPointNumbers, viewer));
1062: }
1063: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */
1065: if (DMPlexStorageVersionGE(version, 2, 1, 0)) {
1066: const char *distribution_name;
1068: PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
1069: PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions"));
1070: PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
1071: PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name));
1072: PetscCall(DMPlexDistributionView_HDF5_Private(dm, globalPointNumbers, viewer));
1073: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */
1074: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */
1075: }
1077: PetscCall(PetscViewerHDF5PopGroup(viewer));
1078: PetscFunctionReturn(PETSC_SUCCESS);
1079: }
1081: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
1082: {
1083: PetscSF sfPoint;
1084: DMLabel cutLabel, cutVertexLabel = NULL;
1085: IS globalVertexNumbers, cutvertices = NULL;
1086: const PetscInt *gcell, *gvertex, *cutverts = NULL;
1087: PetscInt *vertices;
1088: PetscInt conesSize = 0;
1089: PetscInt dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
1091: PetscFunctionBegin;
1092: *numCorners = 0;
1093: PetscCall(DMGetDimension(dm, &dim));
1094: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1095: PetscCall(ISGetIndices(globalCellNumbers, &gcell));
1097: for (cell = cStart; cell < cEnd; ++cell) {
1098: PetscInt *closure = NULL;
1099: PetscInt closureSize, v, Nc = 0;
1101: if (gcell[cell] < 0) continue;
1102: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
1103: for (v = 0; v < closureSize * 2; v += 2) {
1104: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
1105: }
1106: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
1107: conesSize += Nc;
1108: if (!numCornersLocal) numCornersLocal = Nc;
1109: else if (numCornersLocal != Nc) numCornersLocal = 1;
1110: }
1111: PetscCallMPI(MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
1112: PetscCheck(!numCornersLocal || !(numCornersLocal != *numCorners || *numCorners == 0), PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
1113: /* Handle periodic cuts by identifying vertices which should be duplicated */
1114: PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
1115: PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
1116: if (cutVertexLabel) PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices));
1117: if (cutvertices) {
1118: PetscCall(ISGetIndices(cutvertices, &cutverts));
1119: PetscCall(ISGetLocalSize(cutvertices, &vExtra));
1120: }
1121: PetscCall(DMGetPointSF(dm, &sfPoint));
1122: if (cutLabel) {
1123: const PetscInt *ilocal;
1124: const PetscSFNode *iremote;
1125: PetscInt nroots, nleaves;
1127: PetscCall(PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote));
1128: if (nleaves < 0) {
1129: PetscCall(PetscObjectReference((PetscObject)sfPoint));
1130: } else {
1131: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), &sfPoint));
1132: PetscCall(PetscSFSetGraph(sfPoint, nroots + vExtra, nleaves, (PetscInt *)ilocal, PETSC_COPY_VALUES, (PetscSFNode *)iremote, PETSC_COPY_VALUES));
1133: }
1134: } else {
1135: PetscCall(PetscObjectReference((PetscObject)sfPoint));
1136: }
1137: /* Number all vertices */
1138: PetscCall(DMPlexCreateNumbering_Plex(dm, vStart, vEnd + vExtra, 0, NULL, sfPoint, &globalVertexNumbers));
1139: PetscCall(PetscSFDestroy(&sfPoint));
1140: /* Create cones */
1141: PetscCall(ISGetIndices(globalVertexNumbers, &gvertex));
1142: PetscCall(PetscMalloc1(conesSize, &vertices));
1143: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
1144: PetscInt *closure = NULL;
1145: PetscInt closureSize, Nc = 0, p, value = -1;
1146: PetscBool replace;
1148: if (gcell[cell] < 0) continue;
1149: if (cutLabel) PetscCall(DMLabelGetValue(cutLabel, cell, &value));
1150: replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
1151: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
1152: for (p = 0; p < closureSize * 2; p += 2) {
1153: if ((closure[p] >= vStart) && (closure[p] < vEnd)) closure[Nc++] = closure[p];
1154: }
1155: PetscCall(DMPlexReorderCell(dm, cell, closure));
1156: for (p = 0; p < Nc; ++p) {
1157: PetscInt nv, gv = gvertex[closure[p] - vStart];
1159: if (replace) {
1160: PetscCall(PetscFindInt(closure[p], vExtra, cutverts, &nv));
1161: if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
1162: }
1163: vertices[v++] = gv < 0 ? -(gv + 1) : gv;
1164: }
1165: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
1166: }
1167: PetscCall(ISRestoreIndices(globalVertexNumbers, &gvertex));
1168: PetscCall(ISDestroy(&globalVertexNumbers));
1169: PetscCall(ISRestoreIndices(globalCellNumbers, &gcell));
1170: if (cutvertices) PetscCall(ISRestoreIndices(cutvertices, &cutverts));
1171: PetscCall(ISDestroy(&cutvertices));
1172: PetscCall(DMLabelDestroy(&cutVertexLabel));
1173: PetscCheck(v == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %" PetscInt_FMT " != %" PetscInt_FMT, v, conesSize);
1174: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS));
1175: PetscCall(PetscLayoutSetBlockSize((*cellIS)->map, *numCorners));
1176: PetscCall(PetscObjectSetName((PetscObject)*cellIS, "cells"));
1177: PetscFunctionReturn(PETSC_SUCCESS);
1178: }
1180: static PetscErrorCode DMPlexTopologyView_HDF5_XDMF_Private(DM dm, IS globalCellNumbers, PetscViewer viewer)
1181: {
1182: DM cdm;
1183: DMLabel depthLabel, ctLabel;
1184: IS cellIS;
1185: PetscInt dim, depth, cellHeight, c, n = 0;
1187: PetscFunctionBegin;
1188: PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz"));
1189: PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
1191: PetscCall(PetscViewerHDF5PopGroup(viewer));
1192: PetscCall(DMGetDimension(dm, &dim));
1193: PetscCall(DMPlexGetDepth(dm, &depth));
1194: PetscCall(DMGetCoordinateDM(dm, &cdm));
1195: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
1196: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
1197: PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel));
1198: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
1199: const DMPolytopeType ict = (DMPolytopeType)c;
1200: PetscInt pStart, pEnd, dep, numCorners;
1201: PetscBool output = PETSC_FALSE, doOutput;
1203: if (ict == DM_POLYTOPE_FV_GHOST) continue;
1204: PetscCall(DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd));
1205: if (pStart >= 0) {
1206: PetscCall(DMLabelGetValue(depthLabel, pStart, &dep));
1207: if (dep == depth - cellHeight) output = PETSC_TRUE;
1208: }
1209: PetscCallMPI(MPIU_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
1210: if (!doOutput) continue;
1211: PetscCall(CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS));
1212: if (!n) {
1213: PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/topology"));
1214: } else {
1215: char group[PETSC_MAX_PATH_LEN];
1217: PetscCall(PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%" PetscInt_FMT, n));
1218: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1219: }
1220: PetscCall(ISView(cellIS, viewer));
1221: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_corners", PETSC_INT, (void *)&numCorners));
1222: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_dim", PETSC_INT, (void *)&dim));
1223: PetscCall(ISDestroy(&cellIS));
1224: PetscCall(PetscViewerHDF5PopGroup(viewer));
1225: ++n;
1226: }
1227: PetscFunctionReturn(PETSC_SUCCESS);
1228: }
1230: static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
1231: {
1232: DM cdm;
1233: Vec coordinates, newcoords;
1234: PetscReal lengthScale;
1235: PetscInt m, M, bs;
1237: PetscFunctionBegin;
1238: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
1239: PetscCall(DMGetCoordinateDM(dm, &cdm));
1240: PetscCall(DMGetCoordinates(dm, &coordinates));
1241: PetscCall(VecCreate(PetscObjectComm((PetscObject)coordinates), &newcoords));
1242: PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices"));
1243: PetscCall(VecGetSize(coordinates, &M));
1244: PetscCall(VecGetLocalSize(coordinates, &m));
1245: PetscCall(VecSetSizes(newcoords, m, M));
1246: PetscCall(VecGetBlockSize(coordinates, &bs));
1247: PetscCall(VecSetBlockSize(newcoords, bs));
1248: PetscCall(VecSetType(newcoords, VECSTANDARD));
1249: PetscCall(VecCopy(coordinates, newcoords));
1250: PetscCall(VecScale(newcoords, lengthScale));
1251: /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
1252: PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry"));
1253: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
1254: PetscCall(VecView(newcoords, viewer));
1255: PetscCall(PetscViewerPopFormat(viewer));
1256: PetscCall(PetscViewerHDF5PopGroup(viewer));
1257: PetscCall(VecDestroy(&newcoords));
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1261: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
1262: {
1263: DM cdm;
1264: Vec coords, newcoords;
1265: PetscInt m, M, bs;
1266: PetscReal lengthScale;
1267: PetscBool viewSection = PETSC_TRUE, ocompress;
1268: const char *topologydm_name, *coordinatedm_name, *coordinates_name;
1269: PetscViewerFormat format;
1270: DMPlexStorageVersion version;
1272: PetscFunctionBegin;
1273: PetscCall(PetscViewerGetFormat(viewer, &format));
1274: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
1275: if (!DMPlexStorageVersionGE(version, 2, 0, 0) || format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ) {
1276: PetscCall(DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer));
1277: PetscFunctionReturn(PETSC_SUCCESS);
1278: }
1279: /* since 2.0.0 */
1280: if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
1281: PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
1282: PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
1283: }
1284: PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_view_coordinate_section", &viewSection, NULL));
1285: PetscCall(DMGetCoordinateDM(dm, &cdm));
1286: PetscCall(DMGetCoordinates(dm, &coords));
1287: PetscCall(PetscObjectGetName((PetscObject)cdm, &coordinatedm_name));
1288: PetscCall(PetscObjectGetName((PetscObject)coords, &coordinates_name));
1289: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1290: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1291: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1292: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name));
1293: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name));
1294: PetscCall(PetscViewerHDF5PopGroup(viewer));
1295: PetscCall(PetscViewerHDF5PopGroup(viewer));
1296: if (viewSection) PetscCall(DMPlexSectionView(dm, viewer, cdm));
1297: PetscCall(VecCreate(PetscObjectComm((PetscObject)coords), &newcoords));
1298: PetscCall(PetscObjectSetName((PetscObject)newcoords, coordinates_name));
1299: PetscCall(VecGetSize(coords, &M));
1300: PetscCall(VecGetLocalSize(coords, &m));
1301: PetscCall(VecSetSizes(newcoords, m, M));
1302: PetscCall(VecGetBlockSize(coords, &bs));
1303: PetscCall(VecSetBlockSize(newcoords, bs));
1304: PetscCall(VecSetType(newcoords, VECSTANDARD));
1305: PetscCall(VecCopy(coords, newcoords));
1306: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
1307: PetscCall(VecScale(newcoords, lengthScale));
1308: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
1309: PetscCall(DMPlexGlobalVectorView(dm, viewer, cdm, newcoords));
1310: PetscCall(PetscViewerPopFormat(viewer));
1311: PetscCall(VecDestroy(&newcoords));
1312: if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: static PetscErrorCode DMPlexCoordinatesView_HDF5_XDMF_Private(DM dm, PetscViewer viewer)
1317: {
1318: DM cdm;
1319: Vec coordinatesLocal, newcoords;
1320: PetscSection cSection, cGlobalSection;
1321: PetscScalar *coords, *ncoords;
1322: DMLabel cutLabel, cutVertexLabel = NULL;
1323: const PetscReal *L;
1324: PetscReal lengthScale;
1325: PetscInt vStart, vEnd, v, bs, N, coordSize, dof, off, d;
1326: PetscBool localized, embedded;
1328: PetscFunctionBegin;
1329: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1330: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
1331: PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal));
1332: PetscCall(VecGetBlockSize(coordinatesLocal, &bs));
1333: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
1334: if (localized == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1335: PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L));
1336: PetscCall(DMGetCoordinateDM(dm, &cdm));
1337: PetscCall(DMGetLocalSection(cdm, &cSection));
1338: PetscCall(DMGetGlobalSection(cdm, &cGlobalSection));
1339: PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
1340: N = 0;
1342: PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
1343: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &newcoords));
1344: PetscCall(PetscSectionGetDof(cSection, vStart, &dof));
1345: PetscCall(PetscPrintf(PETSC_COMM_SELF, "DOF: %" PetscInt_FMT "\n", dof));
1346: embedded = (PetscBool)(L && dof == 2 && !cutLabel);
1347: if (cutVertexLabel) {
1348: PetscCall(DMLabelGetStratumSize(cutVertexLabel, 1, &v));
1349: N += dof * v;
1350: }
1351: for (v = vStart; v < vEnd; ++v) {
1352: PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof));
1353: if (dof < 0) continue;
1354: if (embedded) N += dof + 1;
1355: else N += dof;
1356: }
1357: if (embedded) PetscCall(VecSetBlockSize(newcoords, bs + 1));
1358: else PetscCall(VecSetBlockSize(newcoords, bs));
1359: PetscCall(VecSetSizes(newcoords, N, PETSC_DETERMINE));
1360: PetscCall(VecSetType(newcoords, VECSTANDARD));
1361: PetscCall(VecGetArray(coordinatesLocal, &coords));
1362: PetscCall(VecGetArray(newcoords, &ncoords));
1363: coordSize = 0;
1364: for (v = vStart; v < vEnd; ++v) {
1365: PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof));
1366: PetscCall(PetscSectionGetOffset(cSection, v, &off));
1367: if (dof < 0) continue;
1368: if (embedded) {
1369: if (L && (L[0] > 0.0) && (L[1] > 0.0)) {
1370: PetscReal theta, phi, r, R;
1371: /* XY-periodic */
1372: /* Suppose its an y-z circle, then
1373: \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
1374: and the circle in that plane is
1375: \hat r cos(phi) + \hat x sin(phi) */
1376: theta = 2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1];
1377: phi = 2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0];
1378: r = L[0] / (2.0 * PETSC_PI * 2.0 * L[1]);
1379: R = L[1] / (2.0 * PETSC_PI);
1380: ncoords[coordSize++] = PetscSinReal(phi) * r;
1381: ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
1382: ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
1383: } else if (L && (L[0] > 0.0)) {
1384: /* X-periodic */
1385: ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
1386: ncoords[coordSize++] = coords[off + 1];
1387: ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
1388: } else if (L && (L[1] > 0.0)) {
1389: /* Y-periodic */
1390: ncoords[coordSize++] = coords[off + 0];
1391: ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
1392: ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
1393: #if 0
1394: } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
1395: PetscReal phi, r, R;
1396: /* Mobius strip */
1397: /* Suppose its an x-z circle, then
1398: \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
1399: and in that plane we rotate by pi as we go around the circle
1400: \hat r cos(phi/2) + \hat y sin(phi/2) */
1401: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
1402: R = L[0];
1403: r = PetscRealPart(coords[off+1]) - L[1]/2.0;
1404: ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
1405: ncoords[coordSize++] = PetscSinReal(phi/2.0) * r;
1406: ncoords[coordSize++] = PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
1407: #endif
1408: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
1409: } else {
1410: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off + d];
1411: }
1412: }
1413: if (cutVertexLabel) {
1414: IS vertices;
1415: const PetscInt *verts;
1416: PetscInt n;
1418: PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &vertices));
1419: if (vertices) {
1420: PetscCall(ISGetIndices(vertices, &verts));
1421: PetscCall(ISGetLocalSize(vertices, &n));
1422: for (v = 0; v < n; ++v) {
1423: PetscCall(PetscSectionGetDof(cSection, verts[v], &dof));
1424: PetscCall(PetscSectionGetOffset(cSection, verts[v], &off));
1425: for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off + d] + ((L[d] > 0.) ? L[d] : 0.0);
1426: }
1427: PetscCall(ISRestoreIndices(vertices, &verts));
1428: PetscCall(ISDestroy(&vertices));
1429: }
1430: }
1431: PetscCheck(coordSize == N, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %" PetscInt_FMT " != %" PetscInt_FMT, coordSize, N);
1432: PetscCall(DMLabelDestroy(&cutVertexLabel));
1433: PetscCall(VecRestoreArray(coordinatesLocal, &coords));
1434: PetscCall(VecRestoreArray(newcoords, &ncoords));
1435: PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices"));
1436: PetscCall(VecScale(newcoords, lengthScale));
1437: PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz"));
1438: PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
1439: PetscCall(PetscViewerHDF5PopGroup(viewer));
1440: PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/geometry"));
1441: PetscCall(VecView(newcoords, viewer));
1442: PetscCall(PetscViewerHDF5PopGroup(viewer));
1443: PetscCall(VecDestroy(&newcoords));
1444: PetscFunctionReturn(PETSC_SUCCESS);
1445: }
1447: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
1448: {
1449: const char *topologydm_name;
1450: const PetscInt *gpoint;
1451: PetscInt numLabels;
1452: PetscBool omitCelltypes = PETSC_FALSE, ocompress;
1453: DMPlexStorageVersion version;
1454: char group[PETSC_MAX_PATH_LEN];
1456: PetscFunctionBegin;
1457: PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_omit_celltypes", &omitCelltypes, NULL));
1458: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
1459: if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
1460: PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
1461: PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
1462: }
1463: PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
1464: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1465: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1466: PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name));
1467: } else {
1468: PetscCall(PetscStrncpy(group, "/labels", sizeof(group)));
1469: }
1470: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1471: PetscCall(DMGetNumLabels(dm, &numLabels));
1472: for (PetscInt l = 0; l < numLabels; ++l) {
1473: DMLabel label;
1474: const char *name;
1475: IS valueIS, pvalueIS, globalValueIS;
1476: const PetscInt *values;
1477: PetscInt numValues, v;
1478: PetscBool isDepth, isCelltype, output;
1480: PetscCall(DMGetLabelByNum(dm, l, &label));
1481: PetscCall(PetscObjectGetName((PetscObject)label, &name));
1482: PetscCall(DMGetLabelOutput(dm, name, &output));
1483: PetscCall(PetscStrncmp(name, "depth", 10, &isDepth));
1484: PetscCall(PetscStrncmp(name, "celltype", 10, &isCelltype));
1485: // TODO Should only filter out celltype if it can be calculated
1486: if (isDepth || (isCelltype && omitCelltypes) || !output) continue;
1487: PetscCall(PetscViewerHDF5PushGroup(viewer, name));
1488: PetscCall(DMLabelGetValueIS(label, &valueIS));
1489: /* Must copy to a new IS on the global comm */
1490: PetscCall(ISGetLocalSize(valueIS, &numValues));
1491: PetscCall(ISGetIndices(valueIS, &values));
1492: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS));
1493: PetscCall(ISRestoreIndices(valueIS, &values));
1494: PetscCall(ISAllGather(pvalueIS, &globalValueIS));
1495: PetscCall(ISDestroy(&pvalueIS));
1496: PetscCall(ISSortRemoveDups(globalValueIS));
1497: PetscCall(ISGetLocalSize(globalValueIS, &numValues));
1498: PetscCall(ISGetIndices(globalValueIS, &values));
1499: for (v = 0; v < numValues; ++v) {
1500: IS stratumIS, globalStratumIS;
1501: const PetscInt *spoints = NULL;
1502: PetscInt *gspoints, n = 0, gn, p;
1503: const char *iname = "indices";
1504: char group[PETSC_MAX_PATH_LEN];
1506: PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, values[v]));
1507: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1508: PetscCall(DMLabelGetStratumIS(label, values[v], &stratumIS));
1510: if (stratumIS) PetscCall(ISGetLocalSize(stratumIS, &n));
1511: if (stratumIS) PetscCall(ISGetIndices(stratumIS, &spoints));
1512: for (gn = 0, p = 0; p < n; ++p)
1513: if (gpoint[spoints[p]] >= 0) ++gn;
1514: PetscCall(PetscMalloc1(gn, &gspoints));
1515: for (gn = 0, p = 0; p < n; ++p)
1516: if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
1517: if (stratumIS) PetscCall(ISRestoreIndices(stratumIS, &spoints));
1518: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS));
1519: PetscCall(PetscObjectSetName((PetscObject)globalStratumIS, iname));
1520: PetscCall(ISView(globalStratumIS, viewer));
1521: PetscCall(ISDestroy(&globalStratumIS));
1522: PetscCall(ISDestroy(&stratumIS));
1523: PetscCall(PetscViewerHDF5PopGroup(viewer));
1524: }
1525: PetscCall(ISRestoreIndices(globalValueIS, &values));
1526: PetscCall(ISDestroy(&globalValueIS));
1527: PetscCall(ISDestroy(&valueIS));
1528: PetscCall(PetscViewerHDF5PopGroup(viewer));
1529: }
1530: PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
1531: PetscCall(PetscViewerHDF5PopGroup(viewer));
1532: if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
1533: PetscFunctionReturn(PETSC_SUCCESS);
1534: }
1536: /* We only write cells and vertices. Does this screw up parallel reading? */
1537: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
1538: {
1539: IS globalPointNumbers;
1540: PetscViewerFormat format;
1541: PetscBool viz_geom = PETSC_FALSE, xdmf_topo = PETSC_FALSE, petsc_topo = PETSC_FALSE;
1543: PetscFunctionBegin;
1544: PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers));
1545: PetscCall(DMPlexCoordinatesView_HDF5_Internal(dm, viewer));
1547: PetscCall(PetscViewerGetFormat(viewer, &format));
1548: switch (format) {
1549: case PETSC_VIEWER_HDF5_VIZ:
1550: viz_geom = PETSC_TRUE;
1551: xdmf_topo = PETSC_TRUE;
1552: break;
1553: case PETSC_VIEWER_HDF5_XDMF:
1554: xdmf_topo = PETSC_TRUE;
1555: break;
1556: case PETSC_VIEWER_HDF5_PETSC:
1557: petsc_topo = PETSC_TRUE;
1558: break;
1559: case PETSC_VIEWER_DEFAULT:
1560: case PETSC_VIEWER_NATIVE:
1561: viz_geom = PETSC_TRUE;
1562: xdmf_topo = PETSC_TRUE;
1563: petsc_topo = PETSC_TRUE;
1564: break;
1565: default:
1566: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
1567: }
1569: if (viz_geom) PetscCall(DMPlexCoordinatesView_HDF5_XDMF_Private(dm, viewer));
1570: if (xdmf_topo) PetscCall(DMPlexTopologyView_HDF5_XDMF_Private(dm, globalPointNumbers, viewer));
1571: if (petsc_topo) {
1572: PetscBool viewLabels = PETSC_TRUE;
1574: PetscCall(DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer));
1575: PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_view_labels", &viewLabels, NULL));
1576: if (viewLabels) PetscCall(DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer));
1577: }
1579: PetscCall(ISDestroy(&globalPointNumbers));
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
1584: {
1585: DMPlexStorageVersion version;
1586: MPI_Comm comm;
1587: const char *topologydm_name;
1588: const char *sectiondm_name;
1589: PetscSection gsection;
1590: PetscBool ocompress;
1592: PetscFunctionBegin;
1593: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
1594: if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
1595: PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
1596: PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
1597: }
1598: PetscCall(PetscObjectGetComm((PetscObject)sectiondm, &comm));
1599: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1600: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
1601: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1602: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1603: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1604: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1605: PetscCall(DMGetGlobalSection(sectiondm, &gsection));
1606: /* Save raw section */
1607: PetscCall(PetscSectionView(gsection, viewer));
1608: /* Save plex wrapper */
1609: {
1610: PetscInt pStart, pEnd, p, n;
1611: IS globalPointNumbers;
1612: const PetscInt *gpoints;
1613: IS orderIS;
1614: PetscInt *order;
1616: PetscCall(PetscSectionGetChart(gsection, &pStart, &pEnd));
1617: PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers));
1618: PetscCall(ISGetIndices(globalPointNumbers, &gpoints));
1619: for (p = pStart, n = 0; p < pEnd; ++p)
1620: if (gpoints[p] >= 0) n++;
1621: /* "order" is an array of global point numbers.
1622: When loading, it is used with topology/order array
1623: to match section points with plex topology points. */
1624: PetscCall(PetscMalloc1(n, &order));
1625: for (p = pStart, n = 0; p < pEnd; ++p)
1626: if (gpoints[p] >= 0) order[n++] = gpoints[p];
1627: PetscCall(ISRestoreIndices(globalPointNumbers, &gpoints));
1628: PetscCall(ISDestroy(&globalPointNumbers));
1629: PetscCall(ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS));
1630: PetscCall(PetscObjectSetName((PetscObject)orderIS, "order"));
1631: PetscCall(ISView(orderIS, viewer));
1632: PetscCall(ISDestroy(&orderIS));
1633: }
1634: PetscCall(PetscViewerHDF5PopGroup(viewer));
1635: PetscCall(PetscViewerHDF5PopGroup(viewer));
1636: PetscCall(PetscViewerHDF5PopGroup(viewer));
1637: PetscCall(PetscViewerHDF5PopGroup(viewer));
1638: if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
1639: PetscFunctionReturn(PETSC_SUCCESS);
1640: }
1642: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1643: {
1644: const char *topologydm_name;
1645: const char *sectiondm_name;
1646: const char *vec_name;
1647: PetscInt bs;
1649: PetscFunctionBegin;
1650: /* Check consistency */
1651: {
1652: PetscSF pointsf, pointsf1;
1654: PetscCall(DMGetPointSF(dm, &pointsf));
1655: PetscCall(DMGetPointSF(sectiondm, &pointsf1));
1656: PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
1657: }
1658: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1659: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
1660: PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
1661: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1662: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1663: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1664: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1665: PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
1666: PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
1667: PetscCall(VecGetBlockSize(vec, &bs));
1668: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs));
1669: PetscCall(VecSetBlockSize(vec, 1));
1670: /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but, */
1671: /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view */
1672: /* is set to VecView_Plex, which would save vec in a predefined location. */
1673: /* To save vec in where we want, we create a new Vec (temp) with */
1674: /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
1675: {
1676: Vec temp;
1677: const PetscScalar *array;
1678: PetscLayout map;
1680: PetscCall(VecCreate(PetscObjectComm((PetscObject)vec), &temp));
1681: PetscCall(PetscObjectSetName((PetscObject)temp, vec_name));
1682: PetscCall(VecGetLayout(vec, &map));
1683: PetscCall(VecSetLayout(temp, map));
1684: PetscCall(VecSetUp(temp));
1685: PetscCall(VecGetArrayRead(vec, &array));
1686: PetscCall(VecPlaceArray(temp, array));
1687: PetscCall(VecView(temp, viewer));
1688: PetscCall(VecResetArray(temp));
1689: PetscCall(VecRestoreArrayRead(vec, &array));
1690: PetscCall(VecDestroy(&temp));
1691: }
1692: PetscCall(VecSetBlockSize(vec, bs));
1693: PetscCall(PetscViewerHDF5PopGroup(viewer));
1694: PetscCall(PetscViewerHDF5PopGroup(viewer));
1695: PetscCall(PetscViewerHDF5PopGroup(viewer));
1696: PetscCall(PetscViewerHDF5PopGroup(viewer));
1697: PetscCall(PetscViewerHDF5PopGroup(viewer));
1698: PetscCall(PetscViewerHDF5PopGroup(viewer));
1699: PetscFunctionReturn(PETSC_SUCCESS);
1700: }
1702: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1703: {
1704: MPI_Comm comm;
1705: const char *topologydm_name;
1706: const char *sectiondm_name;
1707: const char *vec_name;
1708: PetscSection section;
1709: PetscBool includesConstraints;
1710: Vec gvec;
1711: PetscInt m, bs;
1713: PetscFunctionBegin;
1714: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1715: /* Check consistency */
1716: {
1717: PetscSF pointsf, pointsf1;
1719: PetscCall(DMGetPointSF(dm, &pointsf));
1720: PetscCall(DMGetPointSF(sectiondm, &pointsf1));
1721: PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
1722: }
1723: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1724: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
1725: PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
1726: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1727: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1728: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1729: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1730: PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
1731: PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
1732: PetscCall(VecGetBlockSize(vec, &bs));
1733: PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs));
1734: PetscCall(VecCreate(comm, &gvec));
1735: PetscCall(PetscObjectSetName((PetscObject)gvec, vec_name));
1736: PetscCall(DMGetGlobalSection(sectiondm, §ion));
1737: PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints));
1738: if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m));
1739: else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m));
1740: PetscCall(VecSetSizes(gvec, m, PETSC_DECIDE));
1741: PetscCall(VecSetUp(gvec));
1742: PetscCall(DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec));
1743: PetscCall(DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec));
1744: PetscCall(VecView(gvec, viewer));
1745: PetscCall(VecDestroy(&gvec));
1746: PetscCall(PetscViewerHDF5PopGroup(viewer));
1747: PetscCall(PetscViewerHDF5PopGroup(viewer));
1748: PetscCall(PetscViewerHDF5PopGroup(viewer));
1749: PetscCall(PetscViewerHDF5PopGroup(viewer));
1750: PetscCall(PetscViewerHDF5PopGroup(viewer));
1751: PetscCall(PetscViewerHDF5PopGroup(viewer));
1752: PetscFunctionReturn(PETSC_SUCCESS);
1753: }
1755: struct _n_LoadLabelsCtx {
1756: MPI_Comm comm;
1757: PetscMPIInt rank;
1758: DM dm;
1759: PetscViewer viewer;
1760: DMLabel label;
1761: PetscSF sfXC;
1762: PetscLayout layoutX;
1763: };
1764: typedef struct _n_LoadLabelsCtx *LoadLabelsCtx;
1766: static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx)
1767: {
1768: PetscFunctionBegin;
1769: PetscCall(PetscNew(ctx));
1770: PetscCall(PetscObjectReference((PetscObject)((*ctx)->dm = dm)));
1771: PetscCall(PetscObjectReference((PetscObject)((*ctx)->viewer = viewer)));
1772: PetscCall(PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm));
1773: PetscCallMPI(MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank));
1774: (*ctx)->sfXC = sfXC;
1775: if (sfXC) {
1776: PetscInt nX;
1778: PetscCall(PetscObjectReference((PetscObject)sfXC));
1779: PetscCall(PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL));
1780: PetscCall(PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX));
1781: }
1782: PetscFunctionReturn(PETSC_SUCCESS);
1783: }
1785: static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx)
1786: {
1787: PetscFunctionBegin;
1788: if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
1789: PetscCall(DMDestroy(&(*ctx)->dm));
1790: PetscCall(PetscViewerDestroy(&(*ctx)->viewer));
1791: PetscCall(PetscSFDestroy(&(*ctx)->sfXC));
1792: PetscCall(PetscLayoutDestroy(&(*ctx)->layoutX));
1793: PetscCall(PetscFree(*ctx));
1794: PetscFunctionReturn(PETSC_SUCCESS);
1795: }
1797: /*
1798: A: on-disk points
1799: X: global points [0, NX)
1800: C: distributed plex points
1801: */
1802: static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS)
1803: {
1804: MPI_Comm comm = ctx->comm;
1805: PetscSF sfXC = ctx->sfXC;
1806: PetscLayout layoutX = ctx->layoutX;
1807: PetscSF sfXA;
1808: const PetscInt *A_points;
1809: PetscInt nX, nC;
1810: PetscInt n;
1812: PetscFunctionBegin;
1813: PetscCall(PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL));
1814: PetscCall(ISGetLocalSize(stratumIS, &n));
1815: PetscCall(ISGetIndices(stratumIS, &A_points));
1816: PetscCall(PetscSFCreate(comm, &sfXA));
1817: PetscCall(PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points));
1818: PetscCall(ISCreate(comm, newStratumIS));
1819: PetscCall(ISSetType(*newStratumIS, ISGENERAL));
1820: {
1821: PetscInt i;
1822: PetscBool *A_mask, *X_mask, *C_mask;
1824: PetscCall(PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask));
1825: for (i = 0; i < n; i++) A_mask[i] = PETSC_TRUE;
1826: PetscCall(PetscSFReduceBegin(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE));
1827: PetscCall(PetscSFReduceEnd(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE));
1828: PetscCall(PetscSFBcastBegin(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR));
1829: PetscCall(PetscSFBcastEnd(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR));
1830: PetscCall(ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask));
1831: PetscCall(PetscFree3(A_mask, X_mask, C_mask));
1832: }
1833: PetscCall(PetscSFDestroy(&sfXA));
1834: PetscCall(ISRestoreIndices(stratumIS, &A_points));
1835: PetscFunctionReturn(PETSC_SUCCESS);
1836: }
1838: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data)
1839: {
1840: LoadLabelsCtx ctx = (LoadLabelsCtx)op_data;
1841: PetscViewer viewer = ctx->viewer;
1842: DMLabel label = ctx->label;
1843: MPI_Comm comm = ctx->comm;
1844: IS stratumIS;
1845: const PetscInt *ind;
1846: PetscInt value, N, i;
1848: PetscCall(PetscOptionsStringToInt(vname, &value));
1849: PetscCall(ISCreate(comm, &stratumIS));
1850: PetscCall(PetscObjectSetName((PetscObject)stratumIS, "indices"));
1851: PetscCall(PetscViewerHDF5PushGroup(viewer, vname)); /* labels/<lname>/<vname> */
1853: if (!ctx->sfXC) {
1854: /* Force serial load */
1855: PetscCall(PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N));
1856: PetscCall(PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0));
1857: PetscCall(PetscLayoutSetSize(stratumIS->map, N));
1858: }
1859: PetscCall(ISLoad(stratumIS, viewer));
1861: if (ctx->sfXC) {
1862: IS newStratumIS;
1864: PetscCallHDF5(ReadLabelStratumHDF5_Distribute_Private, (stratumIS, ctx, &newStratumIS));
1865: PetscCall(ISDestroy(&stratumIS));
1866: stratumIS = newStratumIS;
1867: }
1869: PetscCall(PetscViewerHDF5PopGroup(viewer));
1870: PetscCall(ISGetLocalSize(stratumIS, &N));
1871: PetscCall(ISGetIndices(stratumIS, &ind));
1872: for (i = 0; i < N; ++i) PetscCall(DMLabelSetValue(label, ind[i], value));
1873: PetscCall(ISRestoreIndices(stratumIS, &ind));
1874: PetscCall(ISDestroy(&stratumIS));
1875: return 0;
1876: }
1878: /* TODO: Fix this code, it is returning PETSc error codes when it should be translating them to herr_t codes */
1879: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data)
1880: {
1881: LoadLabelsCtx ctx = (LoadLabelsCtx)op_data;
1882: DM dm = ctx->dm;
1883: hsize_t idx = 0;
1884: PetscErrorCode ierr;
1885: PetscBool flg;
1886: herr_t err;
1888: PetscCall(DMHasLabel(dm, lname, &flg));
1889: if (flg) PetscCall(DMRemoveLabel(dm, lname, NULL));
1890: ierr = DMCreateLabel(dm, lname);
1891: if (ierr) return (herr_t)ierr;
1892: ierr = DMGetLabel(dm, lname, &ctx->label);
1893: if (ierr) return (herr_t)ierr;
1894: ierr = PetscViewerHDF5PushGroup(ctx->viewer, lname);
1895: if (ierr) return (herr_t)ierr;
1896: /* Iterate over the label's strata */
1897: PetscCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1898: ierr = PetscViewerHDF5PopGroup(ctx->viewer);
1899: if (ierr) return (herr_t)ierr;
1900: return err;
1901: }
1903: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1904: {
1905: const char *topologydm_name;
1906: LoadLabelsCtx ctx;
1907: hsize_t idx = 0;
1908: char group[PETSC_MAX_PATH_LEN];
1909: DMPlexStorageVersion version;
1910: PetscBool distributed, hasGroup;
1912: PetscFunctionBegin;
1913: PetscCall(DMPlexIsDistributed(dm, &distributed));
1914: if (distributed) PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load");
1915: PetscCall(LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx));
1916: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1917: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
1918: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1919: PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name));
1920: } else {
1921: PetscCall(PetscStrncpy(group, "labels", sizeof(group)));
1922: }
1923: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1924: PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup));
1925: if (hasGroup) {
1926: hid_t fileId, groupId;
1928: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &fileId, &groupId));
1929: /* Iterate over labels */
1930: PetscCallHDF5(H5Literate, (groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx));
1931: PetscCallHDF5(H5Gclose, (groupId));
1932: }
1933: PetscCall(PetscViewerHDF5PopGroup(viewer));
1934: PetscCall(LoadLabelsCtxDestroy(&ctx));
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1938: static PetscErrorCode DMPlexDistributionLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF sf, PetscSF *distsf, DM *distdm)
1939: {
1940: MPI_Comm comm;
1941: PetscMPIInt size, rank;
1942: PetscInt dist_size;
1943: const char *distribution_name;
1944: PetscInt p, lsize;
1945: IS chartSizesIS, ownersIS, gpointsIS;
1946: const PetscInt *chartSize, *owners, *gpoints;
1947: PetscLayout layout;
1948: PetscBool has;
1950: PetscFunctionBegin;
1951: *distsf = NULL;
1952: *distdm = NULL;
1953: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1954: PetscCallMPI(MPI_Comm_size(comm, &size));
1955: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1956: PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
1957: if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS);
1958: PetscCall(PetscLogEventBegin(DMPLEX_DistributionLoad, viewer, 0, 0, 0));
1959: PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &has));
1960: if (!has) {
1961: const char *full_group;
1963: PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &full_group));
1964: 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);
1965: }
1966: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "comm_size", PETSC_INT, NULL, (void *)&dist_size));
1967: 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);
1968: PetscCall(ISCreate(comm, &chartSizesIS));
1969: PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes"));
1970: PetscCall(ISCreate(comm, &ownersIS));
1971: PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners"));
1972: PetscCall(ISCreate(comm, &gpointsIS));
1973: PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers"));
1974: PetscCall(PetscLayoutSetLocalSize(chartSizesIS->map, 1));
1975: PetscCall(ISLoad(chartSizesIS, viewer));
1976: PetscCall(ISGetIndices(chartSizesIS, &chartSize));
1977: PetscCall(PetscLayoutSetLocalSize(ownersIS->map, *chartSize));
1978: PetscCall(PetscLayoutSetLocalSize(gpointsIS->map, *chartSize));
1979: PetscCall(ISLoad(ownersIS, viewer));
1980: PetscCall(ISLoad(gpointsIS, viewer));
1981: PetscCall(ISGetIndices(ownersIS, &owners));
1982: PetscCall(ISGetIndices(gpointsIS, &gpoints));
1983: PetscCall(PetscSFCreate(comm, distsf));
1984: PetscCall(PetscSFSetFromOptions(*distsf));
1985: PetscCall(PetscLayoutCreate(comm, &layout));
1986: PetscCall(PetscSFGetGraph(sf, &lsize, NULL, NULL, NULL));
1987: PetscCall(PetscLayoutSetLocalSize(layout, lsize));
1988: PetscCall(PetscLayoutSetBlockSize(layout, 1));
1989: PetscCall(PetscLayoutSetUp(layout));
1990: PetscCall(PetscSFSetGraphLayout(*distsf, layout, *chartSize, NULL, PETSC_OWN_POINTER, gpoints));
1991: PetscCall(PetscLayoutDestroy(&layout));
1992: /* Migrate DM */
1993: {
1994: PetscInt pStart, pEnd;
1995: PetscSFNode *buffer0, *buffer1, *buffer2;
1997: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1998: PetscCall(PetscMalloc2(pEnd - pStart, &buffer0, lsize, &buffer1));
1999: PetscCall(PetscMalloc1(*chartSize, &buffer2));
2000: {
2001: PetscSF workPointSF;
2002: PetscInt workNroots, workNleaves;
2003: const PetscInt *workIlocal;
2004: const PetscSFNode *workIremote;
2006: for (p = pStart; p < pEnd; ++p) {
2007: buffer0[p - pStart].rank = rank;
2008: buffer0[p - pStart].index = p - pStart;
2009: }
2010: PetscCall(DMGetPointSF(dm, &workPointSF));
2011: PetscCall(PetscSFGetGraph(workPointSF, &workNroots, &workNleaves, &workIlocal, &workIremote));
2012: for (p = 0; p < workNleaves; ++p) {
2013: PetscInt workIlocalp = (workIlocal ? workIlocal[p] : p);
2015: buffer0[workIlocalp].rank = -1;
2016: }
2017: }
2018: for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
2019: for (p = 0; p < *chartSize; ++p) buffer2[p].rank = -1;
2020: PetscCall(PetscSFReduceBegin(sf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC));
2021: PetscCall(PetscSFReduceEnd(sf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC));
2022: PetscCall(PetscSFBcastBegin(*distsf, MPIU_SF_NODE, buffer1, buffer2, MPI_REPLACE));
2023: PetscCall(PetscSFBcastEnd(*distsf, MPIU_SF_NODE, buffer1, buffer2, MPI_REPLACE));
2024: if (PetscDefined(USE_DEBUG)) {
2025: for (p = 0; p < *chartSize; ++p) 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);
2026: }
2027: PetscCall(PetscFree2(buffer0, buffer1));
2028: PetscCall(DMCreate(comm, distdm));
2029: PetscCall(DMSetType(*distdm, DMPLEX));
2030: PetscCall(PetscObjectSetName((PetscObject)*distdm, ((PetscObject)dm)->name));
2031: PetscCall(DMPlexDistributionSetName(*distdm, distribution_name));
2032: {
2033: PetscSF migrationSF;
2035: PetscCall(PetscSFCreate(comm, &migrationSF));
2036: PetscCall(PetscSFSetFromOptions(migrationSF));
2037: PetscCall(PetscSFSetGraph(migrationSF, pEnd - pStart, *chartSize, NULL, PETSC_OWN_POINTER, buffer2, PETSC_OWN_POINTER));
2038: PetscCall(PetscSFSetUp(migrationSF));
2039: PetscCall(DMPlexMigrate(dm, migrationSF, *distdm));
2040: PetscCall(PetscSFDestroy(&migrationSF));
2041: }
2042: }
2043: /* Set pointSF */
2044: {
2045: PetscSF pointSF;
2046: PetscInt *ilocal, nleaves, q;
2047: PetscSFNode *iremote, *buffer0, *buffer1;
2049: PetscCall(PetscMalloc2(*chartSize, &buffer0, lsize, &buffer1));
2050: for (p = 0, nleaves = 0; p < *chartSize; ++p) {
2051: if (owners[p] == rank) {
2052: buffer0[p].rank = rank;
2053: } else {
2054: buffer0[p].rank = -1;
2055: nleaves++;
2056: }
2057: buffer0[p].index = p;
2058: }
2059: for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
2060: PetscCall(PetscSFReduceBegin(*distsf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC));
2061: PetscCall(PetscSFReduceEnd(*distsf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC));
2062: for (p = 0; p < *chartSize; ++p) buffer0[p].rank = -1;
2063: PetscCall(PetscSFBcastBegin(*distsf, MPIU_SF_NODE, buffer1, buffer0, MPI_REPLACE));
2064: PetscCall(PetscSFBcastEnd(*distsf, MPIU_SF_NODE, buffer1, buffer0, MPI_REPLACE));
2065: if (PetscDefined(USE_DEBUG)) {
2066: 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);
2067: }
2068: PetscCall(PetscMalloc1(nleaves, &ilocal));
2069: PetscCall(PetscMalloc1(nleaves, &iremote));
2070: for (p = 0, q = 0; p < *chartSize; ++p) {
2071: if (buffer0[p].rank != rank) {
2072: ilocal[q] = p;
2073: iremote[q].rank = buffer0[p].rank;
2074: iremote[q].index = buffer0[p].index;
2075: q++;
2076: }
2077: }
2078: PetscCheck(q == nleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching leaf sizes: %" PetscInt_FMT " != %" PetscInt_FMT, q, nleaves);
2079: PetscCall(PetscFree2(buffer0, buffer1));
2080: PetscCall(PetscSFCreate(comm, &pointSF));
2081: PetscCall(PetscSFSetFromOptions(pointSF));
2082: PetscCall(PetscSFSetGraph(pointSF, *chartSize, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2083: PetscCall(DMSetPointSF(*distdm, pointSF));
2084: {
2085: DM cdm;
2087: PetscCall(DMGetCoordinateDM(*distdm, &cdm));
2088: PetscCall(DMSetPointSF(cdm, pointSF));
2089: }
2090: PetscCall(PetscSFDestroy(&pointSF));
2091: }
2092: PetscCall(ISRestoreIndices(chartSizesIS, &chartSize));
2093: PetscCall(ISRestoreIndices(ownersIS, &owners));
2094: PetscCall(ISRestoreIndices(gpointsIS, &gpoints));
2095: PetscCall(ISDestroy(&chartSizesIS));
2096: PetscCall(ISDestroy(&ownersIS));
2097: PetscCall(ISDestroy(&gpointsIS));
2098: /* Record that overlap has been manually created. */
2099: /* This is to pass `DMPlexCheckPointSF()`, which checks that */
2100: /* pointSF does not contain cells in the leaves if overlap = 0. */
2101: PetscCall(DMPlexSetOverlap_Plex(*distdm, NULL, DMPLEX_OVERLAP_MANUAL));
2102: PetscCall(DMPlexDistributeSetDefault(*distdm, PETSC_FALSE));
2103: PetscCall(DMPlexReorderSetDefault(*distdm, DM_REORDER_DEFAULT_FALSE));
2104: PetscCall(PetscLogEventEnd(DMPLEX_DistributionLoad, viewer, 0, 0, 0));
2105: PetscFunctionReturn(PETSC_SUCCESS);
2106: }
2108: // Serial load of topology
2109: static PetscErrorCode DMPlexTopologyLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer, PetscSF *sf)
2110: {
2111: MPI_Comm comm;
2112: const char *pointsName, *coneSizesName, *conesName, *orientationsName;
2113: IS pointsIS, coneSizesIS, conesIS, orientationsIS;
2114: const PetscInt *points, *coneSizes, *cones, *orientations;
2115: PetscInt *cone, *ornt;
2116: PetscInt dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
2117: PetscMPIInt size, rank;
2119: PetscFunctionBegin;
2120: pointsName = "order";
2121: coneSizesName = "cones";
2122: conesName = "cells";
2123: orientationsName = "orientation";
2124: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2125: PetscCallMPI(MPI_Comm_size(comm, &size));
2126: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2127: PetscCall(ISCreate(comm, &pointsIS));
2128: PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName));
2129: PetscCall(ISCreate(comm, &coneSizesIS));
2130: PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
2131: PetscCall(ISCreate(comm, &conesIS));
2132: PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
2133: PetscCall(ISCreate(comm, &orientationsIS));
2134: PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
2135: PetscCall(PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject)conesIS, "cell_dim", PETSC_INT, NULL, &dim));
2136: PetscCall(DMSetDimension(dm, dim));
2137: {
2138: /* Force serial load */
2139: PetscCall(PetscInfo(dm, "Loading DM %s in serial\n", dm->hdr.name));
2140: PetscCall(PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np));
2141: PetscCall(PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0));
2142: PetscCall(PetscLayoutSetSize(pointsIS->map, Np));
2143: pEnd = rank == 0 ? Np : 0;
2144: PetscCall(PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np));
2145: PetscCall(PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0));
2146: PetscCall(PetscLayoutSetSize(coneSizesIS->map, Np));
2147: PetscCall(PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N));
2148: PetscCall(PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0));
2149: PetscCall(PetscLayoutSetSize(conesIS->map, N));
2150: PetscCall(PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N));
2151: PetscCall(PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0));
2152: PetscCall(PetscLayoutSetSize(orientationsIS->map, N));
2153: }
2154: PetscCall(ISLoad(pointsIS, viewer));
2155: PetscCall(ISLoad(coneSizesIS, viewer));
2156: PetscCall(ISLoad(conesIS, viewer));
2157: PetscCall(ISLoad(orientationsIS, viewer));
2158: /* Create Plex */
2159: PetscCall(DMPlexSetChart(dm, 0, pEnd));
2160: PetscCall(ISGetIndices(pointsIS, &points));
2161: PetscCall(ISGetIndices(coneSizesIS, &coneSizes));
2162: for (p = 0; p < pEnd; ++p) {
2163: PetscCall(DMPlexSetConeSize(dm, points[p], coneSizes[p]));
2164: maxConeSize = PetscMax(maxConeSize, coneSizes[p]);
2165: }
2166: PetscCall(DMSetUp(dm));
2167: PetscCall(ISGetIndices(conesIS, &cones));
2168: PetscCall(ISGetIndices(orientationsIS, &orientations));
2169: PetscCall(PetscMalloc2(maxConeSize, &cone, maxConeSize, &ornt));
2170: for (p = 0, q = 0; p < pEnd; ++p) {
2171: for (c = 0; c < coneSizes[p]; ++c, ++q) {
2172: cone[c] = cones[q];
2173: ornt[c] = orientations[q];
2174: }
2175: PetscCall(DMPlexSetCone(dm, points[p], cone));
2176: PetscCall(DMPlexSetConeOrientation(dm, points[p], ornt));
2177: }
2178: PetscCall(PetscFree2(cone, ornt));
2179: /* Create global section migration SF */
2180: if (sf) {
2181: PetscLayout layout;
2182: PetscInt *globalIndices;
2184: PetscCall(PetscMalloc1(pEnd, &globalIndices));
2185: /* plex point == globalPointNumber in this case */
2186: for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
2187: PetscCall(PetscLayoutCreate(comm, &layout));
2188: PetscCall(PetscLayoutSetSize(layout, Np));
2189: PetscCall(PetscLayoutSetBlockSize(layout, 1));
2190: PetscCall(PetscLayoutSetUp(layout));
2191: PetscCall(PetscSFCreate(comm, sf));
2192: PetscCall(PetscSFSetFromOptions(*sf));
2193: PetscCall(PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices));
2194: PetscCall(PetscLayoutDestroy(&layout));
2195: PetscCall(PetscFree(globalIndices));
2196: }
2197: /* Clean-up */
2198: PetscCall(ISRestoreIndices(pointsIS, &points));
2199: PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes));
2200: PetscCall(ISRestoreIndices(conesIS, &cones));
2201: PetscCall(ISRestoreIndices(orientationsIS, &orientations));
2202: PetscCall(ISDestroy(&pointsIS));
2203: PetscCall(ISDestroy(&coneSizesIS));
2204: PetscCall(ISDestroy(&conesIS));
2205: PetscCall(ISDestroy(&orientationsIS));
2206: /* Fill in the rest of the topology structure */
2207: PetscCall(DMPlexSymmetrize(dm));
2208: PetscCall(DMPlexStratify(dm));
2209: PetscFunctionReturn(PETSC_SUCCESS);
2210: }
2212: /* Representation of two DMPlex strata in 0-based global numbering */
2213: struct _n_PlexLayer {
2214: PetscInt d;
2215: IS conesIS, orientationsIS;
2216: PetscSection coneSizesSection;
2217: PetscLayout vertexLayout;
2218: PetscSF overlapSF, l2gSF; //TODO maybe confusing names (in DMPlex in general)
2219: PetscInt offset, conesOffset, leafOffset;
2220: };
2221: typedef struct _n_PlexLayer *PlexLayer;
2223: static PetscErrorCode PlexLayerDestroy(PlexLayer *layer)
2224: {
2225: PetscFunctionBegin;
2226: if (!*layer) PetscFunctionReturn(PETSC_SUCCESS);
2227: PetscCall(PetscSectionDestroy(&(*layer)->coneSizesSection));
2228: PetscCall(ISDestroy(&(*layer)->conesIS));
2229: PetscCall(ISDestroy(&(*layer)->orientationsIS));
2230: PetscCall(PetscSFDestroy(&(*layer)->overlapSF));
2231: PetscCall(PetscSFDestroy(&(*layer)->l2gSF));
2232: PetscCall(PetscLayoutDestroy(&(*layer)->vertexLayout));
2233: PetscCall(PetscFree(*layer));
2234: PetscFunctionReturn(PETSC_SUCCESS);
2235: }
2237: static PetscErrorCode PlexLayerCreate_Private(PlexLayer *layer)
2238: {
2239: PetscFunctionBegin;
2240: PetscCall(PetscNew(layer));
2241: (*layer)->d = -1;
2242: (*layer)->offset = -1;
2243: (*layer)->conesOffset = -1;
2244: (*layer)->leafOffset = -1;
2245: PetscFunctionReturn(PETSC_SUCCESS);
2246: }
2248: // Parallel load of a depth stratum
2249: static PetscErrorCode PlexLayerLoad_Private(PlexLayer layer, PetscViewer viewer, PetscInt d, PetscLayout pointsLayout)
2250: {
2251: char path[128];
2252: MPI_Comm comm;
2253: const char *coneSizesName, *conesName, *orientationsName;
2254: IS coneSizesIS, conesIS, orientationsIS;
2255: PetscSection coneSizesSection;
2256: PetscLayout vertexLayout = NULL;
2257: PetscInt s;
2259: PetscFunctionBegin;
2260: coneSizesName = "cone_sizes";
2261: conesName = "cones";
2262: orientationsName = "orientations";
2263: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
2265: /* query size of next lower depth stratum (next lower dimension) */
2266: if (d > 0) {
2267: PetscInt NVertices;
2269: PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT "/%s", d - 1, coneSizesName));
2270: PetscCall(PetscViewerHDF5ReadSizes(viewer, path, NULL, &NVertices));
2271: PetscCall(PetscLayoutCreate(comm, &vertexLayout));
2272: PetscCall(PetscLayoutSetSize(vertexLayout, NVertices));
2273: PetscCall(PetscLayoutSetUp(vertexLayout));
2274: }
2276: PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT, d));
2277: PetscCall(PetscViewerHDF5PushGroup(viewer, path));
2279: /* create coneSizesSection from stored IS coneSizes */
2280: {
2281: const PetscInt *coneSizes;
2283: PetscCall(ISCreate(comm, &coneSizesIS));
2284: PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
2285: if (pointsLayout) PetscCall(ISSetLayout(coneSizesIS, pointsLayout));
2286: PetscCall(ISLoad(coneSizesIS, viewer));
2287: if (!pointsLayout) PetscCall(ISGetLayout(coneSizesIS, &pointsLayout));
2288: PetscCall(ISGetIndices(coneSizesIS, &coneSizes));
2289: PetscCall(PetscSectionCreate(comm, &coneSizesSection));
2290: //TODO different start ?
2291: PetscCall(PetscSectionSetChart(coneSizesSection, 0, pointsLayout->n));
2292: for (s = 0; s < pointsLayout->n; ++s) PetscCall(PetscSectionSetDof(coneSizesSection, s, coneSizes[s]));
2293: PetscCall(PetscSectionSetUp(coneSizesSection));
2294: PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes));
2295: {
2296: PetscLayout tmp = NULL;
2298: /* We need to keep the layout until the end of function */
2299: PetscCall(PetscLayoutReference((PetscLayout)pointsLayout, &tmp));
2300: }
2301: PetscCall(ISDestroy(&coneSizesIS));
2302: }
2304: /* use value layout of coneSizesSection as layout of cones and orientations */
2305: {
2306: PetscLayout conesLayout;
2308: PetscCall(PetscSectionGetValueLayout(comm, coneSizesSection, &conesLayout));
2309: PetscCall(ISCreate(comm, &conesIS));
2310: PetscCall(ISCreate(comm, &orientationsIS));
2311: PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
2312: PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
2313: PetscCall(PetscLayoutDuplicate(conesLayout, &conesIS->map));
2314: PetscCall(PetscLayoutDuplicate(conesLayout, &orientationsIS->map));
2315: PetscCall(ISLoad(conesIS, viewer));
2316: PetscCall(ISLoad(orientationsIS, viewer));
2317: PetscCall(PetscLayoutDestroy(&conesLayout));
2318: }
2320: /* check assertion that layout of points is the same as point layout of coneSizesSection */
2321: {
2322: PetscLayout pointsLayout0;
2323: PetscBool flg;
2325: PetscCall(PetscSectionGetPointLayout(comm, coneSizesSection, &pointsLayout0));
2326: PetscCall(PetscLayoutCompare(pointsLayout, pointsLayout0, &flg));
2327: PetscCheck(flg, comm, PETSC_ERR_PLIB, "points layout != coneSizesSection point layout");
2328: PetscCall(PetscLayoutDestroy(&pointsLayout0));
2329: }
2330: PetscCall(PetscViewerHDF5PopGroup(viewer));
2331: PetscCall(PetscLayoutDestroy(&pointsLayout));
2333: layer->d = d;
2334: layer->conesIS = conesIS;
2335: layer->coneSizesSection = coneSizesSection;
2336: layer->orientationsIS = orientationsIS;
2337: layer->vertexLayout = vertexLayout;
2338: PetscFunctionReturn(PETSC_SUCCESS);
2339: }
2341: static PetscErrorCode PlexLayerDistribute_Private(PlexLayer layer, PetscSF cellLocalToGlobalSF)
2342: {
2343: IS newConesIS, newOrientationsIS;
2344: PetscSection newConeSizesSection;
2345: MPI_Comm comm;
2347: PetscFunctionBegin;
2348: PetscCall(PetscObjectGetComm((PetscObject)cellLocalToGlobalSF, &comm));
2349: PetscCall(PetscSectionCreate(comm, &newConeSizesSection));
2350: //TODO rename to something like ISDistribute() with optional PetscSection argument
2351: PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->conesIS, newConeSizesSection, &newConesIS));
2352: PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->orientationsIS, newConeSizesSection, &newOrientationsIS));
2354: PetscCall(PetscObjectSetName((PetscObject)newConeSizesSection, ((PetscObject)layer->coneSizesSection)->name));
2355: PetscCall(PetscObjectSetName((PetscObject)newConesIS, ((PetscObject)layer->conesIS)->name));
2356: PetscCall(PetscObjectSetName((PetscObject)newOrientationsIS, ((PetscObject)layer->orientationsIS)->name));
2357: PetscCall(PetscSectionDestroy(&layer->coneSizesSection));
2358: PetscCall(ISDestroy(&layer->conesIS));
2359: PetscCall(ISDestroy(&layer->orientationsIS));
2360: layer->coneSizesSection = newConeSizesSection;
2361: layer->conesIS = newConesIS;
2362: layer->orientationsIS = newOrientationsIS;
2363: PetscFunctionReturn(PETSC_SUCCESS);
2364: }
2366: //TODO share code with DMPlexBuildFromCellListParallel()
2367: #include <petsc/private/hashseti.h>
2368: static PetscErrorCode PlexLayerCreateSFs_Private(PlexLayer layer, PetscSF *vertexOverlapSF, PetscSF *sfXC)
2369: {
2370: PetscLayout vertexLayout = layer->vertexLayout;
2371: PetscSection coneSection = layer->coneSizesSection;
2372: IS cellVertexData = layer->conesIS;
2373: IS coneOrientations = layer->orientationsIS;
2374: PetscSF vl2gSF, vOverlapSF;
2375: PetscInt *verticesAdj;
2376: PetscInt i, n, numVerticesAdj;
2377: const PetscInt *cvd, *co = NULL;
2378: MPI_Comm comm;
2380: PetscFunctionBegin;
2381: PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm));
2382: PetscCall(PetscSectionGetStorageSize(coneSection, &n));
2383: {
2384: PetscInt n0;
2386: PetscCall(ISGetLocalSize(cellVertexData, &n0));
2387: 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);
2388: PetscCall(ISGetIndices(cellVertexData, &cvd));
2389: }
2390: if (coneOrientations) {
2391: PetscInt n0;
2393: PetscCall(ISGetLocalSize(coneOrientations, &n0));
2394: 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);
2395: PetscCall(ISGetIndices(coneOrientations, &co));
2396: }
2397: /* Get/check global number of vertices */
2398: {
2399: PetscInt NVerticesInCells = PETSC_INT_MIN;
2401: /* NVerticesInCells = max(cellVertexData) + 1 */
2402: for (i = 0; i < n; i++)
2403: if (cvd[i] > NVerticesInCells) NVerticesInCells = cvd[i];
2404: ++NVerticesInCells;
2405: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, comm));
2407: if (vertexLayout->n == PETSC_DECIDE && vertexLayout->N == PETSC_DECIDE) vertexLayout->N = NVerticesInCells;
2408: else
2409: 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,
2410: vertexLayout->N, NVerticesInCells);
2411: PetscCall(PetscLayoutSetUp(vertexLayout));
2412: }
2413: /* Find locally unique vertices in cellVertexData */
2414: {
2415: PetscHSetI vhash;
2416: PetscInt off = 0;
2418: PetscCall(PetscHSetICreate(&vhash));
2419: for (i = 0; i < n; ++i) PetscCall(PetscHSetIAdd(vhash, cvd[i]));
2420: PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj));
2421: PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj));
2422: PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj));
2423: PetscAssert(off == numVerticesAdj, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assertion failed: off == numVerticesAdj (%" PetscInt_FMT " != %" PetscInt_FMT ")", off, numVerticesAdj);
2424: PetscCall(PetscHSetIDestroy(&vhash));
2425: }
2426: /* We must sort vertices to preserve numbering */
2427: PetscCall(PetscSortInt(numVerticesAdj, verticesAdj));
2428: /* Connect locally unique vertices with star forest */
2429: PetscCall(PetscSFCreateByMatchingIndices(vertexLayout, numVerticesAdj, verticesAdj, NULL, 0, numVerticesAdj, verticesAdj, NULL, 0, &vl2gSF, &vOverlapSF));
2430: PetscCall(PetscObjectSetName((PetscObject)vOverlapSF, "overlapSF"));
2431: PetscCall(PetscObjectSetName((PetscObject)vl2gSF, "localToGlobalSF"));
2433: PetscCall(PetscFree(verticesAdj));
2434: *vertexOverlapSF = vOverlapSF;
2435: *sfXC = vl2gSF;
2436: PetscFunctionReturn(PETSC_SUCCESS);
2437: }
2439: static PetscErrorCode PlexLayerCreateCellSFs_Private(PlexLayer layer, PetscSF *cellOverlapSF, PetscSF *cellLocalToGlobalSF)
2440: {
2441: PetscSection coneSection = layer->coneSizesSection;
2442: PetscInt nCells;
2443: MPI_Comm comm;
2445: PetscFunctionBegin;
2446: PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm));
2447: {
2448: PetscInt cStart;
2450: PetscCall(PetscSectionGetChart(coneSection, &cStart, &nCells));
2451: PetscCheck(cStart == 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "coneSection must start at 0");
2452: }
2453: /* Create overlapSF as empty SF with the right number of roots */
2454: PetscCall(PetscSFCreate(comm, cellOverlapSF));
2455: PetscCall(PetscSFSetGraph(*cellOverlapSF, nCells, 0, NULL, PETSC_USE_POINTER, NULL, PETSC_USE_POINTER));
2456: PetscCall(PetscSFSetUp(*cellOverlapSF));
2457: /* Create localToGlobalSF as identity mapping */
2458: {
2459: PetscLayout map;
2461: PetscCall(PetscLayoutCreateFromSizes(comm, nCells, PETSC_DECIDE, 1, &map));
2462: PetscCall(PetscSFCreateFromLayouts(map, map, cellLocalToGlobalSF));
2463: PetscCall(PetscSFSetUp(*cellLocalToGlobalSF));
2464: PetscCall(PetscLayoutDestroy(&map));
2465: }
2466: PetscFunctionReturn(PETSC_SUCCESS);
2467: }
2469: static PetscErrorCode DMPlexTopologyBuildFromLayers_Private(DM dm, PetscInt depth, PlexLayer *layers, IS strataPermutation)
2470: {
2471: const PetscInt *permArr;
2472: PetscInt d, nPoints;
2473: MPI_Comm comm;
2475: PetscFunctionBegin;
2476: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2477: PetscCall(ISGetIndices(strataPermutation, &permArr));
2479: /* Count points, strata offsets and cones offsets (taking strataPermutation into account) */
2480: {
2481: PetscInt stratumOffset = 0;
2482: PetscInt conesOffset = 0;
2484: for (d = 0; d <= depth; d++) {
2485: const PetscInt e = permArr[d];
2486: const PlexLayer l = layers[e];
2487: PetscInt lo, n, size;
2489: PetscCall(PetscSectionGetChart(l->coneSizesSection, &lo, &n));
2490: PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &size));
2491: PetscCheck(lo == 0, comm, PETSC_ERR_PLIB, "starting point should be 0 in coneSizesSection %" PetscInt_FMT, d);
2492: l->offset = stratumOffset;
2493: l->conesOffset = conesOffset;
2494: stratumOffset += n;
2495: conesOffset += size;
2496: }
2497: nPoints = stratumOffset;
2498: }
2500: /* Set interval for all plex points */
2501: //TODO we should store starting point of plex
2502: PetscCall(DMPlexSetChart(dm, 0, nPoints));
2504: /* Set up plex coneSection from layer coneSections */
2505: {
2506: PetscSection coneSection;
2508: PetscCall(DMPlexGetConeSection(dm, &coneSection));
2509: for (d = 0; d <= depth; d++) {
2510: const PlexLayer l = layers[d];
2511: PetscInt n, q;
2513: PetscCall(PetscSectionGetChart(l->coneSizesSection, NULL, &n));
2514: for (q = 0; q < n; q++) {
2515: const PetscInt p = l->offset + q;
2516: PetscInt coneSize;
2518: PetscCall(PetscSectionGetDof(l->coneSizesSection, q, &coneSize));
2519: PetscCall(PetscSectionSetDof(coneSection, p, coneSize));
2520: }
2521: }
2522: }
2523: //TODO this is terrible, DMSetUp_Plex() should be DMPlexSetUpSections() or so
2524: PetscCall(DMSetUp(dm));
2526: /* Renumber cones points from layer-global numbering to plex-local numbering */
2527: {
2528: PetscInt *cones, *ornts;
2530: PetscCall(DMPlexGetCones(dm, &cones));
2531: PetscCall(DMPlexGetConeOrientations(dm, &ornts));
2532: for (d = 1; d <= depth; d++) {
2533: const PlexLayer l = layers[d];
2534: PetscInt i, lConesSize;
2535: PetscInt *lCones;
2536: const PetscInt *lOrnts;
2537: PetscInt *pCones = &cones[l->conesOffset];
2538: PetscInt *pOrnts = &ornts[l->conesOffset];
2540: PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &lConesSize));
2541: /* Get cones in local plex numbering */
2542: {
2543: ISLocalToGlobalMapping l2g;
2544: PetscLayout vertexLayout = l->vertexLayout;
2545: PetscSF vertexSF = layers[d - 1]->l2gSF; /* vertices of this layer are cells of previous layer */
2546: const PetscInt *gCones;
2547: PetscInt lConesSize0;
2549: PetscCall(ISGetLocalSize(l->conesIS, &lConesSize0));
2550: PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(conesIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize);
2551: PetscCall(ISGetLocalSize(l->orientationsIS, &lConesSize0));
2552: PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(orientationsIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize);
2554: PetscCall(PetscMalloc1(lConesSize, &lCones));
2555: PetscCall(ISGetIndices(l->conesIS, &gCones));
2556: PetscCall(ISLocalToGlobalMappingCreateSF(vertexSF, vertexLayout->rstart, &l2g));
2557: PetscCall(ISGlobalToLocalMappingApply(l2g, IS_GTOLM_MASK, lConesSize, gCones, &lConesSize0, lCones));
2558: PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "global to local does not cover all indices (%" PetscInt_FMT " of %" PetscInt_FMT ")", lConesSize0, lConesSize);
2559: PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
2560: PetscCall(ISRestoreIndices(l->conesIS, &gCones));
2561: }
2562: PetscCall(ISGetIndices(l->orientationsIS, &lOrnts));
2563: /* Set cones, need to add stratum offset */
2564: for (i = 0; i < lConesSize; i++) {
2565: pCones[i] = lCones[i] + layers[d - 1]->offset; /* cone points of current layer are points of previous layer */
2566: pOrnts[i] = lOrnts[i];
2567: }
2568: PetscCall(PetscFree(lCones));
2569: PetscCall(ISRestoreIndices(l->orientationsIS, &lOrnts));
2570: }
2571: }
2572: PetscCall(DMPlexSymmetrize(dm));
2573: PetscCall(DMPlexStratify(dm));
2574: PetscCall(ISRestoreIndices(strataPermutation, &permArr));
2575: PetscFunctionReturn(PETSC_SUCCESS);
2576: }
2578: static PetscErrorCode PlexLayerConcatenateSFs_Private(MPI_Comm comm, PetscInt depth, PlexLayer layers[], IS strataPermutation, PetscSF *overlapSF, PetscSF *l2gSF)
2579: {
2580: PetscInt d;
2581: PetscSF *osfs, *lsfs;
2582: PetscInt *leafOffsets;
2583: const PetscInt *permArr;
2585: PetscFunctionBegin;
2586: PetscCall(ISGetIndices(strataPermutation, &permArr));
2587: PetscCall(PetscCalloc3(depth + 1, &osfs, depth + 1, &lsfs, depth + 1, &leafOffsets));
2588: for (d = 0; d <= depth; d++) {
2589: const PetscInt e = permArr[d];
2591: PetscAssert(e == layers[e]->d, PETSC_COMM_SELF, PETSC_ERR_PLIB, "assertion: e == layers[e]->d");
2592: osfs[d] = layers[e]->overlapSF;
2593: lsfs[d] = layers[e]->l2gSF;
2594: leafOffsets[d] = layers[e]->offset;
2595: }
2596: PetscCall(PetscSFConcatenate(comm, depth + 1, osfs, PETSCSF_CONCATENATE_ROOTMODE_LOCAL, leafOffsets, overlapSF));
2597: PetscCall(PetscSFConcatenate(comm, depth + 1, lsfs, PETSCSF_CONCATENATE_ROOTMODE_GLOBAL, leafOffsets, l2gSF));
2598: PetscCall(PetscFree3(osfs, lsfs, leafOffsets));
2599: PetscCall(ISRestoreIndices(strataPermutation, &permArr));
2600: PetscFunctionReturn(PETSC_SUCCESS);
2601: }
2603: // Parallel load of topology
2604: static PetscErrorCode DMPlexTopologyLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF *sfXC)
2605: {
2606: PlexLayer *layers;
2607: IS strataPermutation;
2608: PetscLayout pointsLayout;
2609: PetscInt depth;
2610: PetscInt d;
2611: MPI_Comm comm;
2613: PetscFunctionBegin;
2614: {
2615: PetscInt dim;
2617: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "depth", PETSC_INT, NULL, &depth));
2618: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "cell_dim", PETSC_INT, NULL, &dim));
2619: PetscCall(DMSetDimension(dm, dim));
2620: }
2621: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2623: PetscCall(PetscInfo(dm, "Loading DM %s in parallel\n", dm->hdr.name));
2624: {
2625: IS spOnComm;
2627: PetscCall(ISCreate(comm, &spOnComm));
2628: PetscCall(PetscObjectSetName((PetscObject)spOnComm, "permutation"));
2629: PetscCall(ISLoad(spOnComm, viewer));
2630: /* have the same serial IS on every rank */
2631: PetscCall(ISAllGather(spOnComm, &strataPermutation));
2632: PetscCall(PetscObjectSetName((PetscObject)strataPermutation, ((PetscObject)spOnComm)->name));
2633: PetscCall(ISDestroy(&spOnComm));
2634: }
2636: /* Create layers, load raw data for each layer */
2637: PetscCall(PetscViewerHDF5PushGroup(viewer, "strata"));
2638: PetscCall(PetscMalloc1(depth + 1, &layers));
2639: for (d = depth, pointsLayout = NULL; d >= 0; pointsLayout = layers[d]->vertexLayout, d--) {
2640: PetscCall(PlexLayerCreate_Private(&layers[d]));
2641: PetscCall(PlexLayerLoad_Private(layers[d], viewer, d, pointsLayout));
2642: }
2643: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */
2645: for (d = depth; d >= 0; d--) {
2646: /* Redistribute cells and vertices for each applicable layer */
2647: if (d < depth) PetscCall(PlexLayerDistribute_Private(layers[d], layers[d]->l2gSF));
2648: /* Create vertex overlap SF and vertex localToGlobal SF for each applicable layer */
2649: if (d > 0) PetscCall(PlexLayerCreateSFs_Private(layers[d], &layers[d - 1]->overlapSF, &layers[d - 1]->l2gSF));
2650: }
2651: /* Build trivial SFs for the cell layer as well */
2652: PetscCall(PlexLayerCreateCellSFs_Private(layers[depth], &layers[depth]->overlapSF, &layers[depth]->l2gSF));
2654: /* Build DMPlex topology from the layers */
2655: PetscCall(DMPlexTopologyBuildFromLayers_Private(dm, depth, layers, strataPermutation));
2657: /* Build overall point SF alias overlap SF */
2658: {
2659: PetscSF overlapSF;
2661: PetscCall(PlexLayerConcatenateSFs_Private(comm, depth, layers, strataPermutation, &overlapSF, sfXC));
2662: PetscCall(DMSetPointSF(dm, overlapSF));
2663: PetscCall(PetscSFDestroy(&overlapSF));
2664: }
2666: for (d = depth; d >= 0; d--) PetscCall(PlexLayerDestroy(&layers[d]));
2667: PetscCall(PetscFree(layers));
2668: PetscCall(ISDestroy(&strataPermutation));
2669: PetscFunctionReturn(PETSC_SUCCESS);
2670: }
2672: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sfXC)
2673: {
2674: DMPlexStorageVersion version;
2675: const char *topologydm_name;
2676: char group[PETSC_MAX_PATH_LEN];
2677: PetscSF sfwork = NULL;
2679: PetscFunctionBegin;
2680: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2681: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2682: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
2683: PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name));
2684: } else {
2685: PetscCall(PetscStrncpy(group, "/", sizeof(group)));
2686: }
2687: PetscCall(PetscViewerHDF5PushGroup(viewer, group));
2689: PetscCall(PetscViewerHDF5PushGroup(viewer, "topology"));
2690: if (version->major < 3) {
2691: PetscCall(DMPlexTopologyLoad_HDF5_Legacy_Private(dm, viewer, &sfwork));
2692: } else {
2693: /* since DMPlexStorageVersion 3.0.0 */
2694: PetscCall(DMPlexTopologyLoad_HDF5_Private(dm, viewer, &sfwork));
2695: }
2696: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */
2698: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
2699: DM distdm;
2700: PetscSF distsf;
2701: const char *distribution_name;
2702: PetscBool exists;
2704: PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
2705: /* this group is guaranteed to be present since DMPlexStorageVersion 2.1.0 */
2706: PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions"));
2707: PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &exists));
2708: if (exists) {
2709: PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name));
2710: PetscCall(DMPlexDistributionLoad_HDF5_Private(dm, viewer, sfwork, &distsf, &distdm));
2711: if (distdm) {
2712: PetscCall(DMPlexReplace_Internal(dm, &distdm));
2713: PetscCall(PetscSFDestroy(&sfwork));
2714: sfwork = distsf;
2715: }
2716: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */
2717: }
2718: PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */
2719: }
2720: if (sfXC) {
2721: *sfXC = sfwork;
2722: } else {
2723: PetscCall(PetscSFDestroy(&sfwork));
2724: }
2726: PetscCall(PetscViewerHDF5PopGroup(viewer));
2727: PetscFunctionReturn(PETSC_SUCCESS);
2728: }
2730: /* If the file is old, it not only has different path to the coordinates, but */
2731: /* does not contain coordinateDMs, so must fall back to the old implementation. */
2732: static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
2733: {
2734: PetscSection coordSection;
2735: Vec coordinates;
2736: PetscReal lengthScale;
2737: PetscInt spatialDim, N, numVertices, vStart, vEnd, v;
2738: PetscMPIInt rank;
2740: PetscFunctionBegin;
2741: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
2742: /* Read geometry */
2743: PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry"));
2744: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates));
2745: PetscCall(PetscObjectSetName((PetscObject)coordinates, "vertices"));
2746: {
2747: /* Force serial load */
2748: PetscCall(PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N));
2749: PetscCall(VecSetSizes(coordinates, rank == 0 ? N : 0, N));
2750: PetscCall(VecSetBlockSize(coordinates, spatialDim));
2751: }
2752: PetscCall(VecLoad(coordinates, viewer));
2753: PetscCall(PetscViewerHDF5PopGroup(viewer));
2754: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
2755: PetscCall(VecScale(coordinates, 1.0 / lengthScale));
2756: PetscCall(VecGetLocalSize(coordinates, &numVertices));
2757: PetscCall(VecGetBlockSize(coordinates, &spatialDim));
2758: numVertices /= spatialDim;
2759: /* Create coordinates */
2760: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
2761: 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);
2762: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2763: PetscCall(PetscSectionSetNumFields(coordSection, 1));
2764: PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spatialDim));
2765: PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd));
2766: for (v = vStart; v < vEnd; ++v) {
2767: PetscCall(PetscSectionSetDof(coordSection, v, spatialDim));
2768: PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spatialDim));
2769: }
2770: PetscCall(PetscSectionSetUp(coordSection));
2771: PetscCall(DMSetCoordinates(dm, coordinates));
2772: PetscCall(VecDestroy(&coordinates));
2773: PetscFunctionReturn(PETSC_SUCCESS);
2774: }
2776: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
2777: {
2778: DMPlexStorageVersion version;
2779: DM cdm;
2780: Vec coords;
2781: PetscInt blockSize;
2782: PetscReal lengthScale;
2783: PetscSF lsf;
2784: const char *topologydm_name;
2785: char *coordinatedm_name, *coordinates_name;
2787: PetscFunctionBegin;
2788: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2789: if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
2790: PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer));
2791: PetscFunctionReturn(PETSC_SUCCESS);
2792: }
2793: /* else: since DMPlexStorageVersion 2.0.0 */
2794: PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load");
2795: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2796: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
2797: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
2798: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, NULL, &coordinatedm_name));
2799: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, NULL, &coordinates_name));
2800: PetscCall(PetscViewerHDF5PopGroup(viewer));
2801: PetscCall(PetscViewerHDF5PopGroup(viewer));
2802: PetscCall(DMGetCoordinateDM(dm, &cdm));
2803: PetscCall(PetscObjectSetName((PetscObject)cdm, coordinatedm_name));
2804: PetscCall(PetscFree(coordinatedm_name));
2805: /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
2806: PetscCall(DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf));
2807: PetscCall(DMCreateLocalVector(cdm, &coords));
2808: PetscCall(PetscObjectSetName((PetscObject)coords, coordinates_name));
2809: PetscCall(PetscFree(coordinates_name));
2810: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
2811: PetscCall(DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords));
2812: PetscCall(PetscViewerPopFormat(viewer));
2813: PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
2814: PetscCall(VecScale(coords, 1.0 / lengthScale));
2815: PetscCall(DMSetCoordinatesLocal(dm, coords));
2816: PetscCall(VecGetBlockSize(coords, &blockSize));
2817: PetscCall(DMSetCoordinateDim(dm, blockSize));
2818: PetscCall(VecDestroy(&coords));
2819: PetscCall(PetscSFDestroy(&lsf));
2820: PetscFunctionReturn(PETSC_SUCCESS);
2821: }
2823: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
2824: {
2825: DMPlexStorageVersion version;
2827: PetscFunctionBegin;
2828: PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2829: PetscCall(PetscInfo(dm, "Loading DM %s storage version %d.%d.%d\n", dm->hdr.name, version->major, version->minor, version->subminor));
2830: if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
2831: PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL));
2832: PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL));
2833: PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer));
2834: } else {
2835: PetscSF sfXC;
2837: /* since DMPlexStorageVersion 2.0.0 */
2838: PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC));
2839: PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC));
2840: PetscCall(DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC));
2841: PetscCall(PetscSFDestroy(&sfXC));
2842: }
2843: PetscFunctionReturn(PETSC_SUCCESS);
2844: }
2846: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2847: {
2848: MPI_Comm comm;
2849: PetscInt pStart, pEnd, p, m;
2850: PetscInt *goffs, *ilocal;
2851: PetscBool rootIncludeConstraints, leafIncludeConstraints;
2853: PetscFunctionBegin;
2854: PetscCall(PetscObjectGetComm((PetscObject)leafSection, &comm));
2855: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2856: PetscCall(PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints));
2857: PetscCall(PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints));
2858: if (rootIncludeConstraints && leafIncludeConstraints) PetscCall(PetscSectionGetStorageSize(leafSection, &m));
2859: else PetscCall(PetscSectionGetConstrainedStorageSize(leafSection, &m));
2860: PetscCall(PetscMalloc1(m, &ilocal));
2861: PetscCall(PetscMalloc1(m, &goffs));
2862: /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
2863: /* for the top-level section (not for each field), so one must have */
2864: /* rootSection->pointMajor == PETSC_TRUE. */
2865: PetscCheck(rootSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
2866: /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
2867: PetscCheck(leafSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
2868: for (p = pStart, m = 0; p < pEnd; ++p) {
2869: PetscInt dof, cdof, i, j, off, goff;
2870: const PetscInt *cinds;
2872: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2873: if (dof < 0) continue;
2874: goff = globalOffsets[p - pStart];
2875: PetscCall(PetscSectionGetOffset(leafSection, p, &off));
2876: PetscCall(PetscSectionGetConstraintDof(leafSection, p, &cdof));
2877: PetscCall(PetscSectionGetConstraintIndices(leafSection, p, &cinds));
2878: for (i = 0, j = 0; i < dof; ++i) {
2879: PetscBool constrained = (PetscBool)(j < cdof && i == cinds[j]);
2881: if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {
2882: ilocal[m] = off++;
2883: goffs[m++] = goff++;
2884: } else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
2885: else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff;
2886: if (constrained) ++j;
2887: }
2888: }
2889: PetscCall(PetscSFCreate(comm, sectionSF));
2890: PetscCall(PetscSFSetFromOptions(*sectionSF));
2891: PetscCall(PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs));
2892: PetscCall(PetscFree(goffs));
2893: PetscFunctionReturn(PETSC_SUCCESS);
2894: }
2896: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
2897: {
2898: MPI_Comm comm;
2899: PetscMPIInt size, rank;
2900: const char *topologydm_name;
2901: const char *sectiondm_name;
2902: PetscSection sectionA, sectionB;
2903: PetscBool has;
2904: PetscInt nX, n, i;
2905: PetscSF sfAB;
2907: PetscFunctionBegin;
2908: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2909: PetscCallMPI(MPI_Comm_size(comm, &size));
2910: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2911: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2912: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
2913: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
2914: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
2915: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
2916: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
2917: /* A: on-disk points */
2918: /* X: list of global point numbers, [0, NX) */
2919: /* B: plex points */
2920: /* Load raw section (sectionA) */
2921: PetscCall(PetscSectionCreate(comm, §ionA));
2922: PetscCall(PetscViewerHDF5HasGroup(viewer, "section", &has));
2923: if (has) PetscCall(PetscSectionLoad(sectionA, viewer));
2924: else {
2925: // TODO If section is missing, create the default affine section with dim dofs on each vertex. Use PetscSplitOwnership() to split vertices.
2926: // How do I know the total number of vertices?
2927: PetscInt dim, Nf = 1, Nv, nv = PETSC_DECIDE;
2929: PetscCall(DMGetDimension(dm, &dim));
2930: PetscCall(DMPlexGetDepthStratumGlobalSize(dm, 0, &Nv));
2931: PetscCall(PetscSectionSetNumFields(sectionA, Nf));
2932: PetscCall(PetscSectionSetFieldName(sectionA, 0, "Cartesian"));
2933: PetscCall(PetscSectionSetFieldComponents(sectionA, 0, dim));
2934: for (PetscInt c = 0; c < dim; ++c) {
2935: char axis = 'X' + (char)c;
2937: PetscCall(PetscSectionSetComponentName(sectionA, 0, c, &axis));
2938: }
2939: PetscCall(PetscSplitOwnership(comm, &nv, &Nv));
2940: PetscCall(PetscSectionSetChart(sectionA, 0, nv));
2941: for (PetscInt p = 0; p < nv; ++p) {
2942: PetscCall(PetscSectionSetDof(sectionA, p, dim));
2943: PetscCall(PetscSectionSetFieldDof(sectionA, p, 0, dim));
2944: }
2945: PetscCall(PetscSectionSetUp(sectionA));
2946: }
2947: PetscCall(PetscSectionGetChart(sectionA, NULL, &n));
2948: /* Create sfAB: A -> B */
2949: #if defined(PETSC_USE_DEBUG)
2950: {
2951: PetscInt N, N1;
2953: PetscCall(PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1));
2954: PetscCallMPI(MPIU_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm));
2955: 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);
2956: }
2957: #endif
2958: {
2959: IS orderIS;
2960: const PetscInt *gpoints;
2961: PetscSF sfXA, sfAX;
2962: PetscLayout layout;
2963: PetscSFNode *owners, *buffer;
2964: PetscInt nleaves;
2965: PetscInt *ilocal;
2966: PetscSFNode *iremote;
2968: /* Create sfAX: A -> X */
2969: PetscCall(ISCreate(comm, &orderIS));
2970: PetscCall(PetscObjectSetName((PetscObject)orderIS, "order"));
2971: PetscCall(PetscLayoutSetLocalSize(orderIS->map, n));
2972: PetscCall(ISLoad(orderIS, viewer));
2973: PetscCall(PetscLayoutCreate(comm, &layout));
2974: PetscCall(PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL));
2975: PetscCall(PetscLayoutSetLocalSize(layout, nX));
2976: PetscCall(PetscLayoutSetBlockSize(layout, 1));
2977: PetscCall(PetscLayoutSetUp(layout));
2978: PetscCall(PetscSFCreate(comm, &sfXA));
2979: PetscCall(ISGetIndices(orderIS, &gpoints));
2980: PetscCall(PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints));
2981: PetscCall(ISRestoreIndices(orderIS, &gpoints));
2982: PetscCall(ISDestroy(&orderIS));
2983: PetscCall(PetscLayoutDestroy(&layout));
2984: PetscCall(PetscMalloc1(n, &owners));
2985: PetscCall(PetscMalloc1(nX, &buffer));
2986: for (i = 0; i < n; ++i) {
2987: owners[i].rank = rank;
2988: owners[i].index = i;
2989: }
2990: for (i = 0; i < nX; ++i) {
2991: buffer[i].rank = -1;
2992: buffer[i].index = -1;
2993: }
2994: PetscCall(PetscSFReduceBegin(sfXA, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC));
2995: PetscCall(PetscSFReduceEnd(sfXA, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC));
2996: PetscCall(PetscSFDestroy(&sfXA));
2997: PetscCall(PetscFree(owners));
2998: for (i = 0, nleaves = 0; i < nX; ++i)
2999: if (buffer[i].rank >= 0) nleaves++;
3000: PetscCall(PetscMalloc1(nleaves, &ilocal));
3001: PetscCall(PetscMalloc1(nleaves, &iremote));
3002: for (i = 0, nleaves = 0; i < nX; ++i) {
3003: if (buffer[i].rank >= 0) {
3004: ilocal[nleaves] = i;
3005: iremote[nleaves].rank = buffer[i].rank;
3006: iremote[nleaves].index = buffer[i].index;
3007: nleaves++;
3008: }
3009: }
3010: PetscCall(PetscFree(buffer));
3011: PetscCall(PetscSFCreate(comm, &sfAX));
3012: PetscCall(PetscSFSetFromOptions(sfAX));
3013: PetscCall(PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
3014: PetscCall(PetscSFCompose(sfAX, sfXB, &sfAB));
3015: PetscCall(PetscSFDestroy(&sfAX));
3016: }
3017: PetscCall(PetscViewerHDF5PopGroup(viewer));
3018: PetscCall(PetscViewerHDF5PopGroup(viewer));
3019: PetscCall(PetscViewerHDF5PopGroup(viewer));
3020: PetscCall(PetscViewerHDF5PopGroup(viewer));
3021: /* Create plex section (sectionB) */
3022: PetscCall(DMGetLocalSection(sectiondm, §ionB));
3023: if (lsf || gsf) {
3024: PetscLayout layout;
3025: PetscInt M, m;
3026: PetscInt *offsetsA;
3027: PetscBool includesConstraintsA;
3029: PetscCall(PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB));
3030: PetscCall(PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA));
3031: if (includesConstraintsA) PetscCall(PetscSectionGetStorageSize(sectionA, &m));
3032: else PetscCall(PetscSectionGetConstrainedStorageSize(sectionA, &m));
3033: PetscCallMPI(MPIU_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm));
3034: PetscCall(PetscLayoutCreate(comm, &layout));
3035: PetscCall(PetscLayoutSetSize(layout, M));
3036: PetscCall(PetscLayoutSetUp(layout));
3037: if (lsf) {
3038: PetscSF lsfABdata;
3040: PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata));
3041: *lsf = lsfABdata;
3042: }
3043: if (gsf) {
3044: PetscSection gsectionB, gsectionB1;
3045: PetscBool includesConstraintsB;
3046: PetscSF gsfABdata, pointsf;
3048: PetscCall(DMGetGlobalSection(sectiondm, &gsectionB1));
3049: PetscCall(PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB));
3050: PetscCall(DMGetPointSF(sectiondm, &pointsf));
3051: PetscCall(PetscSectionCreateGlobalSection(sectionB, pointsf, PETSC_TRUE, includesConstraintsB, PETSC_TRUE, &gsectionB));
3052: PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata));
3053: PetscCall(PetscSectionDestroy(&gsectionB));
3054: *gsf = gsfABdata;
3055: }
3056: PetscCall(PetscLayoutDestroy(&layout));
3057: PetscCall(PetscFree(offsetsA));
3058: } else {
3059: PetscCall(PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB));
3060: }
3061: PetscCall(PetscSFDestroy(&sfAB));
3062: PetscCall(PetscSectionDestroy(§ionA));
3063: PetscFunctionReturn(PETSC_SUCCESS);
3064: }
3066: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
3067: {
3068: MPI_Comm comm;
3069: const char *topologydm_name;
3070: const char *sectiondm_name;
3071: const char *vec_name;
3072: Vec vecA;
3073: PetscInt mA, m, bs;
3074: const PetscInt *ilocal;
3075: const PetscScalar *src;
3076: PetscScalar *dest;
3078: PetscFunctionBegin;
3079: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
3080: PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
3081: PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name));
3082: PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
3083: PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
3084: PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
3085: PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
3086: PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
3087: PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
3088: PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
3089: PetscCall(VecCreate(comm, &vecA));
3090: PetscCall(PetscObjectSetName((PetscObject)vecA, vec_name));
3091: PetscCall(PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL));
3092: /* Check consistency */
3093: {
3094: PetscSF pointsf, pointsf1;
3095: PetscInt m1, i, j;
3097: PetscCall(DMGetPointSF(dm, &pointsf));
3098: PetscCall(DMGetPointSF(sectiondm, &pointsf1));
3099: PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
3100: #if defined(PETSC_USE_DEBUG)
3101: {
3102: PetscInt MA, MA1;
3104: PetscCallMPI(MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm));
3105: PetscCall(PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1));
3106: PetscCheck(MA1 == MA, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total SF root size (%" PetscInt_FMT ") != On-disk vector data size (%" PetscInt_FMT ")", MA, MA1);
3107: }
3108: #endif
3109: PetscCall(VecGetLocalSize(vec, &m1));
3110: PetscCheck(m1 >= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Target vector size (%" PetscInt_FMT ") < SF leaf size (%" PetscInt_FMT ")", m1, m);
3111: for (i = 0; i < m; ++i) {
3112: j = ilocal ? ilocal[i] : i;
3113: 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);
3114: }
3115: }
3116: PetscCall(VecSetSizes(vecA, mA, PETSC_DECIDE));
3117: PetscCall(VecLoad(vecA, viewer));
3118: PetscCall(VecGetArrayRead(vecA, &src));
3119: PetscCall(VecGetArray(vec, &dest));
3120: PetscCall(PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE));
3121: PetscCall(PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE));
3122: PetscCall(VecRestoreArray(vec, &dest));
3123: PetscCall(VecRestoreArrayRead(vecA, &src));
3124: PetscCall(VecDestroy(&vecA));
3125: PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *)&bs));
3126: PetscCall(VecSetBlockSize(vec, bs));
3127: PetscCall(PetscViewerHDF5PopGroup(viewer));
3128: PetscCall(PetscViewerHDF5PopGroup(viewer));
3129: PetscCall(PetscViewerHDF5PopGroup(viewer));
3130: PetscCall(PetscViewerHDF5PopGroup(viewer));
3131: PetscCall(PetscViewerHDF5PopGroup(viewer));
3132: PetscCall(PetscViewerHDF5PopGroup(viewer));
3133: PetscFunctionReturn(PETSC_SUCCESS);
3134: }