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