Actual source code: plexvtk.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h>
3: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
5: PetscErrorCode DMPlexVTKGetCellType_Internal(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType)
6: {
7: *cellType = -1;
8: switch (dim) {
9: case 0:
10: switch (corners) {
11: case 1:
12: *cellType = 1; /* VTK_VERTEX */
13: break;
14: default:
15: break;
16: }
17: break;
18: case 1:
19: switch (corners) {
20: case 2:
21: *cellType = 3; /* VTK_LINE */
22: break;
23: case 3:
24: *cellType = 21; /* VTK_QUADRATIC_EDGE */
25: break;
26: default:
27: break;
28: }
29: break;
30: case 2:
31: switch (corners) {
32: case 3:
33: *cellType = 5; /* VTK_TRIANGLE */
34: break;
35: case 4:
36: *cellType = 9; /* VTK_QUAD */
37: break;
38: case 6:
39: *cellType = 22; /* VTK_QUADRATIC_TRIANGLE */
40: break;
41: case 9:
42: *cellType = 23; /* VTK_QUADRATIC_QUAD */
43: break;
44: default:
45: break;
46: }
47: break;
48: case 3:
49: switch (corners) {
50: case 4:
51: *cellType = 10; /* VTK_TETRA */
52: break;
53: case 5:
54: *cellType = 14; /* VTK_PYRAMID */
55: break;
56: case 6:
57: *cellType = 13; /* VTK_WEDGE */
58: break;
59: case 8:
60: *cellType = 12; /* VTK_HEXAHEDRON */
61: break;
62: case 10:
63: *cellType = 24; /* VTK_QUADRATIC_TETRA */
64: break;
65: case 27:
66: *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
67: break;
68: default:
69: break;
70: }
71: }
72: return 0;
73: }
75: static PetscErrorCode DMPlexVTKWriteCells_ASCII(DM dm, FILE *fp, PetscInt *totalCells)
76: {
77: MPI_Comm comm;
78: DMLabel label;
79: IS globalVertexNumbers = NULL;
80: const PetscInt *gvertex;
81: PetscInt dim;
82: PetscInt numCorners = 0, totCorners = 0, maxCorners, *corners;
83: PetscInt numCells = 0, totCells = 0, maxCells, cellHeight;
84: PetscInt numLabelCells, maxLabelCells, cStart, cEnd, c, vStart, vEnd, v;
85: PetscMPIInt size, rank, proc, tag;
87: PetscObjectGetComm((PetscObject)dm, &comm);
88: PetscCommGetNewTag(comm, &tag);
89: MPI_Comm_size(comm, &size);
90: MPI_Comm_rank(comm, &rank);
91: DMGetDimension(dm, &dim);
92: DMPlexGetVTKCellHeight(dm, &cellHeight);
93: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
94: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
95: DMGetLabel(dm, "vtk", &label);
96: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
97: MPIU_Allreduce(&numLabelCells, &maxLabelCells, 1, MPIU_INT, MPI_MAX, comm);
98: if (!maxLabelCells) label = NULL;
99: for (c = cStart; c < cEnd; ++c) {
100: PetscInt *closure = NULL;
101: PetscInt closureSize, value;
103: if (label) {
104: DMLabelGetValue(label, c, &value);
105: if (value != 1) continue;
106: }
107: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
108: for (v = 0; v < closureSize * 2; v += 2) {
109: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
110: }
111: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
112: ++numCells;
113: }
114: maxCells = numCells;
115: MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);
116: MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);
117: MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);
118: MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);
119: DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
120: ISGetIndices(globalVertexNumbers, &gvertex);
121: PetscMalloc1(maxCells, &corners);
122: PetscFPrintf(comm, fp, "CELLS %" PetscInt_FMT " %" PetscInt_FMT "\n", totCells, totCorners + totCells);
123: if (rank == 0) {
124: PetscInt *remoteVertices, *vertices;
126: PetscMalloc1(maxCorners, &vertices);
127: for (c = cStart, numCells = 0; c < cEnd; ++c) {
128: PetscInt *closure = NULL;
129: PetscInt closureSize, value, nC = 0;
131: if (label) {
132: DMLabelGetValue(label, c, &value);
133: if (value != 1) continue;
134: }
135: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
136: for (v = 0; v < closureSize * 2; v += 2) {
137: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
138: const PetscInt gv = gvertex[closure[v] - vStart];
139: vertices[nC++] = gv < 0 ? -(gv + 1) : gv;
140: }
141: }
142: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
143: DMPlexReorderCell(dm, c, vertices);
144: corners[numCells++] = nC;
145: PetscFPrintf(comm, fp, "%" PetscInt_FMT " ", nC);
146: for (v = 0; v < nC; ++v) PetscFPrintf(comm, fp, " %" PetscInt_FMT, vertices[v]);
147: PetscFPrintf(comm, fp, "\n");
148: }
149: if (size > 1) PetscMalloc1(maxCorners + maxCells, &remoteVertices);
150: for (proc = 1; proc < size; ++proc) {
151: MPI_Status status;
153: MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);
154: MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);
155: for (c = 0; c < numCorners;) {
156: PetscInt nC = remoteVertices[c++];
158: for (v = 0; v < nC; ++v, ++c) vertices[v] = remoteVertices[c];
159: PetscFPrintf(comm, fp, "%" PetscInt_FMT " ", nC);
160: for (v = 0; v < nC; ++v) PetscFPrintf(comm, fp, " %" PetscInt_FMT, vertices[v]);
161: PetscFPrintf(comm, fp, "\n");
162: }
163: }
164: if (size > 1) PetscFree(remoteVertices);
165: PetscFree(vertices);
166: } else {
167: PetscInt *localVertices, numSend = numCells + numCorners, k = 0;
169: PetscMalloc1(numSend, &localVertices);
170: for (c = cStart, numCells = 0; c < cEnd; ++c) {
171: PetscInt *closure = NULL;
172: PetscInt closureSize, value, nC = 0;
174: if (label) {
175: DMLabelGetValue(label, c, &value);
176: if (value != 1) continue;
177: }
178: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
179: for (v = 0; v < closureSize * 2; v += 2) {
180: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
181: const PetscInt gv = gvertex[closure[v] - vStart];
182: closure[nC++] = gv < 0 ? -(gv + 1) : gv;
183: }
184: }
185: corners[numCells++] = nC;
186: localVertices[k++] = nC;
187: for (v = 0; v < nC; ++v, ++k) localVertices[k] = closure[v];
188: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
189: DMPlexReorderCell(dm, c, localVertices + k - nC);
190: }
192: MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm);
193: MPI_Send(localVertices, numSend, MPIU_INT, 0, tag, comm);
194: PetscFree(localVertices);
195: }
196: ISRestoreIndices(globalVertexNumbers, &gvertex);
197: PetscFPrintf(comm, fp, "CELL_TYPES %" PetscInt_FMT "\n", totCells);
198: if (rank == 0) {
199: PetscInt cellType;
201: for (c = 0; c < numCells; ++c) {
202: DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);
203: PetscFPrintf(comm, fp, "%" PetscInt_FMT "\n", cellType);
204: }
205: for (proc = 1; proc < size; ++proc) {
206: MPI_Status status;
208: MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
209: MPI_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status);
210: for (c = 0; c < numCells; ++c) {
211: DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);
212: PetscFPrintf(comm, fp, "%" PetscInt_FMT "\n", cellType);
213: }
214: }
215: } else {
216: MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
217: MPI_Send(corners, numCells, MPIU_INT, 0, tag, comm);
218: }
219: PetscFree(corners);
220: *totalCells = totCells;
221: return 0;
222: }
224: static PetscErrorCode DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
225: {
226: MPI_Comm comm;
227: PetscInt numCells = 0, cellHeight;
228: PetscInt numLabelCells, cStart, cEnd, c;
229: PetscMPIInt size, rank, proc, tag;
230: PetscBool hasLabel;
232: PetscObjectGetComm((PetscObject)dm, &comm);
233: PetscCommGetNewTag(comm, &tag);
234: MPI_Comm_size(comm, &size);
235: MPI_Comm_rank(comm, &rank);
236: DMPlexGetVTKCellHeight(dm, &cellHeight);
237: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
238: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
239: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
240: for (c = cStart; c < cEnd; ++c) {
241: if (hasLabel) {
242: PetscInt value;
244: DMGetLabelValue(dm, "vtk", c, &value);
245: if (value != 1) continue;
246: }
247: ++numCells;
248: }
249: if (rank == 0) {
250: for (c = 0; c < numCells; ++c) PetscFPrintf(comm, fp, "%d\n", rank);
251: for (proc = 1; proc < size; ++proc) {
252: MPI_Status status;
254: MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
255: for (c = 0; c < numCells; ++c) PetscFPrintf(comm, fp, "%d\n", proc);
256: }
257: } else {
258: MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
259: }
260: return 0;
261: }
263: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
264: typedef double PetscVTKReal;
265: #elif defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
266: typedef float PetscVTKReal;
267: #else
268: typedef PetscReal PetscVTKReal;
269: #endif
271: static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscInt imag)
272: {
273: MPI_Comm comm;
274: const MPI_Datatype mpiType = MPIU_SCALAR;
275: PetscScalar *array;
276: PetscInt numDof = 0, maxDof;
277: PetscInt numLabelCells, cellHeight, cStart, cEnd, numLabelVertices, vStart, vEnd, pStart, pEnd, p;
278: PetscMPIInt size, rank, proc, tag;
279: PetscBool hasLabel;
281: PetscObjectGetComm((PetscObject)dm, &comm);
284: if (precision < 0) precision = 6;
285: PetscCommGetNewTag(comm, &tag);
286: MPI_Comm_size(comm, &size);
287: MPI_Comm_rank(comm, &rank);
288: PetscSectionGetChart(section, &pStart, &pEnd);
289: /* VTK only wants the values at cells or vertices */
290: DMPlexGetVTKCellHeight(dm, &cellHeight);
291: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
292: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
293: pStart = PetscMax(PetscMin(cStart, vStart), pStart);
294: pEnd = PetscMin(PetscMax(cEnd, vEnd), pEnd);
295: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
296: DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);
297: hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
298: for (p = pStart; p < pEnd; ++p) {
299: /* Reject points not either cells or vertices */
300: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
301: if (hasLabel) {
302: PetscInt value;
304: if (((p >= cStart) && (p < cEnd) && numLabelCells) || ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
305: DMGetLabelValue(dm, "vtk", p, &value);
306: if (value != 1) continue;
307: }
308: }
309: PetscSectionGetDof(section, p, &numDof);
310: if (numDof) break;
311: }
312: MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
313: enforceDof = PetscMax(enforceDof, maxDof);
314: VecGetArray(v, &array);
315: if (rank == 0) {
316: PetscVTKReal dval;
317: PetscScalar val;
318: char formatString[8];
320: PetscSNPrintf(formatString, 8, "%%.%" PetscInt_FMT "e", precision);
321: for (p = pStart; p < pEnd; ++p) {
322: /* Here we lose a way to filter points by keeping them out of the Numbering */
323: PetscInt dof, off, goff, d;
325: /* Reject points not either cells or vertices */
326: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
327: if (hasLabel) {
328: PetscInt value;
330: if (((p >= cStart) && (p < cEnd) && numLabelCells) || ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
331: DMGetLabelValue(dm, "vtk", p, &value);
332: if (value != 1) continue;
333: }
334: }
335: PetscSectionGetDof(section, p, &dof);
336: PetscSectionGetOffset(section, p, &off);
337: PetscSectionGetOffset(globalSection, p, &goff);
338: if (dof && goff >= 0) {
339: for (d = 0; d < dof; d++) {
340: if (d > 0) PetscFPrintf(comm, fp, " ");
341: val = array[off + d];
342: dval = (PetscVTKReal)((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
343: PetscFPrintf(comm, fp, formatString, dval);
344: }
345: for (d = dof; d < enforceDof; d++) PetscFPrintf(comm, fp, " 0.0");
346: PetscFPrintf(comm, fp, "\n");
347: }
348: }
349: for (proc = 1; proc < size; ++proc) {
350: PetscScalar *remoteValues;
351: PetscInt size = 0, d;
352: MPI_Status status;
354: MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
355: PetscMalloc1(size, &remoteValues);
356: MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
357: for (p = 0; p < size / maxDof; ++p) {
358: for (d = 0; d < maxDof; ++d) {
359: if (d > 0) PetscFPrintf(comm, fp, " ");
360: val = remoteValues[p * maxDof + d];
361: dval = (PetscVTKReal)((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
362: PetscFPrintf(comm, fp, formatString, dval);
363: }
364: for (d = maxDof; d < enforceDof; ++d) PetscFPrintf(comm, fp, " 0.0");
365: PetscFPrintf(comm, fp, "\n");
366: }
367: PetscFree(remoteValues);
368: }
369: } else {
370: PetscScalar *localValues;
371: PetscInt size, k = 0;
373: PetscSectionGetStorageSize(section, &size);
374: PetscMalloc1(size, &localValues);
375: for (p = pStart; p < pEnd; ++p) {
376: PetscInt dof, off, goff, d;
378: /* Reject points not either cells or vertices */
379: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
380: if (hasLabel) {
381: PetscInt value;
383: if (((p >= cStart) && (p < cEnd) && numLabelCells) || ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
384: DMGetLabelValue(dm, "vtk", p, &value);
385: if (value != 1) continue;
386: }
387: }
388: PetscSectionGetDof(section, p, &dof);
389: PetscSectionGetOffset(section, p, &off);
390: PetscSectionGetOffset(globalSection, p, &goff);
391: if (goff >= 0) {
392: for (d = 0; d < dof; ++d) localValues[k++] = array[off + d];
393: }
394: }
395: MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
396: MPI_Send(localValues, k, mpiType, 0, tag, comm);
397: PetscFree(localValues);
398: }
399: VecRestoreArray(v, &array);
400: return 0;
401: }
403: static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscBool nameComplex, PetscInt imag)
404: {
405: MPI_Comm comm;
406: PetscInt numDof = 0, maxDof;
407: PetscInt pStart, pEnd, p;
409: PetscObjectGetComm((PetscObject)dm, &comm);
410: PetscSectionGetChart(section, &pStart, &pEnd);
411: for (p = pStart; p < pEnd; ++p) {
412: PetscSectionGetDof(section, p, &numDof);
413: if (numDof) break;
414: }
415: numDof = PetscMax(numDof, enforceDof);
416: MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
417: if (!name) name = "Unknown";
418: if (maxDof == 3) {
419: if (nameComplex) {
420: PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re");
421: } else {
422: PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
423: }
424: } else {
425: if (nameComplex) {
426: PetscFPrintf(comm, fp, "SCALARS %s.%s double %" PetscInt_FMT "\n", name, imag ? "Im" : "Re", maxDof);
427: } else {
428: PetscFPrintf(comm, fp, "SCALARS %s double %" PetscInt_FMT "\n", name, maxDof);
429: }
430: PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
431: }
432: DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag);
433: return 0;
434: }
436: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
437: {
438: MPI_Comm comm;
439: PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data;
440: FILE *fp;
441: PetscViewerVTKObjectLink link;
442: PetscInt totVertices, totCells = 0, loops_per_scalar, l;
443: PetscBool hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
444: const char *dmname;
446: #if defined(PETSC_USE_COMPLEX)
447: loops_per_scalar = 2;
448: writeComplex = PETSC_TRUE;
449: #else
450: loops_per_scalar = 1;
451: writeComplex = PETSC_FALSE;
452: #endif
453: DMGetCoordinatesLocalized(dm, &localized);
454: PetscObjectGetComm((PetscObject)dm, &comm);
456: PetscFOpen(comm, vtk->filename, "wb", &fp);
457: PetscObjectGetName((PetscObject)dm, &dmname);
458: PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
459: PetscFPrintf(comm, fp, "%s\n", dmname);
460: PetscFPrintf(comm, fp, "ASCII\n");
461: PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
462: /* Vertices */
463: {
464: PetscSection section, coordSection, globalCoordSection;
465: Vec coordinates;
466: PetscReal lengthScale;
467: DMLabel label;
468: IS vStratumIS;
469: PetscLayout vLayout;
471: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
472: DMGetCoordinatesLocal(dm, &coordinates);
473: DMPlexGetDepthLabel(dm, &label);
474: DMLabelGetStratumIS(label, 0, &vStratumIS);
475: DMGetCoordinateSection(dm, §ion); /* This section includes all points */
476: PetscSectionCreateSubdomainSection(section, vStratumIS, &coordSection); /* This one includes just vertices */
477: PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);
478: PetscSectionGetPointLayout(comm, globalCoordSection, &vLayout);
479: PetscLayoutGetSize(vLayout, &totVertices);
480: PetscFPrintf(comm, fp, "POINTS %" PetscInt_FMT " double\n", totVertices);
481: DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0);
482: ISDestroy(&vStratumIS);
483: PetscLayoutDestroy(&vLayout);
484: PetscSectionDestroy(&coordSection);
485: PetscSectionDestroy(&globalCoordSection);
486: }
487: /* Cells */
488: DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
489: /* Vertex fields */
490: for (link = vtk->link; link; link = link->next) {
491: if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
492: if ((link->ft == PETSC_VTK_CELL_FIELD) || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD)) hasCell = PETSC_TRUE;
493: }
494: if (hasPoint) {
495: PetscFPrintf(comm, fp, "POINT_DATA %" PetscInt_FMT "\n", totVertices);
496: for (link = vtk->link; link; link = link->next) {
497: Vec X = (Vec)link->vec;
498: PetscSection section = NULL, globalSection, newSection = NULL;
499: char namebuf[256];
500: const char *name;
501: PetscInt enforceDof = PETSC_DETERMINE;
503: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
504: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
505: PetscObjectGetName(link->vec, &name);
506: PetscObjectQuery(link->vec, "section", (PetscObject *)§ion);
507: if (!section) {
508: DM dmX;
510: VecGetDM(X, &dmX);
511: if (dmX) {
512: DMLabel subpointMap, subpointMapX;
513: PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
515: DMGetLocalSection(dmX, §ion);
516: /* Here is where we check whether dmX is a submesh of dm */
517: DMGetDimension(dm, &dim);
518: DMGetDimension(dmX, &dimX);
519: DMPlexGetChart(dm, &pStart, &pEnd);
520: DMPlexGetChart(dmX, &qStart, &qEnd);
521: DMPlexGetSubpointMap(dm, &subpointMap);
522: DMPlexGetSubpointMap(dmX, &subpointMapX);
523: if (((dim != dimX) || ((pEnd - pStart) < (qEnd - qStart))) && subpointMap && !subpointMapX) {
524: const PetscInt *ind = NULL;
525: IS subpointIS;
526: PetscInt n = 0, q;
528: PetscSectionGetChart(section, &qStart, &qEnd);
529: DMPlexGetSubpointIS(dm, &subpointIS);
530: if (subpointIS) {
531: ISGetLocalSize(subpointIS, &n);
532: ISGetIndices(subpointIS, &ind);
533: }
534: PetscSectionCreate(comm, &newSection);
535: PetscSectionSetChart(newSection, pStart, pEnd);
536: for (q = qStart; q < qEnd; ++q) {
537: PetscInt dof, off, p;
539: PetscSectionGetDof(section, q, &dof);
540: if (dof) {
541: PetscFindInt(q, n, ind, &p);
542: if (p >= pStart) {
543: PetscSectionSetDof(newSection, p, dof);
544: PetscSectionGetOffset(section, q, &off);
545: PetscSectionSetOffset(newSection, p, off);
546: }
547: }
548: }
549: if (subpointIS) ISRestoreIndices(subpointIS, &ind);
550: /* No need to setup section */
551: section = newSection;
552: }
553: }
554: }
556: if (link->field >= 0) {
557: const char *fieldname;
559: PetscSectionGetFieldName(section, link->field, &fieldname);
560: PetscSectionGetField(section, link->field, §ion);
561: if (fieldname) {
562: PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);
563: } else {
564: PetscSNPrintf(namebuf, sizeof(namebuf), "%s%" PetscInt_FMT, name, link->field);
565: }
566: } else {
567: PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);
568: }
569: PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));
570: PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
571: for (l = 0; l < loops_per_scalar; l++) DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);
572: PetscSectionDestroy(&globalSection);
573: if (newSection) PetscSectionDestroy(&newSection);
574: }
575: }
576: /* Cell Fields */
577: PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_view_partition", &writePartition, NULL);
578: if (hasCell || writePartition) {
579: PetscFPrintf(comm, fp, "CELL_DATA %" PetscInt_FMT "\n", totCells);
580: for (link = vtk->link; link; link = link->next) {
581: Vec X = (Vec)link->vec;
582: PetscSection section = NULL, globalSection;
583: const char *name = "";
584: char namebuf[256];
585: PetscInt enforceDof = PETSC_DETERMINE;
587: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
588: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
589: PetscObjectGetName(link->vec, &name);
590: PetscObjectQuery(link->vec, "section", (PetscObject *)§ion);
591: if (!section) {
592: DM dmX;
594: VecGetDM(X, &dmX);
595: if (dmX) DMGetLocalSection(dmX, §ion);
596: }
598: if (link->field >= 0) {
599: const char *fieldname;
601: PetscSectionGetFieldName(section, link->field, &fieldname);
602: PetscSectionGetField(section, link->field, §ion);
603: if (fieldname) {
604: PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);
605: } else {
606: PetscSNPrintf(namebuf, sizeof(namebuf), "%s%" PetscInt_FMT, name, link->field);
607: }
608: } else {
609: PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);
610: }
611: PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));
612: PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
613: for (l = 0; l < loops_per_scalar; l++) DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);
614: PetscSectionDestroy(&globalSection);
615: }
616: if (writePartition) {
617: PetscFPrintf(comm, fp, "SCALARS partition int 1\n");
618: PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
619: DMPlexVTKWritePartition_ASCII(dm, fp);
620: }
621: }
622: /* Cleanup */
623: PetscFClose(comm, fp);
624: return 0;
625: }
627: /*@C
628: DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
630: Collective
632: Input Parameters:
633: + odm - The `DMPLEX` specifying the mesh, passed as a `PetscObject`
634: - viewer - viewer of type `PETSCVIEWERVTK`
636: Level: developer
638: Note:
639: This function is a callback used by the VTK viewer to actually write the file.
640: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
641: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
643: .seealso: [](chapter_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `PETSCVIEWERVTK`
644: @*/
645: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
646: {
647: DM dm = (DM)odm;
651: switch (viewer->format) {
652: case PETSC_VIEWER_ASCII_VTK_DEPRECATED:
653: DMPlexVTKWriteAll_ASCII(dm, viewer);
654: break;
655: case PETSC_VIEWER_VTK_VTU:
656: DMPlexVTKWriteAll_VTU(dm, viewer);
657: break;
658: default:
659: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
660: }
661: return 0;
662: }