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