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, §ion);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, §ionGlobal);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