Actual source code: complexvtk.c
petsc-3.3-p7 2013-05-11
1: #define PETSCDM_DLL
2: #include <petsc-private/compleximpl.h> /*I "petscdmcomplex.h" I*/
3: #include <../src/sys/viewer/impls/vtk/vtkvimpl.h>
7: PetscErrorCode DMComplexVTKCellTypeValid(DM dm, PetscInt cellType, PetscBool *valid)
8: {
10: *valid = PETSC_TRUE;
11: return(0);
12: }
16: static PetscErrorCode DMComplexVTKGetCellType(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType) {
18: *cellType = -1;
19: switch(dim) {
20: case 0:
21: switch(corners) {
22: case 1:
23: *cellType = 1; /* VTK_VERTEX */
24: break;
25: default:
26: break;
27: }
28: break;
29: case 1:
30: switch(corners) {
31: case 2:
32: *cellType = 3; /* VTK_LINE */
33: break;
34: case 3:
35: *cellType = 21; /* VTK_QUADRATIC_EDGE */
36: break;
37: default:
38: break;
39: }
40: break;
41: case 2:
42: switch(corners) {
43: case 3:
44: *cellType = 5; /* VTK_TRIANGLE */
45: break;
46: case 4:
47: *cellType = 9; /* VTK_QUAD */
48: break;
49: case 6:
50: *cellType = 22; /* VTK_QUADRATIC_TRIANGLE */
51: break;
52: case 9:
53: *cellType = 23; /* VTK_QUADRATIC_QUAD */
54: break;
55: default:
56: break;
57: }
58: break;
59: case 3:
60: switch(corners) {
61: case 4:
62: *cellType = 10; /* VTK_TETRA */
63: break;
64: case 8:
65: *cellType = 12; /* VTK_HEXAHEDRON */
66: break;
67: case 10:
68: *cellType = 24; /* VTK_QUADRATIC_TETRA */
69: break;
70: case 27:
71: *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
72: break;
73: default:
74: break;
75: }
76: }
77: return(0);
78: }
82: PetscErrorCode DMComplexVTKWriteCells(DM dm, PetscSection globalConeSection, FILE *fp, PetscInt *totalCells)
83: {
84: MPI_Comm comm = ((PetscObject) dm)->comm;
85: PetscInt dim;
86: PetscInt numCorners = 0, maxCorners;
87: PetscInt numCells = 0, totCells, cellType;
88: PetscInt cStart, cEnd, c, vStart, vEnd, v;
89: PetscMPIInt numProcs, rank, proc, tag = 1;
93: MPI_Comm_size(comm, &numProcs);
94: MPI_Comm_rank(comm, &rank);
95: DMComplexGetDimension(dm, &dim);
96: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
97: DMComplexGetDepthStratum(dm, 0, &vStart, &vEnd);
98: for(c = cStart; c < cEnd; ++c) {
99: PetscInt nC;
100: PetscBool valid;
102: /* TODO: Does not work for interpolated meshes */
103: PetscSectionGetDof(globalConeSection, c, &nC);
104: DMComplexVTKGetCellType(dm, dim, nC, &cellType);
105: DMComplexVTKCellTypeValid(dm, cellType, &valid);
106: if (!valid) continue;
107: if (!numCorners) numCorners = nC;
108: ++numCells;
109: }
110: MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);
111: MPI_Allreduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, comm);
112: if (numCorners && numCorners != maxCorners) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "All elements must have %d corners, not %d", maxCorners, numCorners);
113: PetscFPrintf(comm, fp, "CELLS %d %d\n", totCells, totCells*(maxCorners+1));
114: if (!rank) {
115: for(c = cStart; c < cEnd; ++c) {
116: const PetscInt *cone;
117: PetscInt coneSize;
118: PetscBool valid;
120: /* TODO Use closure for interpolated meshes */
121: PetscSectionGetDof(globalConeSection, c, &coneSize);
122: DMComplexVTKGetCellType(dm, dim, coneSize, &cellType);
123: DMComplexVTKCellTypeValid(dm, cellType, &valid);
124: if (!valid) continue;
125: PetscFPrintf(comm, fp, "%d ", maxCorners);
126: DMComplexGetCone(dm, c, &cone);
127: /* TODO Need global vertex numbering */
128: for(v = 0; v < coneSize; ++v) {
129: PetscFPrintf(comm, fp, " %d", cone[v] - vStart);
130: }
131: PetscFPrintf(comm, fp, "\n");
132: }
133: for(proc = 1; proc < numProcs; ++proc) {
134: PetscInt *remoteVertices;
135: MPI_Status status;
137: MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
138: PetscMalloc(numCells*maxCorners * sizeof(PetscInt), &remoteVertices);
139: MPI_Recv(remoteVertices, numCells*maxCorners, MPIU_INT, proc, tag, comm, &status);
140: for(c = 0; c < numCells; ++c) {
141: PetscFPrintf(comm, fp, "%d ", maxCorners);
142: for(v = 0; v < maxCorners; ++v) {
143: PetscFPrintf(comm, fp, " %d", remoteVertices[c*maxCorners+v]);
144: }
145: PetscFPrintf(comm, fp, "\n");
146: }
147: PetscFree(remoteVertices);
148: }
149: } else {
150: PetscInt *localVertices, k = 0;
152: PetscMalloc(numCells*maxCorners * sizeof(PetscInt), &localVertices);
153: for(c = cStart; c < cEnd; ++c) {
154: const PetscInt *cone;
155: PetscInt coneSize;
156: PetscBool valid;
158: /* TODO Use closure for interpolated meshes */
159: PetscSectionGetDof(globalConeSection, c, &coneSize);
160: DMComplexVTKGetCellType(dm, dim, coneSize, &cellType);
161: DMComplexVTKCellTypeValid(dm, cellType, &valid);
162: if (!valid) continue;
163: DMComplexGetCone(dm, c, &cone);
164: /* TODO Need global vertex numbering */
165: for(v = 0; v < coneSize; ++v) {
166: localVertices[k++] = cone[v];
167: }
168: }
169: if (k != numCells*maxCorners) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numCells*maxCorners);
170: MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
171: MPI_Send(localVertices, numCells*maxCorners, MPI_INT, 0, tag, comm);
172: PetscFree(localVertices);
173: }
174: PetscFPrintf(comm, fp, "CELL_TYPES %d\n", totCells);
175: DMComplexVTKGetCellType(dm, dim, maxCorners, &cellType);
176: for(c = 0; c < totCells; ++c) {
177: PetscFPrintf(comm, fp, "%d\n", cellType);
178: }
179: *totalCells = totCells;
180: return(0);
181: }
185: PetscErrorCode DMComplexVTKWriteSection(DM dm, PetscSection section, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision) {
186: MPI_Comm comm = ((PetscObject) v)->comm;
187: const MPI_Datatype mpiType = MPIU_SCALAR;
188: PetscScalar *array;
189: PetscInt numDof = 0, maxDof;
190: PetscInt pStart, pEnd, p;
191: PetscMPIInt numProcs, rank, proc, tag = 1;
192: PetscErrorCode ierr;
195: if (precision < 0) precision = 6;
196: MPI_Comm_size(comm, &numProcs);
197: MPI_Comm_rank(comm, &rank);
198: PetscSectionGetChart(section, &pStart, &pEnd);
199: for(p = pStart; p < pEnd; ++p) {
200: PetscSectionGetDof(section, p, &numDof);
201: if (numDof) {break;}
202: }
203: MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
204: enforceDof = PetscMax(enforceDof, maxDof);
205: VecGetArray(v, &array);
206: if (!rank) {
207: char formatString[8];
209: PetscSNPrintf(formatString, 8, "%%.%de", precision);
210: for(p = pStart; p < pEnd; ++p) {
211: /* Here we lose a way to filter points by keeping them out of the Numbering */
212: PetscInt dof, off, d;
214: PetscSectionGetDof(section, p, &dof);
215: PetscSectionGetOffset(section, p, &off);
216: if (dof && off >= 0) {
217: for(d = 0; d < dof; d++) {
218: if (d > 0) {
219: PetscFPrintf(comm, fp, " ");
220: }
221: PetscFPrintf(comm, fp, formatString, PetscRealPart(array[off+d]));
222: }
223: for(d = dof; d < enforceDof; d++) {
224: PetscFPrintf(comm, fp, " 0.0");
225: }
226: PetscFPrintf(comm, fp, "\n");
227: }
228: }
229: for(proc = 1; proc < numProcs; ++proc) {
230: PetscScalar *remoteValues;
231: PetscInt size, d;
232: MPI_Status status;
234: MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
235: PetscMalloc(size * sizeof(PetscScalar), &remoteValues);
236: MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
237: for(p = 0; p < size/maxDof; ++p) {
238: for(d = 0; d < maxDof; ++d) {
239: if (d > 0) {
240: PetscFPrintf(comm, fp, " ");
241: }
242: PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d]));
243: }
244: for(d = maxDof; d < enforceDof; ++d) {
245: PetscFPrintf(comm, fp, " 0.0");
246: }
247: PetscFPrintf(comm, fp, "\n");
248: }
249: PetscFree(remoteValues);
250: }
251: } else {
252: PetscScalar *localValues;
253: PetscInt size, k = 0;
255: PetscSectionGetStorageSize(section, &size);
256: PetscMalloc(size * sizeof(PetscScalar), &localValues);
257: for(p = pStart; p < pEnd; ++p) {
258: PetscInt dof, off, d;
260: PetscSectionGetDof(section, p, &dof);
261: PetscSectionGetOffset(section, p, &off);
262: if (off >= 0) {
263: for(d = 0; d < dof; ++d) {
264: localValues[k++] = array[off+d];
265: }
266: }
267: }
268: MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
269: MPI_Send(localValues, k, mpiType, 0, tag, comm);
270: PetscFree(localValues);
271: }
272: VecRestoreArray(v, &array);
273: return(0);
274: }
278: PetscErrorCode DMComplexVTKWriteField(DM dm, PetscSection section, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision)
279: {
280: MPI_Comm comm = ((PetscObject) dm)->comm;
281: PetscInt numDof = 0, maxDof;
282: PetscInt pStart, pEnd, p;
286: PetscSectionGetChart(section, &pStart, &pEnd);
287: for(p = pStart; p < pEnd; ++p) {
288: PetscSectionGetDof(section, p, &numDof);
289: if (numDof) {break;}
290: }
291: numDof = PetscMax(numDof, enforceDof);
292: MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, ((PetscObject) dm)->comm);
293: if (!name) name = "Unknown";
294: if (maxDof == 3) {
295: PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
296: } else {
297: PetscFPrintf(comm, fp, "SCALARS %s double %d\n", name, maxDof);
298: PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
299: }
300: DMComplexVTKWriteSection(dm, section, field, fp, enforceDof, precision);
301: return(0);
302: }
306: static PetscErrorCode DMComplexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
307: {
308: MPI_Comm comm = ((PetscObject) dm)->comm;
309: PetscViewer_VTK *vtk = (PetscViewer_VTK *) viewer->data;
310: FILE *fp;
311: PetscViewerVTKObjectLink link;
312: PetscSection coordSection, globalCoordSection;
313: PetscSection coneSection, globalConeSection;
314: PetscLayout vLayout;
315: Vec coordinates;
316: PetscInt totVertices, totCells;
317: PetscBool hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE;
318: PetscErrorCode ierr;
321: PetscFOpen(comm, vtk->filename, "wb", &fp);
322: PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
323: PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");
324: PetscFPrintf(comm, fp, "ASCII\n");
325: PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
326: /* Vertices */
327: /* TODO: Need to check for "coordinates_dimensioned" */
328: DMComplexGetCoordinateSection(dm, &coordSection);
329: PetscSectionCreateGlobalSection(coordSection, dm->sf, &globalCoordSection);
330: DMComplexGetCoordinateVec(dm, &coordinates);
331: PetscSectionGetPointLayout(((PetscObject) dm)->comm, globalCoordSection, &vLayout);
332: PetscLayoutGetSize(vLayout, &totVertices);
333: PetscFPrintf(comm, fp, "POINTS %d double\n", totVertices);
334: DMComplexVTKWriteSection(dm, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE);
335: /* Cells */
336: DMComplexGetConeSection(dm, &coneSection);
337: PetscSectionCreateGlobalSection(coneSection, dm->sf, &globalConeSection);
338: DMComplexVTKWriteCells(dm, globalConeSection, fp, &totCells);
339: /* Vertex fields */
340: for(link = vtk->link; link; link = link->next) {
341: if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
342: if ((link->ft == PETSC_VTK_CELL_FIELD) || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD)) hasCell = PETSC_TRUE;
343: }
344: if (hasPoint) {
345: PetscFPrintf(comm, fp, "POINT_DATA %d\n", totVertices);
346: for(link = vtk->link; link; link = link->next) {
347: Vec X = (Vec) link->vec;
348: PetscContainer c;
349: PetscSection section, globalSection;
350: const char *name;
351: PetscInt enforceDof = PETSC_DETERMINE;
353: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
354: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
355: PetscObjectGetName(link->vec, &name);
356: PetscObjectQuery(link->vec, "section", (PetscObject *) &c);
357: if (!c) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
358: PetscContainerGetPointer(c, (void **) §ion);
359: if (!section) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
360: PetscSectionCreateGlobalSection(section, dm->sf, &globalSection);
361: DMComplexVTKWriteField(dm, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE);
362: PetscSectionDestroy(&globalSection);
363: }
364: }
365: /* Cell Fields */
366: if (hasCell) {
367: PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);
368: for(link = vtk->link; link; link = link->next) {
369: Vec X = (Vec) link->vec;
370: PetscContainer c;
371: PetscSection section, globalSection;
372: const char *name;
373: PetscInt enforceDof = PETSC_DETERMINE;
375: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
376: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
377: PetscObjectGetName(link->vec, &name);
378: PetscObjectQuery(link->vec, "section", (PetscObject *) &c);
379: if (!c) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
380: PetscContainerGetPointer(c, (void **) §ion);
381: if (!section) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
382: PetscSectionCreateGlobalSection(section, dm->sf, &globalSection);
383: DMComplexVTKWriteField(dm, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE);
384: PetscSectionDestroy(&globalSection);
385: }
386: }
387: /* Cleanup */
388: PetscSectionDestroy(&globalCoordSection);
389: PetscLayoutDestroy(&vLayout);
390: PetscSectionDestroy(&globalConeSection);
391: PetscFClose(comm, fp);
392: return(0);
393: }
397: /*@C
398: DMComplexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
400: Collective
402: Input Arguments:
403: + odm - The DMComplex specifying the mesh, passed as a PetscObject
404: - viewer - viewer of type VTK
406: Level: developer
408: Note:
409: This function is a callback used by the VTK viewer to actually write the file.
410: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
411: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
413: .seealso: PETSCVIEWERVTK
414: @*/
415: PetscErrorCode DMComplexVTKWriteAll(PetscObject odm, PetscViewer viewer)
416: {
417: DM dm = (DM) odm;
418: PetscBool isvtk;
419: PetscErrorCode ierr;
424: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
425: if (!isvtk) SETERRQ1(((PetscObject) viewer)->comm, PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
426: switch (viewer->format) {
427: case PETSC_VIEWER_ASCII_VTK:
428: DMComplexVTKWriteAll_ASCII(dm, viewer);
429: break;
430: default: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
431: }
432: return(0);
433: }