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: #if defined(PETSC_HAVE_HDF5)
8: static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *);
9: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer, DMPlexStorageVersion);
10: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer, const char[], DMPlexStorageVersion);
11: static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *);
13: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
15: static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer viewer, const char str[], DMPlexStorageVersion *version)
16: {
17: PetscToken t;
18: char *ts;
19: PetscInt i;
20: PetscInt ti[3];
21: DMPlexStorageVersion v;
23: PetscTokenCreate(str, '.', &t);
24: for (i = 0; i < 3; i++) {
25: PetscTokenFind(t, &ts);
27: PetscOptionsStringToInt(ts, &ti[i]);
28: }
29: PetscTokenFind(t, &ts);
31: PetscTokenDestroy(&t);
32: PetscNew(&v);
33: v->major = ti[0];
34: v->minor = ti[1];
35: v->subminor = ti[2];
36: PetscViewerCheckVersion_Private(viewer, v);
37: *version = v;
38: return 0;
39: }
41: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion v)
42: {
43: PetscContainer cont;
45: PetscContainerCreate(PetscObjectComm((PetscObject)viewer), &cont);
46: PetscContainerSetPointer(cont, v);
47: PetscContainerSetUserDestroy(cont, PetscContainerUserDestroyDefault);
48: PetscObjectCompose((PetscObject)viewer, key, (PetscObject)cont);
49: PetscContainerDestroy(&cont);
50: return 0;
51: }
53: static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion *v)
54: {
55: PetscContainer cont;
57: PetscObjectQuery((PetscObject)viewer, key, (PetscObject *)&cont);
58: *v = NULL;
59: if (cont) PetscContainerGetPointer(cont, (void **)v);
60: return 0;
61: }
63: /*
64: Version log:
65: 1.0.0 legacy version (default if no "dmplex_storage_version" attribute found in file)
66: 2.0.0 introduce versioning and multiple topologies
67: 2.1.0 introduce distributions
68: */
69: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer viewer, DMPlexStorageVersion version)
70: {
71: PetscBool valid = PETSC_FALSE;
73: switch (version->major) {
74: case 1:
75: switch (version->minor) {
76: case 0:
77: switch (version->subminor) {
78: case 0:
79: valid = PETSC_TRUE;
80: break;
81: };
82: break;
83: };
84: break;
85: case 2:
86: switch (version->minor) {
87: case 0:
88: switch (version->subminor) {
89: case 0:
90: valid = PETSC_TRUE;
91: break;
92: };
93: break;
94: case 1:
95: switch (version->subminor) {
96: case 0:
97: valid = PETSC_TRUE;
98: break;
99: };
100: break;
101: };
102: break;
103: }
105: return 0;
106: }
108: static inline PetscBool DMPlexStorageVersionGE(DMPlexStorageVersion version, int major, int minor, int subminor)
109: {
110: return (PetscBool)((version->major == major && version->minor == minor && version->subminor >= subminor) || (version->major == major && version->minor > minor) || (version->major > major));
111: }
113: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer viewer, DMPlexStorageVersion *version)
114: {
115: const char ATTR_NAME[] = "dmplex_storage_version";
116: PetscBool fileHasVersion;
117: char optVersion[16], fileVersion[16];
121: PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, version);
122: if (*version) return 0;
124: PetscStrcpy(fileVersion, DMPLEX_STORAGE_VERSION_STABLE);
125: PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion);
126: if (fileHasVersion) {
127: char *tmp;
129: PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &tmp);
130: PetscStrcpy(fileVersion, tmp);
131: PetscFree(tmp);
132: }
133: PetscStrcpy(optVersion, fileVersion);
134: PetscObjectOptionsBegin((PetscObject)viewer);
135: PetscOptionsString("-dm_plex_view_hdf5_storage_version", "DMPlex HDF5 viewer storage version", NULL, optVersion, optVersion, sizeof(optVersion), NULL);
136: PetscOptionsEnd();
137: if (!fileHasVersion) {
138: PetscViewerHDF5WriteAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, optVersion);
139: } else {
140: PetscBool flg;
142: PetscStrcmp(fileVersion, optVersion, &flg);
144: }
145: PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT);
146: PetscViewerParseVersion_Private(viewer, optVersion, version);
147: PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, *version);
148: return 0;
149: }
151: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer viewer, DMPlexStorageVersion *version)
152: {
153: const char ATTR_NAME[] = "dmplex_storage_version";
154: char *defaultVersion;
155: char *versionString;
159: PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, version);
160: if (*version) return 0;
162: //TODO string HDF5 attribute handling is terrible and should be redesigned
163: PetscStrallocpy("1.0.0", &defaultVersion);
164: PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, &defaultVersion, &versionString);
165: PetscViewerParseVersion_Private(viewer, versionString, version);
166: PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, *version);
167: PetscFree(versionString);
168: PetscFree(defaultVersion);
169: return 0;
170: }
172: static PetscErrorCode DMPlexGetHDF5Name_Private(DM dm, const char *name[])
173: {
174: if (((PetscObject)dm)->name) {
175: PetscObjectGetName((PetscObject)dm, name);
176: } else {
177: *name = "plex";
178: }
179: return 0;
180: }
182: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
183: {
184: Vec stamp;
185: PetscMPIInt rank;
187: if (seqnum < 0) return 0;
188: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
189: VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp);
190: VecSetBlockSize(stamp, 1);
191: PetscObjectSetName((PetscObject)stamp, seqname);
192: if (rank == 0) {
193: PetscReal timeScale;
194: PetscBool istime;
196: PetscStrncmp(seqname, "time", 5, &istime);
197: if (istime) {
198: DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale);
199: value *= timeScale;
200: }
201: VecSetValue(stamp, 0, value, INSERT_VALUES);
202: }
203: VecAssemblyBegin(stamp);
204: VecAssemblyEnd(stamp);
205: PetscViewerHDF5PushGroup(viewer, "/");
206: PetscViewerHDF5PushTimestepping(viewer);
207: PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
208: VecView(stamp, viewer);
209: PetscViewerHDF5PopTimestepping(viewer);
210: PetscViewerHDF5PopGroup(viewer);
211: VecDestroy(&stamp);
212: return 0;
213: }
215: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
216: {
217: Vec stamp;
218: PetscMPIInt rank;
220: if (seqnum < 0) return 0;
221: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
222: VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp);
223: VecSetBlockSize(stamp, 1);
224: PetscObjectSetName((PetscObject)stamp, seqname);
225: PetscViewerHDF5PushGroup(viewer, "/");
226: PetscViewerHDF5PushTimestepping(viewer);
227: PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
228: VecLoad(stamp, viewer);
229: PetscViewerHDF5PopTimestepping(viewer);
230: PetscViewerHDF5PopGroup(viewer);
231: if (rank == 0) {
232: const PetscScalar *a;
233: PetscReal timeScale;
234: PetscBool istime;
236: VecGetArrayRead(stamp, &a);
237: *value = a[0];
238: VecRestoreArrayRead(stamp, &a);
239: PetscStrncmp(seqname, "time", 5, &istime);
240: if (istime) {
241: DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale);
242: *value /= timeScale;
243: }
244: }
245: VecDestroy(&stamp);
246: return 0;
247: }
249: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
250: {
251: IS cutcells = NULL;
252: const PetscInt *cutc;
253: PetscInt cellHeight, vStart, vEnd, cStart, cEnd, c;
255: if (!cutLabel) return 0;
256: DMPlexGetVTKCellHeight(dm, &cellHeight);
257: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
258: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
259: /* Label vertices that should be duplicated */
260: DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel);
261: DMLabelGetStratumIS(cutLabel, 2, &cutcells);
262: if (cutcells) {
263: PetscInt n;
265: ISGetIndices(cutcells, &cutc);
266: ISGetLocalSize(cutcells, &n);
267: for (c = 0; c < n; ++c) {
268: if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
269: PetscInt *closure = NULL;
270: PetscInt closureSize, cl, value;
272: DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
273: for (cl = 0; cl < closureSize * 2; cl += 2) {
274: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
275: DMLabelGetValue(cutLabel, closure[cl], &value);
276: if (value == 1) DMLabelSetValue(*cutVertexLabel, closure[cl], 1);
277: }
278: }
279: DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
280: }
281: }
282: ISRestoreIndices(cutcells, &cutc);
283: }
284: ISDestroy(&cutcells);
285: return 0;
286: }
288: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
289: {
290: DM dm;
291: DM dmBC;
292: PetscSection section, sectionGlobal;
293: Vec gv;
294: const char *name;
295: PetscViewerVTKFieldType ft;
296: PetscViewerFormat format;
297: PetscInt seqnum;
298: PetscReal seqval;
299: PetscBool isseq;
301: PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq);
302: VecGetDM(v, &dm);
303: DMGetLocalSection(dm, §ion);
304: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
305: DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar)seqval, viewer);
306: if (seqnum >= 0) {
307: PetscViewerHDF5PushTimestepping(viewer);
308: PetscViewerHDF5SetTimestep(viewer, seqnum);
309: }
310: PetscViewerGetFormat(viewer, &format);
311: DMGetOutputDM(dm, &dmBC);
312: DMGetGlobalSection(dmBC, §ionGlobal);
313: DMGetGlobalVector(dmBC, &gv);
314: PetscObjectGetName((PetscObject)v, &name);
315: PetscObjectSetName((PetscObject)gv, name);
316: DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
317: DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
318: PetscObjectTypeCompare((PetscObject)gv, VECSEQ, &isseq);
319: if (isseq) VecView_Seq(gv, viewer);
320: else VecView_MPI(gv, viewer);
321: if (format == PETSC_VIEWER_HDF5_VIZ) {
322: /* Output visualization representation */
323: PetscInt numFields, f;
324: DMLabel cutLabel, cutVertexLabel = NULL;
326: PetscSectionGetNumFields(section, &numFields);
327: DMGetLabel(dm, "periodic_cut", &cutLabel);
328: for (f = 0; f < numFields; ++f) {
329: Vec subv;
330: IS is;
331: const char *fname, *fgroup, *componentName;
332: char subname[PETSC_MAX_PATH_LEN];
333: PetscInt pStart, pEnd, Nc, c;
335: DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
336: fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
337: PetscSectionGetFieldName(section, f, &fname);
338: if (!fname || ft == PETSC_VTK_INVALID) continue;
339: PetscViewerHDF5PushGroup(viewer, fgroup);
340: if (cutLabel) {
341: const PetscScalar *ga;
342: PetscScalar *suba;
343: PetscInt gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;
345: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
346: PetscSectionGetFieldComponents(section, f, &Nc);
347: for (p = pStart; p < pEnd; ++p) {
348: PetscInt gdof, fdof = 0, val;
350: PetscSectionGetDof(sectionGlobal, p, &gdof);
351: if (gdof > 0) PetscSectionGetFieldDof(section, p, f, &fdof);
352: subSize += fdof;
353: DMLabelGetValue(cutVertexLabel, p, &val);
354: if (val == 1) extSize += fdof;
355: }
356: VecCreate(PetscObjectComm((PetscObject)gv), &subv);
357: VecSetSizes(subv, subSize + extSize, PETSC_DETERMINE);
358: VecSetBlockSize(subv, Nc);
359: VecSetType(subv, VECSTANDARD);
360: VecGetOwnershipRange(gv, &gstart, NULL);
361: VecGetArrayRead(gv, &ga);
362: VecGetArray(subv, &suba);
363: for (p = pStart; p < pEnd; ++p) {
364: PetscInt gdof, goff, val;
366: PetscSectionGetDof(sectionGlobal, p, &gdof);
367: if (gdof > 0) {
368: PetscInt fdof, fc, f2, poff = 0;
370: PetscSectionGetOffset(sectionGlobal, p, &goff);
371: /* Can get rid of this loop by storing field information in the global section */
372: for (f2 = 0; f2 < f; ++f2) {
373: PetscSectionGetFieldDof(section, p, f2, &fdof);
374: poff += fdof;
375: }
376: PetscSectionGetFieldDof(section, p, f, &fdof);
377: for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff + poff + fc - gstart];
378: DMLabelGetValue(cutVertexLabel, p, &val);
379: if (val == 1) {
380: for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize + newOff] = ga[goff + poff + fc - gstart];
381: }
382: }
383: }
384: VecRestoreArrayRead(gv, &ga);
385: VecRestoreArray(subv, &suba);
386: DMLabelDestroy(&cutVertexLabel);
387: } else {
388: PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
389: }
390: PetscStrncpy(subname, name, sizeof(subname));
391: PetscStrlcat(subname, "_", sizeof(subname));
392: PetscStrlcat(subname, fname, sizeof(subname));
393: PetscObjectSetName((PetscObject)subv, subname);
394: if (isseq) VecView_Seq(subv, viewer);
395: else VecView_MPI(subv, viewer);
396: if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
397: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "vector");
398: } else {
399: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "scalar");
400: }
402: /* Output the component names in the field if available */
403: PetscSectionGetFieldComponents(section, f, &Nc);
404: for (c = 0; c < Nc; ++c) {
405: char componentNameLabel[PETSC_MAX_PATH_LEN];
406: PetscSectionGetComponentName(section, f, c, &componentName);
407: PetscSNPrintf(componentNameLabel, sizeof(componentNameLabel), "componentName%" PetscInt_FMT, c);
408: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, componentNameLabel, PETSC_STRING, componentName);
409: }
411: if (cutLabel) VecDestroy(&subv);
412: else PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
413: PetscViewerHDF5PopGroup(viewer);
414: }
415: }
416: if (seqnum >= 0) PetscViewerHDF5PopTimestepping(viewer);
417: DMRestoreGlobalVector(dmBC, &gv);
418: return 0;
419: }
421: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
422: {
423: DM dm;
424: Vec locv;
425: PetscObject isZero;
426: const char *name;
427: PetscReal time;
429: VecGetDM(v, &dm);
430: DMGetLocalVector(dm, &locv);
431: PetscObjectGetName((PetscObject)v, &name);
432: PetscObjectSetName((PetscObject)locv, name);
433: PetscObjectQuery((PetscObject)v, "__Vec_bc_zero__", &isZero);
434: PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", isZero);
435: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
436: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
437: DMGetOutputSequenceNumber(dm, NULL, &time);
438: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
439: PetscViewerHDF5PushGroup(viewer, "/fields");
440: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
441: VecView_Plex_Local_HDF5_Internal(locv, viewer);
442: PetscViewerPopFormat(viewer);
443: PetscViewerHDF5PopGroup(viewer);
444: PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", NULL);
445: DMRestoreLocalVector(dm, &locv);
446: return 0;
447: }
449: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
450: {
451: PetscBool isseq;
453: PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq);
454: PetscViewerHDF5PushGroup(viewer, "/fields");
455: if (isseq) VecView_Seq(v, viewer);
456: else VecView_MPI(v, viewer);
457: PetscViewerHDF5PopGroup(viewer);
458: return 0;
459: }
461: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
462: {
463: DM dm;
464: Vec locv;
465: const char *name;
466: PetscInt seqnum;
468: VecGetDM(v, &dm);
469: DMGetLocalVector(dm, &locv);
470: PetscObjectGetName((PetscObject)v, &name);
471: PetscObjectSetName((PetscObject)locv, name);
472: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
473: PetscViewerHDF5PushGroup(viewer, "/fields");
474: if (seqnum >= 0) {
475: PetscViewerHDF5PushTimestepping(viewer);
476: PetscViewerHDF5SetTimestep(viewer, seqnum);
477: }
478: VecLoad_Plex_Local(locv, viewer);
479: if (seqnum >= 0) PetscViewerHDF5PopTimestepping(viewer);
480: PetscViewerHDF5PopGroup(viewer);
481: DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
482: DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
483: DMRestoreLocalVector(dm, &locv);
484: return 0;
485: }
487: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
488: {
489: DM dm;
490: PetscInt seqnum;
492: VecGetDM(v, &dm);
493: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
494: PetscViewerHDF5PushGroup(viewer, "/fields");
495: if (seqnum >= 0) {
496: PetscViewerHDF5PushTimestepping(viewer);
497: PetscViewerHDF5SetTimestep(viewer, seqnum);
498: }
499: VecLoad_Default(v, viewer);
500: if (seqnum >= 0) PetscViewerHDF5PopTimestepping(viewer);
501: PetscViewerHDF5PopGroup(viewer);
502: return 0;
503: }
505: static PetscErrorCode DMPlexDistributionView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
506: {
507: MPI_Comm comm;
508: PetscMPIInt size, rank;
509: const char *topologydm_name, *distribution_name;
510: const PetscInt *gpoint;
511: PetscInt pStart, pEnd, p;
512: PetscSF pointSF;
513: PetscInt nroots, nleaves;
514: const PetscInt *ilocal;
515: const PetscSFNode *iremote;
516: IS chartSizesIS, ownersIS, gpointsIS;
517: PetscInt *chartSize, *owners, *gpoints;
519: PetscObjectGetComm((PetscObject)dm, &comm);
520: MPI_Comm_size(comm, &size);
521: MPI_Comm_rank(comm, &rank);
522: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
523: DMPlexDistributionGetName(dm, &distribution_name);
524: if (!distribution_name) return 0;
525: PetscViewerHDF5WriteAttribute(viewer, NULL, "comm_size", PETSC_INT, (void *)&size);
526: ISGetIndices(globalPointNumbers, &gpoint);
527: DMPlexGetChart(dm, &pStart, &pEnd);
528: PetscMalloc1(1, &chartSize);
529: *chartSize = pEnd - pStart;
530: PetscMalloc1(*chartSize, &owners);
531: PetscMalloc1(*chartSize, &gpoints);
532: DMGetPointSF(dm, &pointSF);
533: PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote);
534: for (p = pStart; p < pEnd; ++p) {
535: PetscInt gp = gpoint[p - pStart];
537: owners[p - pStart] = rank;
538: gpoints[p - pStart] = (gp < 0 ? -(gp + 1) : gp);
539: }
540: for (p = 0; p < nleaves; ++p) {
541: PetscInt ilocalp = (ilocal ? ilocal[p] : p);
543: owners[ilocalp] = iremote[p].rank;
544: }
545: ISCreateGeneral(comm, 1, chartSize, PETSC_OWN_POINTER, &chartSizesIS);
546: ISCreateGeneral(comm, *chartSize, owners, PETSC_OWN_POINTER, &ownersIS);
547: ISCreateGeneral(comm, *chartSize, gpoints, PETSC_OWN_POINTER, &gpointsIS);
548: PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes");
549: PetscObjectSetName((PetscObject)ownersIS, "owners");
550: PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers");
551: ISView(chartSizesIS, viewer);
552: ISView(ownersIS, viewer);
553: ISView(gpointsIS, viewer);
554: ISDestroy(&chartSizesIS);
555: ISDestroy(&ownersIS);
556: ISDestroy(&gpointsIS);
557: ISRestoreIndices(globalPointNumbers, &gpoint);
558: return 0;
559: }
561: 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[])
562: {
563: IS coneSizesIS, conesIS, orientationsIS;
564: PetscInt *coneSizes, *cones, *orientations;
565: const PetscInt *gpoint;
566: PetscInt nPoints = 0, conesSize = 0;
567: PetscInt p, c, s;
568: MPI_Comm comm;
570: PetscObjectGetComm((PetscObject)dm, &comm);
571: ISGetIndices(globalPointNumbers, &gpoint);
572: for (p = pStart; p < pEnd; ++p) {
573: if (gpoint[p] >= 0) {
574: PetscInt coneSize;
576: DMPlexGetConeSize(dm, p, &coneSize);
577: nPoints += 1;
578: conesSize += coneSize;
579: }
580: }
581: PetscMalloc1(nPoints, &coneSizes);
582: PetscMalloc1(conesSize, &cones);
583: PetscMalloc1(conesSize, &orientations);
584: for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
585: if (gpoint[p] >= 0) {
586: const PetscInt *cone, *ornt;
587: PetscInt coneSize, cp;
589: DMPlexGetConeSize(dm, p, &coneSize);
590: DMPlexGetCone(dm, p, &cone);
591: DMPlexGetConeOrientation(dm, p, &ornt);
592: coneSizes[s] = coneSize;
593: for (cp = 0; cp < coneSize; ++cp, ++c) {
594: cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]] + 1) : gpoint[cone[cp]];
595: orientations[c] = ornt[cp];
596: }
597: ++s;
598: }
599: }
602: ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS);
603: ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS);
604: ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS);
605: PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName);
606: PetscObjectSetName((PetscObject)conesIS, conesName);
607: PetscObjectSetName((PetscObject)orientationsIS, orientationsName);
608: ISView(coneSizesIS, viewer);
609: ISView(conesIS, viewer);
610: ISView(orientationsIS, viewer);
611: ISDestroy(&coneSizesIS);
612: ISDestroy(&conesIS);
613: ISDestroy(&orientationsIS);
614: if (pointsName) {
615: IS pointsIS;
616: PetscInt *points;
618: PetscMalloc1(nPoints, &points);
619: for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
620: if (gpoint[p] >= 0) {
621: points[s] = gpoint[p];
622: ++s;
623: }
624: }
625: ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS);
626: PetscObjectSetName((PetscObject)pointsIS, pointsName);
627: ISView(pointsIS, viewer);
628: ISDestroy(&pointsIS);
629: }
630: ISRestoreIndices(globalPointNumbers, &gpoint);
631: return 0;
632: }
634: static PetscErrorCode DMPlexTopologyView_HDF5_Legacy_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
635: {
636: const char *pointsName, *coneSizesName, *conesName, *orientationsName;
637: PetscInt pStart, pEnd, dim;
639: pointsName = "order";
640: coneSizesName = "cones";
641: conesName = "cells";
642: orientationsName = "orientation";
643: DMPlexGetChart(dm, &pStart, &pEnd);
644: DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers, viewer, pStart, pEnd, pointsName, coneSizesName, conesName, orientationsName);
645: DMGetDimension(dm, &dim);
646: PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *)&dim);
647: return 0;
648: }
650: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
651: {
652: DMPlexStorageVersion version;
653: const char *topologydm_name;
654: char group[PETSC_MAX_PATH_LEN];
656: PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version);
657: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
658: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
659: PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name);
660: } else {
661: PetscStrcpy(group, "/");
662: }
663: PetscViewerHDF5PushGroup(viewer, group);
665: PetscViewerHDF5PushGroup(viewer, "topology");
666: DMPlexTopologyView_HDF5_Legacy_Private(dm, globalPointNumbers, viewer);
667: PetscViewerHDF5PopGroup(viewer); /* "topology" */
669: if (DMPlexStorageVersionGE(version, 2, 1, 0)) {
670: const char *distribution_name;
672: DMPlexDistributionGetName(dm, &distribution_name);
673: PetscViewerHDF5PushGroup(viewer, "distributions");
674: PetscViewerHDF5WriteGroup(viewer, NULL);
675: PetscViewerHDF5PushGroup(viewer, distribution_name);
676: DMPlexDistributionView_HDF5_Private(dm, globalPointNumbers, viewer);
677: PetscViewerHDF5PopGroup(viewer); /* distribution_name */
678: PetscViewerHDF5PopGroup(viewer); /* "distributions" */
679: }
681: PetscViewerHDF5PopGroup(viewer);
682: return 0;
683: }
685: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
686: {
687: PetscSF sfPoint;
688: DMLabel cutLabel, cutVertexLabel = NULL;
689: IS globalVertexNumbers, cutvertices = NULL;
690: const PetscInt *gcell, *gvertex, *cutverts = NULL;
691: PetscInt *vertices;
692: PetscInt conesSize = 0;
693: PetscInt dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
695: *numCorners = 0;
696: DMGetDimension(dm, &dim);
697: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
698: ISGetIndices(globalCellNumbers, &gcell);
700: for (cell = cStart; cell < cEnd; ++cell) {
701: PetscInt *closure = NULL;
702: PetscInt closureSize, v, Nc = 0;
704: if (gcell[cell] < 0) continue;
705: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
706: for (v = 0; v < closureSize * 2; v += 2) {
707: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
708: }
709: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
710: conesSize += Nc;
711: if (!numCornersLocal) numCornersLocal = Nc;
712: else if (numCornersLocal != Nc) numCornersLocal = 1;
713: }
714: MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
716: /* Handle periodic cuts by identifying vertices which should be duplicated */
717: DMGetLabel(dm, "periodic_cut", &cutLabel);
718: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
719: if (cutVertexLabel) DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);
720: if (cutvertices) {
721: ISGetIndices(cutvertices, &cutverts);
722: ISGetLocalSize(cutvertices, &vExtra);
723: }
724: DMGetPointSF(dm, &sfPoint);
725: if (cutLabel) {
726: const PetscInt *ilocal;
727: const PetscSFNode *iremote;
728: PetscInt nroots, nleaves;
730: PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
731: if (nleaves < 0) {
732: PetscObjectReference((PetscObject)sfPoint);
733: } else {
734: PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), &sfPoint);
735: PetscSFSetGraph(sfPoint, nroots + vExtra, nleaves, (PetscInt *)ilocal, PETSC_COPY_VALUES, (PetscSFNode *)iremote, PETSC_COPY_VALUES);
736: }
737: } else {
738: PetscObjectReference((PetscObject)sfPoint);
739: }
740: /* Number all vertices */
741: DMPlexCreateNumbering_Plex(dm, vStart, vEnd + vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
742: PetscSFDestroy(&sfPoint);
743: /* Create cones */
744: ISGetIndices(globalVertexNumbers, &gvertex);
745: PetscMalloc1(conesSize, &vertices);
746: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
747: PetscInt *closure = NULL;
748: PetscInt closureSize, Nc = 0, p, value = -1;
749: PetscBool replace;
751: if (gcell[cell] < 0) continue;
752: if (cutLabel) DMLabelGetValue(cutLabel, cell, &value);
753: replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
754: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
755: for (p = 0; p < closureSize * 2; p += 2) {
756: if ((closure[p] >= vStart) && (closure[p] < vEnd)) closure[Nc++] = closure[p];
757: }
758: DMPlexReorderCell(dm, cell, closure);
759: for (p = 0; p < Nc; ++p) {
760: PetscInt nv, gv = gvertex[closure[p] - vStart];
762: if (replace) {
763: PetscFindInt(closure[p], vExtra, cutverts, &nv);
764: if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
765: }
766: vertices[v++] = gv < 0 ? -(gv + 1) : gv;
767: }
768: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
769: }
770: ISRestoreIndices(globalVertexNumbers, &gvertex);
771: ISDestroy(&globalVertexNumbers);
772: ISRestoreIndices(globalCellNumbers, &gcell);
773: if (cutvertices) ISRestoreIndices(cutvertices, &cutverts);
774: ISDestroy(&cutvertices);
775: DMLabelDestroy(&cutVertexLabel);
777: ISCreateGeneral(PetscObjectComm((PetscObject)dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS);
778: PetscLayoutSetBlockSize((*cellIS)->map, *numCorners);
779: PetscObjectSetName((PetscObject)*cellIS, "cells");
780: return 0;
781: }
783: static PetscErrorCode DMPlexTopologyView_HDF5_XDMF_Private(DM dm, IS globalCellNumbers, PetscViewer viewer)
784: {
785: DM cdm;
786: DMLabel depthLabel, ctLabel;
787: IS cellIS;
788: PetscInt dim, depth, cellHeight, c, n = 0;
790: PetscViewerHDF5PushGroup(viewer, "/viz");
791: PetscViewerHDF5WriteGroup(viewer, NULL);
793: PetscViewerHDF5PopGroup(viewer);
794: DMGetDimension(dm, &dim);
795: DMPlexGetDepth(dm, &depth);
796: DMGetCoordinateDM(dm, &cdm);
797: DMPlexGetVTKCellHeight(dm, &cellHeight);
798: DMPlexGetDepthLabel(dm, &depthLabel);
799: DMPlexGetCellTypeLabel(dm, &ctLabel);
800: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
801: const DMPolytopeType ict = (DMPolytopeType)c;
802: PetscInt pStart, pEnd, dep, numCorners;
803: PetscBool output = PETSC_FALSE, doOutput;
805: if (ict == DM_POLYTOPE_FV_GHOST) continue;
806: DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
807: if (pStart >= 0) {
808: DMLabelGetValue(depthLabel, pStart, &dep);
809: if (dep == depth - cellHeight) output = PETSC_TRUE;
810: }
811: MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm));
812: if (!doOutput) continue;
813: CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS);
814: if (!n) {
815: PetscViewerHDF5PushGroup(viewer, "/viz/topology");
816: } else {
817: char group[PETSC_MAX_PATH_LEN];
819: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%" PetscInt_FMT, n);
820: PetscViewerHDF5PushGroup(viewer, group);
821: }
822: ISView(cellIS, viewer);
823: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_corners", PETSC_INT, (void *)&numCorners);
824: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_dim", PETSC_INT, (void *)&dim);
825: ISDestroy(&cellIS);
826: PetscViewerHDF5PopGroup(viewer);
827: ++n;
828: }
829: return 0;
830: }
832: static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
833: {
834: DM cdm;
835: Vec coordinates, newcoords;
836: PetscReal lengthScale;
837: PetscInt m, M, bs;
839: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
840: DMGetCoordinateDM(dm, &cdm);
841: DMGetCoordinates(dm, &coordinates);
842: VecCreate(PetscObjectComm((PetscObject)coordinates), &newcoords);
843: PetscObjectSetName((PetscObject)newcoords, "vertices");
844: VecGetSize(coordinates, &M);
845: VecGetLocalSize(coordinates, &m);
846: VecSetSizes(newcoords, m, M);
847: VecGetBlockSize(coordinates, &bs);
848: VecSetBlockSize(newcoords, bs);
849: VecSetType(newcoords, VECSTANDARD);
850: VecCopy(coordinates, newcoords);
851: VecScale(newcoords, lengthScale);
852: /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
853: PetscViewerHDF5PushGroup(viewer, "/geometry");
854: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
855: VecView(newcoords, viewer);
856: PetscViewerPopFormat(viewer);
857: PetscViewerHDF5PopGroup(viewer);
858: VecDestroy(&newcoords);
859: return 0;
860: }
862: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
863: {
864: DM cdm;
865: Vec coords, newcoords;
866: PetscInt m, M, bs;
867: PetscReal lengthScale;
868: const char *topologydm_name, *coordinatedm_name, *coordinates_name;
870: {
871: PetscViewerFormat format;
872: DMPlexStorageVersion version;
874: PetscViewerGetFormat(viewer, &format);
875: PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version);
876: if (!DMPlexStorageVersionGE(version, 2, 0, 0) || format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ) {
877: DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer);
878: return 0;
879: }
880: }
881: /* since 2.0.0 */
882: DMGetCoordinateDM(dm, &cdm);
883: DMGetCoordinates(dm, &coords);
884: PetscObjectGetName((PetscObject)cdm, &coordinatedm_name);
885: PetscObjectGetName((PetscObject)coords, &coordinates_name);
886: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
887: PetscViewerHDF5PushGroup(viewer, "topologies");
888: PetscViewerHDF5PushGroup(viewer, topologydm_name);
889: PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name);
890: PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name);
891: PetscViewerHDF5PopGroup(viewer);
892: PetscViewerHDF5PopGroup(viewer);
893: DMPlexSectionView(dm, viewer, cdm);
894: VecCreate(PetscObjectComm((PetscObject)coords), &newcoords);
895: PetscObjectSetName((PetscObject)newcoords, coordinates_name);
896: VecGetSize(coords, &M);
897: VecGetLocalSize(coords, &m);
898: VecSetSizes(newcoords, m, M);
899: VecGetBlockSize(coords, &bs);
900: VecSetBlockSize(newcoords, bs);
901: VecSetType(newcoords, VECSTANDARD);
902: VecCopy(coords, newcoords);
903: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
904: VecScale(newcoords, lengthScale);
905: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
906: DMPlexGlobalVectorView(dm, viewer, cdm, newcoords);
907: PetscViewerPopFormat(viewer);
908: VecDestroy(&newcoords);
909: return 0;
910: }
912: static PetscErrorCode DMPlexCoordinatesView_HDF5_XDMF_Private(DM dm, PetscViewer viewer)
913: {
914: DM cdm;
915: Vec coordinatesLocal, newcoords;
916: PetscSection cSection, cGlobalSection;
917: PetscScalar *coords, *ncoords;
918: DMLabel cutLabel, cutVertexLabel = NULL;
919: const PetscReal *L;
920: PetscReal lengthScale;
921: PetscInt vStart, vEnd, v, bs, N, coordSize, dof, off, d;
922: PetscBool localized, embedded;
924: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
925: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
926: DMGetCoordinatesLocal(dm, &coordinatesLocal);
927: VecGetBlockSize(coordinatesLocal, &bs);
928: DMGetCoordinatesLocalized(dm, &localized);
929: if (localized == PETSC_FALSE) return 0;
930: DMGetPeriodicity(dm, NULL, NULL, &L);
931: DMGetCoordinateDM(dm, &cdm);
932: DMGetLocalSection(cdm, &cSection);
933: DMGetGlobalSection(cdm, &cGlobalSection);
934: DMGetLabel(dm, "periodic_cut", &cutLabel);
935: N = 0;
937: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
938: VecCreate(PetscObjectComm((PetscObject)dm), &newcoords);
939: PetscSectionGetDof(cSection, vStart, &dof);
940: PetscPrintf(PETSC_COMM_SELF, "DOF: %" PetscInt_FMT "\n", dof);
941: embedded = (PetscBool)(L && dof == 2 && !cutLabel);
942: if (cutVertexLabel) {
943: DMLabelGetStratumSize(cutVertexLabel, 1, &v);
944: N += dof * v;
945: }
946: for (v = vStart; v < vEnd; ++v) {
947: PetscSectionGetDof(cGlobalSection, v, &dof);
948: if (dof < 0) continue;
949: if (embedded) N += dof + 1;
950: else N += dof;
951: }
952: if (embedded) VecSetBlockSize(newcoords, bs + 1);
953: else VecSetBlockSize(newcoords, bs);
954: VecSetSizes(newcoords, N, PETSC_DETERMINE);
955: VecSetType(newcoords, VECSTANDARD);
956: VecGetArray(coordinatesLocal, &coords);
957: VecGetArray(newcoords, &ncoords);
958: coordSize = 0;
959: for (v = vStart; v < vEnd; ++v) {
960: PetscSectionGetDof(cGlobalSection, v, &dof);
961: PetscSectionGetOffset(cSection, v, &off);
962: if (dof < 0) continue;
963: if (embedded) {
964: if (L && (L[0] > 0.0) && (L[1] > 0.0)) {
965: PetscReal theta, phi, r, R;
966: /* XY-periodic */
967: /* Suppose its an y-z circle, then
968: \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
969: and the circle in that plane is
970: \hat r cos(phi) + \hat x sin(phi) */
971: theta = 2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1];
972: phi = 2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0];
973: r = L[0] / (2.0 * PETSC_PI * 2.0 * L[1]);
974: R = L[1] / (2.0 * PETSC_PI);
975: ncoords[coordSize++] = PetscSinReal(phi) * r;
976: ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
977: ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
978: } else if (L && (L[0] > 0.0)) {
979: /* X-periodic */
980: ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
981: ncoords[coordSize++] = coords[off + 1];
982: ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
983: } else if (L && (L[1] > 0.0)) {
984: /* Y-periodic */
985: ncoords[coordSize++] = coords[off + 0];
986: ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
987: ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
988: #if 0
989: } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
990: PetscReal phi, r, R;
991: /* Mobius strip */
992: /* Suppose its an x-z circle, then
993: \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
994: and in that plane we rotate by pi as we go around the circle
995: \hat r cos(phi/2) + \hat y sin(phi/2) */
996: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
997: R = L[0];
998: r = PetscRealPart(coords[off+1]) - L[1]/2.0;
999: ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
1000: ncoords[coordSize++] = PetscSinReal(phi/2.0) * r;
1001: ncoords[coordSize++] = PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
1002: #endif
1003: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
1004: } else {
1005: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off + d];
1006: }
1007: }
1008: if (cutVertexLabel) {
1009: IS vertices;
1010: const PetscInt *verts;
1011: PetscInt n;
1013: DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
1014: if (vertices) {
1015: ISGetIndices(vertices, &verts);
1016: ISGetLocalSize(vertices, &n);
1017: for (v = 0; v < n; ++v) {
1018: PetscSectionGetDof(cSection, verts[v], &dof);
1019: PetscSectionGetOffset(cSection, verts[v], &off);
1020: for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off + d] + ((L[d] > 0.) ? L[d] : 0.0);
1021: }
1022: ISRestoreIndices(vertices, &verts);
1023: ISDestroy(&vertices);
1024: }
1025: }
1027: DMLabelDestroy(&cutVertexLabel);
1028: VecRestoreArray(coordinatesLocal, &coords);
1029: VecRestoreArray(newcoords, &ncoords);
1030: PetscObjectSetName((PetscObject)newcoords, "vertices");
1031: VecScale(newcoords, lengthScale);
1032: PetscViewerHDF5PushGroup(viewer, "/viz");
1033: PetscViewerHDF5WriteGroup(viewer, NULL);
1034: PetscViewerHDF5PopGroup(viewer);
1035: PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
1036: VecView(newcoords, viewer);
1037: PetscViewerHDF5PopGroup(viewer);
1038: VecDestroy(&newcoords);
1039: return 0;
1040: }
1042: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
1043: {
1044: const char *topologydm_name;
1045: const PetscInt *gpoint;
1046: PetscInt numLabels, l;
1047: DMPlexStorageVersion version;
1048: char group[PETSC_MAX_PATH_LEN];
1050: PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version);
1051: ISGetIndices(globalPointNumbers, &gpoint);
1052: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1053: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1054: PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name);
1055: } else {
1056: PetscStrcpy(group, "/labels");
1057: }
1058: PetscViewerHDF5PushGroup(viewer, group);
1059: DMGetNumLabels(dm, &numLabels);
1060: for (l = 0; l < numLabels; ++l) {
1061: DMLabel label;
1062: const char *name;
1063: IS valueIS, pvalueIS, globalValueIS;
1064: const PetscInt *values;
1065: PetscInt numValues, v;
1066: PetscBool isDepth, output;
1068: DMGetLabelByNum(dm, l, &label);
1069: PetscObjectGetName((PetscObject)label, &name);
1070: DMGetLabelOutput(dm, name, &output);
1071: PetscStrncmp(name, "depth", 10, &isDepth);
1072: if (isDepth || !output) continue;
1073: PetscViewerHDF5PushGroup(viewer, name);
1074: DMLabelGetValueIS(label, &valueIS);
1075: /* Must copy to a new IS on the global comm */
1076: ISGetLocalSize(valueIS, &numValues);
1077: ISGetIndices(valueIS, &values);
1078: ISCreateGeneral(PetscObjectComm((PetscObject)dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
1079: ISRestoreIndices(valueIS, &values);
1080: ISAllGather(pvalueIS, &globalValueIS);
1081: ISDestroy(&pvalueIS);
1082: ISSortRemoveDups(globalValueIS);
1083: ISGetLocalSize(globalValueIS, &numValues);
1084: ISGetIndices(globalValueIS, &values);
1085: for (v = 0; v < numValues; ++v) {
1086: IS stratumIS, globalStratumIS;
1087: const PetscInt *spoints = NULL;
1088: PetscInt *gspoints, n = 0, gn, p;
1089: const char *iname = "indices";
1090: char group[PETSC_MAX_PATH_LEN];
1092: PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, values[v]);
1093: PetscViewerHDF5PushGroup(viewer, group);
1094: DMLabelGetStratumIS(label, values[v], &stratumIS);
1096: if (stratumIS) ISGetLocalSize(stratumIS, &n);
1097: if (stratumIS) ISGetIndices(stratumIS, &spoints);
1098: for (gn = 0, p = 0; p < n; ++p)
1099: if (gpoint[spoints[p]] >= 0) ++gn;
1100: PetscMalloc1(gn, &gspoints);
1101: for (gn = 0, p = 0; p < n; ++p)
1102: if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
1103: if (stratumIS) ISRestoreIndices(stratumIS, &spoints);
1104: ISCreateGeneral(PetscObjectComm((PetscObject)dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
1105: PetscObjectSetName((PetscObject)globalStratumIS, iname);
1107: ISView(globalStratumIS, viewer);
1108: ISDestroy(&globalStratumIS);
1109: ISDestroy(&stratumIS);
1110: PetscViewerHDF5PopGroup(viewer);
1111: }
1112: ISRestoreIndices(globalValueIS, &values);
1113: ISDestroy(&globalValueIS);
1114: ISDestroy(&valueIS);
1115: PetscViewerHDF5PopGroup(viewer);
1116: }
1117: ISRestoreIndices(globalPointNumbers, &gpoint);
1118: PetscViewerHDF5PopGroup(viewer);
1119: return 0;
1120: }
1122: /* We only write cells and vertices. Does this screw up parallel reading? */
1123: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
1124: {
1125: IS globalPointNumbers;
1126: PetscViewerFormat format;
1127: PetscBool viz_geom = PETSC_FALSE, xdmf_topo = PETSC_FALSE, petsc_topo = PETSC_FALSE;
1129: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
1130: DMPlexCoordinatesView_HDF5_Internal(dm, viewer);
1132: PetscViewerGetFormat(viewer, &format);
1133: switch (format) {
1134: case PETSC_VIEWER_HDF5_VIZ:
1135: viz_geom = PETSC_TRUE;
1136: xdmf_topo = PETSC_TRUE;
1137: break;
1138: case PETSC_VIEWER_HDF5_XDMF:
1139: xdmf_topo = PETSC_TRUE;
1140: break;
1141: case PETSC_VIEWER_HDF5_PETSC:
1142: petsc_topo = PETSC_TRUE;
1143: break;
1144: case PETSC_VIEWER_DEFAULT:
1145: case PETSC_VIEWER_NATIVE:
1146: viz_geom = PETSC_TRUE;
1147: xdmf_topo = PETSC_TRUE;
1148: petsc_topo = PETSC_TRUE;
1149: break;
1150: default:
1151: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
1152: }
1154: if (viz_geom) DMPlexCoordinatesView_HDF5_XDMF_Private(dm, viewer);
1155: if (xdmf_topo) DMPlexTopologyView_HDF5_XDMF_Private(dm, globalPointNumbers, viewer);
1156: if (petsc_topo) {
1157: DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer);
1158: DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer);
1159: }
1161: ISDestroy(&globalPointNumbers);
1162: return 0;
1163: }
1165: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
1166: {
1167: MPI_Comm comm;
1168: const char *topologydm_name;
1169: const char *sectiondm_name;
1170: PetscSection gsection;
1172: PetscObjectGetComm((PetscObject)sectiondm, &comm);
1173: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1174: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1175: PetscViewerHDF5PushGroup(viewer, "topologies");
1176: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1177: PetscViewerHDF5PushGroup(viewer, "dms");
1178: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1179: DMGetGlobalSection(sectiondm, &gsection);
1180: /* Save raw section */
1181: PetscSectionView(gsection, viewer);
1182: /* Save plex wrapper */
1183: {
1184: PetscInt pStart, pEnd, p, n;
1185: IS globalPointNumbers;
1186: const PetscInt *gpoints;
1187: IS orderIS;
1188: PetscInt *order;
1190: PetscSectionGetChart(gsection, &pStart, &pEnd);
1191: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
1192: ISGetIndices(globalPointNumbers, &gpoints);
1193: for (p = pStart, n = 0; p < pEnd; ++p)
1194: if (gpoints[p] >= 0) n++;
1195: /* "order" is an array of global point numbers.
1196: When loading, it is used with topology/order array
1197: to match section points with plex topology points. */
1198: PetscMalloc1(n, &order);
1199: for (p = pStart, n = 0; p < pEnd; ++p)
1200: if (gpoints[p] >= 0) order[n++] = gpoints[p];
1201: ISRestoreIndices(globalPointNumbers, &gpoints);
1202: ISDestroy(&globalPointNumbers);
1203: ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS);
1204: PetscObjectSetName((PetscObject)orderIS, "order");
1205: ISView(orderIS, viewer);
1206: ISDestroy(&orderIS);
1207: }
1208: PetscViewerHDF5PopGroup(viewer);
1209: PetscViewerHDF5PopGroup(viewer);
1210: PetscViewerHDF5PopGroup(viewer);
1211: PetscViewerHDF5PopGroup(viewer);
1212: return 0;
1213: }
1215: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1216: {
1217: const char *topologydm_name;
1218: const char *sectiondm_name;
1219: const char *vec_name;
1220: PetscInt bs;
1222: /* Check consistency */
1223: {
1224: PetscSF pointsf, pointsf1;
1226: DMGetPointSF(dm, &pointsf);
1227: DMGetPointSF(sectiondm, &pointsf1);
1229: }
1230: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1231: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1232: PetscObjectGetName((PetscObject)vec, &vec_name);
1233: PetscViewerHDF5PushGroup(viewer, "topologies");
1234: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1235: PetscViewerHDF5PushGroup(viewer, "dms");
1236: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1237: PetscViewerHDF5PushGroup(viewer, "vecs");
1238: PetscViewerHDF5PushGroup(viewer, vec_name);
1239: VecGetBlockSize(vec, &bs);
1240: PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs);
1241: VecSetBlockSize(vec, 1);
1242: /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but, */
1243: /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view */
1244: /* is set to VecView_Plex, which would save vec in a predefined location. */
1245: /* To save vec in where we want, we create a new Vec (temp) with */
1246: /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
1247: {
1248: Vec temp;
1249: const PetscScalar *array;
1250: PetscLayout map;
1252: VecCreate(PetscObjectComm((PetscObject)vec), &temp);
1253: PetscObjectSetName((PetscObject)temp, vec_name);
1254: VecGetLayout(vec, &map);
1255: VecSetLayout(temp, map);
1256: VecSetUp(temp);
1257: VecGetArrayRead(vec, &array);
1258: VecPlaceArray(temp, array);
1259: VecView(temp, viewer);
1260: VecResetArray(temp);
1261: VecRestoreArrayRead(vec, &array);
1262: VecDestroy(&temp);
1263: }
1264: VecSetBlockSize(vec, bs);
1265: PetscViewerHDF5PopGroup(viewer);
1266: PetscViewerHDF5PopGroup(viewer);
1267: PetscViewerHDF5PopGroup(viewer);
1268: PetscViewerHDF5PopGroup(viewer);
1269: PetscViewerHDF5PopGroup(viewer);
1270: PetscViewerHDF5PopGroup(viewer);
1271: return 0;
1272: }
1274: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1275: {
1276: MPI_Comm comm;
1277: const char *topologydm_name;
1278: const char *sectiondm_name;
1279: const char *vec_name;
1280: PetscSection section;
1281: PetscBool includesConstraints;
1282: Vec gvec;
1283: PetscInt m, bs;
1285: PetscObjectGetComm((PetscObject)dm, &comm);
1286: /* Check consistency */
1287: {
1288: PetscSF pointsf, pointsf1;
1290: DMGetPointSF(dm, &pointsf);
1291: DMGetPointSF(sectiondm, &pointsf1);
1293: }
1294: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1295: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1296: PetscObjectGetName((PetscObject)vec, &vec_name);
1297: PetscViewerHDF5PushGroup(viewer, "topologies");
1298: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1299: PetscViewerHDF5PushGroup(viewer, "dms");
1300: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1301: PetscViewerHDF5PushGroup(viewer, "vecs");
1302: PetscViewerHDF5PushGroup(viewer, vec_name);
1303: VecGetBlockSize(vec, &bs);
1304: PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs);
1305: VecCreate(comm, &gvec);
1306: PetscObjectSetName((PetscObject)gvec, vec_name);
1307: DMGetGlobalSection(sectiondm, §ion);
1308: PetscSectionGetIncludesConstraints(section, &includesConstraints);
1309: if (includesConstraints) PetscSectionGetStorageSize(section, &m);
1310: else PetscSectionGetConstrainedStorageSize(section, &m);
1311: VecSetSizes(gvec, m, PETSC_DECIDE);
1312: VecSetUp(gvec);
1313: DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec);
1314: DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec);
1315: VecView(gvec, viewer);
1316: VecDestroy(&gvec);
1317: PetscViewerHDF5PopGroup(viewer);
1318: PetscViewerHDF5PopGroup(viewer);
1319: PetscViewerHDF5PopGroup(viewer);
1320: PetscViewerHDF5PopGroup(viewer);
1321: PetscViewerHDF5PopGroup(viewer);
1322: PetscViewerHDF5PopGroup(viewer);
1323: return 0;
1324: }
1326: struct _n_LoadLabelsCtx {
1327: MPI_Comm comm;
1328: PetscMPIInt rank;
1329: DM dm;
1330: PetscViewer viewer;
1331: DMLabel label;
1332: PetscSF sfXC;
1333: PetscLayout layoutX;
1334: };
1335: typedef struct _n_LoadLabelsCtx *LoadLabelsCtx;
1337: static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx)
1338: {
1339: PetscNew(ctx);
1340: PetscObjectReference((PetscObject)((*ctx)->dm = dm));
1341: PetscObjectReference((PetscObject)((*ctx)->viewer = viewer));
1342: PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm);
1343: MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank);
1344: (*ctx)->sfXC = sfXC;
1345: if (sfXC) {
1346: PetscInt nX;
1348: PetscObjectReference((PetscObject)sfXC);
1349: PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL);
1350: PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX);
1351: }
1352: return 0;
1353: }
1355: static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx)
1356: {
1357: if (!*ctx) return 0;
1358: DMDestroy(&(*ctx)->dm);
1359: PetscViewerDestroy(&(*ctx)->viewer);
1360: PetscSFDestroy(&(*ctx)->sfXC);
1361: PetscLayoutDestroy(&(*ctx)->layoutX);
1362: PetscFree(*ctx);
1363: return 0;
1364: }
1366: /*
1367: A: on-disk points
1368: X: global points [0, NX)
1369: C: distributed plex points
1370: */
1371: static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS)
1372: {
1373: MPI_Comm comm = ctx->comm;
1374: PetscSF sfXC = ctx->sfXC;
1375: PetscLayout layoutX = ctx->layoutX;
1376: PetscSF sfXA;
1377: const PetscInt *A_points;
1378: PetscInt nX, nC;
1379: PetscInt n;
1381: PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL);
1382: ISGetLocalSize(stratumIS, &n);
1383: ISGetIndices(stratumIS, &A_points);
1384: PetscSFCreate(comm, &sfXA);
1385: PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points);
1386: ISCreate(comm, newStratumIS);
1387: ISSetType(*newStratumIS, ISGENERAL);
1388: {
1389: PetscInt i;
1390: PetscBool *A_mask, *X_mask, *C_mask;
1392: PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask);
1393: for (i = 0; i < n; i++) A_mask[i] = PETSC_TRUE;
1394: PetscSFReduceBegin(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE);
1395: PetscSFReduceEnd(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE);
1396: PetscSFBcastBegin(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR);
1397: PetscSFBcastEnd(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR);
1398: ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask);
1399: PetscFree3(A_mask, X_mask, C_mask);
1400: }
1401: PetscSFDestroy(&sfXA);
1402: ISRestoreIndices(stratumIS, &A_points);
1403: return 0;
1404: }
1406: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data)
1407: {
1408: LoadLabelsCtx ctx = (LoadLabelsCtx)op_data;
1409: PetscViewer viewer = ctx->viewer;
1410: DMLabel label = ctx->label;
1411: MPI_Comm comm = ctx->comm;
1412: IS stratumIS;
1413: const PetscInt *ind;
1414: PetscInt value, N, i;
1416: PetscOptionsStringToInt(vname, &value);
1417: ISCreate(comm, &stratumIS);
1418: PetscObjectSetName((PetscObject)stratumIS, "indices");
1419: PetscViewerHDF5PushGroup(viewer, vname); /* labels/<lname>/<vname> */
1421: if (!ctx->sfXC) {
1422: /* Force serial load */
1423: PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
1424: PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0);
1425: PetscLayoutSetSize(stratumIS->map, N);
1426: }
1427: ISLoad(stratumIS, viewer);
1429: if (ctx->sfXC) {
1430: IS newStratumIS;
1432: ReadLabelStratumHDF5_Distribute_Private(stratumIS, ctx, &newStratumIS);
1433: ISDestroy(&stratumIS);
1434: stratumIS = newStratumIS;
1435: }
1437: PetscViewerHDF5PopGroup(viewer);
1438: ISGetLocalSize(stratumIS, &N);
1439: ISGetIndices(stratumIS, &ind);
1440: for (i = 0; i < N; ++i) DMLabelSetValue(label, ind[i], value);
1441: ISRestoreIndices(stratumIS, &ind);
1442: ISDestroy(&stratumIS);
1443: return 0;
1444: }
1446: /* TODO: Fix this code, it is returning PETSc error codes when it should be translating them to herr_t codes */
1447: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data)
1448: {
1449: LoadLabelsCtx ctx = (LoadLabelsCtx)op_data;
1450: DM dm = ctx->dm;
1451: hsize_t idx = 0;
1453: PetscBool flg;
1454: herr_t err;
1456: DMHasLabel(dm, lname, &flg);
1457: if (flg) DMRemoveLabel(dm, lname, NULL);
1458: DMCreateLabel(dm, lname);
1459: if (ierr) return (herr_t)ierr;
1460: DMGetLabel(dm, lname, &ctx->label);
1461: if (ierr) return (herr_t)ierr;
1462: PetscViewerHDF5PushGroup(ctx->viewer, lname);
1463: if (ierr) return (herr_t)ierr;
1464: /* Iterate over the label's strata */
1465: PetscCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1466: PetscViewerHDF5PopGroup(ctx->viewer);
1467: if (ierr) return (herr_t)ierr;
1468: return err;
1469: }
1471: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1472: {
1473: const char *topologydm_name;
1474: LoadLabelsCtx ctx;
1475: hsize_t idx = 0;
1476: char group[PETSC_MAX_PATH_LEN];
1477: DMPlexStorageVersion version;
1478: PetscBool distributed, hasGroup;
1480: DMPlexIsDistributed(dm, &distributed);
1482: LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx);
1483: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1484: PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version);
1485: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1486: PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name);
1487: } else {
1488: PetscStrcpy(group, "labels");
1489: }
1490: PetscViewerHDF5PushGroup(viewer, group);
1491: PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup);
1492: if (hasGroup) {
1493: hid_t fileId, groupId;
1495: PetscViewerHDF5OpenGroup(viewer, NULL, &fileId, &groupId);
1496: /* Iterate over labels */
1497: PetscCallHDF5(H5Literate, (groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx));
1498: PetscCallHDF5(H5Gclose, (groupId));
1499: }
1500: PetscViewerHDF5PopGroup(viewer);
1501: LoadLabelsCtxDestroy(&ctx);
1502: return 0;
1503: }
1505: static PetscErrorCode DMPlexDistributionLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF sf, PetscSF *distsf, DM *distdm)
1506: {
1507: MPI_Comm comm;
1508: PetscMPIInt size, rank, dist_size;
1509: const char *distribution_name;
1510: PetscInt p, lsize;
1511: IS chartSizesIS, ownersIS, gpointsIS;
1512: const PetscInt *chartSize, *owners, *gpoints;
1513: PetscLayout layout;
1514: PetscBool has;
1516: *distsf = NULL;
1517: *distdm = NULL;
1518: PetscObjectGetComm((PetscObject)dm, &comm);
1519: MPI_Comm_size(comm, &size);
1520: MPI_Comm_rank(comm, &rank);
1521: DMPlexDistributionGetName(dm, &distribution_name);
1522: if (!distribution_name) return 0;
1523: PetscViewerHDF5HasGroup(viewer, NULL, &has);
1524: if (!has) {
1525: char *full_group;
1527: PetscViewerHDF5GetGroup(viewer, NULL, &full_group);
1529: }
1530: PetscViewerHDF5ReadAttribute(viewer, NULL, "comm_size", PETSC_INT, NULL, (void *)&dist_size);
1532: ISCreate(comm, &chartSizesIS);
1533: PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes");
1534: ISCreate(comm, &ownersIS);
1535: PetscObjectSetName((PetscObject)ownersIS, "owners");
1536: ISCreate(comm, &gpointsIS);
1537: PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers");
1538: PetscLayoutSetLocalSize(chartSizesIS->map, 1);
1539: ISLoad(chartSizesIS, viewer);
1540: ISGetIndices(chartSizesIS, &chartSize);
1541: PetscLayoutSetLocalSize(ownersIS->map, *chartSize);
1542: PetscLayoutSetLocalSize(gpointsIS->map, *chartSize);
1543: ISLoad(ownersIS, viewer);
1544: ISLoad(gpointsIS, viewer);
1545: ISGetIndices(ownersIS, &owners);
1546: ISGetIndices(gpointsIS, &gpoints);
1547: PetscSFCreate(comm, distsf);
1548: PetscSFSetFromOptions(*distsf);
1549: PetscLayoutCreate(comm, &layout);
1550: PetscSFGetGraph(sf, &lsize, NULL, NULL, NULL);
1551: PetscLayoutSetLocalSize(layout, lsize);
1552: PetscLayoutSetBlockSize(layout, 1);
1553: PetscLayoutSetUp(layout);
1554: PetscSFSetGraphLayout(*distsf, layout, *chartSize, NULL, PETSC_OWN_POINTER, gpoints);
1555: PetscLayoutDestroy(&layout);
1556: /* Migrate DM */
1557: {
1558: PetscInt pStart, pEnd;
1559: PetscSFNode *buffer0, *buffer1, *buffer2;
1561: DMPlexGetChart(dm, &pStart, &pEnd);
1562: PetscMalloc2(pEnd - pStart, &buffer0, lsize, &buffer1);
1563: PetscMalloc1(*chartSize, &buffer2);
1564: {
1565: PetscSF workPointSF;
1566: PetscInt workNroots, workNleaves;
1567: const PetscInt *workIlocal;
1568: const PetscSFNode *workIremote;
1570: for (p = pStart; p < pEnd; ++p) {
1571: buffer0[p - pStart].rank = rank;
1572: buffer0[p - pStart].index = p - pStart;
1573: }
1574: DMGetPointSF(dm, &workPointSF);
1575: PetscSFGetGraph(workPointSF, &workNroots, &workNleaves, &workIlocal, &workIremote);
1576: for (p = 0; p < workNleaves; ++p) {
1577: PetscInt workIlocalp = (workIlocal ? workIlocal[p] : p);
1579: buffer0[workIlocalp].rank = -1;
1580: }
1581: }
1582: for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
1583: for (p = 0; p < *chartSize; ++p) buffer2[p].rank = -1;
1584: PetscSFReduceBegin(sf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC);
1585: PetscSFReduceEnd(sf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC);
1586: PetscSFBcastBegin(*distsf, MPIU_2INT, buffer1, buffer2, MPI_REPLACE);
1587: PetscSFBcastEnd(*distsf, MPIU_2INT, buffer1, buffer2, MPI_REPLACE);
1588: if (PetscDefined(USE_DEBUG)) {
1589: for (p = 0; p < *chartSize; ++p) {
1591: }
1592: }
1593: PetscFree2(buffer0, buffer1);
1594: DMCreate(comm, distdm);
1595: DMSetType(*distdm, DMPLEX);
1596: PetscObjectSetName((PetscObject)(*distdm), ((PetscObject)dm)->name);
1597: DMPlexDistributionSetName(*distdm, distribution_name);
1598: {
1599: PetscSF migrationSF;
1601: PetscSFCreate(comm, &migrationSF);
1602: PetscSFSetFromOptions(migrationSF);
1603: PetscSFSetGraph(migrationSF, pEnd - pStart, *chartSize, NULL, PETSC_OWN_POINTER, buffer2, PETSC_OWN_POINTER);
1604: PetscSFSetUp(migrationSF);
1605: DMPlexMigrate(dm, migrationSF, *distdm);
1606: PetscSFDestroy(&migrationSF);
1607: }
1608: }
1609: /* Set pointSF */
1610: {
1611: PetscSF pointSF;
1612: PetscInt *ilocal, nleaves, q;
1613: PetscSFNode *iremote, *buffer0, *buffer1;
1615: PetscMalloc2(*chartSize, &buffer0, lsize, &buffer1);
1616: for (p = 0, nleaves = 0; p < *chartSize; ++p) {
1617: if (owners[p] == rank) {
1618: buffer0[p].rank = rank;
1619: } else {
1620: buffer0[p].rank = -1;
1621: nleaves++;
1622: }
1623: buffer0[p].index = p;
1624: }
1625: for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
1626: PetscSFReduceBegin(*distsf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC);
1627: PetscSFReduceEnd(*distsf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC);
1628: for (p = 0; p < *chartSize; ++p) buffer0[p].rank = -1;
1629: PetscSFBcastBegin(*distsf, MPIU_2INT, buffer1, buffer0, MPI_REPLACE);
1630: PetscSFBcastEnd(*distsf, MPIU_2INT, buffer1, buffer0, MPI_REPLACE);
1631: if (PetscDefined(USE_DEBUG)) {
1633: }
1634: PetscMalloc1(nleaves, &ilocal);
1635: PetscMalloc1(nleaves, &iremote);
1636: for (p = 0, q = 0; p < *chartSize; ++p) {
1637: if (buffer0[p].rank != rank) {
1638: ilocal[q] = p;
1639: iremote[q].rank = buffer0[p].rank;
1640: iremote[q].index = buffer0[p].index;
1641: q++;
1642: }
1643: }
1645: PetscFree2(buffer0, buffer1);
1646: PetscSFCreate(comm, &pointSF);
1647: PetscSFSetFromOptions(pointSF);
1648: PetscSFSetGraph(pointSF, *chartSize, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1649: DMSetPointSF(*distdm, pointSF);
1650: {
1651: DM cdm;
1653: DMGetCoordinateDM(*distdm, &cdm);
1654: DMSetPointSF(cdm, pointSF);
1655: }
1656: PetscSFDestroy(&pointSF);
1657: }
1658: ISRestoreIndices(chartSizesIS, &chartSize);
1659: ISRestoreIndices(ownersIS, &owners);
1660: ISRestoreIndices(gpointsIS, &gpoints);
1661: ISDestroy(&chartSizesIS);
1662: ISDestroy(&ownersIS);
1663: ISDestroy(&gpointsIS);
1664: /* Record that overlap has been manually created. */
1665: /* This is to pass `DMPlexCheckPointSF()`, which checks that */
1666: /* pointSF does not contain cells in the leaves if overlap = 0. */
1667: DMPlexSetOverlap_Plex(*distdm, NULL, DMPLEX_OVERLAP_MANUAL);
1668: DMPlexDistributeSetDefault(*distdm, PETSC_FALSE);
1669: DMPlexReorderSetDefault(*distdm, DMPLEX_REORDER_DEFAULT_FALSE);
1670: return 0;
1671: }
1673: static PetscErrorCode DMPlexTopologyLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer, PetscSF *sf)
1674: {
1675: MPI_Comm comm;
1676: const char *pointsName, *coneSizesName, *conesName, *orientationsName;
1677: IS pointsIS, coneSizesIS, conesIS, orientationsIS;
1678: const PetscInt *points, *coneSizes, *cones, *orientations;
1679: PetscInt *cone, *ornt;
1680: PetscInt dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
1681: PetscMPIInt size, rank;
1683: pointsName = "order";
1684: coneSizesName = "cones";
1685: conesName = "cells";
1686: orientationsName = "orientation";
1687: PetscObjectGetComm((PetscObject)dm, &comm);
1688: MPI_Comm_size(comm, &size);
1689: MPI_Comm_rank(comm, &rank);
1690: ISCreate(comm, &pointsIS);
1691: PetscObjectSetName((PetscObject)pointsIS, pointsName);
1692: ISCreate(comm, &coneSizesIS);
1693: PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName);
1694: ISCreate(comm, &conesIS);
1695: PetscObjectSetName((PetscObject)conesIS, conesName);
1696: ISCreate(comm, &orientationsIS);
1697: PetscObjectSetName((PetscObject)orientationsIS, orientationsName);
1698: PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject)conesIS, "cell_dim", PETSC_INT, NULL, &dim);
1699: DMSetDimension(dm, dim);
1700: {
1701: /* Force serial load */
1702: PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np);
1703: PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0);
1704: PetscLayoutSetSize(pointsIS->map, Np);
1705: pEnd = rank == 0 ? Np : 0;
1706: PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np);
1707: PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0);
1708: PetscLayoutSetSize(coneSizesIS->map, Np);
1709: PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N);
1710: PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0);
1711: PetscLayoutSetSize(conesIS->map, N);
1712: PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N);
1713: PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0);
1714: PetscLayoutSetSize(orientationsIS->map, N);
1715: }
1716: ISLoad(pointsIS, viewer);
1717: ISLoad(coneSizesIS, viewer);
1718: ISLoad(conesIS, viewer);
1719: ISLoad(orientationsIS, viewer);
1720: /* Create Plex */
1721: DMPlexSetChart(dm, 0, pEnd);
1722: ISGetIndices(pointsIS, &points);
1723: ISGetIndices(coneSizesIS, &coneSizes);
1724: for (p = 0; p < pEnd; ++p) {
1725: DMPlexSetConeSize(dm, points[p], coneSizes[p]);
1726: maxConeSize = PetscMax(maxConeSize, coneSizes[p]);
1727: }
1728: DMSetUp(dm);
1729: ISGetIndices(conesIS, &cones);
1730: ISGetIndices(orientationsIS, &orientations);
1731: PetscMalloc2(maxConeSize, &cone, maxConeSize, &ornt);
1732: for (p = 0, q = 0; p < pEnd; ++p) {
1733: for (c = 0; c < coneSizes[p]; ++c, ++q) {
1734: cone[c] = cones[q];
1735: ornt[c] = orientations[q];
1736: }
1737: DMPlexSetCone(dm, points[p], cone);
1738: DMPlexSetConeOrientation(dm, points[p], ornt);
1739: }
1740: PetscFree2(cone, ornt);
1741: /* Create global section migration SF */
1742: if (sf) {
1743: PetscLayout layout;
1744: PetscInt *globalIndices;
1746: PetscMalloc1(pEnd, &globalIndices);
1747: /* plex point == globalPointNumber in this case */
1748: for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
1749: PetscLayoutCreate(comm, &layout);
1750: PetscLayoutSetSize(layout, Np);
1751: PetscLayoutSetBlockSize(layout, 1);
1752: PetscLayoutSetUp(layout);
1753: PetscSFCreate(comm, sf);
1754: PetscSFSetFromOptions(*sf);
1755: PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices);
1756: PetscLayoutDestroy(&layout);
1757: PetscFree(globalIndices);
1758: }
1759: /* Clean-up */
1760: ISRestoreIndices(pointsIS, &points);
1761: ISRestoreIndices(coneSizesIS, &coneSizes);
1762: ISRestoreIndices(conesIS, &cones);
1763: ISRestoreIndices(orientationsIS, &orientations);
1764: ISDestroy(&pointsIS);
1765: ISDestroy(&coneSizesIS);
1766: ISDestroy(&conesIS);
1767: ISDestroy(&orientationsIS);
1768: /* Fill in the rest of the topology structure */
1769: DMPlexSymmetrize(dm);
1770: DMPlexStratify(dm);
1771: return 0;
1772: }
1774: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sfXC)
1775: {
1776: DMPlexStorageVersion version;
1777: const char *topologydm_name;
1778: char group[PETSC_MAX_PATH_LEN];
1779: PetscSF sfwork = NULL;
1781: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1782: PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version);
1783: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1784: PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name);
1785: } else {
1786: PetscStrcpy(group, "/");
1787: }
1788: PetscViewerHDF5PushGroup(viewer, group);
1790: PetscViewerHDF5PushGroup(viewer, "topology");
1791: DMPlexTopologyLoad_HDF5_Legacy_Private(dm, viewer, &sfwork);
1792: PetscViewerHDF5PopGroup(viewer); /* "topology" */
1794: if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1795: DM distdm;
1796: PetscSF distsf;
1797: const char *distribution_name;
1798: PetscBool exists;
1800: DMPlexDistributionGetName(dm, &distribution_name);
1801: /* this group is guaranteed to be present since DMPlexStorageVersion 2.1.0 */
1802: PetscViewerHDF5PushGroup(viewer, "distributions");
1803: PetscViewerHDF5HasGroup(viewer, NULL, &exists);
1804: if (exists) {
1805: PetscViewerHDF5PushGroup(viewer, distribution_name);
1806: DMPlexDistributionLoad_HDF5_Private(dm, viewer, sfwork, &distsf, &distdm);
1807: if (distdm) {
1808: DMPlexReplace_Internal(dm, &distdm);
1809: PetscSFDestroy(&sfwork);
1810: sfwork = distsf;
1811: }
1812: PetscViewerHDF5PopGroup(viewer); /* distribution_name */
1813: }
1814: PetscViewerHDF5PopGroup(viewer); /* "distributions" */
1815: }
1816: if (sfXC) {
1817: *sfXC = sfwork;
1818: } else {
1819: PetscSFDestroy(&sfwork);
1820: }
1822: PetscViewerHDF5PopGroup(viewer);
1823: return 0;
1824: }
1826: /* If the file is old, it not only has different path to the coordinates, but */
1827: /* does not contain coordinateDMs, so must fall back to the old implementation. */
1828: static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
1829: {
1830: PetscSection coordSection;
1831: Vec coordinates;
1832: PetscReal lengthScale;
1833: PetscInt spatialDim, N, numVertices, vStart, vEnd, v;
1834: PetscMPIInt rank;
1836: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1837: /* Read geometry */
1838: PetscViewerHDF5PushGroup(viewer, "/geometry");
1839: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1840: PetscObjectSetName((PetscObject)coordinates, "vertices");
1841: {
1842: /* Force serial load */
1843: PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
1844: VecSetSizes(coordinates, rank == 0 ? N : 0, N);
1845: VecSetBlockSize(coordinates, spatialDim);
1846: }
1847: VecLoad(coordinates, viewer);
1848: PetscViewerHDF5PopGroup(viewer);
1849: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1850: VecScale(coordinates, 1.0 / lengthScale);
1851: VecGetLocalSize(coordinates, &numVertices);
1852: VecGetBlockSize(coordinates, &spatialDim);
1853: numVertices /= spatialDim;
1854: /* Create coordinates */
1855: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1857: DMGetCoordinateSection(dm, &coordSection);
1858: PetscSectionSetNumFields(coordSection, 1);
1859: PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
1860: PetscSectionSetChart(coordSection, vStart, vEnd);
1861: for (v = vStart; v < vEnd; ++v) {
1862: PetscSectionSetDof(coordSection, v, spatialDim);
1863: PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
1864: }
1865: PetscSectionSetUp(coordSection);
1866: DMSetCoordinates(dm, coordinates);
1867: VecDestroy(&coordinates);
1868: return 0;
1869: }
1871: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1872: {
1873: DMPlexStorageVersion version;
1874: DM cdm;
1875: Vec coords;
1876: PetscInt blockSize;
1877: PetscReal lengthScale;
1878: PetscSF lsf;
1879: const char *topologydm_name;
1880: char *coordinatedm_name, *coordinates_name;
1882: PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version);
1883: if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
1884: DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer);
1885: return 0;
1886: }
1887: /* else: since DMPlexStorageVersion 2.0.0 */
1889: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1890: PetscViewerHDF5PushGroup(viewer, "topologies");
1891: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1892: PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, NULL, &coordinatedm_name);
1893: PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, NULL, &coordinates_name);
1894: PetscViewerHDF5PopGroup(viewer);
1895: PetscViewerHDF5PopGroup(viewer);
1896: DMGetCoordinateDM(dm, &cdm);
1897: PetscObjectSetName((PetscObject)cdm, coordinatedm_name);
1898: PetscFree(coordinatedm_name);
1899: /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
1900: DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf);
1901: DMCreateLocalVector(cdm, &coords);
1902: PetscObjectSetName((PetscObject)coords, coordinates_name);
1903: PetscFree(coordinates_name);
1904: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
1905: DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords);
1906: PetscViewerPopFormat(viewer);
1907: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1908: VecScale(coords, 1.0 / lengthScale);
1909: DMSetCoordinatesLocal(dm, coords);
1910: VecGetBlockSize(coords, &blockSize);
1911: DMSetCoordinateDim(dm, blockSize);
1912: VecDestroy(&coords);
1913: PetscSFDestroy(&lsf);
1914: return 0;
1915: }
1917: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1918: {
1919: DMPlexStorageVersion version;
1921: PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version);
1922: if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
1923: DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL);
1924: DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL);
1925: DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer);
1926: } else {
1927: PetscSF sfXC;
1929: /* since DMPlexStorageVersion 2.0.0 */
1930: DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC);
1931: DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC);
1932: DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC);
1933: PetscSFDestroy(&sfXC);
1934: }
1935: return 0;
1936: }
1938: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1939: {
1940: MPI_Comm comm;
1941: PetscInt pStart, pEnd, p, m;
1942: PetscInt *goffs, *ilocal;
1943: PetscBool rootIncludeConstraints, leafIncludeConstraints;
1945: PetscObjectGetComm((PetscObject)leafSection, &comm);
1946: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1947: PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints);
1948: PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints);
1949: if (rootIncludeConstraints && leafIncludeConstraints) PetscSectionGetStorageSize(leafSection, &m);
1950: else PetscSectionGetConstrainedStorageSize(leafSection, &m);
1951: PetscMalloc1(m, &ilocal);
1952: PetscMalloc1(m, &goffs);
1953: /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
1954: /* for the top-level section (not for each field), so one must have */
1955: /* rootSection->pointMajor == PETSC_TRUE. */
1957: /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
1959: for (p = pStart, m = 0; p < pEnd; ++p) {
1960: PetscInt dof, cdof, i, j, off, goff;
1961: const PetscInt *cinds;
1963: PetscSectionGetDof(leafSection, p, &dof);
1964: if (dof < 0) continue;
1965: goff = globalOffsets[p - pStart];
1966: PetscSectionGetOffset(leafSection, p, &off);
1967: PetscSectionGetConstraintDof(leafSection, p, &cdof);
1968: PetscSectionGetConstraintIndices(leafSection, p, &cinds);
1969: for (i = 0, j = 0; i < dof; ++i) {
1970: PetscBool constrained = (PetscBool)(j < cdof && i == cinds[j]);
1972: if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {
1973: ilocal[m] = off++;
1974: goffs[m++] = goff++;
1975: } else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
1976: else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff;
1977: if (constrained) ++j;
1978: }
1979: }
1980: PetscSFCreate(comm, sectionSF);
1981: PetscSFSetFromOptions(*sectionSF);
1982: PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs);
1983: PetscFree(goffs);
1984: return 0;
1985: }
1987: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
1988: {
1989: MPI_Comm comm;
1990: PetscMPIInt size, rank;
1991: const char *topologydm_name;
1992: const char *sectiondm_name;
1993: PetscSection sectionA, sectionB;
1994: PetscInt nX, n, i;
1995: PetscSF sfAB;
1997: PetscObjectGetComm((PetscObject)dm, &comm);
1998: MPI_Comm_size(comm, &size);
1999: MPI_Comm_rank(comm, &rank);
2000: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
2001: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
2002: PetscViewerHDF5PushGroup(viewer, "topologies");
2003: PetscViewerHDF5PushGroup(viewer, topologydm_name);
2004: PetscViewerHDF5PushGroup(viewer, "dms");
2005: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
2006: /* A: on-disk points */
2007: /* X: list of global point numbers, [0, NX) */
2008: /* B: plex points */
2009: /* Load raw section (sectionA) */
2010: PetscSectionCreate(comm, §ionA);
2011: PetscSectionLoad(sectionA, viewer);
2012: PetscSectionGetChart(sectionA, NULL, &n);
2013: /* Create sfAB: A -> B */
2014: #if defined(PETSC_USE_DEBUG)
2015: {
2016: PetscInt N, N1;
2018: PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1);
2019: MPI_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm);
2021: }
2022: #endif
2023: {
2024: IS orderIS;
2025: const PetscInt *gpoints;
2026: PetscSF sfXA, sfAX;
2027: PetscLayout layout;
2028: PetscSFNode *owners, *buffer;
2029: PetscInt nleaves;
2030: PetscInt *ilocal;
2031: PetscSFNode *iremote;
2033: /* Create sfAX: A -> X */
2034: ISCreate(comm, &orderIS);
2035: PetscObjectSetName((PetscObject)orderIS, "order");
2036: PetscLayoutSetLocalSize(orderIS->map, n);
2037: ISLoad(orderIS, viewer);
2038: PetscLayoutCreate(comm, &layout);
2039: PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL);
2040: PetscLayoutSetLocalSize(layout, nX);
2041: PetscLayoutSetBlockSize(layout, 1);
2042: PetscLayoutSetUp(layout);
2043: PetscSFCreate(comm, &sfXA);
2044: ISGetIndices(orderIS, &gpoints);
2045: PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints);
2046: ISRestoreIndices(orderIS, &gpoints);
2047: ISDestroy(&orderIS);
2048: PetscLayoutDestroy(&layout);
2049: PetscMalloc1(n, &owners);
2050: PetscMalloc1(nX, &buffer);
2051: for (i = 0; i < n; ++i) {
2052: owners[i].rank = rank;
2053: owners[i].index = i;
2054: }
2055: for (i = 0; i < nX; ++i) {
2056: buffer[i].rank = -1;
2057: buffer[i].index = -1;
2058: }
2059: PetscSFReduceBegin(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
2060: PetscSFReduceEnd(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
2061: PetscSFDestroy(&sfXA);
2062: PetscFree(owners);
2063: for (i = 0, nleaves = 0; i < nX; ++i)
2064: if (buffer[i].rank >= 0) nleaves++;
2065: PetscMalloc1(nleaves, &ilocal);
2066: PetscMalloc1(nleaves, &iremote);
2067: for (i = 0, nleaves = 0; i < nX; ++i) {
2068: if (buffer[i].rank >= 0) {
2069: ilocal[nleaves] = i;
2070: iremote[nleaves].rank = buffer[i].rank;
2071: iremote[nleaves].index = buffer[i].index;
2072: nleaves++;
2073: }
2074: }
2075: PetscFree(buffer);
2076: PetscSFCreate(comm, &sfAX);
2077: PetscSFSetFromOptions(sfAX);
2078: PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
2079: PetscSFCompose(sfAX, sfXB, &sfAB);
2080: PetscSFDestroy(&sfAX);
2081: }
2082: PetscViewerHDF5PopGroup(viewer);
2083: PetscViewerHDF5PopGroup(viewer);
2084: PetscViewerHDF5PopGroup(viewer);
2085: PetscViewerHDF5PopGroup(viewer);
2086: /* Create plex section (sectionB) */
2087: DMGetLocalSection(sectiondm, §ionB);
2088: if (lsf || gsf) {
2089: PetscLayout layout;
2090: PetscInt M, m;
2091: PetscInt *offsetsA;
2092: PetscBool includesConstraintsA;
2094: PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB);
2095: PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA);
2096: if (includesConstraintsA) PetscSectionGetStorageSize(sectionA, &m);
2097: else PetscSectionGetConstrainedStorageSize(sectionA, &m);
2098: MPI_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm);
2099: PetscLayoutCreate(comm, &layout);
2100: PetscLayoutSetSize(layout, M);
2101: PetscLayoutSetUp(layout);
2102: if (lsf) {
2103: PetscSF lsfABdata;
2105: DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata);
2106: *lsf = lsfABdata;
2107: }
2108: if (gsf) {
2109: PetscSection gsectionB, gsectionB1;
2110: PetscBool includesConstraintsB;
2111: PetscSF gsfABdata, pointsf;
2113: DMGetGlobalSection(sectiondm, &gsectionB1);
2114: PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB);
2115: DMGetPointSF(sectiondm, &pointsf);
2116: PetscSectionCreateGlobalSection(sectionB, pointsf, includesConstraintsB, PETSC_TRUE, &gsectionB);
2117: DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata);
2118: PetscSectionDestroy(&gsectionB);
2119: *gsf = gsfABdata;
2120: }
2121: PetscLayoutDestroy(&layout);
2122: PetscFree(offsetsA);
2123: } else {
2124: PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB);
2125: }
2126: PetscSFDestroy(&sfAB);
2127: PetscSectionDestroy(§ionA);
2128: return 0;
2129: }
2131: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
2132: {
2133: MPI_Comm comm;
2134: const char *topologydm_name;
2135: const char *sectiondm_name;
2136: const char *vec_name;
2137: Vec vecA;
2138: PetscInt mA, m, bs;
2139: const PetscInt *ilocal;
2140: const PetscScalar *src;
2141: PetscScalar *dest;
2143: PetscObjectGetComm((PetscObject)dm, &comm);
2144: DMPlexGetHDF5Name_Private(dm, &topologydm_name);
2145: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
2146: PetscObjectGetName((PetscObject)vec, &vec_name);
2147: PetscViewerHDF5PushGroup(viewer, "topologies");
2148: PetscViewerHDF5PushGroup(viewer, topologydm_name);
2149: PetscViewerHDF5PushGroup(viewer, "dms");
2150: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
2151: PetscViewerHDF5PushGroup(viewer, "vecs");
2152: PetscViewerHDF5PushGroup(viewer, vec_name);
2153: VecCreate(comm, &vecA);
2154: PetscObjectSetName((PetscObject)vecA, vec_name);
2155: PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL);
2156: /* Check consistency */
2157: {
2158: PetscSF pointsf, pointsf1;
2159: PetscInt m1, i, j;
2161: DMGetPointSF(dm, &pointsf);
2162: DMGetPointSF(sectiondm, &pointsf1);
2164: #if defined(PETSC_USE_DEBUG)
2165: {
2166: PetscInt MA, MA1;
2168: MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm);
2169: PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1);
2171: }
2172: #endif
2173: VecGetLocalSize(vec, &m1);
2175: for (i = 0; i < m; ++i) {
2176: j = ilocal ? ilocal[i] : i;
2178: }
2179: }
2180: VecSetSizes(vecA, mA, PETSC_DECIDE);
2181: VecLoad(vecA, viewer);
2182: VecGetArrayRead(vecA, &src);
2183: VecGetArray(vec, &dest);
2184: PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
2185: PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
2186: VecRestoreArray(vec, &dest);
2187: VecRestoreArrayRead(vecA, &src);
2188: VecDestroy(&vecA);
2189: PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *)&bs);
2190: VecSetBlockSize(vec, bs);
2191: PetscViewerHDF5PopGroup(viewer);
2192: PetscViewerHDF5PopGroup(viewer);
2193: PetscViewerHDF5PopGroup(viewer);
2194: PetscViewerHDF5PopGroup(viewer);
2195: PetscViewerHDF5PopGroup(viewer);
2196: PetscViewerHDF5PopGroup(viewer);
2197: return 0;
2198: }
2199: #endif