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