Actual source code: plexvtu.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>

  4: typedef struct {
  5:   PetscInt nvertices;
  6:   PetscInt ncells;
  7:   PetscInt nconn; /* number of entries in cell->vertex connectivity array */
  8: } PieceInfo;

 10: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
 11: /* output in float if single or half precision in memory */
 12: static const char precision[] = "Float32";
 13: typedef float     PetscVTUReal;
 14:   #define MPIU_VTUREAL MPI_FLOAT
 15: #elif defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
 16: /* output in double if double or quad precision in memory */
 17: static const char precision[] = "Float64";
 18: typedef double    PetscVTUReal;
 19:   #define MPIU_VTUREAL MPI_DOUBLE
 20: #else
 21: static const char precision[] = "UnknownPrecision";
 22: typedef PetscReal PetscVTUReal;
 23:   #define MPIU_VTUREAL MPIU_REAL
 24: #endif

 26: static PetscErrorCode TransferWrite(MPI_Comm comm, PetscViewer viewer, FILE *fp, PetscMPIInt srank, PetscMPIInt root, const void *send, void *recv, PetscMPIInt count, MPI_Datatype mpidatatype, PetscMPIInt tag)
 27: {
 28:   PetscMPIInt rank;

 30:   MPI_Comm_rank(comm, &rank);
 31:   if (rank == srank && rank != root) {
 32:     MPI_Send((void *)send, count, mpidatatype, root, tag, comm);
 33:   } else if (rank == root) {
 34:     const void *buffer;
 35:     if (root == srank) { /* self */
 36:       buffer = send;
 37:     } else {
 38:       MPI_Status  status;
 39:       PetscMPIInt nrecv;
 40:       MPI_Recv(recv, count, mpidatatype, srank, tag, comm, &status);
 41:       MPI_Get_count(&status, mpidatatype, &nrecv);
 43:       buffer = recv;
 44:     }
 45:     PetscViewerVTKFWrite(viewer, fp, buffer, count, mpidatatype);
 46:   }
 47:   return 0;
 48: }

 50: static PetscErrorCode DMPlexGetVTKConnectivity(DM dm, PetscBool localized, PieceInfo *piece, PetscVTKInt **oconn, PetscVTKInt **ooffsets, PetscVTKType **otypes)
 51: {
 52:   PetscSection  coordSection, cellCoordSection;
 53:   PetscVTKInt  *conn, *offsets;
 54:   PetscVTKType *types;
 55:   PetscInt      dim, vStart, vEnd, cStart, cEnd, pStart, pEnd, cellHeight, numLabelCells, hasLabel, c, v, countcell, countconn;

 57:   PetscMalloc3(piece->nconn, &conn, piece->ncells, &offsets, piece->ncells, &types);
 58:   DMGetCoordinateSection(dm, &coordSection);
 59:   DMGetCellCoordinateSection(dm, &cellCoordSection);
 60:   DMGetDimension(dm, &dim);
 61:   DMPlexGetChart(dm, &pStart, &pEnd);
 62:   DMPlexGetVTKCellHeight(dm, &cellHeight);
 63:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
 64:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 65:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
 66:   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;

 68:   countcell = 0;
 69:   countconn = 0;
 70:   for (c = cStart; c < cEnd; ++c) {
 71:     PetscInt nverts, dof = 0, celltype, startoffset, nC = 0;

 73:     if (hasLabel) {
 74:       PetscInt value;

 76:       DMGetLabelValue(dm, "vtk", c, &value);
 77:       if (value != 1) continue;
 78:     }
 79:     startoffset = countconn;
 80:     if (localized) PetscSectionGetDof(cellCoordSection, c, &dof);
 81:     if (!dof) {
 82:       PetscInt *closure = NULL;
 83:       PetscInt  closureSize;

 85:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
 86:       for (v = 0; v < closureSize * 2; v += 2) {
 87:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
 88:           if (!localized) conn[countconn++] = closure[v] - vStart;
 89:           else conn[countconn++] = startoffset + nC;
 90:           ++nC;
 91:         }
 92:       }
 93:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
 94:     } else {
 95:       for (nC = 0; nC < dof / dim; nC++) conn[countconn++] = startoffset + nC;
 96:     }

 98:     {
 99:       PetscInt n = PetscMin(nC, 8), s = countconn - nC, i, cone[8];
100:       for (i = 0; i < n; ++i) cone[i] = conn[s + i];
101:       DMPlexReorderCell(dm, c, cone);
102:       for (i = 0; i < n; ++i) conn[s + i] = (int)cone[i];
103:     }

105:     offsets[countcell] = countconn;

107:     nverts = countconn - startoffset;
108:     DMPlexVTKGetCellType_Internal(dm, dim, nverts, &celltype);

110:     types[countcell] = celltype;
111:     countcell++;
112:   }
115:   *oconn    = conn;
116:   *ooffsets = offsets;
117:   *otypes   = types;
118:   return 0;
119: }

