Actual source code: plexhdf5.c
petsc-3.10.5 2019-03-28
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/isimpl.h>
3: #include <petsc/private/vecimpl.h>
4: #include <petscviewerhdf5.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: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
71: {
72: DM dm;
73: DM dmBC;
74: PetscSection section, sectionGlobal;
75: Vec gv;
76: const char *name;
77: PetscViewerVTKFieldType ft;
78: PetscViewerFormat format;
79: PetscInt seqnum;
80: PetscReal seqval;
81: PetscBool isseq;
82: PetscErrorCode ierr;
85: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
86: VecGetDM(v, &dm);
87: DMGetSection(dm, §ion);
88: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
89: PetscViewerHDF5SetTimestep(viewer, seqnum);
90: DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
91: PetscViewerGetFormat(viewer, &format);
92: DMGetOutputDM(dm, &dmBC);
93: DMGetGlobalSection(dmBC, §ionGlobal);
94: DMGetGlobalVector(dmBC, &gv);
95: PetscObjectGetName((PetscObject) v, &name);
96: PetscObjectSetName((PetscObject) gv, name);
97: DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
98: DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
99: PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
100: if (isseq) {VecView_Seq(gv, viewer);}
101: else {VecView_MPI(gv, viewer);}
102: if (format == PETSC_VIEWER_HDF5_VIZ) {
103: /* Output visualization representation */
104: PetscInt numFields, f;
105: DMLabel cutLabel;
107: PetscSectionGetNumFields(section, &numFields);
108: DMGetLabel(dm, "periodic_cut", &cutLabel);
109: for (f = 0; f < numFields; ++f) {
110: Vec subv;
111: IS is;
112: const char *fname, *fgroup;
113: char subname[PETSC_MAX_PATH_LEN];
114: char group[PETSC_MAX_PATH_LEN];
115: PetscInt pStart, pEnd;
116: PetscBool flag;
118: DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
119: fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
120: PetscSectionGetFieldName(section, f, &fname);
121: if (!fname) continue;
122: PetscViewerHDF5PushGroup(viewer, fgroup);
123: if (cutLabel) {
124: const PetscScalar *ga;
125: PetscScalar *suba;
126: PetscInt Nc, gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;
128: PetscSectionGetFieldComponents(section, f, &Nc);
129: for (p = pStart; p < pEnd; ++p) {
130: PetscInt gdof, fdof = 0, val;
132: PetscSectionGetDof(sectionGlobal, p, &gdof);
133: if (gdof > 0) {PetscSectionGetFieldDof(section, p, f, &fdof);}
134: subSize += fdof;
135: DMLabelGetValue(cutLabel, p, &val);
136: if (val == 1) extSize += fdof;
137: }
138: VecCreate(PetscObjectComm((PetscObject) gv), &subv);
139: VecSetSizes(subv, subSize+extSize, PETSC_DETERMINE);
140: VecSetBlockSize(subv, Nc);
141: VecSetType(subv, VECSTANDARD);
142: VecGetOwnershipRange(gv, &gstart, NULL);
143: VecGetArrayRead(gv, &ga);
144: VecGetArray(subv, &suba);
145: for (p = pStart; p < pEnd; ++p) {
146: PetscInt gdof, goff, val;
148: PetscSectionGetDof(sectionGlobal, p, &gdof);
149: if (gdof > 0) {
150: PetscInt fdof, fc, f2, poff = 0;
152: PetscSectionGetOffset(sectionGlobal, p, &goff);
153: /* Can get rid of this loop by storing field information in the global section */
154: for (f2 = 0; f2 < f; ++f2) {
155: PetscSectionGetFieldDof(section, p, f2, &fdof);
156: poff += fdof;
157: }
158: PetscSectionGetFieldDof(section, p, f, &fdof);
159: for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff+poff+fc - gstart];
160: DMLabelGetValue(cutLabel, p, &val);
161: if (val == 1) {
162: for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize+newOff] = ga[goff+poff+fc - gstart];
163: }
164: }
165: }
166: VecRestoreArrayRead(gv, &ga);
167: VecRestoreArray(subv, &suba);
168: } else {
169: PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
170: }
171: PetscStrncpy(subname, name,sizeof(subname));
172: PetscStrlcat(subname, "_",sizeof(subname));
173: PetscStrlcat(subname, fname,sizeof(subname));
174: PetscObjectSetName((PetscObject) subv, subname);
175: if (isseq) {VecView_Seq(subv, viewer);}
176: else {VecView_MPI(subv, viewer);}
177: if (cutLabel) {
178: VecDestroy(&subv);
179: } else {
180: PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
181: }
182: PetscViewerHDF5PopGroup(viewer);
183: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "%s/%s", fgroup, subname);
184: PetscViewerHDF5HasAttribute(viewer, group, "vector_field_type", &flag);
185: if (!flag) {
186: if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
187: PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "vector");
188: } else {
189: PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "scalar");
190: }
191: }
192: }
193: }
194: DMRestoreGlobalVector(dmBC, &gv);
195: return(0);
196: }
198: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
199: {
200: DM dm;
201: Vec locv;
202: const char *name;
203: PetscReal time;
207: VecGetDM(v, &dm);
208: DMGetLocalVector(dm, &locv);
209: PetscObjectGetName((PetscObject) v, &name);
210: PetscObjectSetName((PetscObject) locv, name);
211: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
212: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
213: DMGetOutputSequenceNumber(dm, NULL, &time);
214: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
215: PetscViewerHDF5PushGroup(viewer, "/fields");
216: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
217: VecView_Plex_Local_HDF5_Internal(locv, viewer);
218: PetscViewerPopFormat(viewer);
219: PetscViewerHDF5PopGroup(viewer);
220: DMRestoreLocalVector(dm, &locv);
221: return(0);
222: }
224: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
225: {
226: PetscBool isseq;
230: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
231: PetscViewerHDF5PushGroup(viewer, "/fields");
232: if (isseq) {VecView_Seq(v, viewer);}
233: else {VecView_MPI(v, viewer);}
234: PetscViewerHDF5PopGroup(viewer);
235: return(0);
236: }
238: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
239: {
240: DM dm;
241: Vec locv;
242: const char *name;
243: PetscInt seqnum;
247: VecGetDM(v, &dm);
248: DMGetLocalVector(dm, &locv);
249: PetscObjectGetName((PetscObject) v, &name);
250: PetscObjectSetName((PetscObject) locv, name);
251: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
252: PetscViewerHDF5SetTimestep(viewer, seqnum);
253: PetscViewerHDF5PushGroup(viewer, "/fields");
254: VecLoad_Plex_Local(locv, viewer);
255: PetscViewerHDF5PopGroup(viewer);
256: DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
257: DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
258: DMRestoreLocalVector(dm, &locv);
259: return(0);
260: }
262: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
263: {
264: DM dm;
265: PetscInt seqnum;
269: VecGetDM(v, &dm);
270: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
271: PetscViewerHDF5SetTimestep(viewer, seqnum);
272: PetscViewerHDF5PushGroup(viewer, "/fields");
273: VecLoad_Default(v, viewer);
274: PetscViewerHDF5PopGroup(viewer);
275: return(0);
276: }
278: static PetscErrorCode DMPlexWriteTopology_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
279: {
280: IS orderIS, conesIS, cellsIS, orntsIS;
281: const PetscInt *gpoint;
282: PetscInt *order, *sizes, *cones, *ornts;
283: PetscInt dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
284: PetscErrorCode ierr;
287: ISGetIndices(globalPointNumbers, &gpoint);
288: DMGetDimension(dm, &dim);
289: DMPlexGetChart(dm, &pStart, &pEnd);
290: for (p = pStart; p < pEnd; ++p) {
291: if (gpoint[p] >= 0) {
292: PetscInt coneSize;
294: DMPlexGetConeSize(dm, p, &coneSize);
295: conesSize += 1;
296: cellsSize += coneSize;
297: }
298: }
299: PetscMalloc1(conesSize, &order);
300: PetscMalloc1(conesSize, &sizes);
301: PetscMalloc1(cellsSize, &cones);
302: PetscMalloc1(cellsSize, &ornts);
303: for (p = pStart; p < pEnd; ++p) {
304: if (gpoint[p] >= 0) {
305: const PetscInt *cone, *ornt;
306: PetscInt coneSize, cp;
308: DMPlexGetConeSize(dm, p, &coneSize);
309: DMPlexGetCone(dm, p, &cone);
310: DMPlexGetConeOrientation(dm, p, &ornt);
311: order[s] = gpoint[p];
312: sizes[s++] = coneSize;
313: for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
314: }
315: }
316: if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %d != %d", s, conesSize);
317: if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %d != %d", c, cellsSize);
318: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
319: PetscObjectSetName((PetscObject) orderIS, "order");
320: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
321: PetscObjectSetName((PetscObject) conesIS, "cones");
322: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
323: PetscObjectSetName((PetscObject) cellsIS, "cells");
324: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
325: PetscObjectSetName((PetscObject) orntsIS, "orientation");
326: PetscViewerHDF5PushGroup(viewer, "/topology");
327: ISView(orderIS, viewer);
328: ISView(conesIS, viewer);
329: ISView(cellsIS, viewer);
330: ISView(orntsIS, viewer);
331: PetscViewerHDF5PopGroup(viewer);
332: ISDestroy(&orderIS);
333: ISDestroy(&conesIS);
334: ISDestroy(&cellsIS);
335: ISDestroy(&orntsIS);
336: ISRestoreIndices(globalPointNumbers, &gpoint);
338: PetscViewerHDF5WriteAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
339: return(0);
340: }
342: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
343: {
344: DMLabel cutLabel;
345: PetscSection cSection;
346: IS cellIS, globalVertexNumbers;
347: const PetscInt *gvertex;
348: PetscInt *vertices;
349: IS cutvertices;
350: const PetscInt *cutverts;
351: PetscInt dim, depth, vStart, vEnd, vExtra = 0, v, cellHeight, cStart, cMax, cEnd, cell, conesSize = 0, numCornersLocal = 0, numCorners;
352: hid_t fileId, groupId;
353: PetscErrorCode ierr;
356: DMGetDimension(dm, &dim);
357: DMPlexGetDepth(dm, &depth);
358: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
359: DMGetCoordinateSection(dm, &cSection);
360: DMPlexGetVTKCellHeight(dm, &cellHeight);
361: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
362: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
363: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
364: for (cell = cStart; cell < cEnd; ++cell) {
365: PetscInt *closure = NULL;
366: PetscInt closureSize, v, Nc = 0;
368: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
369: for (v = 0; v < closureSize*2; v += 2) {
370: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
371: }
372: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
373: conesSize += Nc;
374: if (!numCornersLocal) numCornersLocal = Nc;
375: else if (numCornersLocal != Nc) numCornersLocal = 1;
376: }
377: MPIU_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
378: if (numCornersLocal && (numCornersLocal != numCorners || numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
380: PetscViewerHDF5PushGroup(viewer, "/viz");
381: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
382: PetscStackCallHDF5(H5Gclose,(groupId));
383: PetscViewerHDF5PopGroup(viewer);
385: DMGetLabel(dm, "periodic_cut", &cutLabel);
386: if (cutLabel) {
387: DMLabelGetStratumIS(cutLabel, 1, &cutvertices);
388: ISGetIndices(cutvertices, &cutverts);
389: ISGetLocalSize(cutvertices, &vExtra);
390: }
391: DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
392: ISGetIndices(globalVertexNumbers, &gvertex);
393: PetscMalloc1(conesSize, &vertices);
394: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
395: PetscInt *closure = NULL;
396: PetscInt closureSize, Nc = 0, p;
397: PetscBool replace = PETSC_FALSE;
399: if (cutLabel) {
400: PetscInt value;
401: DMLabelGetValue(cutLabel, cell, &value);
402: if (value == 2) replace = PETSC_TRUE;
403: }
404: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
405: for (p = 0; p < closureSize*2; p += 2) {
406: if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
407: closure[Nc++] = closure[p];
408: }
409: }
410: DMPlexInvertCell_Internal(dim, Nc, closure);
411: for (p = 0; p < Nc; ++p) {
412: const PetscInt gv = gvertex[closure[p] - vStart];
413: vertices[v++] = gv < 0 ? -(gv+1) : gv;
414: if (replace) {
415: PetscInt newv;
416: PetscFindInt(closure[p], vExtra, cutverts, &newv);
417: if (newv >= 0) vertices[v-1] = vEnd - vStart + newv;
418: }
419: }
420: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
421: }
422: if (cutLabel) {
423: ISRestoreIndices(cutvertices, &cutverts);
424: ISDestroy(&cutvertices);
425: }
426: if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %d != %d", v, conesSize);
427: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, &cellIS);
428: PetscLayoutSetBlockSize(cellIS->map, numCorners);
429: PetscObjectSetName((PetscObject) cellIS, "cells");
430: PetscViewerHDF5PushGroup(viewer, "/viz/topology");
431: ISView(cellIS, viewer);
432: PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
433: ISDestroy(&cellIS);
435: PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
436: return(0);
437: }
439: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
440: {
441: DM cdm;
442: Vec coordinates, newcoords;
443: PetscReal lengthScale;
444: PetscInt m, M, bs;
448: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
449: DMGetCoordinateDM(dm, &cdm);
450: DMGetCoordinates(dm, &coordinates);
451: VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
452: PetscObjectSetName((PetscObject) newcoords, "vertices");
453: VecGetSize(coordinates, &M);
454: VecGetLocalSize(coordinates, &m);
455: VecSetSizes(newcoords, m, M);
456: VecGetBlockSize(coordinates, &bs);
457: VecSetBlockSize(newcoords, bs);
458: VecSetType(newcoords,VECSTANDARD);
459: VecCopy(coordinates, newcoords);
460: VecScale(newcoords, lengthScale);
461: /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
462: PetscViewerHDF5PushGroup(viewer, "/geometry");
463: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
464: VecView(newcoords, viewer);
465: PetscViewerPopFormat(viewer);
466: PetscViewerHDF5PopGroup(viewer);
467: VecDestroy(&newcoords);
468: return(0);
469: }
471: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
472: {
473: Vec coordinates, newcoords;
474: PetscSection cSection;
475: PetscScalar *coords, *ncoords;
476: DMLabel cutLabel;
477: const PetscReal *L;
478: const DMBoundaryType *bd;
479: PetscReal lengthScale;
480: PetscInt vStart, vEnd, vExtra = 0, v, bs, coordSize, dof, off, d;
481: PetscBool localized, embedded;
482: hid_t fileId, groupId;
483: PetscErrorCode ierr;
486: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
487: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
488: DMGetCoordinatesLocal(dm, &coordinates);
489: VecGetBlockSize(coordinates, &bs);
490: VecGetLocalSize(coordinates, &coordSize);
491: DMGetCoordinatesLocalized(dm,&localized);
492: if (localized == PETSC_FALSE) return(0);
493: DMGetLabel(dm, "periodic_cut", &cutLabel);
494: if (cutLabel) {
495: IS vertices;
496: const PetscInt *verts;
497: PetscInt n;
499: DMLabelGetStratumIS(cutLabel, 1, &vertices);
500: ISGetIndices(vertices, &verts);
501: ISGetLocalSize(vertices, &n);
502: for (v = 0; v < n; ++v) {
503: if ((verts[v] >= vStart) && (verts[v] < vEnd)) ++vExtra;
504: }
505: ISRestoreIndices(vertices, &verts);
506: ISDestroy(&vertices);
507: }
508: DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
509: DMGetCoordinateSection(dm, &cSection);
510: VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
511: PetscSectionGetDof(cSection, vStart, &dof);
512: embedded = (PetscBool) (L && dof == 2 && !cutLabel);
513: coordSize = 0;
514: coordSize += dof*vExtra;
515: for (v = vStart; v < vEnd; ++v) {
516: PetscSectionGetDof(cSection, v, &dof);
517: if (embedded) coordSize += dof+1;
518: else coordSize += dof;
519: }
520: if (embedded) {VecSetBlockSize(newcoords, bs+1);}
521: else {VecSetBlockSize(newcoords, bs);}
522: VecSetSizes(newcoords, coordSize, PETSC_DETERMINE);
523: VecSetType(newcoords, VECSTANDARD);
524: VecGetArray(coordinates, &coords);
525: VecGetArray(newcoords, &ncoords);
526: coordSize = 0;
527: for (v = vStart; v < vEnd; ++v) {
528: PetscSectionGetDof(cSection, v, &dof);
529: PetscSectionGetOffset(cSection, v, &off);
530: if (embedded) {
531: if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
532: PetscReal theta, phi, r, R;
533: /* XY-periodic */
534: /* Suppose its an y-z circle, then
535: \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
536: and the circle in that plane is
537: \hat r cos(phi) + \hat x sin(phi) */
538: theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
539: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
540: r = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
541: R = L[1]/(2.0*PETSC_PI);
542: ncoords[coordSize++] = PetscSinReal(phi) * r;
543: ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
544: ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
545: } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
546: /* X-periodic */
547: ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
548: ncoords[coordSize++] = coords[off+1];
549: ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
550: } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
551: /* Y-periodic */
552: ncoords[coordSize++] = coords[off+0];
553: ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
554: ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
555: } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
556: PetscReal phi, r, R;
557: /* Mobius strip */
558: /* Suppose its an x-z circle, then
559: \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
560: and in that plane we rotate by pi as we go around the circle
561: \hat r cos(phi/2) + \hat y sin(phi/2) */
562: phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
563: R = L[0];
564: r = PetscRealPart(coords[off+1]) - L[1]/2.0;
565: ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
566: ncoords[coordSize++] = PetscSinReal(phi/2.0) * r;
567: ncoords[coordSize++] = PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
568: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
569: } else {
570: if (cutLabel) {
571: DMLocalizeCoordinate(dm, &coords[off], PETSC_TRUE, &ncoords[coordSize]);
572: coordSize += dof;
573: } else {
574: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
575: }
576: }
577: }
578: if (cutLabel) {
579: IS vertices;
580: const PetscInt *verts;
581: PetscInt n;
583: DMLabelGetStratumIS(cutLabel, 1, &vertices);
584: ISGetIndices(vertices, &verts);
585: ISGetLocalSize(vertices, &n);
586: for (v = 0; v < n; ++v) {
587: if ((verts[v] < vStart) || (verts[v] >= vEnd)) continue;
588: PetscSectionGetDof(cSection, verts[v], &dof);
589: PetscSectionGetOffset(cSection, verts[v], &off);
590: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
591: }
592: ISRestoreIndices(vertices, &verts);
593: ISDestroy(&vertices);
594: }
595: VecRestoreArray(coordinates, &coords);
596: VecRestoreArray(newcoords, &ncoords);
597: PetscObjectSetName((PetscObject) newcoords, "vertices");
598: VecScale(newcoords, lengthScale);
599: PetscViewerHDF5PushGroup(viewer, "/viz");
600: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
601: PetscStackCallHDF5(H5Gclose,(groupId));
602: PetscViewerHDF5PopGroup(viewer);
603: PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
604: VecView(newcoords, viewer);
605: PetscViewerHDF5PopGroup(viewer);
606: VecDestroy(&newcoords);
607: return(0);
608: }
611: static PetscErrorCode DMPlexWriteLabels_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
612: {
613: const PetscInt *gpoint;
614: PetscInt numLabels, l;
615: hid_t fileId, groupId;
616: PetscErrorCode ierr;
619: ISGetIndices(globalPointNumbers, &gpoint);
620: PetscViewerHDF5PushGroup(viewer, "/labels");
621: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
622: if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
623: PetscViewerHDF5PopGroup(viewer);
624: DMGetNumLabels(dm, &numLabels);
625: for (l = 0; l < numLabels; ++l) {
626: DMLabel label;
627: const char *name;
628: IS valueIS, pvalueIS, globalValueIS;
629: const PetscInt *values;
630: PetscInt numValues, v;
631: PetscBool isDepth, output;
632: char group[PETSC_MAX_PATH_LEN];
634: DMGetLabelName(dm, l, &name);
635: DMGetLabelOutput(dm, name, &output);
636: PetscStrncmp(name, "depth", 10, &isDepth);
637: if (isDepth || !output) continue;
638: DMGetLabel(dm, name, &label);
639: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
640: PetscViewerHDF5PushGroup(viewer, group);
641: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
642: if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
643: PetscViewerHDF5PopGroup(viewer);
644: DMLabelGetValueIS(label, &valueIS);
645: /* Must copy to a new IS on the global comm */
646: ISGetLocalSize(valueIS, &numValues);
647: ISGetIndices(valueIS, &values);
648: ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
649: ISRestoreIndices(valueIS, &values);
650: ISAllGather(pvalueIS, &globalValueIS);
651: ISDestroy(&pvalueIS);
652: ISSortRemoveDups(globalValueIS);
653: ISGetLocalSize(globalValueIS, &numValues);
654: ISGetIndices(globalValueIS, &values);
655: for (v = 0; v < numValues; ++v) {
656: IS stratumIS, globalStratumIS;
657: const PetscInt *spoints = NULL;
658: PetscInt *gspoints, n = 0, gn, p;
659: const char *iname = "indices";
661: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%d", name, values[v]);
662: DMLabelGetStratumIS(label, values[v], &stratumIS);
664: if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
665: if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
666: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
667: PetscMalloc1(gn,&gspoints);
668: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
669: if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
670: ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
671: if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
672: PetscObjectSetName((PetscObject) globalStratumIS, iname);
674: PetscViewerHDF5PushGroup(viewer, group);
675: ISView(globalStratumIS, viewer);
676: PetscViewerHDF5PopGroup(viewer);
677: ISDestroy(&globalStratumIS);
678: ISDestroy(&stratumIS);
679: }
680: ISRestoreIndices(globalValueIS, &values);
681: ISDestroy(&globalValueIS);
682: ISDestroy(&valueIS);
683: }
684: ISRestoreIndices(globalPointNumbers, &gpoint);
685: return(0);
686: }
688: /* We only write cells and vertices. Does this screw up parallel reading? */
689: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
690: {
691: IS globalPointNumbers;
692: PetscViewerFormat format;
693: PetscBool viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
694: PetscErrorCode ierr;
697: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
698: DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
699: DMPlexWriteLabels_HDF5_Static(dm, globalPointNumbers, viewer);
701: PetscViewerGetFormat(viewer, &format);
702: switch (format) {
703: case PETSC_VIEWER_HDF5_VIZ:
704: viz_geom = PETSC_TRUE;
705: xdmf_topo = PETSC_TRUE;
706: break;
707: case PETSC_VIEWER_HDF5_XDMF:
708: xdmf_topo = PETSC_TRUE;
709: break;
710: case PETSC_VIEWER_HDF5_PETSC:
711: petsc_topo = PETSC_TRUE;
712: break;
713: case PETSC_VIEWER_DEFAULT:
714: case PETSC_VIEWER_NATIVE:
715: viz_geom = PETSC_TRUE;
716: xdmf_topo = PETSC_TRUE;
717: petsc_topo = PETSC_TRUE;
718: break;
719: default:
720: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
721: }
723: if (viz_geom) {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
724: if (xdmf_topo) {DMPlexWriteTopology_Vertices_HDF5_Static(dm, viewer);}
725: if (petsc_topo) {DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);}
727: ISDestroy(&globalPointNumbers);
728: return(0);
729: }
731: typedef struct {
732: PetscMPIInt rank;
733: DM dm;
734: PetscViewer viewer;
735: DMLabel label;
736: } LabelCtx;
738: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
739: {
740: PetscViewer viewer = ((LabelCtx *) op_data)->viewer;
741: DMLabel label = ((LabelCtx *) op_data)->label;
742: IS stratumIS;
743: const PetscInt *ind;
744: PetscInt value, N, i;
745: const char *lname;
746: char group[PETSC_MAX_PATH_LEN];
747: PetscErrorCode ierr;
749: PetscOptionsStringToInt(name, &value);
750: ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
751: PetscObjectSetName((PetscObject) stratumIS, "indices");
752: DMLabelGetName(label, &lname);
753: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
754: PetscViewerHDF5PushGroup(viewer, group);
755: {
756: /* Force serial load */
757: PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
758: PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
759: PetscLayoutSetSize(stratumIS->map, N);
760: }
761: ISLoad(stratumIS, viewer);
762: PetscViewerHDF5PopGroup(viewer);
763: ISGetLocalSize(stratumIS, &N);
764: ISGetIndices(stratumIS, &ind);
765: for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
766: ISRestoreIndices(stratumIS, &ind);
767: ISDestroy(&stratumIS);
768: return 0;
769: }
771: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
772: {
773: DM dm = ((LabelCtx *) op_data)->dm;
774: hsize_t idx = 0;
776: herr_t err;
778: DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
779: DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
780: PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
781: return err;
782: }
784: PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM dm, PetscViewer viewer)
785: {
786: LabelCtx ctx;
787: hid_t fileId, groupId;
788: hsize_t idx = 0;
789: PetscErrorCode ierr;
792: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &ctx.rank);
793: ctx.dm = dm;
794: ctx.viewer = viewer;
795: PetscViewerHDF5PushGroup(viewer, "/labels");
796: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
797: PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
798: PetscStackCallHDF5(H5Gclose,(groupId));
799: PetscViewerHDF5PopGroup(viewer);
800: return(0);
801: }
803: /* The first version will read everything onto proc 0, letting the user distribute
804: The next will create a naive partition, and then rebalance after reading
805: */
806: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
807: {
808: PetscSection coordSection;
809: Vec coordinates;
810: IS orderIS, conesIS, cellsIS, orntsIS;
811: const PetscInt *order, *cones, *cells, *ornts;
812: PetscReal lengthScale;
813: PetscInt *cone, *ornt;
814: PetscInt dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
815: PetscMPIInt rank;
816: PetscErrorCode ierr;
819: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
820: /* Read toplogy */
821: PetscViewerHDF5ReadAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
822: DMSetDimension(dm, dim);
823: PetscViewerHDF5PushGroup(viewer, "/topology");
825: ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
826: PetscObjectSetName((PetscObject) orderIS, "order");
827: ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
828: PetscObjectSetName((PetscObject) conesIS, "cones");
829: ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
830: PetscObjectSetName((PetscObject) cellsIS, "cells");
831: ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
832: PetscObjectSetName((PetscObject) orntsIS, "orientation");
833: {
834: /* Force serial load */
835: PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
836: PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
837: PetscLayoutSetSize(orderIS->map, pEnd);
838: PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
839: PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
840: PetscLayoutSetSize(conesIS->map, pEnd);
841: pEnd = !rank ? pEnd : 0;
842: PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
843: PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
844: PetscLayoutSetSize(cellsIS->map, N);
845: PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
846: PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
847: PetscLayoutSetSize(orntsIS->map, N);
848: }
849: ISLoad(orderIS, viewer);
850: ISLoad(conesIS, viewer);
851: ISLoad(cellsIS, viewer);
852: ISLoad(orntsIS, viewer);
853: PetscViewerHDF5PopGroup(viewer);
854: /* Read geometry */
855: PetscViewerHDF5PushGroup(viewer, "/geometry");
856: VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
857: PetscObjectSetName((PetscObject) coordinates, "vertices");
858: {
859: /* Force serial load */
860: PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
861: VecSetSizes(coordinates, !rank ? N : 0, N);
862: VecSetBlockSize(coordinates, spatialDim);
863: }
864: VecLoad(coordinates, viewer);
865: PetscViewerHDF5PopGroup(viewer);
866: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
867: VecScale(coordinates, 1.0/lengthScale);
868: VecGetLocalSize(coordinates, &numVertices);
869: VecGetBlockSize(coordinates, &spatialDim);
870: numVertices /= spatialDim;
871: /* Create Plex */
872: DMPlexSetChart(dm, 0, pEnd);
873: ISGetIndices(orderIS, &order);
874: ISGetIndices(conesIS, &cones);
875: for (p = 0; p < pEnd; ++p) {
876: DMPlexSetConeSize(dm, order[p], cones[p]);
877: maxConeSize = PetscMax(maxConeSize, cones[p]);
878: }
879: DMSetUp(dm);
880: ISGetIndices(cellsIS, &cells);
881: ISGetIndices(orntsIS, &ornts);
882: PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
883: for (p = 0, q = 0; p < pEnd; ++p) {
884: for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
885: DMPlexSetCone(dm, order[p], cone);
886: DMPlexSetConeOrientation(dm, order[p], ornt);
887: }
888: PetscFree2(cone,ornt);
889: ISRestoreIndices(orderIS, &order);
890: ISRestoreIndices(conesIS, &cones);
891: ISRestoreIndices(cellsIS, &cells);
892: ISRestoreIndices(orntsIS, &ornts);
893: ISDestroy(&orderIS);
894: ISDestroy(&conesIS);
895: ISDestroy(&cellsIS);
896: ISDestroy(&orntsIS);
897: DMPlexSymmetrize(dm);
898: DMPlexStratify(dm);
899: /* Create coordinates */
900: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
901: 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);
902: DMGetCoordinateSection(dm, &coordSection);
903: PetscSectionSetNumFields(coordSection, 1);
904: PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
905: PetscSectionSetChart(coordSection, vStart, vEnd);
906: for (v = vStart; v < vEnd; ++v) {
907: PetscSectionSetDof(coordSection, v, spatialDim);
908: PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
909: }
910: PetscSectionSetUp(coordSection);
911: DMSetCoordinates(dm, coordinates);
912: VecDestroy(&coordinates);
913: /* Read Labels */
914: DMPlexLoadLabels_HDF5_Internal(dm, viewer);
915: return(0);
916: }
917: #endif