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