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