121: /*
122:   Write all fields that have been provided to the viewer
123:   Multi-block XML format with binary appended data.
124: */
125: PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm, PetscViewer viewer)
126: {
127:   MPI_Comm                 comm;
128:   PetscSection             coordSection, cellCoordSection;
129:   PetscViewer_VTK         *vtk = (PetscViewer_VTK *)viewer->data;
130:   PetscViewerVTKObjectLink link;
131:   FILE                    *fp;
132:   PetscMPIInt              rank, size, tag;
133:   PetscInt                 dimEmbed, cellHeight, cStart, cEnd, vStart, vEnd, numLabelCells, hasLabel, c, v, r, i;
134:   PetscBool                localized;
135:   PieceInfo                piece, *gpiece = NULL;
136:   void                    *buffer     = NULL;
137:   const char              *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
138:   PetscInt                 loops_per_scalar;

140:   PetscObjectGetComm((PetscObject)dm, &comm);
141: #if defined(PETSC_USE_COMPLEX)
142:   loops_per_scalar = 2;
143: #else
144:   loops_per_scalar = 1;
145: #endif
146:   MPI_Comm_size(comm, &size);
147:   MPI_Comm_rank(comm, &rank);
148:   PetscCommGetNewTag(comm, &tag);

150:   PetscFOpen(comm, vtk->filename, "wb", &fp);
151:   PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n");
152:   PetscFPrintf(comm, fp, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order);
153:   PetscFPrintf(comm, fp, "  <UnstructuredGrid>\n");

155:   DMGetCoordinateDim(dm, &dimEmbed);
156:   DMPlexGetVTKCellHeight(dm, &cellHeight);
157:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
158:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
159:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
160:   DMGetCoordinatesLocalized(dm, &localized);
161:   DMGetCoordinateSection(dm, &coordSection);
162:   DMGetCellCoordinateSection(dm, &cellCoordSection);

164:   hasLabel        = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
165:   piece.nvertices = 0;
166:   piece.ncells    = 0;
167:   piece.nconn     = 0;
168:   if (!localized) piece.nvertices = vEnd - vStart;
169:   for (c = cStart; c < cEnd; ++c) {
170:     PetscInt dof = 0;
171:     if (hasLabel) {
172:       PetscInt value;

174:       DMGetLabelValue(dm, "vtk", c, &value);
175:       if (value != 1) continue;
176:     }
177:     if (localized) PetscSectionGetDof(cellCoordSection, c, &dof);
178:     if (!dof) {
179:       PetscInt *closure = NULL;
180:       PetscInt  closureSize;

182:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
183:       for (v = 0; v < closureSize * 2; v += 2) {
184:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
185:           piece.nconn++;
186:           if (localized) piece.nvertices++;
187:         }
188:       }
189:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
190:     } else {
191:       piece.nvertices += dof / dimEmbed;
192:       piece.nconn += dof / dimEmbed;
193:     }
194:     piece.ncells++;
195:   }
196:   if (rank == 0) PetscMalloc1(size, &gpiece);
197:   MPI_Gather((PetscInt *)&piece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, (PetscInt *)gpiece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, 0, comm);

