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) {
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: PetscViewerHDF5SetTimestep(viewer, seqnum);
34: VecView(stamp, viewer);
35: PetscViewerHDF5PopGroup(viewer);
36: VecDestroy(&stamp);
37: return(0);
38: }
40: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
41: {
42: Vec stamp;
43: PetscMPIInt rank;
47: if (seqnum < 0) return(0);
48: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
49: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
50: VecSetBlockSize(stamp, 1);
51: PetscObjectSetName((PetscObject) stamp, seqname);
52: PetscViewerHDF5PushGroup(viewer, "/");
53: PetscViewerHDF5SetTimestep(viewer, seqnum);
54: VecLoad(stamp, viewer);
55: PetscViewerHDF5PopGroup(viewer);
56: if (!rank) {
57: const PetscScalar *a;
58: PetscReal timeScale;
59: PetscBool istime;
61: VecGetArrayRead(stamp, &a);
62: *value = a[0];
63: VecRestoreArrayRead(stamp, &a);
64: PetscStrncmp(seqname, "time", 5, &istime);
65: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); *value /= timeScale;}
66: }
67: VecDestroy(&stamp);
68: return(0);
69: }
71: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
72: {
73: IS cutcells = NULL;
74: const PetscInt *cutc;
75: PetscInt cellHeight, vStart, vEnd, cStart, cEnd, c;
76: PetscErrorCode ierr;
79: if (!cutLabel) return(0);
80: DMPlexGetVTKCellHeight(dm, &cellHeight);
81: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
82: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
83: /* Label vertices that should be duplicated */
84: DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel);
85: DMLabelGetStratumIS(cutLabel, 2, &cutcells);
86: if (cutcells) {
87: PetscInt n;
89: ISGetIndices(cutcells, &cutc);
90: ISGetLocalSize(cutcells, &n);
91: for (c = 0; c < n; ++c) {
92: if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
93: PetscInt *closure = NULL;
94: PetscInt closureSize, cl, value;
96: DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
97: for (cl = 0; cl < closureSize*2; cl += 2) {
98: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
99: DMLabelGetValue(cutLabel, closure[cl], &value);
100: if (value == 1) {
101: DMLabelSetValue(*cutVertexLabel, closure[cl], 1);
102: }
103: }
104: }
105: DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
106: }
107: }
108: ISRestoreIndices(cutcells, &cutc);
109: }
110: ISDestroy(&cutcells);
111: return(0);
112: }
114: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
115: {
116: DM dm;
117: DM dmBC;
118: PetscSection section, sectionGlobal;
119: Vec gv;
120: const char *name;
121: PetscViewerVTKFieldType ft;
122: PetscViewerFormat format;
123: PetscInt seqnum;
124: PetscReal seqval;
125: PetscBool isseq;
126: PetscErrorCode ierr;
129: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
130: VecGetDM(v, &dm);
131: DMGetLocalSection(dm, §ion);
132: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
133: PetscViewerHDF5SetTimestep(viewer, seqnum);
134: DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
135: PetscViewerGetFormat(viewer, &format);
136: DMGetOutputDM(dm, &dmBC);
137: DMGetGlobalSection(dmBC, §ionGlobal);
138: DMGetGlobalVector(dmBC, &gv);
139: PetscObjectGetName((PetscObject) v, &name);
140: PetscObjectSetName((PetscObject) gv, name);
141: DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
142: DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
143: PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
144: if (isseq) {VecView_Seq(gv, viewer);}
145: else {VecView_MPI(gv, viewer);}
146: if (format == PETSC_VIEWER_HDF5_VIZ) {
147: /* Output visualization representation */
148: PetscInt numFields, f;
149: DMLabel cutLabel, cutVertexLabel = NULL;
151: PetscSectionGetNumFields(section, &numFields);
152: DMGetLabel(dm, "periodic_cut", &cutLabel);
153: for (f = 0; f < numFields; ++f) {
154: Vec subv;
155: IS is;
156: const char *fname, *fgroup;
157: char subname[PETSC_MAX_PATH_LEN];
158: PetscInt pStart, pEnd;
160: DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
161: fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
162: PetscSectionGetFieldName(section, f, &fname);
163: if (!fname) continue;
164: PetscViewerHDF5PushGroup(viewer, fgroup);
165: if (cutLabel) {
166: const PetscScalar *ga;
167: PetscScalar *suba;
168: PetscInt Nc, gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;
170: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
171: PetscSectionGetFieldComponents(section, f, &Nc);
172: for (p = pStart; p < pEnd; ++p) {
173: PetscInt gdof, fdof = 0, val;
175: PetscSectionGetDof(sectionGlobal, p, &gdof);
176: if (gdof > 0) {PetscSectionGetFieldDof(section, p, f, &fdof);}
177: subSize += fdof;
178: DMLabelGetValue(cutVertexLabel, p, &val);
179: if (val == 1) extSize += fdof;
180: }
181: VecCreate(PetscObjectComm((PetscObject) gv), &subv);
182: VecSetSizes(subv, subSize+extSize, PETSC_DETERMINE);
183: VecSetBlockSize(subv, Nc);
184: VecSetType(subv, VECSTANDARD);
185: VecGetOwnershipRange(gv, &gstart, NULL);
186: VecGetArrayRead(gv, &ga);
187: VecGetArray(subv, &suba);
188: for (p = pStart; p < pEnd; ++p) {
189: PetscInt gdof, goff, val;
191: PetscSectionGetDof(sectionGlobal, p, &gdof);
192: if (gdof > 0) {
193: PetscInt fdof, fc, f2, poff = 0;
195: PetscSectionGetOffset(sectionGlobal, p, &goff);
196: /* Can get rid of this loop by storing field information in the global section */
197: for (f2 = 0; f2 < f; ++f2) {
198: PetscSectionGetFieldDof(section, p, f2, &fdof);
199: poff += fdof;
200: }
201: PetscSectionGetFieldDof(section, p, f, &fdof);
202: for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff+poff+fc - gstart];
203: DMLabelGetValue(cutVertexLabel, p, &val);
204: if (val == 1) {
205: for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize+newOff] = ga[goff+poff+fc - gstart];
206: }
207: }
208: }
209: VecRestoreArrayRead(gv, &ga);
210: VecRestoreArray(subv, &suba);
211: DMLabelDestroy(&cutVertexLabel);
212: } else {
213: PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
214: }
215: PetscStrncpy(subname, name,sizeof(subname));
216: PetscStrlcat(subname, "_",sizeof(subname));
217: PetscStrlcat(subname, fname,sizeof(subname));
218: PetscObjectSetName((PetscObject) subv, subname);
219: if (isseq) {VecView_Seq(subv, viewer);}
220: else {VecView_MPI(subv, viewer);}
221: if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
222: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "vector");
223: } else {
224: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "scalar");
225: }
226: if (cutLabel) {VecDestroy(&subv);}
227: else {PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);}
228: PetscViewerHDF5PopGroup(viewer);
229: }
230: }
231: DMRestoreGlobalVector(dmBC, &gv);
232: return(0);
233: }
235: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
236: {
237: DM dm;
238: Vec locv;
239: PetscObject isZero;
240: const char *name;
241: PetscReal time;
245: VecGetDM(v, &dm);
246: DMGetLocalVector(dm, &locv);
247: PetscObjectGetName((PetscObject) v, &name);
248: PetscObjectSetName((PetscObject) locv, name);
249: PetscObjectQuery((PetscObject) v, "__Vec_bc_zero__", &isZero);
250: PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", isZero);
251: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
252: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
253: DMGetOutputSequenceNumber(dm, NULL, &time);
254: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
255: PetscViewerHDF5PushGroup(viewer, "/fields");
256: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
257: VecView_Plex_Local_HDF5_Internal(locv, viewer);
258: PetscViewerPopFormat(viewer);
259: PetscViewerHDF5PopGroup(viewer);
260: PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", NULL);
261: DMRestoreLocalVector(dm, &locv);
262: return(0);
263: }
265: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
266: {
267: PetscBool isseq;
271: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
272: PetscViewerHDF5PushGroup(viewer, "/fields");
273: if (isseq) {VecView_Seq(v, viewer);}
274: else {VecView_MPI(v, viewer);}
275: PetscViewerHDF5PopGroup(viewer);
276: return(0);
277: }
279: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
280: {
281: DM dm;
282: Vec locv;
283: const char *name;
284: PetscInt seqnum;
288: VecGetDM(v, &dm);
289: DMGetLocalVector(dm, &locv);
290: PetscObjectGetName((PetscObject) v, &name);
291: PetscObjectSetName((PetscObject) locv, name);
292: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
293: PetscViewerHDF5SetTimestep(viewer, seqnum);
294: PetscViewerHDF5PushGroup(viewer, "/fields");
295: VecLoad_Plex_Local(locv, viewer);
296: PetscViewerHDF5PopGroup(viewer);
297: DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
298: DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
299: DMRestoreLocalVector(dm, &locv);
300: return(0);
301: }
303: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
304: {
305: DM dm;
306: PetscInt seqnum;
310: VecGetDM(v, &dm);
311: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
312: PetscViewerHDF5SetTimestep(viewer, seqnum);
313: PetscViewerHDF5PushGroup(viewer, "/fields");
314: VecLoad_Default(v, viewer);
315: PetscViewerHDF5PopGroup(viewer);
316: return(0);
317: }
319: static PetscErrorCode DMPlexWriteTopology_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
320: {
321: IS orderIS, conesIS, cellsIS, orntsIS;
322: const PetscInt *gpoint;
323: PetscInt *order, *sizes, *cones, *ornts;
324: PetscInt dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
325: PetscErrorCode ierr;
328: ISGetIndices(globalPointNumbers, &gpoint);
329: DMGetDimension(dm, &dim);
330: DMPlexGetChart(dm, &pStart, &pEnd);
331: for (p = pStart; p < pEnd; ++p) {
332: if (gpoint[p] >= 0) {
333: PetscInt coneSize;
335: DMPlexGetConeSize(dm, p, &coneSize);
336: conesSize += 1;
337: cellsSize += coneSize;
338: }
339: }
340: PetscMalloc1(conesSize, &order);
341: PetscMalloc1(conesSize, &sizes);
342: PetscMalloc1(cellsSize, &cones);
343: PetscMalloc1(cellsSize, &ornts);
344: for (p = pStart; p < pEnd; ++p) {
345: if (gpoint[p] >= 0) {
346: const PetscInt *cone, *ornt;
347: PetscInt coneSize, cp;
349: DMPlexGetConeSize(dm, p, &coneSize);
350: DMPlexGetCone(dm, p, &cone);
351: DMPlexGetConeOrientation(dm, p, &ornt);
352: order[s] = gpoint[p];
353: sizes[s++] = coneSize;
354: for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
355: }
356: }
357: if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %D != %D", s, conesSize);
358: if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %D != %D", c, cellsSize);
359: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
360: PetscObjectSetName((PetscObject) orderIS, "order");
361: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
362: PetscObjectSetName((PetscObject) conesIS, "cones");
363: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
364: PetscObjectSetName((PetscObject) cellsIS, "cells");
365: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
366: PetscObjectSetName((PetscObject) orntsIS, "orientation");
367: PetscViewerHDF5PushGroup(viewer, "/topology");
368: ISView(orderIS, viewer);
369: ISView(conesIS, viewer);
370: ISView(cellsIS, viewer);
371: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, (void *) &dim);
372: ISView(orntsIS, viewer);
373: PetscViewerHDF5PopGroup(viewer);
374: ISDestroy(&orderIS);
375: ISDestroy(&conesIS);
376: ISDestroy(&cellsIS);
377: ISDestroy(&orntsIS);
378: ISRestoreIndices(globalPointNumbers, &gpoint);
379: return(0);
380: }
382: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
383: {
384: PetscSF sfPoint;
385: DMLabel cutLabel, cutVertexLabel = NULL;
386: IS globalVertexNumbers, cutvertices = NULL;
387: const PetscInt *gcell, *gvertex, *cutverts = NULL;
388: PetscInt *vertices;
389: PetscInt conesSize = 0;
390: PetscInt dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
391: PetscErrorCode ierr;
394: *numCorners = 0;
395: DMGetDimension(dm, &dim);
396: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
397: ISGetIndices(globalCellNumbers, &gcell);
399: for (cell = cStart; cell < cEnd; ++cell) {
400: PetscInt *closure = NULL;
401: PetscInt closureSize, v, Nc = 0;
403: if (gcell[cell] < 0) continue;
404: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
405: for (v = 0; v < closureSize*2; v += 2) {
406: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
407: }
408: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
409: conesSize += Nc;
410: if (!numCornersLocal) numCornersLocal = Nc;
411: else if (numCornersLocal != Nc) numCornersLocal = 1;
412: }
413: MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
414: if (numCornersLocal && (numCornersLocal != *numCorners || *numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
415: /* Handle periodic cuts by identifying vertices which should be duplicated */
416: DMGetLabel(dm, "periodic_cut", &cutLabel);
417: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
418: if (cutVertexLabel) {DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);}
419: if (cutvertices) {
420: ISGetIndices(cutvertices, &cutverts);
421: ISGetLocalSize(cutvertices, &vExtra);
422: }
423: DMGetPointSF(dm, &sfPoint);
424: if (cutLabel) {
425: const PetscInt *ilocal;
426: const PetscSFNode *iremote;
427: PetscInt nroots, nleaves;
429: PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
430: if (nleaves < 0) {
431: PetscObjectReference((PetscObject) sfPoint);
432: } else {
433: PetscSFCreate(PetscObjectComm((PetscObject) sfPoint), &sfPoint);
434: PetscSFSetGraph(sfPoint, nroots+vExtra, nleaves, ilocal, PETSC_USE_POINTER, iremote, PETSC_USE_POINTER);
435: }
436: } else {
437: PetscObjectReference((PetscObject) sfPoint);
438: }
439: /* Number all vertices */
440: DMPlexCreateNumbering_Plex(dm, vStart, vEnd+vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
441: PetscSFDestroy(&sfPoint);
442: /* Create cones */
443: ISGetIndices(globalVertexNumbers, &gvertex);
444: PetscMalloc1(conesSize, &vertices);
445: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
446: PetscInt *closure = NULL;
447: PetscInt closureSize, Nc = 0, p, value = -1;
448: PetscBool replace;
450: if (gcell[cell] < 0) continue;
451: if (cutLabel) {DMLabelGetValue(cutLabel, cell, &value);}
452: replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
453: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
454: for (p = 0; p < closureSize*2; p += 2) {
455: if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
456: closure[Nc++] = closure[p];
457: }
458: }
459: DMPlexReorderCell(dm, cell, closure);
460: for (p = 0; p < Nc; ++p) {
461: PetscInt nv, gv = gvertex[closure[p] - vStart];
463: if (replace) {
464: PetscFindInt(closure[p], vExtra, cutverts, &nv);
465: if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
466: }
467: vertices[v++] = gv < 0 ? -(gv+1) : gv;
468: }
469: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
470: }
471: ISRestoreIndices(globalVertexNumbers, &gvertex);
472: ISDestroy(&globalVertexNumbers);
473: ISRestoreIndices(globalCellNumbers, &gcell);
474: if (cutvertices) {ISRestoreIndices(cutvertices, &cutverts);}
475: ISDestroy(&cutvertices);
476: DMLabelDestroy(&cutVertexLabel);
477: if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %D != %D", v, conesSize);
478: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS);
479: PetscLayoutSetBlockSize((*cellIS)->map, *numCorners);
480: PetscObjectSetName((PetscObject) *cellIS, "cells");
481: return(0);
482: }
484: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, IS globalCellNumbers, PetscViewer viewer)
485: {
486: DM cdm;
487: DMLabel depthLabel, ctLabel;
488: IS cellIS;
489: PetscInt dim, depth, cellHeight, c;
490: hid_t fileId, groupId;
491: PetscErrorCode ierr;
494: PetscViewerHDF5PushGroup(viewer, "/viz");
495: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
496: PetscStackCallHDF5(H5Gclose,(groupId));
498: PetscViewerHDF5PopGroup(viewer);
499: DMGetDimension(dm, &dim);
500: DMPlexGetDepth(dm, &depth);
501: DMGetCoordinateDM(dm, &cdm);
502: DMPlexGetVTKCellHeight(dm, &cellHeight);
503: DMPlexGetDepthLabel(dm, &depthLabel);
504: DMPlexGetCellTypeLabel(dm, &ctLabel);
505: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
506: const DMPolytopeType ict = (DMPolytopeType) c;
507: PetscInt pStart, pEnd, dep, numCorners, n = 0;
508: PetscBool output = PETSC_FALSE, doOutput;
510: if (ict == DM_POLYTOPE_FV_GHOST) continue;
511: DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
512: if (pStart >= 0) {
513: DMLabelGetValue(depthLabel, pStart, &dep);
514: if (dep == depth - cellHeight) output = PETSC_TRUE;
515: }
516: MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
517: if (!doOutput) continue;
518: CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS);
519: if (!n) {
520: PetscViewerHDF5PushGroup(viewer, "/viz/topology");
521: } else {
522: char group[PETSC_MAX_PATH_LEN];
524: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%D", n);
525: PetscViewerHDF5PushGroup(viewer, group);
526: }
527: ISView(cellIS, viewer);
528: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_corners", PETSC_INT, (void *) &numCorners);
529: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_dim", PETSC_INT, (void *) &dim);
530: ISDestroy(&cellIS);
531: PetscViewerHDF5PopGroup(viewer);
532: ++n;
533: }
534: return(0);
535: }
537: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
538: {
539: DM cdm;
540: Vec coordinates, newcoords;
541: PetscReal lengthScale;
542: PetscInt m, M, bs;
546: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
547: DMGetCoordinateDM(dm, &cdm);
548: DMGetCoordinates(dm, &coordinates);
549: VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
550: PetscObjectSetName((PetscObject) newcoords, "vertices");
551: VecGetSize(coordinates, &M);
552: VecGetLocalSize(coordinates, &m);
553: VecSetSizes(newcoords, m, M);
554: VecGetBlockSize(coordinates, &bs);
555: VecSetBlockSize(newcoords, bs);
556: VecSetType(newcoords,VECSTANDARD);
557: VecCopy(coordinates, newcoords);
558: VecScale(newcoords, lengthScale);
559: /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
560: PetscViewerHDF5PushGroup(viewer, "/geometry");
561: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
562: VecView(newcoords, viewer);
563: PetscViewerPopFormat(viewer);
564: PetscViewerHDF5PopGroup(viewer);
565: VecDestroy(&newcoords);
566: return(0);
567: }
569: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
570: {
571: DM cdm;
572: Vec coordinatesLocal, newcoords;
573: PetscSection cSection, cGlobalSection;
574: PetscScalar *coords, *ncoords;
575: DMLabel cutLabel, cutVertexLabel = NULL;
576: const PetscReal *L;
577: const DMBoundaryType *bd;
578: PetscReal lengthScale;
579: PetscInt vStart, vEnd, v, bs, N, coordSize, dof, off, d;
580: PetscBool localized, embedded;
581: hid_t fileId, groupId;
582: PetscErrorCode ierr;
585: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
586: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
587: DMGetCoordinatesLocal(dm, &coordinatesLocal);
588: VecGetBlockSize(coordinatesLocal, &bs);
589: DMGetCoordinatesLocalized(dm, &localized);
590: if (localized == PETSC_FALSE) return(0);
591: DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
592: DMGetCoordinateDM(dm, &cdm);
593: DMGetLocalSection(cdm, &cSection);
594: DMGetGlobalSection(cdm, &cGlobalSection);
595: DMGetLabel(dm, "periodic_cut", &cutLabel);
596: N = 0;
598: DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
599: VecCreate(PetscObjectComm((PetscObject) dm), &newcoords);
600: PetscSectionGetDof(cSection, vStart, &dof);
601: PetscPrintf(PETSC_COMM_SELF, "DOF: %D\n", dof);
602: embedded = (PetscBool) (L && dof == 2 && !cutLabel);
603: if (cutVertexLabel) {
604: DMLabelGetStratumSize(cutVertexLabel, 1, &v);
605: N += dof*v;
606: }
607: for (v = vStart; v < vEnd; ++v) {
608: PetscSectionGetDof(cGlobalSection, v, &dof);
609: if (dof < 0) continue;
610: if (embedded) N += dof+1;
611: else N += dof;
612: }
613: if (embedded) {VecSetBlockSize(newcoords, bs+1);}
614: else {VecSetBlockSize(newcoords, bs);}
615: VecSetSizes(newcoords, N, PETSC_DETERMINE);
616: VecSetType(newcoords, VECSTANDARD);
617: VecGetArray(coordinatesLocal, &coords);
618: VecGetArray(newcoords, &ncoords);
619: coordSize = 0;
620: for (v = vStart; v < vEnd; ++v) {
621: PetscSectionGetDof(cGlobalSection, v, &dof);
622: PetscSectionGetOffset(cSection, v, &off);
623: if (dof < 0) continue;
624: if (embedded) {
625: if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
626: PetscReal theta, phi, r, R;
627: /* XY-periodic */
628: /* Suppose its an y-z circle, then
629: \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
630: and the circle in that plane is
631: \hat r cos(phi) + \hat x sin(phi) */
632: theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
633: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
634: r = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
635: R = L[1]/(2.0*PETSC_PI);
636: ncoords[coordSize++] = PetscSinReal(phi) * r;
637: ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
638: ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
639: } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
640: /* X-periodic */
641: ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
642: ncoords[coordSize++] = coords[off+1];
643: ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
644: } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
645: /* Y-periodic */
646: ncoords[coordSize++] = coords[off+0];
647: ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
648: ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
649: } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
650: PetscReal phi, r, R;
651: /* Mobius strip */
652: /* Suppose its an x-z circle, then
653: \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
654: and in that plane we rotate by pi as we go around the circle
655: \hat r cos(phi/2) + \hat y sin(phi/2) */
656: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
657: R = L[0];
658: r = PetscRealPart(coords[off+1]) - L[1]/2.0;
659: ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
660: ncoords[coordSize++] = PetscSinReal(phi/2.0) * r;
661: ncoords[coordSize++] = PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
662: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
663: } else {
664: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
665: }
666: }
667: if (cutVertexLabel) {
668: IS vertices;
669: const PetscInt *verts;
670: PetscInt n;
672: DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
673: if (vertices) {
674: ISGetIndices(vertices, &verts);
675: ISGetLocalSize(vertices, &n);
676: for (v = 0; v < n; ++v) {
677: PetscSectionGetDof(cSection, verts[v], &dof);
678: PetscSectionGetOffset(cSection, verts[v], &off);
679: for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
680: }
681: ISRestoreIndices(vertices, &verts);
682: ISDestroy(&vertices);
683: }
684: }
685: if (coordSize != N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %D != %D", coordSize, N);
686: DMLabelDestroy(&cutVertexLabel);
687: VecRestoreArray(coordinatesLocal, &coords);
688: VecRestoreArray(newcoords, &ncoords);
689: PetscObjectSetName((PetscObject) newcoords, "vertices");
690: VecScale(newcoords, lengthScale);
691: PetscViewerHDF5PushGroup(viewer, "/viz");
692: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
693: PetscStackCallHDF5(H5Gclose,(groupId));
694: PetscViewerHDF5PopGroup(viewer);
695: PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
696: VecView(newcoords, viewer);
697: PetscViewerHDF5PopGroup(viewer);
698: VecDestroy(&newcoords);
699: return(0);
700: }
703: static PetscErrorCode DMPlexWriteLabels_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
704: {
705: const PetscInt *gpoint;
706: PetscInt numLabels, l;
707: hid_t fileId, groupId;
708: PetscErrorCode ierr;
711: ISGetIndices(globalPointNumbers, &gpoint);
712: PetscViewerHDF5PushGroup(viewer, "/labels");
713: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
714: if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
715: PetscViewerHDF5PopGroup(viewer);
716: DMGetNumLabels(dm, &numLabels);
717: for (l = 0; l < numLabels; ++l) {
718: DMLabel label;
719: const char *name;
720: IS valueIS, pvalueIS, globalValueIS;
721: const PetscInt *values;
722: PetscInt numValues, v;
723: PetscBool isDepth, output;
724: char group[PETSC_MAX_PATH_LEN];
726: DMGetLabelName(dm, l, &name);
727: DMGetLabelOutput(dm, name, &output);
728: PetscStrncmp(name, "depth", 10, &isDepth);
729: if (isDepth || !output) continue;
730: DMGetLabel(dm, name, &label);
731: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
732: PetscViewerHDF5PushGroup(viewer, group);
733: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
734: if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
735: PetscViewerHDF5PopGroup(viewer);
736: DMLabelGetValueIS(label, &valueIS);
737: /* Must copy to a new IS on the global comm */
738: ISGetLocalSize(valueIS, &numValues);
739: ISGetIndices(valueIS, &values);
740: ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
741: ISRestoreIndices(valueIS, &values);
742: ISAllGather(pvalueIS, &globalValueIS);
743: ISDestroy(&pvalueIS);
744: ISSortRemoveDups(globalValueIS);
745: ISGetLocalSize(globalValueIS, &numValues);
746: ISGetIndices(globalValueIS, &values);
747: for (v = 0; v < numValues; ++v) {
748: IS stratumIS, globalStratumIS;
749: const PetscInt *spoints = NULL;
750: PetscInt *gspoints, n = 0, gn, p;
751: const char *iname = "indices";
753: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%D", name, values[v]);
754: DMLabelGetStratumIS(label, values[v], &stratumIS);
756: if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
757: if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
758: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
759: PetscMalloc1(gn,&gspoints);
760: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
761: if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
762: ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
763: if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
764: PetscObjectSetName((PetscObject) globalStratumIS, iname);
766: PetscViewerHDF5PushGroup(viewer, group);
767: ISView(globalStratumIS, viewer);
768: PetscViewerHDF5PopGroup(viewer);
769: ISDestroy(&globalStratumIS);
770: ISDestroy(&stratumIS);
771: }
772: ISRestoreIndices(globalValueIS, &values);
773: ISDestroy(&globalValueIS);
774: ISDestroy(&valueIS);
775: }
776: ISRestoreIndices(globalPointNumbers, &gpoint);
777: return(0);
778: }
780: /* We only write cells and vertices. Does this screw up parallel reading? */
781: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
782: {
783: IS globalPointNumbers;
784: PetscViewerFormat format;
785: PetscBool viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
786: PetscErrorCode ierr;
789: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
790: DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
791: DMPlexWriteLabels_HDF5_Static(dm, globalPointNumbers, viewer);
793: PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT);
795: PetscViewerGetFormat(viewer, &format);
796: switch (format) {
797: case PETSC_VIEWER_HDF5_VIZ:
798: viz_geom = PETSC_TRUE;
799: xdmf_topo = PETSC_TRUE;
800: break;
801: case PETSC_VIEWER_HDF5_XDMF:
802: xdmf_topo = PETSC_TRUE;
803: break;
804: case PETSC_VIEWER_HDF5_PETSC:
805: petsc_topo = PETSC_TRUE;
806: break;
807: case PETSC_VIEWER_DEFAULT:
808: case PETSC_VIEWER_NATIVE:
809: viz_geom = PETSC_TRUE;
810: xdmf_topo = PETSC_TRUE;
811: petsc_topo = PETSC_TRUE;
812: break;
813: default:
814: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
815: }
817: if (viz_geom) {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
818: if (xdmf_topo) {DMPlexWriteTopology_Vertices_HDF5_Static(dm, globalPointNumbers, viewer);}
819: if (petsc_topo) {DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);}
821: ISDestroy(&globalPointNumbers);
822: return(0);
823: }
825: typedef struct {
826: PetscMPIInt rank;
827: DM dm;
828: PetscViewer viewer;
829: DMLabel label;
830: } LabelCtx;
832: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
833: {
834: PetscViewer viewer = ((LabelCtx *) op_data)->viewer;
835: DMLabel label = ((LabelCtx *) op_data)->label;
836: IS stratumIS;
837: const PetscInt *ind;
838: PetscInt value, N, i;
839: const char *lname;
840: char group[PETSC_MAX_PATH_LEN];
841: PetscErrorCode ierr;
843: PetscOptionsStringToInt(name, &value);
844: ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
845: PetscObjectSetName((PetscObject) stratumIS, "indices");
846: PetscObjectGetName((PetscObject) label, &lname);
847: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
848: PetscViewerHDF5PushGroup(viewer, group);
849: {
850: /* Force serial load */
851: PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
852: PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
853: PetscLayoutSetSize(stratumIS->map, N);
854: }
855: ISLoad(stratumIS, viewer);
856: PetscViewerHDF5PopGroup(viewer);
857: ISGetLocalSize(stratumIS, &N);
858: ISGetIndices(stratumIS, &ind);
859: for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
860: ISRestoreIndices(stratumIS, &ind);
861: ISDestroy(&stratumIS);
862: return 0;
863: }
865: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
866: {
867: DM dm = ((LabelCtx *) op_data)->dm;
868: hsize_t idx = 0;
870: herr_t err;
872: DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
873: DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
874: PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
875: return err;
876: }
878: PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM dm, PetscViewer viewer)
879: {
880: LabelCtx ctx;
881: hid_t fileId, groupId;
882: hsize_t idx = 0;
883: PetscErrorCode ierr;
886: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &ctx.rank);
887: ctx.dm = dm;
888: ctx.viewer = viewer;
889: PetscViewerHDF5PushGroup(viewer, "/labels");
890: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
891: PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
892: PetscStackCallHDF5(H5Gclose,(groupId));
893: PetscViewerHDF5PopGroup(viewer);
894: return(0);
895: }
897: /* The first version will read everything onto proc 0, letting the user distribute
898: The next will create a naive partition, and then rebalance after reading
899: */
900: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
901: {
902: PetscSection coordSection;
903: Vec coordinates;
904: IS orderIS, conesIS, cellsIS, orntsIS;
905: const PetscInt *order, *cones, *cells, *ornts;
906: PetscReal lengthScale;
907: PetscInt *cone, *ornt;
908: PetscInt dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
909: PetscMPIInt rank;
910: PetscErrorCode ierr;
913: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
914: /* Read toplogy */
915: PetscViewerHDF5PushGroup(viewer, "/topology");
916: ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
917: PetscObjectSetName((PetscObject) orderIS, "order");
918: ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
919: PetscObjectSetName((PetscObject) conesIS, "cones");
920: ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
921: PetscObjectSetName((PetscObject) cellsIS, "cells");
922: ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
923: PetscObjectSetName((PetscObject) orntsIS, "orientation");
924: PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, (void *) &dim);
925: DMSetDimension(dm, dim);
926: {
927: /* Force serial load */
928: PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
929: PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
930: PetscLayoutSetSize(orderIS->map, pEnd);
931: PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
932: PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
933: PetscLayoutSetSize(conesIS->map, pEnd);
934: pEnd = !rank ? pEnd : 0;
935: PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
936: PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
937: PetscLayoutSetSize(cellsIS->map, N);
938: PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
939: PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
940: PetscLayoutSetSize(orntsIS->map, N);
941: }
942: ISLoad(orderIS, viewer);
943: ISLoad(conesIS, viewer);
944: ISLoad(cellsIS, viewer);
945: ISLoad(orntsIS, viewer);
946: PetscViewerHDF5PopGroup(viewer);
947: /* Read geometry */
948: PetscViewerHDF5PushGroup(viewer, "/geometry");
949: VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
950: PetscObjectSetName((PetscObject) coordinates, "vertices");
951: {
952: /* Force serial load */
953: PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
954: VecSetSizes(coordinates, !rank ? N : 0, N);
955: VecSetBlockSize(coordinates, spatialDim);
956: }
957: VecLoad(coordinates, viewer);
958: PetscViewerHDF5PopGroup(viewer);
959: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
960: VecScale(coordinates, 1.0/lengthScale);
961: VecGetLocalSize(coordinates, &numVertices);
962: VecGetBlockSize(coordinates, &spatialDim);
963: numVertices /= spatialDim;
964: /* Create Plex */
965: DMPlexSetChart(dm, 0, pEnd);
966: ISGetIndices(orderIS, &order);
967: ISGetIndices(conesIS, &cones);
968: for (p = 0; p < pEnd; ++p) {
969: DMPlexSetConeSize(dm, order[p], cones[p]);
970: maxConeSize = PetscMax(maxConeSize, cones[p]);
971: }
972: DMSetUp(dm);
973: ISGetIndices(cellsIS, &cells);
974: ISGetIndices(orntsIS, &ornts);
975: PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
976: for (p = 0, q = 0; p < pEnd; ++p) {
977: for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
978: DMPlexSetCone(dm, order[p], cone);
979: DMPlexSetConeOrientation(dm, order[p], ornt);
980: }
981: PetscFree2(cone,ornt);
982: ISRestoreIndices(orderIS, &order);
983: ISRestoreIndices(conesIS, &cones);
984: ISRestoreIndices(cellsIS, &cells);
985: ISRestoreIndices(orntsIS, &ornts);
986: ISDestroy(&orderIS);
987: ISDestroy(&conesIS);
988: ISDestroy(&cellsIS);
989: ISDestroy(&orntsIS);
990: DMPlexSymmetrize(dm);
991: DMPlexStratify(dm);
992: /* Create coordinates */
993: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
994: 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);
995: DMGetCoordinateSection(dm, &coordSection);
996: PetscSectionSetNumFields(coordSection, 1);
997: PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
998: PetscSectionSetChart(coordSection, vStart, vEnd);
999: for (v = vStart; v < vEnd; ++v) {
1000: PetscSectionSetDof(coordSection, v, spatialDim);
1001: PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
1002: }
1003: PetscSectionSetUp(coordSection);
1004: DMSetCoordinates(dm, coordinates);
1005: VecDestroy(&coordinates);
1006: /* Read Labels */
1007: DMPlexLoadLabels_HDF5_Internal(dm, viewer);
1008: return(0);
1009: }
1010: #endif