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 *)§ion);
242: if (!section) DMGetLocalSection(dmX, §ion);
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 *)§ion);
339: if (!section) DMGetLocalSection(dmX, §ion);
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 *)§ion);
513: if (!section) DMGetLocalSection(dmX, §ion);
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 *)§ion);
614: if (!section) DMGetLocalSection(dmX, §ion);
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 *)§ion);
743: if (!section) DMGetLocalSection(dmX, §ion);
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 *)§ion);
795: if (!section) DMGetLocalSection(dmX, §ion);
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: }