199:   /*
200:    * Write file header
201:    */
202:   if (rank == 0) {
203:     PetscInt64 boffset = 0;

205:     for (r = 0; r < size; r++) {
206:       PetscFPrintf(PETSC_COMM_SELF, fp, "    <Piece NumberOfPoints=\"%" PetscInt_FMT "\" NumberOfCells=\"%" PetscInt_FMT "\">\n", gpiece[r].nvertices, gpiece[r].ncells);
207:       /* Coordinate positions */
208:       PetscFPrintf(PETSC_COMM_SELF, fp, "      <Points>\n");
209:       PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset);
210:       boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + sizeof(int);
211:       PetscFPrintf(PETSC_COMM_SELF, fp, "      </Points>\n");
212:       /* Cell connectivity */
213:       PetscFPrintf(PETSC_COMM_SELF, fp, "      <Cells>\n");
214:       PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset);
215:       boffset += gpiece[r].nconn * sizeof(PetscInt) + sizeof(int);
216:       PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"offsets\"      NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset);
217:       boffset += gpiece[r].ncells * sizeof(PetscInt) + sizeof(int);
218:       PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"UInt8\" Name=\"types\"        NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset);
219:       boffset += gpiece[r].ncells * sizeof(unsigned char) + sizeof(int);
220:       PetscFPrintf(PETSC_COMM_SELF, fp, "      </Cells>\n");

222:       /*
223:        * Cell Data headers
224:        */
225:       PetscFPrintf(PETSC_COMM_SELF, fp, "      <CellData>\n");
226:       PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset);
227:       boffset += gpiece[r].ncells * sizeof(int) + sizeof(int);
228:       /* all the vectors */
229:       for (link = vtk->link; link; link = link->next) {
230:         Vec          X       = (Vec)link->vec;
231:         DM           dmX     = NULL;
232:         PetscInt     bs      = 1, nfields, field;
233:         const char  *vecname = "";
234:         PetscSection section;
235:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
236:         if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */
237:           PetscObjectGetName((PetscObject)X, &vecname);
238:         }
239:         VecGetDM(X, &dmX);
240:         if (!dmX) dmX = dm;
241:         PetscObjectQuery(link->vec, "section", (PetscObject *)&section);
242:         if (!section) DMGetLocalSection(dmX, &section);
243:         if (cEnd > cStart) PetscSectionGetDof(section, cStart, &bs);
244:         PetscSectionGetNumFields(section, &nfields);
245:         field = 0;
246:         if (link->field >= 0) {
247:           field   = link->field;
248:           nfields = field + 1;
249:         }
250:         for (i = 0; field < (nfields ? nfields : 1); field++) {
251:           PetscInt     fbs, j;
252:           PetscFV      fv = NULL;
253:           PetscObject  f;
254:           PetscClassId fClass;
255:           const char  *fieldname = NULL;
256:           char         buf[256];
257:           PetscBool    vector;
258:           if (nfields) { /* We have user-defined fields/components */
259:             PetscSectionGetFieldDof(section, cStart, field, &fbs);
260:             PetscSectionGetFieldName(section, field, &fieldname);
261:           } else fbs = bs; /* Say we have one field with 'bs' components */
262:           DMGetField(dmX, field, NULL, &f);
263:           PetscObjectGetClassId(f, &fClass);
264:           if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
265:           if (nfields && !fieldname) {
266:             PetscSNPrintf(buf, sizeof(buf), "CellField%" PetscInt_FMT, field);
267:             fieldname = buf;
268:           }
269:           vector = PETSC_FALSE;
270:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
271:             vector = PETSC_TRUE;
273:             for (j = 0; j < fbs; j++) {
274:               const char *compName = NULL;
275:               if (fv) {
276:                 PetscFVGetComponentName(fv, j, &compName);
277:                 if (compName) break;
278:               }
279:             }
280:             if (j < fbs) vector = PETSC_FALSE;
281:           }
282:           if (vector) {
283: #if defined(PETSC_USE_COMPLEX)
284:             PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset);
285:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + sizeof(int);
286:             PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset);
287:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + sizeof(int);
288: #else
289:             PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset);
290:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + sizeof(int);
291: #endif
292:             i += fbs;
293:           } else {
294:             for (j = 0; j < fbs; j++) {
295:               const char *compName = NULL;
296:               char        finalname[256];
297:               if (fv) PetscFVGetComponentName(fv, j, &compName);
298:               if (compName) {
299:                 PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName);
300:               } else if (fbs > 1) {
301:                 PetscSNPrintf(finalname, 255, "%s%s.%" PetscInt_FMT, vecname, fieldname, j);
302:               } else {
303:                 PetscSNPrintf(finalname, 255, "%s%s", vecname, fieldname);
304:               }
305: #if defined(PETSC_USE_COMPLEX)
306:               PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset);
307:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + sizeof(int);
308:               PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset);
309:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + sizeof(int);
310: #else
311:               PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset);
312:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + sizeof(int);
313: #endif
314:               i++;
315:             }
316:           }
317:         }
319:       }
320:       PetscFPrintf(PETSC_COMM_SELF, fp, "      </CellData>\n");

