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: PetscFunctionBegin;
31: PetscCallMPI(MPI_Comm_rank(comm, &rank));
32: if (rank == srank && rank != root) {
33: PetscCallMPI(MPI_Send((void *)send, count, mpidatatype, root, tag, comm));
34: } else if (rank == root) {
35: const void *buffer;
36: if (root == srank) { /* self */
37: buffer = send;
38: } else {
39: MPI_Status status;
40: PetscMPIInt nrecv;
41: PetscCallMPI(MPI_Recv(recv, count, mpidatatype, srank, tag, comm, &status));
42: PetscCallMPI(MPI_Get_count(&status, mpidatatype, &nrecv));
43: PetscCheck(count == nrecv, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch");
44: buffer = recv;
45: }
46: PetscCall(PetscViewerVTKFWrite(viewer, fp, buffer, count, mpidatatype));
47: }
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode DMPlexGetVTKConnectivity(DM dm, PetscBool localized, PieceInfo *piece, PetscVTKInt **oconn, PetscVTKInt **ooffsets, PetscVTKType **otypes)
52: {
53: PetscSection coordSection, cellCoordSection;
54: PetscVTKInt *conn, *offsets;
55: PetscVTKType *types;
56: PetscInt dim, vStart, vEnd, cStart, cEnd, pStart, pEnd, cellHeight, numLabelCells, hasLabel, c, v, countcell, countconn;
58: PetscFunctionBegin;
59: PetscCall(PetscMalloc3(piece->nconn, &conn, piece->ncells, &offsets, piece->ncells, &types));
60: PetscCall(DMGetCoordinateSection(dm, &coordSection));
61: PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection));
62: PetscCall(DMGetDimension(dm, &dim));
63: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
64: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
65: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
66: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
67: PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells));
68: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
70: countcell = 0;
71: countconn = 0;
72: for (c = cStart; c < cEnd; ++c) {
73: PetscInt nverts, dof = 0, celltype, startoffset, nC = 0;
75: if (hasLabel) {
76: PetscInt value;
78: PetscCall(DMGetLabelValue(dm, "vtk", c, &value));
79: if (value != 1) continue;
80: }
81: startoffset = countconn;
82: if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
83: if (!dof) {
84: PetscInt *closure = NULL;
85: PetscInt closureSize;
87: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
88: for (v = 0; v < closureSize * 2; v += 2) {
89: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
90: if (!localized) conn[countconn++] = closure[v] - vStart;
91: else conn[countconn++] = startoffset + nC;
92: ++nC;
93: }
94: }
95: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
96: } else {
97: for (nC = 0; nC < dof / dim; nC++) conn[countconn++] = startoffset + nC;
98: }
100: {
101: PetscInt n = PetscMin(nC, 8), s = countconn - nC, i, cone[8];
102: for (i = 0; i < n; ++i) cone[i] = conn[s + i];
103: PetscCall(DMPlexReorderCell(dm, c, cone));
104: for (i = 0; i < n; ++i) conn[s + i] = (int)cone[i];
105: }
107: offsets[countcell] = countconn;
109: nverts = countconn - startoffset;
110: PetscCall(DMPlexVTKGetCellType_Internal(dm, dim, nverts, &celltype));
112: types[countcell] = celltype;
113: countcell++;
114: }
115: PetscCheck(countcell == piece->ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent cell count");
116: PetscCheck(countconn == piece->nconn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent connectivity count");
117: *oconn = conn;
118: *ooffsets = offsets;
119: *otypes = types;
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: /*
124: Write all fields that have been provided to the viewer
125: Multi-block XML format with binary appended data.
126: */
127: PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm, PetscViewer viewer)
128: {
129: MPI_Comm comm;
130: PetscSection coordSection, cellCoordSection;
131: PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data;
132: PetscViewerVTKObjectLink link;
133: FILE *fp;
134: PetscMPIInt rank, size, tag;
135: PetscInt dimEmbed, cellHeight, cStart, cEnd, vStart, vEnd, numLabelCells, hasLabel, c, v, r, i;
136: PetscBool localized;
137: PieceInfo piece, *gpiece = NULL;
138: void *buffer = NULL;
139: const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
140: PetscInt loops_per_scalar;
142: PetscFunctionBegin;
143: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
144: #if defined(PETSC_USE_COMPLEX)
145: loops_per_scalar = 2;
146: #else
147: loops_per_scalar = 1;
148: #endif
149: PetscCallMPI(MPI_Comm_size(comm, &size));
150: PetscCallMPI(MPI_Comm_rank(comm, &rank));
151: PetscCall(PetscCommGetNewTag(comm, &tag));
153: PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
154: PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n"));
155: PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order));
156: PetscCall(PetscFPrintf(comm, fp, " <UnstructuredGrid>\n"));
158: PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
159: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
160: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
161: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
162: PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells));
163: PetscCall(DMGetCoordinatesLocalized(dm, &localized));
164: PetscCall(DMGetCoordinateSection(dm, &coordSection));
165: PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection));
167: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
168: piece.nvertices = 0;
169: piece.ncells = 0;
170: piece.nconn = 0;
171: if (!localized) piece.nvertices = vEnd - vStart;
172: for (c = cStart; c < cEnd; ++c) {
173: PetscInt dof = 0;
174: if (hasLabel) {
175: PetscInt value;
177: PetscCall(DMGetLabelValue(dm, "vtk", c, &value));
178: if (value != 1) continue;
179: }
180: if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
181: if (!dof) {
182: PetscInt *closure = NULL;
183: PetscInt closureSize;
185: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
186: for (v = 0; v < closureSize * 2; v += 2) {
187: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
188: piece.nconn++;
189: if (localized) piece.nvertices++;
190: }
191: }
192: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
193: } else {
194: piece.nvertices += dof / dimEmbed;
195: piece.nconn += dof / dimEmbed;
196: }
197: piece.ncells++;
198: }
199: if (rank == 0) PetscCall(PetscMalloc1(size, &gpiece));
200: PetscCallMPI(MPI_Gather((PetscInt *)&piece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, (PetscInt *)gpiece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, 0, comm));
202: /*
203: * Write file header
204: */
205: if (rank == 0) {
206: PetscInt64 boffset = 0;
208: for (r = 0; r < size; r++) {
209: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Piece NumberOfPoints=\"%" PetscInt_FMT "\" NumberOfCells=\"%" PetscInt_FMT "\">\n", gpiece[r].nvertices, gpiece[r].ncells));
210: /* Coordinate positions */
211: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Points>\n"));
212: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
213: boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
214: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Points>\n"));
215: /* Cell connectivity */
216: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Cells>\n"));
217: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
218: boffset += gpiece[r].nconn * sizeof(PetscVTKInt) + (gpiece[r].nconn ? sizeof(PetscInt64) : 0);
219: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
220: boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
221: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
222: boffset += gpiece[r].ncells * sizeof(unsigned char) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
223: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Cells>\n"));
225: /*
226: * Cell Data headers
227: */
228: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <CellData>\n"));
229: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
230: boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
231: /* all the vectors */
232: for (link = vtk->link; link; link = link->next) {
233: Vec X = (Vec)link->vec;
234: DM dmX = NULL;
235: PetscInt bs = 1, nfields, field;
236: const char *vecname = "";
237: PetscSection section;
238: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
239: 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. */
240: PetscCall(PetscObjectGetName((PetscObject)X, &vecname));
241: }
242: PetscCall(VecGetDM(X, &dmX));
243: if (!dmX) dmX = dm;
244: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
245: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
246: if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
247: PetscCall(PetscSectionGetNumFields(section, &nfields));
248: field = 0;
249: if (link->field >= 0) {
250: field = link->field;
251: nfields = field + 1;
252: }
253: for (i = 0; field < (nfields ? nfields : 1); field++) {
254: PetscInt fbs, j;
255: PetscFV fv = NULL;
256: PetscObject f;
257: PetscClassId fClass;
258: const char *fieldname = NULL;
259: char buf[256];
260: PetscBool vector;
261: if (nfields) { /* We have user-defined fields/components */
262: PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
263: PetscCall(PetscSectionGetFieldName(section, field, &fieldname));
264: } else fbs = bs; /* Say we have one field with 'bs' components */
265: PetscCall(DMGetField(dmX, field, NULL, &f));
266: PetscCall(PetscObjectGetClassId(f, &fClass));
267: if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
268: if (nfields && !fieldname) {
269: PetscCall(PetscSNPrintf(buf, sizeof(buf), "CellField%" PetscInt_FMT, field));
270: fieldname = buf;
271: }
272: vector = PETSC_FALSE;
273: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
274: vector = PETSC_TRUE;
275: PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Cell vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs);
276: for (j = 0; j < fbs; j++) {
277: const char *compName = NULL;
278: if (fv) {
279: PetscCall(PetscFVGetComponentName(fv, j, &compName));
280: if (compName) break;
281: }
282: }
283: if (j < fbs) vector = PETSC_FALSE;
284: }
285: if (vector) {
286: #if defined(PETSC_USE_COMPLEX)
287: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
288: boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
289: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
290: boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
291: #else
292: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
293: boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
294: #endif
295: i += fbs;
296: } else {
297: for (j = 0; j < fbs; j++) {
298: const char *compName = NULL;
299: char finalname[256];
300: if (fv) PetscCall(PetscFVGetComponentName(fv, j, &compName));
301: if (compName) {
302: PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName));
303: } else if (fbs > 1) {
304: PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%" PetscInt_FMT, vecname, fieldname, j));
305: } else {
306: PetscCall(PetscSNPrintf(finalname, 255, "%s%s", vecname, fieldname));
307: }
308: #if defined(PETSC_USE_COMPLEX)
309: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
310: boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
311: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
312: boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
313: #else
314: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
315: boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
316: #endif
317: i++;
318: }
319: }
320: }
321: //PetscCheck(i == bs,comm,PETSC_ERR_PLIB,"Total number of field components %" PetscInt_FMT " != block size %" PetscInt_FMT,i,bs);
322: }
323: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </CellData>\n"));
325: /*
326: * Point Data headers
327: */
328: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <PointData>\n"));
329: for (link = vtk->link; link; link = link->next) {
330: Vec X = (Vec)link->vec;
331: DM dmX;
332: PetscInt bs = 1, nfields, field;
333: const char *vecname = "";
334: PetscSection section;
335: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
336: 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. */
337: PetscCall(PetscObjectGetName((PetscObject)X, &vecname));
338: }
339: PetscCall(VecGetDM(X, &dmX));
340: if (!dmX) dmX = dm;
341: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
342: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
343: if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
344: PetscCall(PetscSectionGetNumFields(section, &nfields));
345: field = 0;
346: if (link->field >= 0) {
347: field = link->field;
348: nfields = field + 1;
349: }
350: for (i = 0; field < (nfields ? nfields : 1); field++) {
351: PetscInt fbs, j;
352: const char *fieldname = NULL;
353: char buf[256];
354: if (nfields) { /* We have user-defined fields/components */
355: PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
356: PetscCall(PetscSectionGetFieldName(section, field, &fieldname));
357: } else fbs = bs; /* Say we have one field with 'bs' components */
358: if (nfields && !fieldname) {
359: PetscCall(PetscSNPrintf(buf, sizeof(buf), "PointField%" PetscInt_FMT, field));
360: fieldname = buf;
361: }
362: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
363: PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Point vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs);
364: #if defined(PETSC_USE_COMPLEX)
365: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
366: boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
367: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
368: boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
369: #else
370: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
371: boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
372: #endif
373: } else {
374: for (j = 0; j < fbs; j++) {
375: const char *compName = NULL;
376: char finalname[256];
377: PetscCall(PetscSectionGetComponentName(section, field, j, &compName));
378: PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName));
379: #if defined(PETSC_USE_COMPLEX)
380: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
381: boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
382: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
383: boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
384: #else
385: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
386: boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
387: #endif
388: }
389: }
390: }
391: }
392: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </PointData>\n"));
393: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Piece>\n"));
394: }
395: }
397: PetscCall(PetscFPrintf(comm, fp, " </UnstructuredGrid>\n"));
398: PetscCall(PetscFPrintf(comm, fp, " <AppendedData encoding=\"raw\">\n"));
399: PetscCall(PetscFPrintf(comm, fp, "_"));
401: if (rank == 0) {
402: PetscInt maxsize = 0;
403: for (r = 0; r < size; r++) {
404: maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nvertices * 3 * sizeof(PetscVTUReal)));
405: maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].ncells * 3 * sizeof(PetscVTUReal)));
406: maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nconn * sizeof(PetscVTKInt)));
407: }
408: PetscCall(PetscMalloc(maxsize, &buffer));
409: }
410: for (r = 0; r < size; r++) {
411: if (r == rank) {
412: PetscInt nsend;
413: { /* Position */
414: const PetscScalar *x, *cx = NULL;
415: PetscVTUReal *y = NULL;
416: Vec coords, cellCoords;
417: PetscBool copy;
419: PetscCall(DMGetCoordinatesLocal(dm, &coords));
420: PetscCall(VecGetArrayRead(coords, &x));
421: PetscCall(DMGetCellCoordinatesLocal(dm, &cellCoords));
422: if (cellCoords) PetscCall(VecGetArrayRead(cellCoords, &cx));
423: #if defined(PETSC_USE_COMPLEX)
424: copy = PETSC_TRUE;
425: #else
426: copy = (PetscBool)(dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal)));
427: #endif
428: if (copy) {
429: PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
430: if (localized) {
431: PetscInt cnt;
432: for (c = cStart, cnt = 0; c < cEnd; c++) {
433: PetscInt off, dof;
435: PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
436: if (!dof) {
437: PetscInt *closure = NULL;
438: PetscInt closureSize;
440: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
441: for (v = 0; v < closureSize * 2; v += 2) {
442: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
443: PetscCall(PetscSectionGetOffset(coordSection, closure[v], &off));
444: if (dimEmbed != 3) {
445: y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]);
446: y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[off + 1]) : 0.0);
447: y[cnt * 3 + 2] = (PetscVTUReal)0.0;
448: } else {
449: y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]);
450: y[cnt * 3 + 1] = (PetscVTUReal)PetscRealPart(x[off + 1]);
451: y[cnt * 3 + 2] = (PetscVTUReal)PetscRealPart(x[off + 2]);
452: }
453: cnt++;
454: }
455: }
456: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
457: } else {
458: PetscCall(PetscSectionGetOffset(cellCoordSection, c, &off));
459: if (dimEmbed != 3) {
460: for (i = 0; i < dof / dimEmbed; i++) {
461: y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(cx[off + i * dimEmbed + 0]);
462: y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(cx[off + i * dimEmbed + 1]) : 0.0);
463: y[cnt * 3 + 2] = (PetscVTUReal)0.0;
464: cnt++;
465: }
466: } else {
467: for (i = 0; i < dof; i++) y[cnt * 3 + i] = (PetscVTUReal)PetscRealPart(cx[off + i]);
468: cnt += dof / dimEmbed;
469: }
470: }
471: }
472: PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
473: } else {
474: for (i = 0; i < piece.nvertices; i++) {
475: y[i * 3 + 0] = (PetscVTUReal)PetscRealPart(x[i * dimEmbed + 0]);
476: y[i * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[i * dimEmbed + 1]) : 0.);
477: y[i * 3 + 2] = (PetscVTUReal)((dimEmbed > 2) ? PetscRealPart(x[i * dimEmbed + 2]) : 0.);
478: }
479: }
480: }
481: nsend = piece.nvertices * 3;
482: PetscCall(TransferWrite(comm, viewer, fp, r, 0, copy ? (const void *)y : (const void *)x, buffer, nsend, MPIU_VTUREAL, tag));
483: PetscCall(PetscFree(y));
484: PetscCall(VecRestoreArrayRead(coords, &x));
485: if (cellCoords) PetscCall(VecRestoreArrayRead(cellCoords, &cx));
486: }
487: { /* Connectivity, offsets, types */
488: PetscVTKInt *connectivity = NULL, *offsets = NULL;
489: PetscVTKType *types = NULL;
490: PetscCall(DMPlexGetVTKConnectivity(dm, localized, &piece, &connectivity, &offsets, &types));
491: PetscCall(TransferWrite(comm, viewer, fp, r, 0, connectivity, buffer, piece.nconn, MPI_INT, tag));
492: PetscCall(TransferWrite(comm, viewer, fp, r, 0, offsets, buffer, piece.ncells, MPI_INT, tag));
493: PetscCall(TransferWrite(comm, viewer, fp, r, 0, types, buffer, piece.ncells, MPI_CHAR, tag));
494: PetscCall(PetscFree3(connectivity, offsets, types));
495: }
496: { /* Owners (cell data) */
497: PetscVTKInt *owners;
498: PetscCall(PetscMalloc1(piece.ncells, &owners));
499: for (i = 0; i < piece.ncells; i++) owners[i] = rank;
500: PetscCall(TransferWrite(comm, viewer, fp, r, 0, owners, buffer, piece.ncells, MPI_INT, tag));
501: PetscCall(PetscFree(owners));
502: }
503: /* Cell data */
504: for (link = vtk->link; link; link = link->next) {
505: Vec X = (Vec)link->vec;
506: DM dmX;
507: const PetscScalar *x;
508: PetscVTUReal *y;
509: PetscInt bs = 1, nfields, field;
510: PetscSection section = NULL;
512: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
513: PetscCall(VecGetDM(X, &dmX));
514: if (!dmX) dmX = dm;
515: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
516: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
517: if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
518: PetscCall(PetscSectionGetNumFields(section, &nfields));
519: field = 0;
520: if (link->field >= 0) {
521: field = link->field;
522: nfields = field + 1;
523: }
524: PetscCall(VecGetArrayRead(X, &x));
525: PetscCall(PetscMalloc1(piece.ncells * 3, &y));
526: for (i = 0; field < (nfields ? nfields : 1); field++) {
527: PetscInt fbs, j;
528: PetscFV fv = NULL;
529: PetscObject f;
530: PetscClassId fClass;
531: PetscBool vector;
532: if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
533: PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
534: } else fbs = bs; /* Say we have one field with 'bs' components */
535: PetscCall(DMGetField(dmX, field, NULL, &f));
536: PetscCall(PetscObjectGetClassId(f, &fClass));
537: if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
538: vector = PETSC_FALSE;
539: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
540: vector = PETSC_TRUE;
541: for (j = 0; j < fbs; j++) {
542: const char *compName = NULL;
543: if (fv) {
544: PetscCall(PetscFVGetComponentName(fv, j, &compName));
545: if (compName) break;
546: }
547: }
548: if (j < fbs) vector = PETSC_FALSE;
549: }
550: if (vector) {
551: PetscInt cnt, l;
552: for (l = 0; l < loops_per_scalar; l++) {
553: for (c = cStart, cnt = 0; c < cEnd; c++) {
554: const PetscScalar *xpoint;
555: PetscInt off, j;
557: if (hasLabel) { /* Ignore some cells */
558: PetscInt value;
559: PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
560: if (value != 1) continue;
561: }
562: if (nfields) {
563: PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
564: } else {
565: PetscCall(PetscSectionGetOffset(section, c, &off));
566: }
567: xpoint = &x[off];
568: for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
569: for (; j < 3; j++) y[cnt++] = 0.;
570: }
571: PetscCheck(cnt == piece.ncells * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
572: PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells * 3, MPIU_VTUREAL, tag));
573: }
574: } else {
575: for (i = 0; i < fbs; i++) {
576: PetscInt cnt, l;
577: for (l = 0; l < loops_per_scalar; l++) {
578: for (c = cStart, cnt = 0; c < cEnd; c++) {
579: const PetscScalar *xpoint;
580: PetscInt off;
582: if (hasLabel) { /* Ignore some cells */
583: PetscInt value;
584: PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
585: if (value != 1) continue;
586: }
587: if (nfields) {
588: PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
589: } else {
590: PetscCall(PetscSectionGetOffset(section, c, &off));
591: }
592: xpoint = &x[off];
593: y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
594: }
595: PetscCheck(cnt == piece.ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
596: PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells, MPIU_VTUREAL, tag));
597: }
598: }
599: }
600: }
601: PetscCall(PetscFree(y));
602: PetscCall(VecRestoreArrayRead(X, &x));
603: }
604: /* point data */
605: for (link = vtk->link; link; link = link->next) {
606: Vec X = (Vec)link->vec;
607: DM dmX;
608: const PetscScalar *x;
609: PetscVTUReal *y;
610: PetscInt bs = 1, nfields, field;
611: PetscSection section = NULL;
613: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
614: PetscCall(VecGetDM(X, &dmX));
615: if (!dmX) dmX = dm;
616: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
617: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
618: if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
619: PetscCall(PetscSectionGetNumFields(section, &nfields));
620: field = 0;
621: if (link->field >= 0) {
622: field = link->field;
623: nfields = field + 1;
624: }
625: PetscCall(VecGetArrayRead(X, &x));
626: PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
627: for (i = 0; field < (nfields ? nfields : 1); field++) {
628: PetscInt fbs, j;
629: if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
630: PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
631: } else fbs = bs; /* Say we have one field with 'bs' components */
632: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
633: PetscInt cnt, l;
634: for (l = 0; l < loops_per_scalar; l++) {
635: if (!localized) {
636: for (v = vStart, cnt = 0; v < vEnd; v++) {
637: PetscInt off;
638: const PetscScalar *xpoint;
640: if (nfields) {
641: PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
642: } else {
643: PetscCall(PetscSectionGetOffset(section, v, &off));
644: }
645: xpoint = &x[off];
646: for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
647: for (; j < 3; j++) y[cnt++] = 0.;
648: }
649: } else {
650: for (c = cStart, cnt = 0; c < cEnd; c++) {
651: PetscInt *closure = NULL;
652: PetscInt closureSize, off;
654: PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
655: for (v = 0, off = 0; v < closureSize * 2; v += 2) {
656: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
657: PetscInt voff;
658: const PetscScalar *xpoint;
660: if (nfields) {
661: PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
662: } else {
663: PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
664: }
665: xpoint = &x[voff];
666: for (j = 0; j < fbs; j++) y[cnt + off++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
667: for (; j < 3; j++) y[cnt + off++] = 0.;
668: }
669: }
670: cnt += off;
671: PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
672: }
673: }
674: PetscCheck(cnt == piece.nvertices * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
675: PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices * 3, MPIU_VTUREAL, tag));
676: }
677: } else {
678: for (i = 0; i < fbs; i++) {
679: PetscInt cnt, l;
680: for (l = 0; l < loops_per_scalar; l++) {
681: if (!localized) {
682: for (v = vStart, cnt = 0; v < vEnd; v++) {
683: PetscInt off;
684: const PetscScalar *xpoint;
686: if (nfields) {
687: PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
688: } else {
689: PetscCall(PetscSectionGetOffset(section, v, &off));
690: }
691: xpoint = &x[off];
692: y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
693: }
694: } else {
695: for (c = cStart, cnt = 0; c < cEnd; c++) {
696: PetscInt *closure = NULL;
697: PetscInt closureSize, off;
699: PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
700: for (v = 0, off = 0; v < closureSize * 2; v += 2) {
701: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
702: PetscInt voff;
703: const PetscScalar *xpoint;
705: if (nfields) {
706: PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
707: } else {
708: PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
709: }
710: xpoint = &x[voff];
711: y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
712: }
713: }
714: cnt += off;
715: PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
716: }
717: }
718: PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
719: PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices, MPIU_VTUREAL, tag));
720: }
721: }
722: }
723: }
724: PetscCall(PetscFree(y));
725: PetscCall(VecRestoreArrayRead(X, &x));
726: }
727: } else if (rank == 0) {
728: PetscInt l;
730: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag)); /* positions */
731: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nconn, MPI_INT, tag)); /* connectivity */
732: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag)); /* offsets */
733: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_CHAR, tag)); /* types */
734: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag)); /* owner rank (cells) */
735: /* all cell data */
736: for (link = vtk->link; link; link = link->next) {
737: Vec X = (Vec)link->vec;
738: PetscInt bs = 1, nfields, field;
739: DM dmX;
740: PetscSection section = NULL;
742: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
743: PetscCall(VecGetDM(X, &dmX));
744: if (!dmX) dmX = dm;
745: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
746: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
747: if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
748: PetscCall(PetscSectionGetNumFields(section, &nfields));
749: field = 0;
750: if (link->field >= 0) {
751: field = link->field;
752: nfields = field + 1;
753: }
754: for (i = 0; field < (nfields ? nfields : 1); field++) {
755: PetscInt fbs, j;
756: PetscFV fv = NULL;
757: PetscObject f;
758: PetscClassId fClass;
759: PetscBool vector;
760: if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
761: PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
762: } else fbs = bs; /* Say we have one field with 'bs' components */
763: PetscCall(DMGetField(dmX, field, NULL, &f));
764: PetscCall(PetscObjectGetClassId(f, &fClass));
765: if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
766: vector = PETSC_FALSE;
767: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
768: vector = PETSC_TRUE;
769: for (j = 0; j < fbs; j++) {
770: const char *compName = NULL;
771: if (fv) {
772: PetscCall(PetscFVGetComponentName(fv, j, &compName));
773: if (compName) break;
774: }
775: }
776: if (j < fbs) vector = PETSC_FALSE;
777: }
778: if (vector) {
779: for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells * 3, MPIU_VTUREAL, tag));
780: } else {
781: for (i = 0; i < fbs; i++) {
782: for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPIU_VTUREAL, tag));
783: }
784: }
785: }
786: }
787: /* all point data */
788: for (link = vtk->link; link; link = link->next) {
789: Vec X = (Vec)link->vec;
790: DM dmX;
791: PetscInt bs = 1, nfields, field;
792: PetscSection section = NULL;
794: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
795: PetscCall(VecGetDM(X, &dmX));
796: if (!dmX) dmX = dm;
797: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
798: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
799: if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
800: PetscCall(PetscSectionGetNumFields(section, &nfields));
801: field = 0;
802: if (link->field >= 0) {
803: field = link->field;
804: nfields = field + 1;
805: }
806: for (i = 0; field < (nfields ? nfields : 1); field++) {
807: PetscInt fbs;
808: if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
809: PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
810: } else fbs = bs; /* Say we have one field with 'bs' components */
811: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
812: for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag));
813: } else {
814: for (i = 0; i < fbs; i++) {
815: for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices, MPIU_VTUREAL, tag));
816: }
817: }
818: }
819: }
820: }
821: }
822: PetscCall(PetscFree(gpiece));
823: PetscCall(PetscFree(buffer));
824: PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n"));
825: PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
826: PetscCall(PetscFClose(comm, fp));
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }