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