322:       /*
323:        * Point Data headers
324:        */
325:       PetscFPrintf(PETSC_COMM_SELF, fp, "      <PointData>\n");
326:       for (link = vtk->link; link; link = link->next) {
327:         Vec          X = (Vec)link->vec;
328:         DM           dmX;
329:         PetscInt     bs      = 1, nfields, field;
330:         const char  *vecname = "";
331:         PetscSection section;
332:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
333:         if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */
334:           PetscObjectGetName((PetscObject)X, &vecname);
335:         }
336:         VecGetDM(X, &dmX);
337:         if (!dmX) dmX = dm;
338:         PetscObjectQuery(link->vec, "section", (PetscObject *)&section);
339:         if (!section) DMGetLocalSection(dmX, &section);
340:         if (vEnd > vStart) PetscSectionGetDof(section, vStart, &bs);
341:         PetscSectionGetNumFields(section, &nfields);
342:         field = 0;
343:         if (link->field >= 0) {
344:           field   = link->field;
345:           nfields = field + 1;
346:         }
347:         for (i = 0; field < (nfields ? nfields : 1); field++) {
348:           PetscInt    fbs, j;
349:           const char *fieldname = NULL;
350:           char        buf[256];
351:           if (nfields) { /* We have user-defined fields/components */
352:             PetscSectionGetFieldDof(section, vStart, field, &fbs);
353:             PetscSectionGetFieldName(section, field, &fieldname);
354:           } else fbs = bs; /* Say we have one field with 'bs' components */
355:           if (nfields && !fieldname) {
356:             PetscSNPrintf(buf, sizeof(buf), "PointField%" PetscInt_FMT, field);
357:             fieldname = buf;
358:           }
359:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
361: #if defined(PETSC_USE_COMPLEX)
362:             PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset);
363:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + sizeof(int);
364:             PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset);
365:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + sizeof(int);
366: #else
367:             PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset);
368:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + sizeof(int);
369: #endif
370:           } else {
371:             for (j = 0; j < fbs; j++) {
372:               const char *compName = NULL;
373:               char        finalname[256];
374:               PetscSectionGetComponentName(section, field, j, &compName);
375:               PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName);
376: #if defined(PETSC_USE_COMPLEX)
377:               PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset);
378:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + sizeof(int);
379:               PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset);
380:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + sizeof(int);
381: #else
382:               PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset);
383:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + sizeof(int);
384: #endif
385:             }
386:           }
387:         }
388:       }
389:       PetscFPrintf(PETSC_COMM_SELF, fp, "      </PointData>\n");
390:       PetscFPrintf(PETSC_COMM_SELF, fp, "    </Piece>\n");
391:     }
392:   }

394:   PetscFPrintf(comm, fp, "  </UnstructuredGrid>\n");
395:   PetscFPrintf(comm, fp, "  <AppendedData encoding=\"raw\">\n");
396:   PetscFPrintf(comm, fp, "_");

