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