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