398:   if (rank == 0) {
399:     PetscInt maxsize = 0;
400:     for (r = 0; r < size; r++) {
401:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nvertices * 3 * sizeof(PetscVTUReal)));
402:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].ncells * 3 * sizeof(PetscVTUReal)));
403:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nconn * sizeof(PetscVTKInt)));
404:     }
405:     PetscMalloc(maxsize, &buffer);
406:   }
407:   for (r = 0; r < size; r++) {
408:     if (r == rank) {
409:       PetscInt nsend;
410:       { /* Position */
411:         const PetscScalar *x, *cx = NULL;
412:         PetscVTUReal      *y = NULL;
413:         Vec                coords, cellCoords;
414:         PetscBool          copy;

416:         DMGetCoordinatesLocal(dm, &coords);
417:         VecGetArrayRead(coords, &x);
418:         DMGetCellCoordinatesLocal(dm, &cellCoords);
419:         if (cellCoords) VecGetArrayRead(cellCoords, &cx);
420: #if defined(PETSC_USE_COMPLEX)
421:         copy = PETSC_TRUE;
422: #else
423:         copy = (PetscBool)(dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal)));
424: #endif
425:         if (copy) {
426:           PetscMalloc1(piece.nvertices * 3, &y);
427:           if (localized) {
428:             PetscInt cnt;
429:             for (c = cStart, cnt = 0; c < cEnd; c++) {
430:               PetscInt off, dof;

432:               PetscSectionGetDof(cellCoordSection, c, &dof);
433:               if (!dof) {
434:                 PetscInt *closure = NULL;
435:                 PetscInt  closureSize;

437:                 DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
438:                 for (v = 0; v < closureSize * 2; v += 2) {
439:                   if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
440:                     PetscSectionGetOffset(coordSection, closure[v], &off);
441:                     if (dimEmbed != 3) {
442:                       y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]);
443:                       y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[off + 1]) : 0.0);
444:                       y[cnt * 3 + 2] = (PetscVTUReal)0.0;
445:                     } else {
446:                       y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]);
447:                       y[cnt * 3 + 1] = (PetscVTUReal)PetscRealPart(x[off + 1]);
448:                       y[cnt * 3 + 2] = (PetscVTUReal)PetscRealPart(x[off + 2]);
449:                     }
450:                     cnt++;
451:                   }
452:                 }
453:                 DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
454:               } else {
455:                 PetscSectionGetOffset(cellCoordSection, c, &off);
456:                 if (dimEmbed != 3) {
457:                   for (i = 0; i < dof / dimEmbed; i++) {
458:                     y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(cx[off + i * dimEmbed + 0]);
459:                     y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(cx[off + i * dimEmbed + 1]) : 0.0);
460:                     y[cnt * 3 + 2] = (PetscVTUReal)0.0;
461:                     cnt++;
462:                   }
463:                 } else {
464:                   for (i = 0; i < dof; i++) y[cnt * 3 + i] = (PetscVTUReal)PetscRealPart(cx[off + i]);
465:                   cnt += dof / dimEmbed;
466:                 }
467:               }
468:             }
470:           } else {
471:             for (i = 0; i < piece.nvertices; i++) {
472:               y[i * 3 + 0] = (PetscVTUReal)PetscRealPart(x[i * dimEmbed + 0]);
473:               y[i * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[i * dimEmbed + 1]) : 0.);
474:               y[i * 3 + 2] = (PetscVTUReal)((dimEmbed > 2) ? PetscRealPart(x[i * dimEmbed + 2]) : 0.);
475:             }
476:           }
477:         }
478:         nsend = piece.nvertices * 3;
479:         TransferWrite(comm, viewer, fp, r, 0, copy ? (const void *)y : (const void *)x, buffer, nsend, MPIU_VTUREAL, tag);
480:         PetscFree(y);
481:         VecRestoreArrayRead(coords, &x);
482:         if (cellCoords) VecRestoreArrayRead(cellCoords, &cx);
483:       }
484:       { /* Connectivity, offsets, types */
485:         PetscVTKInt  *connectivity = NULL, *offsets = NULL;
486:         PetscVTKType *types = NULL;
487:         DMPlexGetVTKConnectivity(dm, localized, &piece, &connectivity, &offsets, &types);
488:         TransferWrite(comm, viewer, fp, r, 0, connectivity, buffer, piece.nconn, MPI_INT, tag);
489:         TransferWrite(comm, viewer, fp, r, 0, offsets, buffer, piece.ncells, MPI_INT, tag);
490:         TransferWrite(comm, viewer, fp, r, 0, types, buffer, piece.ncells, MPI_CHAR, tag);
491:         PetscFree3(connectivity, offsets, types);
492:       }
493:       { /* Owners (cell data) */
494:         PetscVTKInt *owners;
495:         PetscMalloc1(piece.ncells, &owners);
496:         for (i = 0; i < piece.ncells; i++) owners[i] = rank;
497:         TransferWrite(comm, viewer, fp, r, 0, owners, buffer, piece.ncells, MPI_INT, tag);
498:         PetscFree(owners);
499:       }
500:       /* Cell data */
501:       for (link = vtk->link; link; link = link->next) {
502:         Vec                X = (Vec)link->vec;
503:         DM                 dmX;
504:         const PetscScalar *x;
505:         PetscVTUReal      *y;
506:         PetscInt           bs      = 1, nfields, field;
507:         PetscSection       section = NULL;

509:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
510:         VecGetDM(X, &dmX);
511:         if (!dmX) dmX = dm;
512:         PetscObjectQuery(link->vec, "section", (PetscObject *)&section);
513:         if (!section) DMGetLocalSection(dmX, &section);
514:         if (cEnd > cStart) PetscSectionGetDof(section, cStart, &bs);
515:         PetscSectionGetNumFields(section, &nfields);
516:         field = 0;
517:         if (link->field >= 0) {
518:           field   = link->field;
519:           nfields = field + 1;
520:         }
521:         VecGetArrayRead(X, &x);
522:         PetscMalloc1(piece.ncells * 3, &y);
523:         for (i = 0; field < (nfields ? nfields : 1); field++) {
524:           PetscInt     fbs, j;
525:           PetscFV      fv = NULL;
526:           PetscObject  f;
527:           PetscClassId fClass;
528:           PetscBool    vector;
529:           if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
530:             PetscSectionGetFieldDof(section, cStart, field, &fbs);
531:           } else fbs = bs; /* Say we have one field with 'bs' components */
532:           DMGetField(dmX, field, NULL, &f);
533:           PetscObjectGetClassId(f, &fClass);
534:           if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
535:           vector = PETSC_FALSE;
536:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
537:             vector = PETSC_TRUE;
538:             for (j = 0; j < fbs; j++) {
539:               const char *compName = NULL;
540:               if (fv) {
541:                 PetscFVGetComponentName(fv, j, &compName);
542:                 if (compName) break;
543:               }
544:             }
545:             if (j < fbs) vector = PETSC_FALSE;
546:           }
547:           if (vector) {
548:             PetscInt cnt, l;
549:             for (l = 0; l < loops_per_scalar; l++) {
550:               for (c = cStart, cnt = 0; c < cEnd; c++) {
551:                 const PetscScalar *xpoint;
552:                 PetscInt           off, j;

554:                 if (hasLabel) { /* Ignore some cells */
555:                   PetscInt value;
556:                   DMGetLabelValue(dmX, "vtk", c, &value);
557:                   if (value != 1) continue;
558:                 }
559:                 if (nfields) {
560:                   PetscSectionGetFieldOffset(section, c, field, &off);
561:                 } else {
562:                   PetscSectionGetOffset(section, c, &off);
563:                 }
564:                 xpoint = &x[off];
565:                 for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
566:                 for (; j < 3; j++) y[cnt++] = 0.;
567:               }
569:               TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells * 3, MPIU_VTUREAL, tag);
570:             }
571:           } else {
572:             for (i = 0; i < fbs; i++) {
573:               PetscInt cnt, l;
574:               for (l = 0; l < loops_per_scalar; l++) {
575:                 for (c = cStart, cnt = 0; c < cEnd; c++) {
576:                   const PetscScalar *xpoint;
577:                   PetscInt           off;

579:                   if (hasLabel) { /* Ignore some cells */
580:                     PetscInt value;
581:                     DMGetLabelValue(dmX, "vtk", c, &value);
582:                     if (value != 1) continue;
583:                   }
584:                   if (nfields) {
585:                     PetscSectionGetFieldOffset(section, c, field, &off);
586:                   } else {
587:                     PetscSectionGetOffset(section, c, &off);
588:                   }
589:                   xpoint   = &x[off];
590:                   y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
591:                 }
593:                 TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells, MPIU_VTUREAL, tag);
594:               }
595:             }
596:           }
597:         }
598:         PetscFree(y);
599:         VecRestoreArrayRead(X, &x);
600:       }
601:       /* point data */
602:       for (link = vtk->link; link; link = link->next) {
603:         Vec                X = (Vec)link->vec;
604:         DM                 dmX;
605:         const PetscScalar *x;
606:         PetscVTUReal      *y;
607:         PetscInt           bs      = 1, nfields, field;
608:         PetscSection       section = NULL;

610:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
611:         VecGetDM(X, &dmX);
612:         if (!dmX) dmX = dm;
613:         PetscObjectQuery(link->vec, "section", (PetscObject *)&section);
614:         if (!section) DMGetLocalSection(dmX, &section);
615:         if (vEnd > vStart) PetscSectionGetDof(section, vStart, &bs);
616:         PetscSectionGetNumFields(section, &nfields);
617:         field = 0;
618:         if (link->field >= 0) {
619:           field   = link->field;
620:           nfields = field + 1;
621:         }
622:         VecGetArrayRead(X, &x);
623:         PetscMalloc1(piece.nvertices * 3, &y);
624:         for (i = 0; field < (nfields ? nfields : 1); field++) {
625:           PetscInt fbs, j;
626:           if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
627:             PetscSectionGetFieldDof(section, vStart, field, &fbs);
628:           } else fbs = bs; /* Say we have one field with 'bs' components */
629:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
630:             PetscInt cnt, l;
631:             for (l = 0; l < loops_per_scalar; l++) {
632:               if (!localized) {
633:                 for (v = vStart, cnt = 0; v < vEnd; v++) {
634:                   PetscInt           off;
635:                   const PetscScalar *xpoint;

637:                   if (nfields) {
638:                     PetscSectionGetFieldOffset(section, v, field, &off);
639:                   } else {
640:                     PetscSectionGetOffset(section, v, &off);
641:                   }
642:                   xpoint = &x[off];
643:                   for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
644:                   for (; j < 3; j++) y[cnt++] = 0.;
645:                 }
646:               } else {
647:                 for (c = cStart, cnt = 0; c < cEnd; c++) {
648:                   PetscInt *closure = NULL;
649:                   PetscInt  closureSize, off;

651:                   DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
652:                   for (v = 0, off = 0; v < closureSize * 2; v += 2) {
653:                     if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
654:                       PetscInt           voff;
655:                       const PetscScalar *xpoint;

657:                       if (nfields) {
658:                         PetscSectionGetFieldOffset(section, closure[v], field, &voff);
659:                       } else {
660:                         PetscSectionGetOffset(section, closure[v], &voff);
661:                       }
662:                       xpoint = &x[voff];
663:                       for (j = 0; j < fbs; j++) y[cnt + off++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
664:                       for (; j < 3; j++) y[cnt + off++] = 0.;
665:                     }
666:                   }
667:                   cnt += off;
668:                   DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
669:                 }
670:               }
672:               TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices * 3, MPIU_VTUREAL, tag);
673:             }
674:           } else {
675:             for (i = 0; i < fbs; i++) {
676:               PetscInt cnt, l;
677:               for (l = 0; l < loops_per_scalar; l++) {
678:                 if (!localized) {
679:                   for (v = vStart, cnt = 0; v < vEnd; v++) {
680:                     PetscInt           off;
681:                     const PetscScalar *xpoint;

683:                     if (nfields) {
684:                       PetscSectionGetFieldOffset(section, v, field, &off);
685:                     } else {
686:                       PetscSectionGetOffset(section, v, &off);
687:                     }
688:                     xpoint   = &x[off];
689:                     y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
690:                   }
691:                 } else {
692:                   for (c = cStart, cnt = 0; c < cEnd; c++) {
693:                     PetscInt *closure = NULL;
694:                     PetscInt  closureSize, off;

696:                     DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
697:                     for (v = 0, off = 0; v < closureSize * 2; v += 2) {
698:                       if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
699:                         PetscInt           voff;
700:                         const PetscScalar *xpoint;

702:                         if (nfields) {
703:                           PetscSectionGetFieldOffset(section, closure[v], field, &voff);
704:                         } else {
705:                           PetscSectionGetOffset(section, closure[v], &voff);
706:                         }
707:                         xpoint         = &x[voff];
708:                         y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
709:                       }
710:                     }
711:                     cnt += off;
712:                     DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
713:                   }
714:                 }
716:                 TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices, MPIU_VTUREAL, tag);
717:               }
718:             }
719:           }
720:         }
721:         PetscFree(y);
722:         VecRestoreArrayRead(X, &x);
723:       }
724:     } else if (rank == 0) {
725:       PetscInt l;

727:       TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag); /* positions */
728:       TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nconn, MPI_INT, tag);              /* connectivity */
729:       TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag);             /* offsets */
730:       TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_CHAR, tag);            /* types */
731:       TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag);             /* owner rank (cells) */
732:       /* all cell data */
733:       for (link = vtk->link; link; link = link->next) {
734:         Vec          X  = (Vec)link->vec;
735:         PetscInt     bs = 1, nfields, field;
736:         DM           dmX;
737:         PetscSection section = NULL;

739:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
740:         VecGetDM(X, &dmX);
741:         if (!dmX) dmX = dm;
742:         PetscObjectQuery(link->vec, "section", (PetscObject *)&section);
743:         if (!section) DMGetLocalSection(dmX, &section);
744:         if (cEnd > cStart) PetscSectionGetDof(section, cStart, &bs);
745:         PetscSectionGetNumFields(section, &nfields);
746:         field = 0;
747:         if (link->field >= 0) {
748:           field   = link->field;
749:           nfields = field + 1;
750:         }
751:         for (i = 0; field < (nfields ? nfields : 1); field++) {
752:           PetscInt     fbs, j;
753:           PetscFV      fv = NULL;
754:           PetscObject  f;
755:           PetscClassId fClass;
756:           PetscBool    vector;
757:           if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
758:             PetscSectionGetFieldDof(section, cStart, field, &fbs);
759:           } else fbs = bs; /* Say we have one field with 'bs' components */
760:           DMGetField(dmX, field, NULL, &f);
761:           PetscObjectGetClassId(f, &fClass);
762:           if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
763:           vector = PETSC_FALSE;
764:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
765:             vector = PETSC_TRUE;
766:             for (j = 0; j < fbs; j++) {
767:               const char *compName = NULL;
768:               if (fv) {
769:                 PetscFVGetComponentName(fv, j, &compName);
770:                 if (compName) break;
771:               }
772:             }
773:             if (j < fbs) vector = PETSC_FALSE;
774:           }
775:           if (vector) {
776:             for (l = 0; l < loops_per_scalar; l++) TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells * 3, MPIU_VTUREAL, tag);
777:           } else {
778:             for (i = 0; i < fbs; i++) {
779:               for (l = 0; l < loops_per_scalar; l++) TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPIU_VTUREAL, tag);
780:             }
781:           }
782:         }
783:       }
784:       /* all point data */
785:       for (link = vtk->link; link; link = link->next) {
786:         Vec          X = (Vec)link->vec;
787:         DM           dmX;
788:         PetscInt     bs      = 1, nfields, field;
789:         PetscSection section = NULL;

791:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
792:         VecGetDM(X, &dmX);
793:         if (!dmX) dmX = dm;
794:         PetscObjectQuery(link->vec, "section", (PetscObject *)&section);
795:         if (!section) DMGetLocalSection(dmX, &section);
796:         if (vEnd > vStart) PetscSectionGetDof(section, vStart, &bs);
797:         PetscSectionGetNumFields(section, &nfields);
798:         field = 0;
799:         if (link->field >= 0) {
800:           field   = link->field;
801:           nfields = field + 1;
802:         }
803:         for (i = 0; field < (nfields ? nfields : 1); field++) {
804:           PetscInt fbs;
805:           if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
806:             PetscSectionGetFieldDof(section, vStart, field, &fbs);
807:           } else fbs = bs; /* Say we have one field with 'bs' components */
808:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
809:             for (l = 0; l < loops_per_scalar; l++) TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag);
810:           } else {
811:             for (i = 0; i < fbs; i++) {
812:               for (l = 0; l < loops_per_scalar; l++) TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices, MPIU_VTUREAL, tag);
813:             }
814:           }
815:         }
816:       }
817:     }
818:   }
819:   PetscFree(gpiece);
820:   PetscFree(buffer);
821:   PetscFPrintf(comm, fp, "\n  </AppendedData>\n");
822:   PetscFPrintf(comm, fp, "</VTKFile>\n");
823:   PetscFClose(comm, fp);
824:   return 0;
825: }