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