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: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
9: #if defined(PETSC_HAVE_HDF5)
10: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
11: {
12: Vec stamp;
13: PetscMPIInt rank;
17: if (seqnum < 0) return(0);
18: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
19: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
20: VecSetBlockSize(stamp, 1);
21: PetscObjectSetName((PetscObject) stamp, seqname);
22: if (rank == 0) {
23: PetscReal timeScale;
24: PetscBool istime;
26: PetscStrncmp(seqname, "time", 5, &istime);
27: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); value *= timeScale;}
28: VecSetValue(stamp, 0, value, INSERT_VALUES);
29: }
30: VecAssemblyBegin(stamp);
31: VecAssemblyEnd(stamp);
32: PetscViewerHDF5PushGroup(viewer, "/");
33: PetscViewerHDF5PushTimestepping(viewer);
34: PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
35: VecView(stamp, viewer);
36: PetscViewerHDF5PopTimestepping(viewer);
37: PetscViewerHDF5PopGroup(viewer);
38: VecDestroy(&stamp);
39: return(0);
40: }
42: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
43: {
44: Vec stamp;
45: PetscMPIInt rank;
49: if (seqnum < 0) return(0);
50: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
51: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
52: VecSetBlockSize(stamp, 1);
53: PetscObjectSetName((PetscObject) stamp, seqname);
54: PetscViewerHDF5PushGroup(viewer, "/");
55: PetscViewerHDF5PushTimestepping(viewer);
56: PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
57: VecLoad(stamp, viewer);
58: PetscViewerHDF5PopTimestepping(viewer);
59: PetscViewerHDF5PopGroup(viewer);
60: if (rank == 0) {
61: const PetscScalar *a;
62: PetscReal timeScale;
63: PetscBool istime;
65: VecGetArrayRead(stamp, &a);
66: *value = a[0];
67: VecRestoreArrayRead(stamp, &a);
68: PetscStrncmp(seqname, "time", 5, &istime);
69: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); *value /= timeScale;}
70: }
71: VecDestroy(&stamp);
72: return(0);
73: }
75: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
76: {
77: IS cutcells = NULL;
78: const PetscInt *cutc;
79: PetscInt cellHeight, vStart, vEnd, cStart, cEnd, c;
80: PetscErrorCode ierr;
83: if (!cutLabel) return(0);
84: DMPlexGetVTKCellHeight(dm, &cellHeight);
85: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
86: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
87: /* Label vertices that should be duplicated */
88: DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel);
89: DMLabelGetStratumIS(cutLabel, 2, &cutcells);
90: if (cutcells) {
91: PetscInt n;
93: ISGetIndices(cutcells, &cutc);
94: ISGetLocalSize(cutcells, &n);
95: for (c = 0; c < n; ++c) {
96: if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
97: PetscInt *closure = NULL;
98: PetscInt closureSize, cl, value;
100: DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
101: for (cl = 0; cl < closureSize*2; cl += 2) {
102: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
103: DMLabelGetValue(cutLabel, closure[cl], &value);
104: if (value == 1) {
105: DMLabelSetValue(*cutVertexLabel, closure[cl], 1);
106: }
107: }
108: }
109: DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
110: }
111: }
112: ISRestoreIndices(cutcells, &cutc);
113: }
114: ISDestroy(&cutcells);
115: return(0);
116: }
118: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
119: {
120: DM dm;
121: DM dmBC;
122: PetscSection section, sectionGlobal;
123: Vec gv;
124: const char *name;
125: PetscViewerVTKFieldType ft;
126: PetscViewerFormat format;
127: PetscInt seqnum;
128: PetscReal seqval;
129: PetscBool isseq;
130: PetscErrorCode ierr;
133: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
134: VecGetDM(v, &dm);
135: DMGetLocalSection(dm, §ion);
136: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
137: DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
138: if (seqnum >= 0) {
139: PetscViewerHDF5PushTimestepping(viewer);
140: PetscViewerHDF5SetTimestep(viewer, seqnum);
141: }
142: PetscViewerGetFormat(viewer, &format);
143: DMGetOutputDM(dm, &dmBC);
144: DMGetGlobalSection(dmBC, §ionGlobal);
145: DMGetGlobalVector(dmBC, &gv);
146: PetscObjectGetName((PetscObject) v, &name);
147: PetscObjectSetName((PetscObject) gv, name);
148: DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
149: DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
150: PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
151: if (isseq) {VecView_Seq(gv, viewer);}
152: else {VecView_MPI(gv, viewer);}
153: if (format == PETSC_VIEWER_HDF5_VIZ) {
154: /* Output visualization representation */
155: PetscInt numFields, f;
156: DMLabel cutLabel, cutVertexLabel = NULL;
158: PetscSectionGetNumFields(section, &numFields);
159: DMGetLabel(dm, "periodic_cut", &cutLabel);
160: for (f = 0; f < numFields; ++f) {
161: Vec subv;
162: IS is;
163: const char *fname, *fgroup;
164: char subname[PETSC_MAX_PATH_LEN];
165: PetscInt pStart, pEnd;
167: DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
168: fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
169: PetscSectionGetFieldName(section, f, &fname);
170: if (!fname) continue;
171: PetscViewerHDF5PushGroup(viewer, fgroup);
172: if (cutLabel) {
173: const PetscScalar *ga;
174: PetscScalar *suba;
175: PetscInt Nc, gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;
177: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
178: PetscSectionGetFieldComponents(section, f, &Nc);
179: for (p = pStart; p < pEnd; ++p) {
180: PetscInt gdof, fdof = 0, val;
182: PetscSectionGetDof(sectionGlobal, p, &gdof);
183: if (gdof > 0) {PetscSectionGetFieldDof(section, p, f, &fdof);}
184: subSize += fdof;
185: DMLabelGetValue(cutVertexLabel, p, &val);
186: if (val == 1) extSize += fdof;
187: }
188: VecCreate(PetscObjectComm((PetscObject) gv), &subv);
189: VecSetSizes(subv, subSize+extSize, PETSC_DETERMINE);
190: VecSetBlockSize(subv, Nc);
191: VecSetType(subv, VECSTANDARD);
192: VecGetOwnershipRange(gv, &gstart, NULL);
193: VecGetArrayRead(gv, &ga);
194: VecGetArray(subv, &suba);
195: for (p = pStart; p < pEnd; ++p) {
196: PetscInt gdof, goff, val;
198: PetscSectionGetDof(sectionGlobal, p, &gdof);
199: if (gdof > 0) {
200: PetscInt fdof, fc, f2, poff = 0;
202: PetscSectionGetOffset(sectionGlobal, p, &goff);
203: /* Can get rid of this loop by storing field information in the global section */
204: for (f2 = 0; f2 < f; ++f2) {
205: PetscSectionGetFieldDof(section, p, f2, &fdof);
206: poff += fdof;
207: }
208: PetscSectionGetFieldDof(section, p, f, &fdof);
209: for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff+poff+fc - gstart];
210: DMLabelGetValue(cutVertexLabel, p, &val);
211: if (val == 1) {
212: for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize+newOff] = ga[goff+poff+fc - gstart];
213: }
214: }
215: }
216: VecRestoreArrayRead(gv, &ga);
217: VecRestoreArray(subv, &suba);
218: DMLabelDestroy(&cutVertexLabel);
219: } else {
220: PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
221: }
222: PetscStrncpy(subname, name,sizeof(subname));
223: PetscStrlcat(subname, "_",sizeof(subname));
224: PetscStrlcat(subname, fname,sizeof(subname));
225: PetscObjectSetName((PetscObject) subv, subname);
226: if (isseq) {VecView_Seq(subv, viewer);}
227: else {VecView_MPI(subv, viewer);}
228: if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
229: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "vector");
230: } else {
231: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "scalar");
232: }
233: if (cutLabel) {VecDestroy(&subv);}
234: else {PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);}
235: PetscViewerHDF5PopGroup(viewer);
236: }
237: }
238: if (seqnum >= 0) {
239: PetscViewerHDF5PopTimestepping(viewer);
240: }
241: DMRestoreGlobalVector(dmBC, &gv);
242: return(0);
243: }
245: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
246: {
247: DM dm;
248: Vec locv;
249: PetscObject isZero;
250: const char *name;
251: PetscReal time;
255: VecGetDM(v, &dm);
256: DMGetLocalVector(dm, &locv);
257: PetscObjectGetName((PetscObject) v, &name);
258: PetscObjectSetName((PetscObject) locv, name);
259: PetscObjectQuery((PetscObject) v, "__Vec_bc_zero__", &isZero);
260: PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", isZero);
261: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
262: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
263: DMGetOutputSequenceNumber(dm, NULL, &time);
264: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
265: PetscViewerHDF5PushGroup(viewer, "/fields");
266: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
267: VecView_Plex_Local_HDF5_Internal(locv, viewer);
268: PetscViewerPopFormat(viewer);
269: PetscViewerHDF5PopGroup(viewer);
270: PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", NULL);
271: DMRestoreLocalVector(dm, &locv);
272: return(0);
273: }
275: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
276: {
277: PetscBool isseq;
281: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
282: PetscViewerHDF5PushGroup(viewer, "/fields");
283: if (isseq) {VecView_Seq(v, viewer);}
284: else {VecView_MPI(v, viewer);}
285: PetscViewerHDF5PopGroup(viewer);
286: return(0);
287: }
289: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
290: {
291: DM dm;
292: Vec locv;
293: const char *name;
294: PetscInt seqnum;
298: VecGetDM(v, &dm);
299: DMGetLocalVector(dm, &locv);
300: PetscObjectGetName((PetscObject) v, &name);
301: PetscObjectSetName((PetscObject) locv, name);
302: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
303: PetscViewerHDF5PushGroup(viewer, "/fields");
304: if (seqnum >= 0) {
305: PetscViewerHDF5PushTimestepping(viewer);
306: PetscViewerHDF5SetTimestep(viewer, seqnum);
307: }
308: VecLoad_Plex_Local(locv, viewer);
309: if (seqnum >= 0) {
310: PetscViewerHDF5PopTimestepping(viewer);
311: }
312: PetscViewerHDF5PopGroup(viewer);
313: DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
314: DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
315: DMRestoreLocalVector(dm, &locv);
316: return(0);
317: }
319: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
320: {
321: DM dm;
322: PetscInt seqnum;
326: VecGetDM(v, &dm);
327: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
328: PetscViewerHDF5PushGroup(viewer, "/fields");
329: if (seqnum >= 0) {
330: PetscViewerHDF5PushTimestepping(viewer);
331: PetscViewerHDF5SetTimestep(viewer, seqnum);
332: }
333: VecLoad_Default(v, viewer);
334: if (seqnum >= 0) {
335: PetscViewerHDF5PopTimestepping(viewer);
336: }
337: PetscViewerHDF5PopGroup(viewer);
338: return(0);
339: }
341: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
342: {
343: IS orderIS, conesIS, cellsIS, orntsIS;
344: const PetscInt *gpoint;
345: PetscInt *order, *sizes, *cones, *ornts;
346: PetscInt dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
347: PetscErrorCode ierr;
350: ISGetIndices(globalPointNumbers, &gpoint);
351: DMGetDimension(dm, &dim);
352: DMPlexGetChart(dm, &pStart, &pEnd);
353: for (p = pStart; p < pEnd; ++p) {
354: if (gpoint[p] >= 0) {
355: PetscInt coneSize;
357: DMPlexGetConeSize(dm, p, &coneSize);
358: conesSize += 1;
359: cellsSize += coneSize;
360: }
361: }
362: PetscMalloc1(conesSize, &order);
363: PetscMalloc1(conesSize, &sizes);
364: PetscMalloc1(cellsSize, &cones);
365: PetscMalloc1(cellsSize, &ornts);
366: for (p = pStart; p < pEnd; ++p) {
367: if (gpoint[p] >= 0) {
368: const PetscInt *cone, *ornt;
369: PetscInt coneSize, cp;
371: DMPlexGetConeSize(dm, p, &coneSize);
372: DMPlexGetCone(dm, p, &cone);
373: DMPlexGetConeOrientation(dm, p, &ornt);
374: order[s] = gpoint[p];
375: sizes[s++] = coneSize;
376: for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
377: }
378: }
379: if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %D != %D", s, conesSize);
380: if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %D != %D", c, cellsSize);
381: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
382: PetscObjectSetName((PetscObject) orderIS, "order");
383: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
384: PetscObjectSetName((PetscObject) conesIS, "cones");
385: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
386: PetscObjectSetName((PetscObject) cellsIS, "cells");
387: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
388: PetscObjectSetName((PetscObject) orntsIS, "orientation");
389: PetscViewerHDF5PushGroup(viewer, "/topology");
390: ISView(orderIS, viewer);
391: ISView(conesIS, viewer);
392: ISView(cellsIS, viewer);
393: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, (void *) &dim);
394: ISView(orntsIS, viewer);
395: PetscViewerHDF5PopGroup(viewer);
396: ISDestroy(&orderIS);
397: ISDestroy(&conesIS);
398: ISDestroy(&cellsIS);
399: ISDestroy(&orntsIS);
400: ISRestoreIndices(globalPointNumbers, &gpoint);
401: return(0);
402: }
404: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
405: {
406: PetscSF sfPoint;
407: DMLabel cutLabel, cutVertexLabel = NULL;
408: IS globalVertexNumbers, cutvertices = NULL;
409: const PetscInt *gcell, *gvertex, *cutverts = NULL;
410: PetscInt *vertices;
411: PetscInt conesSize = 0;
412: PetscInt dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
413: PetscErrorCode ierr;
416: *numCorners = 0;
417: DMGetDimension(dm, &dim);
418: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
419: ISGetIndices(globalCellNumbers, &gcell);
421: for (cell = cStart; cell < cEnd; ++cell) {
422: PetscInt *closure = NULL;
423: PetscInt closureSize, v, Nc = 0;
425: if (gcell[cell] < 0) continue;
426: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
427: for (v = 0; v < closureSize*2; v += 2) {
428: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
429: }
430: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
431: conesSize += Nc;
432: if (!numCornersLocal) numCornersLocal = Nc;
433: else if (numCornersLocal != Nc) numCornersLocal = 1;
434: }
435: MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
436: if (numCornersLocal && (numCornersLocal != *numCorners || *numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
437: /* Handle periodic cuts by identifying vertices which should be duplicated */
438: DMGetLabel(dm, "periodic_cut", &cutLabel);
439: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
440: if (cutVertexLabel) {DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);}
441: if (cutvertices) {
442: ISGetIndices(cutvertices, &cutverts);
443: ISGetLocalSize(cutvertices, &vExtra);
444: }
445: DMGetPointSF(dm, &sfPoint);
446: if (cutLabel) {
447: const PetscInt *ilocal;
448: const PetscSFNode *iremote;
449: PetscInt nroots, nleaves;
451: PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
452: if (nleaves < 0) {
453: PetscObjectReference((PetscObject) sfPoint);
454: } else {
455: PetscSFCreate(PetscObjectComm((PetscObject) sfPoint), &sfPoint);
456: PetscSFSetGraph(sfPoint, nroots+vExtra, nleaves, ilocal, PETSC_USE_POINTER, iremote, PETSC_USE_POINTER);
457: }
458: } else {
459: PetscObjectReference((PetscObject) sfPoint);
460: }
461: /* Number all vertices */
462: DMPlexCreateNumbering_Plex(dm, vStart, vEnd+vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
463: PetscSFDestroy(&sfPoint);
464: /* Create cones */
465: ISGetIndices(globalVertexNumbers, &gvertex);
466: PetscMalloc1(conesSize, &vertices);
467: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
468: PetscInt *closure = NULL;
469: PetscInt closureSize, Nc = 0, p, value = -1;
470: PetscBool replace;
472: if (gcell[cell] < 0) continue;
473: if (cutLabel) {DMLabelGetValue(cutLabel, cell, &value);}
474: replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
475: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
476: for (p = 0; p < closureSize*2; p += 2) {
477: if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
478: closure[Nc++] = closure[p];
479: }
480: }
481: DMPlexReorderCell(dm, cell, closure);
482: for (p = 0; p < Nc; ++p) {
483: PetscInt nv, gv = gvertex[closure[p] - vStart];
485: if (replace) {
486: PetscFindInt(closure[p], vExtra, cutverts, &nv);
487: if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
488: }
489: vertices[v++] = gv < 0 ? -(gv+1) : gv;
490: }
491: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
492: }
493: ISRestoreIndices(globalVertexNumbers, &gvertex);
494: ISDestroy(&globalVertexNumbers);
495: ISRestoreIndices(globalCellNumbers, &gcell);
496: if (cutvertices) {ISRestoreIndices(cutvertices, &cutverts);}
497: ISDestroy(&cutvertices);
498: DMLabelDestroy(&cutVertexLabel);
499: if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %D != %D", v, conesSize);
500: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS);
501: PetscLayoutSetBlockSize((*cellIS)->map, *numCorners);
502: PetscObjectSetName((PetscObject) *cellIS, "cells");
503: return(0);
504: }
506: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, IS globalCellNumbers, PetscViewer viewer)
507: {
508: DM cdm;
509: DMLabel depthLabel, ctLabel;
510: IS cellIS;
511: PetscInt dim, depth, cellHeight, c;
512: hid_t fileId, groupId;
513: PetscErrorCode ierr;
516: PetscViewerHDF5PushGroup(viewer, "/viz");
517: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
518: PetscStackCallHDF5(H5Gclose,(groupId));
520: PetscViewerHDF5PopGroup(viewer);
521: DMGetDimension(dm, &dim);
522: DMPlexGetDepth(dm, &depth);
523: DMGetCoordinateDM(dm, &cdm);
524: DMPlexGetVTKCellHeight(dm, &cellHeight);
525: DMPlexGetDepthLabel(dm, &depthLabel);
526: DMPlexGetCellTypeLabel(dm, &ctLabel);
527: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
528: const DMPolytopeType ict = (DMPolytopeType) c;
529: PetscInt pStart, pEnd, dep, numCorners, n = 0;
530: PetscBool output = PETSC_FALSE, doOutput;
532: if (ict == DM_POLYTOPE_FV_GHOST) continue;
533: DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
534: if (pStart >= 0) {
535: DMLabelGetValue(depthLabel, pStart, &dep);
536: if (dep == depth - cellHeight) output = PETSC_TRUE;
537: }
538: MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
539: if (!doOutput) continue;
540: CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS);
541: if (!n) {
542: PetscViewerHDF5PushGroup(viewer, "/viz/topology");
543: } else {
544: char group[PETSC_MAX_PATH_LEN];
546: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%D", n);
547: PetscViewerHDF5PushGroup(viewer, group);
548: }
549: ISView(cellIS, viewer);
550: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_corners", PETSC_INT, (void *) &numCorners);
551: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_dim", PETSC_INT, (void *) &dim);
552: ISDestroy(&cellIS);
553: PetscViewerHDF5PopGroup(viewer);
554: ++n;
555: }
556: return(0);
557: }
559: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
560: {
561: DM cdm;
562: Vec coordinates, newcoords;
563: PetscReal lengthScale;
564: PetscInt m, M, bs;
568: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
569: DMGetCoordinateDM(dm, &cdm);
570: DMGetCoordinates(dm, &coordinates);
571: VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
572: PetscObjectSetName((PetscObject) newcoords, "vertices");
573: VecGetSize(coordinates, &M);
574: VecGetLocalSize(coordinates, &m);
575: VecSetSizes(newcoords, m, M);
576: VecGetBlockSize(coordinates, &bs);
577: VecSetBlockSize(newcoords, bs);
578: VecSetType(newcoords,VECSTANDARD);
579: VecCopy(coordinates, newcoords);
580: VecScale(newcoords, lengthScale);
581: /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
582: PetscViewerHDF5PushGroup(viewer, "/geometry");
583: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
584: VecView(newcoords, viewer);
585: PetscViewerPopFormat(viewer);
586: PetscViewerHDF5PopGroup(viewer);
587: VecDestroy(&newcoords);
588: return(0);
589: }
591: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
592: {
593: DM cdm;
594: Vec coordinatesLocal, newcoords;
595: PetscSection cSection, cGlobalSection;
596: PetscScalar *coords, *ncoords;
597: DMLabel cutLabel, cutVertexLabel = NULL;
598: const PetscReal *L;
599: const DMBoundaryType *bd;
600: PetscReal lengthScale;
601: PetscInt vStart, vEnd, v, bs, N, coordSize, dof, off, d;
602: PetscBool localized, embedded;
603: hid_t fileId, groupId;
604: PetscErrorCode ierr;
607: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
608: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
609: DMGetCoordinatesLocal(dm, &coordinatesLocal);
610: VecGetBlockSize(coordinatesLocal, &bs);
611: DMGetCoordinatesLocalized(dm, &localized);
612: if (localized == PETSC_FALSE) return(0);
613: DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
614: DMGetCoordinateDM(dm, &cdm);
615: DMGetLocalSection(cdm, &cSection);
616: DMGetGlobalSection(cdm, &cGlobalSection);
617: DMGetLabel(dm, "periodic_cut", &cutLabel);
618: N = 0;
620: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
621: VecCreate(PetscObjectComm((PetscObject) dm), &newcoords);
622: PetscSectionGetDof(cSection, vStart, &dof);
623: PetscPrintf(PETSC_COMM_SELF, "DOF: %D\n", dof);
624: embedded = (PetscBool) (L && dof == 2 && !cutLabel);
625: if (cutVertexLabel) {
626: DMLabelGetStratumSize(cutVertexLabel, 1, &v);
627: N += dof*v;
628: }
629: for (v = vStart; v < vEnd; ++v) {
630: PetscSectionGetDof(cGlobalSection, v, &dof);
631: if (dof < 0) continue;
632: if (embedded) N += dof+1;
633: else N += dof;
634: }
635: if (embedded) {VecSetBlockSize(newcoords, bs+1);}
636: else {VecSetBlockSize(newcoords, bs);}
637: VecSetSizes(newcoords, N, PETSC_DETERMINE);
638: VecSetType(newcoords, VECSTANDARD);
639: VecGetArray(coordinatesLocal, &coords);
640: VecGetArray(newcoords, &ncoords);
641: coordSize = 0;
642: for (v = vStart; v < vEnd; ++v) {
643: PetscSectionGetDof(cGlobalSection, v, &dof);
644: PetscSectionGetOffset(cSection, v, &off);
645: if (dof < 0) continue;
646: if (embedded) {
647: if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
648: PetscReal theta, phi, r, R;
649: /* XY-periodic */
650: /* Suppose its an y-z circle, then
651: \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
652: and the circle in that plane is
653: \hat r cos(phi) + \hat x sin(phi) */
654: theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
655: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
656: r = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
657: R = L[1]/(2.0*PETSC_PI);
658: ncoords[coordSize++] = PetscSinReal(phi) * r;
659: ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
660: ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
661: } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
662: /* X-periodic */
663: ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
664: ncoords[coordSize++] = coords[off+1];
665: ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
666: } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
667: /* Y-periodic */
668: ncoords[coordSize++] = coords[off+0];
669: ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
670: ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
671: } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
672: PetscReal phi, r, R;
673: /* Mobius strip */
674: /* Suppose its an x-z circle, then
675: \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
676: and in that plane we rotate by pi as we go around the circle
677: \hat r cos(phi/2) + \hat y sin(phi/2) */
678: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
679: R = L[0];
680: r = PetscRealPart(coords[off+1]) - L[1]/2.0;
681: ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
682: ncoords[coordSize++] = PetscSinReal(phi/2.0) * r;
683: ncoords[coordSize++] = PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
684: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
685: } else {
686: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
687: }
688: }
689: if (cutVertexLabel) {
690: IS vertices;
691: const PetscInt *verts;
692: PetscInt n;
694: DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
695: if (vertices) {
696: ISGetIndices(vertices, &verts);
697: ISGetLocalSize(vertices, &n);
698: for (v = 0; v < n; ++v) {
699: PetscSectionGetDof(cSection, verts[v], &dof);
700: PetscSectionGetOffset(cSection, verts[v], &off);
701: for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
702: }
703: ISRestoreIndices(vertices, &verts);
704: ISDestroy(&vertices);
705: }
706: }
707: if (coordSize != N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %D != %D", coordSize, N);
708: DMLabelDestroy(&cutVertexLabel);
709: VecRestoreArray(coordinatesLocal, &coords);
710: VecRestoreArray(newcoords, &ncoords);
711: PetscObjectSetName((PetscObject) newcoords, "vertices");
712: VecScale(newcoords, lengthScale);
713: PetscViewerHDF5PushGroup(viewer, "/viz");
714: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
715: PetscStackCallHDF5(H5Gclose,(groupId));
716: PetscViewerHDF5PopGroup(viewer);
717: PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
718: VecView(newcoords, viewer);
719: PetscViewerHDF5PopGroup(viewer);
720: VecDestroy(&newcoords);
721: return(0);
722: }
724: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
725: {
726: const PetscInt *gpoint;
727: PetscInt numLabels, l;
728: hid_t fileId, groupId;
729: PetscErrorCode ierr;
732: ISGetIndices(globalPointNumbers, &gpoint);
733: PetscViewerHDF5PushGroup(viewer, "/labels");
734: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
735: if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
736: PetscViewerHDF5PopGroup(viewer);
737: DMGetNumLabels(dm, &numLabels);
738: for (l = 0; l < numLabels; ++l) {
739: DMLabel label;
740: const char *name;
741: IS valueIS, pvalueIS, globalValueIS;
742: const PetscInt *values;
743: PetscInt numValues, v;
744: PetscBool isDepth, output;
745: char group[PETSC_MAX_PATH_LEN];
747: DMGetLabelName(dm, l, &name);
748: DMGetLabelOutput(dm, name, &output);
749: PetscStrncmp(name, "depth", 10, &isDepth);
750: if (isDepth || !output) continue;
751: DMGetLabel(dm, name, &label);
752: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
753: PetscViewerHDF5PushGroup(viewer, group);
754: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
755: if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
756: PetscViewerHDF5PopGroup(viewer);
757: DMLabelGetValueIS(label, &valueIS);
758: /* Must copy to a new IS on the global comm */
759: ISGetLocalSize(valueIS, &numValues);
760: ISGetIndices(valueIS, &values);
761: ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
762: ISRestoreIndices(valueIS, &values);
763: ISAllGather(pvalueIS, &globalValueIS);
764: ISDestroy(&pvalueIS);
765: ISSortRemoveDups(globalValueIS);
766: ISGetLocalSize(globalValueIS, &numValues);
767: ISGetIndices(globalValueIS, &values);
768: for (v = 0; v < numValues; ++v) {
769: IS stratumIS, globalStratumIS;
770: const PetscInt *spoints = NULL;
771: PetscInt *gspoints, n = 0, gn, p;
772: const char *iname = "indices";
774: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%D", name, values[v]);
775: DMLabelGetStratumIS(label, values[v], &stratumIS);
777: if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
778: if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
779: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
780: PetscMalloc1(gn,&gspoints);
781: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
782: if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
783: ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
784: if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
785: PetscObjectSetName((PetscObject) globalStratumIS, iname);
787: PetscViewerHDF5PushGroup(viewer, group);
788: ISView(globalStratumIS, viewer);
789: PetscViewerHDF5PopGroup(viewer);
790: ISDestroy(&globalStratumIS);
791: ISDestroy(&stratumIS);
792: }
793: ISRestoreIndices(globalValueIS, &values);
794: ISDestroy(&globalValueIS);
795: ISDestroy(&valueIS);
796: }
797: ISRestoreIndices(globalPointNumbers, &gpoint);
798: return(0);
799: }
801: /* We only write cells and vertices. Does this screw up parallel reading? */
802: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
803: {
804: IS globalPointNumbers;
805: PetscViewerFormat format;
806: PetscBool viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
807: PetscErrorCode ierr;
810: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
811: DMPlexCoordinatesView_HDF5_Internal(dm, viewer);
812: DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer);
814: PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT);
816: PetscViewerGetFormat(viewer, &format);
817: switch (format) {
818: case PETSC_VIEWER_HDF5_VIZ:
819: viz_geom = PETSC_TRUE;
820: xdmf_topo = PETSC_TRUE;
821: break;
822: case PETSC_VIEWER_HDF5_XDMF:
823: xdmf_topo = PETSC_TRUE;
824: break;
825: case PETSC_VIEWER_HDF5_PETSC:
826: petsc_topo = PETSC_TRUE;
827: break;
828: case PETSC_VIEWER_DEFAULT:
829: case PETSC_VIEWER_NATIVE:
830: viz_geom = PETSC_TRUE;
831: xdmf_topo = PETSC_TRUE;
832: petsc_topo = PETSC_TRUE;
833: break;
834: default:
835: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
836: }
838: if (viz_geom) {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
839: if (xdmf_topo) {DMPlexWriteTopology_Vertices_HDF5_Static(dm, globalPointNumbers, viewer);}
840: if (petsc_topo) {DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer);}
842: ISDestroy(&globalPointNumbers);
843: return(0);
844: }
846: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
847: {
848: MPI_Comm comm;
849: const char *topologydm_name;
850: const char *sectiondm_name;
851: PetscSection gsection;
855: PetscObjectGetComm((PetscObject)sectiondm, &comm);
856: PetscObjectGetName((PetscObject)dm, &topologydm_name);
857: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
858: PetscViewerHDF5PushGroup(viewer, "topologies");
859: PetscViewerHDF5PushGroup(viewer, topologydm_name);
860: PetscViewerHDF5PushGroup(viewer, "dms");
861: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
862: DMGetGlobalSection(sectiondm, &gsection);
863: /* Save raw section */
864: PetscSectionView(gsection, viewer);
865: /* Save plex wrapper */
866: {
867: PetscInt pStart, pEnd, p, n;
868: IS globalPointNumbers;
869: const PetscInt *gpoints;
870: IS orderIS;
871: PetscInt *order;
873: PetscSectionGetChart(gsection, &pStart, &pEnd);
874: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
875: ISGetIndices(globalPointNumbers, &gpoints);
876: for (p = pStart, n = 0; p < pEnd; ++p) if (gpoints[p] >= 0) n++;
877: /* "order" is an array of global point numbers.
878: When loading, it is used with topology/order array
879: to match section points with plex topology points. */
880: PetscMalloc1(n, &order);
881: for (p = pStart, n = 0; p < pEnd; ++p) if (gpoints[p] >= 0) order[n++] = gpoints[p];
882: ISRestoreIndices(globalPointNumbers, &gpoints);
883: ISDestroy(&globalPointNumbers);
884: ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS);
885: PetscObjectSetName((PetscObject)orderIS, "order");
886: ISView(orderIS, viewer);
887: ISDestroy(&orderIS);
888: }
889: PetscViewerHDF5PopGroup(viewer);
890: PetscViewerHDF5PopGroup(viewer);
891: PetscViewerHDF5PopGroup(viewer);
892: PetscViewerHDF5PopGroup(viewer);
893: return(0);
894: }
896: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
897: {
898: const char *topologydm_name;
899: const char *sectiondm_name;
900: const char *vec_name;
901: PetscInt bs;
902: PetscErrorCode ierr;
905: /* Check consistency */
906: {
907: PetscSF pointsf, pointsf1;
909: DMGetPointSF(dm, &pointsf);
910: DMGetPointSF(sectiondm, &pointsf1);
911: if (pointsf1 != pointsf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
912: }
913: PetscObjectGetName((PetscObject)dm, &topologydm_name);
914: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
915: PetscObjectGetName((PetscObject)vec, &vec_name);
916: PetscViewerHDF5PushGroup(viewer, "topologies");
917: PetscViewerHDF5PushGroup(viewer, topologydm_name);
918: PetscViewerHDF5PushGroup(viewer, "dms");
919: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
920: PetscViewerHDF5PushGroup(viewer, "vecs");
921: PetscViewerHDF5PushGroup(viewer, vec_name);
922: VecGetBlockSize(vec, &bs);
923: PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *) &bs);
924: VecSetBlockSize(vec, 1);
925: /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but, */
926: /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view */
927: /* is set to VecView_Plex, which would save vec in a predefined location. */
928: /* To save vec in where we want, we create a new Vec (temp) with */
929: /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
930: {
931: Vec temp;
932: const PetscScalar *array;
933: PetscLayout map;
935: VecCreate(PetscObjectComm((PetscObject)vec), &temp);
936: PetscObjectSetName((PetscObject)temp, vec_name);
937: VecGetLayout(vec, &map);
938: VecSetLayout(temp, map);
939: VecSetUp(temp);
940: VecGetArrayRead(vec, &array);
941: VecPlaceArray(temp, array);
942: VecView(temp, viewer);
943: VecResetArray(temp);
944: VecRestoreArrayRead(vec, &array);
945: VecDestroy(&temp);
946: }
947: VecSetBlockSize(vec, bs);
948: PetscViewerHDF5PopGroup(viewer);
949: PetscViewerHDF5PopGroup(viewer);
950: PetscViewerHDF5PopGroup(viewer);
951: PetscViewerHDF5PopGroup(viewer);
952: PetscViewerHDF5PopGroup(viewer);
953: PetscViewerHDF5PopGroup(viewer);
954: return(0);
955: }
957: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
958: {
959: MPI_Comm comm;
960: const char *topologydm_name;
961: const char *sectiondm_name;
962: const char *vec_name;
963: PetscSection section;
964: PetscBool includesConstraints;
965: Vec gvec;
966: PetscInt m, bs;
967: PetscErrorCode ierr;
970: PetscObjectGetComm((PetscObject)dm, &comm);
971: /* Check consistency */
972: {
973: PetscSF pointsf, pointsf1;
975: DMGetPointSF(dm, &pointsf);
976: DMGetPointSF(sectiondm, &pointsf1);
977: if (pointsf1 != pointsf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
978: }
979: PetscObjectGetName((PetscObject)dm, &topologydm_name);
980: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
981: PetscObjectGetName((PetscObject)vec, &vec_name);
982: PetscViewerHDF5PushGroup(viewer, "topologies");
983: PetscViewerHDF5PushGroup(viewer, topologydm_name);
984: PetscViewerHDF5PushGroup(viewer, "dms");
985: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
986: PetscViewerHDF5PushGroup(viewer, "vecs");
987: PetscViewerHDF5PushGroup(viewer, vec_name);
988: VecGetBlockSize(vec, &bs);
989: PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *) &bs);
990: VecCreate(comm, &gvec);
991: PetscObjectSetName((PetscObject)gvec, vec_name);
992: DMGetGlobalSection(sectiondm, §ion);
993: PetscSectionGetIncludesConstraints(section, &includesConstraints);
994: if (includesConstraints) {PetscSectionGetStorageSize(section, &m);}
995: else {PetscSectionGetConstrainedStorageSize(section, &m);}
996: VecSetSizes(gvec, m, PETSC_DECIDE);
997: VecSetUp(gvec);
998: DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec);
999: DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec);
1000: VecView(gvec, viewer);
1001: VecDestroy(&gvec);
1002: PetscViewerHDF5PopGroup(viewer);
1003: PetscViewerHDF5PopGroup(viewer);
1004: PetscViewerHDF5PopGroup(viewer);
1005: PetscViewerHDF5PopGroup(viewer);
1006: PetscViewerHDF5PopGroup(viewer);
1007: PetscViewerHDF5PopGroup(viewer);
1008: return(0);
1009: }
1011: typedef struct {
1012: PetscMPIInt rank;
1013: DM dm;
1014: PetscViewer viewer;
1015: DMLabel label;
1016: } LabelCtx;
1018: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
1019: {
1020: PetscViewer viewer = ((LabelCtx *) op_data)->viewer;
1021: DMLabel label = ((LabelCtx *) op_data)->label;
1022: IS stratumIS;
1023: const PetscInt *ind;
1024: PetscInt value, N, i;
1025: const char *lname;
1026: char group[PETSC_MAX_PATH_LEN];
1027: PetscErrorCode ierr;
1029: PetscOptionsStringToInt(name, &value);
1030: ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
1031: PetscObjectSetName((PetscObject) stratumIS, "indices");
1032: PetscObjectGetName((PetscObject) label, &lname);
1033: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
1034: PetscViewerHDF5PushGroup(viewer, group);
1035: {
1036: /* Force serial load */
1037: PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
1038: PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
1039: PetscLayoutSetSize(stratumIS->map, N);
1040: }
1041: ISLoad(stratumIS, viewer);
1042: PetscViewerHDF5PopGroup(viewer);
1043: ISGetLocalSize(stratumIS, &N);
1044: ISGetIndices(stratumIS, &ind);
1045: for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
1046: ISRestoreIndices(stratumIS, &ind);
1047: ISDestroy(&stratumIS);
1048: return 0;
1049: }
1051: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
1052: {
1053: DM dm = ((LabelCtx *) op_data)->dm;
1054: hsize_t idx = 0;
1056: herr_t err;
1058: DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
1059: DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
1060: PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1061: return err;
1062: }
1064: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1065: {
1066: LabelCtx ctx;
1067: hid_t fileId, groupId;
1068: hsize_t idx = 0;
1069: PetscErrorCode ierr;
1072: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &ctx.rank);
1073: ctx.dm = dm;
1074: ctx.viewer = viewer;
1075: PetscViewerHDF5PushGroup(viewer, "/labels");
1076: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
1077: PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
1078: PetscStackCallHDF5(H5Gclose,(groupId));
1079: PetscViewerHDF5PopGroup(viewer);
1080: return(0);
1081: }
1083: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sf)
1084: {
1085: MPI_Comm comm;
1086: IS orderIS, conesIS, cellsIS, orntsIS;
1087: const PetscInt *order, *cones, *cells, *ornts;
1088: PetscInt *cone, *ornt;
1089: PetscInt dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
1090: PetscMPIInt size, rank;
1091: PetscErrorCode ierr;
1094: PetscObjectGetComm((PetscObject)dm, &comm);
1095: MPI_Comm_size(comm, &size);
1096: MPI_Comm_rank(comm, &rank);
1097: /* Read toplogy */
1098: PetscViewerHDF5PushGroup(viewer, "/topology");
1099: ISCreate(comm, &orderIS);
1100: PetscObjectSetName((PetscObject) orderIS, "order");
1101: ISCreate(comm, &conesIS);
1102: PetscObjectSetName((PetscObject) conesIS, "cones");
1103: ISCreate(comm, &cellsIS);
1104: PetscObjectSetName((PetscObject) cellsIS, "cells");
1105: ISCreate(comm, &orntsIS);
1106: PetscObjectSetName((PetscObject) orntsIS, "orientation");
1107: PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, NULL, &dim);
1108: DMSetDimension(dm, dim);
1109: {
1110: /* Force serial load */
1111: PetscViewerHDF5ReadSizes(viewer, "order", NULL, &Np);
1112: PetscLayoutSetLocalSize(orderIS->map, rank == 0 ? Np : 0);
1113: PetscLayoutSetSize(orderIS->map, Np);
1114: PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &Np);
1115: PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? Np : 0);
1116: PetscLayoutSetSize(conesIS->map, Np);
1117: pEnd = rank == 0 ? Np : 0;
1118: PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
1119: PetscLayoutSetLocalSize(cellsIS->map, rank == 0 ? N : 0);
1120: PetscLayoutSetSize(cellsIS->map, N);
1121: PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
1122: PetscLayoutSetLocalSize(orntsIS->map, rank == 0 ? N : 0);
1123: PetscLayoutSetSize(orntsIS->map, N);
1124: }
1125: ISLoad(orderIS, viewer);
1126: ISLoad(conesIS, viewer);
1127: ISLoad(cellsIS, viewer);
1128: ISLoad(orntsIS, viewer);
1129: PetscViewerHDF5PopGroup(viewer);
1130: /* Create Plex */
1131: DMPlexSetChart(dm, 0, pEnd);
1132: ISGetIndices(orderIS, &order);
1133: ISGetIndices(conesIS, &cones);
1134: for (p = 0; p < pEnd; ++p) {
1135: DMPlexSetConeSize(dm, order[p], cones[p]);
1136: maxConeSize = PetscMax(maxConeSize, cones[p]);
1137: }
1138: DMSetUp(dm);
1139: ISGetIndices(cellsIS, &cells);
1140: ISGetIndices(orntsIS, &ornts);
1141: PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
1142: for (p = 0, q = 0; p < pEnd; ++p) {
1143: for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
1144: DMPlexSetCone(dm, order[p], cone);
1145: DMPlexSetConeOrientation(dm, order[p], ornt);
1146: }
1147: PetscFree2(cone,ornt);
1148: /* Create global section migration SF */
1149: if (sf) {
1150: PetscLayout layout;
1151: PetscInt *globalIndices;
1153: PetscMalloc1(pEnd, &globalIndices);
1154: /* plex point == globalPointNumber in this case */
1155: for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
1156: PetscLayoutCreate(comm, &layout);
1157: PetscLayoutSetSize(layout, Np);
1158: PetscLayoutSetBlockSize(layout, 1);
1159: PetscLayoutSetUp(layout);
1160: PetscSFCreate(comm, sf);
1161: PetscSFSetFromOptions(*sf);
1162: PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices);
1163: PetscLayoutDestroy(&layout);
1164: PetscFree(globalIndices);
1165: }
1166: /* */
1167: ISRestoreIndices(orderIS, &order);
1168: ISRestoreIndices(conesIS, &cones);
1169: ISRestoreIndices(cellsIS, &cells);
1170: ISRestoreIndices(orntsIS, &ornts);
1171: ISDestroy(&orderIS);
1172: ISDestroy(&conesIS);
1173: ISDestroy(&cellsIS);
1174: ISDestroy(&orntsIS);
1175: DMPlexSymmetrize(dm);
1176: DMPlexStratify(dm);
1177: return(0);
1178: }
1180: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1181: {
1182: PetscSection coordSection;
1183: Vec coordinates;
1184: PetscReal lengthScale;
1185: PetscInt spatialDim, N, numVertices, vStart, vEnd, v;
1186: PetscMPIInt rank;
1187: PetscErrorCode ierr;
1190: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1191: /* Read geometry */
1192: PetscViewerHDF5PushGroup(viewer, "/geometry");
1193: VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
1194: PetscObjectSetName((PetscObject) coordinates, "vertices");
1195: {
1196: /* Force serial load */
1197: PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
1198: VecSetSizes(coordinates, !rank ? N : 0, N);
1199: VecSetBlockSize(coordinates, spatialDim);
1200: }
1201: VecLoad(coordinates, viewer);
1202: PetscViewerHDF5PopGroup(viewer);
1203: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1204: VecScale(coordinates, 1.0/lengthScale);
1205: VecGetLocalSize(coordinates, &numVertices);
1206: VecGetBlockSize(coordinates, &spatialDim);
1207: numVertices /= spatialDim;
1208: /* Create coordinates */
1209: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1210: if (numVertices != vEnd - vStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of coordinates loaded %d does not match number of vertices %d", numVertices, vEnd - vStart);
1211: DMGetCoordinateSection(dm, &coordSection);
1212: PetscSectionSetNumFields(coordSection, 1);
1213: PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
1214: PetscSectionSetChart(coordSection, vStart, vEnd);
1215: for (v = vStart; v < vEnd; ++v) {
1216: PetscSectionSetDof(coordSection, v, spatialDim);
1217: PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
1218: }
1219: PetscSectionSetUp(coordSection);
1220: DMSetCoordinates(dm, coordinates);
1221: VecDestroy(&coordinates);
1222: return(0);
1223: }
1225: /* The first version will read everything onto proc 0, letting the user distribute
1226: The next will create a naive partition, and then rebalance after reading
1227: */
1228: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1229: {
1230: PetscErrorCode ierr;
1233: DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL);
1234: DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer);
1235: DMPlexLabelsLoad_HDF5_Internal(dm, viewer);
1236: return(0);
1237: }
1239: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1240: {
1241: MPI_Comm comm;
1242: PetscInt pStart, pEnd, p, m;
1243: PetscInt *goffs, *ilocal;
1244: PetscBool rootIncludeConstraints, leafIncludeConstraints;
1245: PetscErrorCode ierr;
1248: PetscObjectGetComm((PetscObject)leafSection, &comm);
1249: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1250: PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints);
1251: PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints);
1252: if (rootIncludeConstraints && leafIncludeConstraints) {PetscSectionGetStorageSize(leafSection, &m);}
1253: else {PetscSectionGetConstrainedStorageSize(leafSection, &m);}
1254: PetscMalloc1(m, &ilocal);
1255: PetscMalloc1(m, &goffs);
1256: /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
1257: /* for the top-level section (not for each field), so one must have */
1258: /* rootSection->pointMajor == PETSC_TRUE. */
1259: if (!rootSection->pointMajor) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for field major ordering");
1260: /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
1261: if (!leafSection->pointMajor) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for field major ordering");
1262: for (p = pStart, m = 0; p < pEnd; ++p) {
1263: PetscInt dof, cdof, i, j, off, goff;
1264: const PetscInt *cinds;
1266: PetscSectionGetDof(leafSection, p, &dof);
1267: if (dof < 0) continue;
1268: goff = globalOffsets[p-pStart];
1269: PetscSectionGetOffset(leafSection, p, &off);
1270: PetscSectionGetConstraintDof(leafSection, p, &cdof);
1271: PetscSectionGetConstraintIndices(leafSection, p, &cinds);
1272: for (i = 0, j = 0; i < dof; ++i) {
1273: PetscBool constrained = (PetscBool) (j < cdof && i == cinds[j]);
1275: if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {ilocal[m] = off++; goffs[m++] = goff++;}
1276: else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
1277: else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff;
1278: if (constrained) ++j;
1279: }
1280: }
1281: PetscSFCreate(comm, sectionSF);
1282: PetscSFSetFromOptions(*sectionSF);
1283: PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs);
1284: PetscFree(goffs);
1285: return(0);
1286: }
1288: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
1289: {
1290: MPI_Comm comm;
1291: PetscMPIInt size, rank;
1292: const char *topologydm_name;
1293: const char *sectiondm_name;
1294: PetscSection sectionA, sectionB;
1295: PetscInt NX, nX, n, i;
1296: PetscSF sfAB;
1300: PetscObjectGetComm((PetscObject)dm, &comm);
1301: MPI_Comm_size(comm, &size);
1302: MPI_Comm_rank(comm, &rank);
1303: PetscObjectGetName((PetscObject)dm, &topologydm_name);
1304: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1305: PetscViewerHDF5PushGroup(viewer, "topologies");
1306: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1307: PetscViewerHDF5ReadSizes(viewer, "/topology/order", NULL, &NX);
1308: PetscViewerHDF5PushGroup(viewer, "dms");
1309: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1310: /* A: on-disk points */
1311: /* X: list of global point numbers, [0, NX) */
1312: /* B: plex points */
1313: /* Load raw section (sectionA) */
1314: PetscSectionCreate(comm, §ionA);
1315: PetscSectionLoad(sectionA, viewer);
1316: PetscSectionGetChart(sectionA, NULL, &n);
1317: /* Create sfAB: A -> B */
1318: #if defined(PETSC_USE_DEBUG)
1319: {
1320: PetscInt N, N1;
1322: PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1);
1323: MPI_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm);
1324: if (N1 != N) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Mismatching sizes: on-disk order array size (%D) != number of loaded section points (%D)", N1, N);
1325: }
1326: #endif
1327: {
1328: IS orderIS;
1329: const PetscInt *gpoints;
1330: PetscSF sfXA, sfAX;
1331: PetscLayout layout;
1332: PetscSFNode *owners, *buffer;
1333: PetscInt nleaves;
1334: PetscInt *ilocal;
1335: PetscSFNode *iremote;
1337: /* Create sfAX: A -> X */
1338: ISCreate(comm, &orderIS);
1339: PetscObjectSetName((PetscObject)orderIS, "order");
1340: PetscLayoutSetLocalSize(orderIS->map, n);
1341: ISLoad(orderIS, viewer);
1342: PetscLayoutCreate(comm, &layout);
1343: PetscLayoutSetSize(layout, NX);
1344: PetscLayoutSetBlockSize(layout, 1);
1345: PetscLayoutSetUp(layout);
1346: PetscSFCreate(comm, &sfXA);
1347: ISGetIndices(orderIS, &gpoints);
1348: PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints);
1349: ISRestoreIndices(orderIS, &gpoints);
1350: ISDestroy(&orderIS);
1351: PetscLayoutDestroy(&layout);
1352: PetscSFGetGraph(sfXA, &nX, NULL, NULL, NULL);
1353: PetscMalloc1(n, &owners);
1354: PetscMalloc1(nX, &buffer);
1355: for (i = 0; i < n; ++i) {owners[i].rank = rank; owners[i].index = i;}
1356: for (i = 0; i < nX; ++i) {buffer[i].rank = -1; buffer[i].index = -1;}
1357: PetscSFReduceBegin(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
1358: PetscSFReduceEnd(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
1359: PetscSFDestroy(&sfXA);
1360: PetscFree(owners);
1361: for (i = 0, nleaves = 0; i < nX; ++i) if (buffer[i].rank >= 0) nleaves++;
1362: PetscMalloc1(nleaves, &ilocal);
1363: PetscMalloc1(nleaves, &iremote);
1364: for (i = 0, nleaves = 0; i < nX; ++i) {
1365: if (buffer[i].rank >= 0) {
1366: ilocal[nleaves] = i;
1367: iremote[nleaves].rank = buffer[i].rank;
1368: iremote[nleaves].index = buffer[i].index;
1369: nleaves++;
1370: }
1371: }
1372: PetscSFCreate(comm, &sfAX);
1373: PetscSFSetFromOptions(sfAX);
1374: PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1375: /* Fix PetscSFCompose() and replace the code-block below with: */
1376: /* PetscSFCompose(sfAX, sfXB, &sfAB); */
1377: /* which currently causes segmentation fault due to sparse map. */
1378: {
1379: PetscInt npoints;
1380: PetscInt mleaves;
1381: PetscInt *jlocal;
1382: PetscSFNode *jremote;
1384: PetscSFGetGraph(sfXB, NULL, &npoints, NULL, NULL);
1385: PetscMalloc1(npoints, &owners);
1386: for (i = 0; i < npoints; ++i) {owners[i].rank = -1; owners[i].index = -1;}
1387: PetscSFBcastBegin(sfXB, MPIU_2INT, buffer, owners, MPI_REPLACE);
1388: PetscSFBcastEnd(sfXB, MPIU_2INT, buffer, owners, MPI_REPLACE);
1389: for (i = 0, mleaves = 0; i < npoints; ++i) if (owners[i].rank >= 0) mleaves++;
1390: PetscMalloc1(mleaves, &jlocal);
1391: PetscMalloc1(mleaves, &jremote);
1392: for (i = 0, mleaves = 0; i < npoints; ++i) {
1393: if (owners[i].rank >= 0) {
1394: jlocal[mleaves] = i;
1395: jremote[mleaves].rank = owners[i].rank;
1396: jremote[mleaves].index = owners[i].index;
1397: mleaves++;
1398: }
1399: }
1400: PetscSFCreate(comm, &sfAB);
1401: PetscSFSetFromOptions(sfAB);
1402: PetscSFSetGraph(sfAB, n, mleaves, jlocal, PETSC_OWN_POINTER, jremote, PETSC_OWN_POINTER);
1403: PetscFree(owners);
1404: }
1405: PetscFree(buffer);
1406: PetscSFDestroy(&sfAX);
1407: }
1408: PetscViewerHDF5PopGroup(viewer);
1409: PetscViewerHDF5PopGroup(viewer);
1410: PetscViewerHDF5PopGroup(viewer);
1411: PetscViewerHDF5PopGroup(viewer);
1412: /* Create plex section (sectionB) */
1413: DMGetLocalSection(sectiondm, §ionB);
1414: if (lsf || gsf) {
1415: PetscLayout layout;
1416: PetscInt M, m;
1417: PetscInt *offsetsA;
1418: PetscBool includesConstraintsA;
1420: PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB);
1421: PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA);
1422: if (includesConstraintsA) {PetscSectionGetStorageSize(sectionA, &m);}
1423: else {PetscSectionGetConstrainedStorageSize(sectionA, &m);}
1424: MPI_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm);
1425: PetscLayoutCreate(comm, &layout);
1426: PetscLayoutSetSize(layout, M);
1427: PetscLayoutSetUp(layout);
1428: if (lsf) {
1429: PetscSF lsfABdata;
1431: DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata);
1432: *lsf = lsfABdata;
1433: }
1434: if (gsf) {
1435: PetscSection gsectionB, gsectionB1;
1436: PetscBool includesConstraintsB;
1437: PetscSF gsfABdata, pointsf;
1439: DMGetGlobalSection(sectiondm, &gsectionB1);
1440: PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB);
1441: DMGetPointSF(sectiondm, &pointsf);
1442: PetscSectionCreateGlobalSection(sectionB, pointsf, includesConstraintsB, PETSC_TRUE, &gsectionB);
1443: DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata);
1444: PetscSectionDestroy(&gsectionB);
1445: *gsf = gsfABdata;
1446: }
1447: PetscLayoutDestroy(&layout);
1448: PetscFree(offsetsA);
1449: } else {
1450: PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB);
1451: }
1452: PetscSFDestroy(&sfAB);
1453: PetscSectionDestroy(§ionA);
1454: return(0);
1455: }
1457: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
1458: {
1459: MPI_Comm comm;
1460: const char *topologydm_name;
1461: const char *sectiondm_name;
1462: const char *vec_name;
1463: Vec vecA;
1464: PetscInt mA, m, bs;
1465: const PetscInt *ilocal;
1466: const PetscScalar *src;
1467: PetscScalar *dest;
1468: PetscErrorCode ierr;
1471: PetscObjectGetComm((PetscObject)dm, &comm);
1472: PetscObjectGetName((PetscObject)dm, &topologydm_name);
1473: PetscObjectGetName((PetscObject)sectiondm, §iondm_name);
1474: PetscObjectGetName((PetscObject)vec, &vec_name);
1475: PetscViewerHDF5PushGroup(viewer, "topologies");
1476: PetscViewerHDF5PushGroup(viewer, topologydm_name);
1477: PetscViewerHDF5PushGroup(viewer, "dms");
1478: PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1479: PetscViewerHDF5PushGroup(viewer, "vecs");
1480: PetscViewerHDF5PushGroup(viewer, vec_name);
1481: VecCreate(comm, &vecA);
1482: PetscObjectSetName((PetscObject)vecA, vec_name);
1483: PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL);
1484: /* Check consistency */
1485: {
1486: PetscSF pointsf, pointsf1;
1487: PetscInt m1, i, j;
1489: DMGetPointSF(dm, &pointsf);
1490: DMGetPointSF(sectiondm, &pointsf1);
1491: if (pointsf1 != pointsf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
1492: #if defined(PETSC_USE_DEBUG)
1493: {
1494: PetscInt MA, MA1;
1496: MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm);
1497: PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1);
1498: if (MA1 != MA) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total SF root size (%D) != On-disk vector data size (%D)", MA, MA1);
1499: }
1500: #endif
1501: VecGetLocalSize(vec, &m1);
1502: if (m1 < m) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Target vector size (%D) < SF leaf size (%D)", m1, m);
1503: for (i = 0; i < m; ++i) {
1504: j = ilocal ? ilocal[i] : i;
1505: if (j < 0 || j >= m1) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Leaf's %D-th index, %D, not in [%D, %D)", i, j, 0, m1);
1506: }
1507: }
1508: VecSetSizes(vecA, mA, PETSC_DECIDE);
1509: VecLoad(vecA, viewer);
1510: VecGetArrayRead(vecA, &src);
1511: VecGetArray(vec, &dest);
1512: PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
1513: PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
1514: VecRestoreArray(vec, &dest);
1515: VecRestoreArrayRead(vecA, &src);
1516: VecDestroy(&vecA);
1517: PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *) &bs);
1518: VecSetBlockSize(vec, bs);
1519: PetscViewerHDF5PopGroup(viewer);
1520: PetscViewerHDF5PopGroup(viewer);
1521: PetscViewerHDF5PopGroup(viewer);
1522: PetscViewerHDF5PopGroup(viewer);
1523: PetscViewerHDF5PopGroup(viewer);
1524: PetscViewerHDF5PopGroup(viewer);
1525: return(0);
1526: }
1527: #endif