1 #include <petsc/private/dmpleximpl.h> 2 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h> 3 4 typedef struct { 5 PetscInt nvertices; 6 PetscInt ncells; 7 PetscInt nconn; /* number of entries in cell->vertex connectivity array */ 8 } PieceInfo; 9 10 #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16) 11 /* output in float if single or half precision in memory */ 12 static const char precision[] = "Float32"; 13 typedef float PetscVTUReal; 14 #define MPIU_VTUREAL MPI_FLOAT 15 #elif defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 16 /* output in double if double or quad precision in memory */ 17 static const char precision[] = "Float64"; 18 typedef double PetscVTUReal; 19 #define MPIU_VTUREAL MPI_DOUBLE 20 #else 21 static const char precision[] = "UnknownPrecision"; 22 typedef PetscReal PetscVTUReal; 23 #define MPIU_VTUREAL MPIU_REAL 24 #endif 25 26 static PetscErrorCode TransferWrite(MPI_Comm comm, PetscViewer viewer, FILE *fp, PetscMPIInt srank, PetscMPIInt root, const void *send, void *recv, PetscCount count, MPI_Datatype mpidatatype, PetscMPIInt tag) 27 { 28 PetscMPIInt rank; 29 30 PetscFunctionBegin; 31 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 32 if (rank == srank && rank != root) { 33 PetscCallMPI(MPIU_Send((void *)send, count, mpidatatype, root, tag, comm)); 34 } else if (rank == root) { 35 const void *buffer; 36 if (root == srank) { /* self */ 37 buffer = send; 38 } else { 39 MPI_Status status; 40 PetscMPIInt nrecv; 41 PetscCallMPI(MPIU_Recv(recv, count, mpidatatype, srank, tag, comm, &status)); 42 PetscCallMPI(MPI_Get_count(&status, mpidatatype, &nrecv)); 43 PetscCheck(count == nrecv, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch"); 44 buffer = recv; 45 } 46 PetscCall(PetscViewerVTKFWrite(viewer, fp, buffer, count, mpidatatype)); 47 } 48 PetscFunctionReturn(PETSC_SUCCESS); 49 } 50 51 static PetscErrorCode DMPlexGetVTKConnectivity(DM dm, PetscBool localized, PieceInfo *piece, PetscVTKInt **oconn, PetscVTKInt **ooffsets, PetscVTKType **otypes) 52 { 53 PetscSection coordSection, cellCoordSection; 54 PetscVTKInt *conn, *offsets; 55 PetscVTKType *types; 56 PetscInt dim, vStart, vEnd, cStart, cEnd, pStart, pEnd, cellHeight, numLabelCells, hasLabel, c, v, countcell, countconn; 57 58 PetscFunctionBegin; 59 PetscCall(PetscMalloc3(piece->nconn, &conn, piece->ncells, &offsets, piece->ncells, &types)); 60 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 61 PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection)); 62 PetscCall(DMGetDimension(dm, &dim)); 63 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 64 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 65 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 66 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 67 PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells)); 68 hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE; 69 70 countcell = 0; 71 countconn = 0; 72 for (c = cStart; c < cEnd; ++c) { 73 PetscInt nverts, dof = 0, celltype, startoffset, nC = 0; 74 75 if (hasLabel) { 76 PetscInt value; 77 78 PetscCall(DMGetLabelValue(dm, "vtk", c, &value)); 79 if (value != 1) continue; 80 } 81 startoffset = countconn; 82 if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof)); 83 if (!dof) { 84 PetscInt *closure = NULL; 85 PetscInt closureSize; 86 87 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 88 for (v = 0; v < closureSize * 2; v += 2) { 89 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 90 if (!localized) PetscCall(PetscVTKIntCast(closure[v] - vStart, &conn[countconn++])); 91 else PetscCall(PetscVTKIntCast(startoffset + nC, &conn[countconn++])); 92 ++nC; 93 } 94 } 95 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 96 } else { 97 for (nC = 0; nC < dof / dim; nC++) PetscCall(PetscVTKIntCast(startoffset + nC, &conn[countconn++])); 98 } 99 100 { 101 PetscInt n = PetscMin(nC, 8), s = countconn - nC, i, cone[8]; 102 for (i = 0; i < n; ++i) cone[i] = conn[s + i]; 103 PetscCall(DMPlexReorderCell(dm, c, cone)); 104 for (i = 0; i < n; ++i) PetscCall(PetscVTKIntCast(cone[i], &conn[s + i])); 105 } 106 PetscCall(PetscVTKIntCast(countconn, &offsets[countcell])); 107 108 nverts = countconn - startoffset; 109 PetscCall(DMPlexVTKGetCellType_Internal(dm, dim, nverts, &celltype)); 110 111 types[countcell] = (PetscVTKType)celltype; 112 countcell++; 113 } 114 PetscCheck(countcell == piece->ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent cell count"); 115 PetscCheck(countconn == piece->nconn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent connectivity count"); 116 *oconn = conn; 117 *ooffsets = offsets; 118 *otypes = types; 119 PetscFunctionReturn(PETSC_SUCCESS); 120 } 121 122 PETSC_INTERN PetscErrorCode DMPlexGetNonEmptyComm_Private(DM dm, MPI_Comm *comm) 123 { 124 DM_Plex *mesh = (DM_Plex *)dm->data; 125 126 PetscFunctionBegin; 127 if (mesh->nonempty_comm == MPI_COMM_SELF) { /* Not yet setup */ 128 PetscInt cStart, cEnd, cellHeight; 129 MPI_Comm dmcomm = PetscObjectComm((PetscObject)dm); 130 PetscMPIInt color, rank; 131 132 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 133 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 134 color = (cStart < cEnd) ? 0 : 1; 135 PetscCallMPI(MPI_Comm_rank(dmcomm, &rank)); 136 PetscCallMPI(MPI_Comm_split(dmcomm, color, rank, &mesh->nonempty_comm)); 137 if (color == 1) { 138 PetscCallMPI(MPI_Comm_free(&mesh->nonempty_comm)); 139 mesh->nonempty_comm = MPI_COMM_NULL; 140 } 141 } 142 *comm = mesh->nonempty_comm; 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 /* 147 Write all fields that have been provided to the viewer 148 Multi-block XML format with binary appended data. 149 */ 150 PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm, PetscViewer viewer) 151 { 152 MPI_Comm comm; 153 PetscSection coordSection, cellCoordSection; 154 PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data; 155 PetscViewerVTKObjectLink link; 156 FILE *fp; 157 PetscMPIInt rank, size, tag; 158 PetscInt dimEmbed, cellHeight, cStart, cEnd, vStart, vEnd, numLabelCells, hasLabel, c, v, i; 159 PetscBool localized; 160 PieceInfo piece, *gpiece = NULL; 161 void *buffer = NULL; 162 const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian"; 163 PetscInt loops_per_scalar; 164 165 PetscFunctionBegin; 166 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 167 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 168 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 169 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 170 PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells)); 171 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 172 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 173 PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection)); 174 PetscCall(PetscCommGetNewTag(PetscObjectComm((PetscObject)dm), &tag)); 175 176 PetscCall(DMPlexGetNonEmptyComm_Private(dm, &comm)); 177 #if defined(PETSC_USE_COMPLEX) 178 loops_per_scalar = 2; 179 #else 180 loops_per_scalar = 1; 181 #endif 182 if (comm == MPI_COMM_NULL) goto finalize; 183 PetscCallMPI(MPI_Comm_size(comm, &size)); 184 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 185 186 PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp)); 187 PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n")); 188 PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order)); 189 PetscCall(PetscFPrintf(comm, fp, " <UnstructuredGrid>\n")); 190 191 hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE; 192 piece.nvertices = 0; 193 piece.ncells = 0; 194 piece.nconn = 0; 195 if (!localized) piece.nvertices = vEnd - vStart; 196 for (c = cStart; c < cEnd; ++c) { 197 PetscInt dof = 0; 198 if (hasLabel) { 199 PetscInt value; 200 201 PetscCall(DMGetLabelValue(dm, "vtk", c, &value)); 202 if (value != 1) continue; 203 } 204 if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof)); 205 if (!dof) { 206 PetscInt *closure = NULL; 207 PetscInt closureSize; 208 209 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 210 for (v = 0; v < closureSize * 2; v += 2) { 211 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 212 piece.nconn++; 213 if (localized) piece.nvertices++; 214 } 215 } 216 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 217 } else { 218 piece.nvertices += dof / dimEmbed; 219 piece.nconn += dof / dimEmbed; 220 } 221 piece.ncells++; 222 } 223 if (rank == 0) PetscCall(PetscMalloc1(size, &gpiece)); 224 PetscCallMPI(MPI_Gather((PetscInt *)&piece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, (PetscInt *)gpiece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, 0, comm)); 225 226 /* 227 * Write file header 228 */ 229 if (rank == 0) { 230 PetscInt64 boffset = 0; 231 232 for (PetscMPIInt r = 0; r < size; r++) { 233 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Piece NumberOfPoints=\"%" PetscInt_FMT "\" NumberOfCells=\"%" PetscInt_FMT "\">\n", gpiece[r].nvertices, gpiece[r].ncells)); 234 /* Coordinate positions */ 235 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Points>\n")); 236 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset)); 237 boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 238 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Points>\n")); 239 /* Cell connectivity */ 240 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Cells>\n")); 241 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset)); 242 boffset += gpiece[r].nconn * sizeof(PetscVTKInt) + (gpiece[r].nconn ? sizeof(PetscInt64) : 0); 243 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset)); 244 boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 245 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset)); 246 boffset += gpiece[r].ncells * sizeof(unsigned char) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 247 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Cells>\n")); 248 249 /* 250 * Cell Data headers 251 */ 252 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <CellData>\n")); 253 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset)); 254 boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 255 /* all the vectors */ 256 for (link = vtk->link; link; link = link->next) { 257 Vec X = (Vec)link->vec; 258 DM dmX = NULL; 259 PetscInt bs = 1, nfields, field; 260 const char *vecname = ""; 261 PetscSection section; 262 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 263 if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */ 264 PetscCall(PetscObjectGetName((PetscObject)X, &vecname)); 265 } 266 PetscCall(VecGetDM(X, &dmX)); 267 if (!dmX) dmX = dm; 268 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 269 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 270 if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs)); 271 PetscCall(PetscSectionGetNumFields(section, &nfields)); 272 field = 0; 273 if (link->field >= 0) { 274 field = link->field; 275 nfields = field + 1; 276 } 277 for (i = 0; field < (nfields ? nfields : 1); field++) { 278 PetscInt fbs, j; 279 PetscFV fv = NULL; 280 PetscObject f; 281 PetscClassId fClass; 282 const char *fieldname = NULL; 283 char buf[256]; 284 PetscBool vector; 285 if (nfields) { /* We have user-defined fields/components */ 286 PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs)); 287 PetscCall(PetscSectionGetFieldName(section, field, &fieldname)); 288 } else fbs = bs; /* Say we have one field with 'bs' components */ 289 PetscCall(DMGetField(dmX, field, NULL, &f)); 290 PetscCall(PetscObjectGetClassId(f, &fClass)); 291 if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f; 292 if (nfields && !fieldname) { 293 PetscCall(PetscSNPrintf(buf, sizeof(buf), "CellField%" PetscInt_FMT, field)); 294 fieldname = buf; 295 } 296 vector = PETSC_FALSE; 297 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 298 vector = PETSC_TRUE; 299 PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Cell vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs); 300 for (j = 0; j < fbs; j++) { 301 const char *compName = NULL; 302 if (fv) { 303 PetscCall(PetscFVGetComponentName(fv, j, &compName)); 304 if (compName) break; 305 } 306 } 307 if (j < fbs) vector = PETSC_FALSE; 308 } 309 if (vector) { 310 #if defined(PETSC_USE_COMPLEX) 311 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 312 boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 313 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 314 boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 315 #else 316 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 317 boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 318 #endif 319 i += fbs; 320 } else { 321 for (j = 0; j < fbs; j++) { 322 const char *compName = NULL; 323 char finalname[256]; 324 if (fv) PetscCall(PetscFVGetComponentName(fv, j, &compName)); 325 if (compName) { 326 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName)); 327 } else if (fbs > 1) { 328 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%" PetscInt_FMT, vecname, fieldname, j)); 329 } else { 330 PetscCall(PetscSNPrintf(finalname, 255, "%s%s", vecname, fieldname)); 331 } 332 #if defined(PETSC_USE_COMPLEX) 333 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 334 boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 335 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 336 boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 337 #else 338 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 339 boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 340 #endif 341 i++; 342 } 343 } 344 } 345 //PetscCheck(i == bs,comm,PETSC_ERR_PLIB,"Total number of field components %" PetscInt_FMT " != block size %" PetscInt_FMT,i,bs); 346 } 347 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </CellData>\n")); 348 349 /* 350 * Point Data headers 351 */ 352 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <PointData>\n")); 353 for (link = vtk->link; link; link = link->next) { 354 Vec X = (Vec)link->vec; 355 DM dmX; 356 PetscInt bs = 1, nfields, field; 357 const char *vecname = ""; 358 PetscSection section; 359 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 360 if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */ 361 PetscCall(PetscObjectGetName((PetscObject)X, &vecname)); 362 } 363 PetscCall(VecGetDM(X, &dmX)); 364 if (!dmX) dmX = dm; 365 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 366 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 367 if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs)); 368 PetscCall(PetscSectionGetNumFields(section, &nfields)); 369 field = 0; 370 if (link->field >= 0) { 371 field = link->field; 372 nfields = field + 1; 373 } 374 for (i = 0; field < (nfields ? nfields : 1); field++) { 375 PetscInt fbs, j; 376 const char *fieldname = NULL; 377 char buf[256]; 378 if (nfields) { /* We have user-defined fields/components */ 379 PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs)); 380 PetscCall(PetscSectionGetFieldName(section, field, &fieldname)); 381 } else fbs = bs; /* Say we have one field with 'bs' components */ 382 if (nfields && !fieldname) { 383 PetscCall(PetscSNPrintf(buf, sizeof(buf), "PointField%" PetscInt_FMT, field)); 384 fieldname = buf; 385 } 386 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 387 PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Point vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs); 388 #if defined(PETSC_USE_COMPLEX) 389 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 390 boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 391 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 392 boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 393 #else 394 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 395 boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 396 #endif 397 } else { 398 for (j = 0; j < fbs; j++) { 399 const char *compName = NULL; 400 char finalname[256]; 401 PetscCall(PetscSectionGetComponentName(section, field, j, &compName)); 402 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName)); 403 #if defined(PETSC_USE_COMPLEX) 404 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 405 boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 406 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 407 boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 408 #else 409 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 410 boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 411 #endif 412 } 413 } 414 } 415 } 416 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </PointData>\n")); 417 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Piece>\n")); 418 } 419 } 420 421 PetscCall(PetscFPrintf(comm, fp, " </UnstructuredGrid>\n")); 422 PetscCall(PetscFPrintf(comm, fp, " <AppendedData encoding=\"raw\">\n")); 423 PetscCall(PetscFPrintf(comm, fp, "_")); 424 425 if (rank == 0) { 426 PetscInt maxsize = 0; 427 for (PetscMPIInt r = 0; r < size; r++) { 428 maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nvertices * 3 * sizeof(PetscVTUReal))); 429 maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].ncells * 3 * sizeof(PetscVTUReal))); 430 maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nconn * sizeof(PetscVTKInt))); 431 } 432 PetscCall(PetscMalloc(maxsize, &buffer)); 433 } 434 for (PetscMPIInt r = 0; r < size; r++) { 435 if (r == rank) { 436 PetscInt nsend; 437 { /* Position */ 438 const PetscScalar *x, *cx = NULL; 439 PetscVTUReal *y = NULL; 440 Vec coords, cellCoords; 441 PetscBool copy; 442 443 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 444 PetscCall(VecGetArrayRead(coords, &x)); 445 PetscCall(DMGetCellCoordinatesLocal(dm, &cellCoords)); 446 if (cellCoords) PetscCall(VecGetArrayRead(cellCoords, &cx)); 447 #if defined(PETSC_USE_COMPLEX) 448 copy = PETSC_TRUE; 449 #else 450 copy = (PetscBool)(dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal))); 451 #endif 452 if (copy) { 453 PetscCall(PetscMalloc1(piece.nvertices * 3, &y)); 454 if (localized) { 455 PetscInt cnt; 456 for (c = cStart, cnt = 0; c < cEnd; c++) { 457 PetscInt off, dof; 458 459 PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof)); 460 if (!dof) { 461 PetscInt *closure = NULL; 462 PetscInt closureSize; 463 464 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 465 for (v = 0; v < closureSize * 2; v += 2) { 466 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 467 PetscCall(PetscSectionGetOffset(coordSection, closure[v], &off)); 468 if (dimEmbed != 3) { 469 y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]); 470 y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[off + 1]) : 0.0); 471 y[cnt * 3 + 2] = (PetscVTUReal)0.0; 472 } else { 473 y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]); 474 y[cnt * 3 + 1] = (PetscVTUReal)PetscRealPart(x[off + 1]); 475 y[cnt * 3 + 2] = (PetscVTUReal)PetscRealPart(x[off + 2]); 476 } 477 cnt++; 478 } 479 } 480 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 481 } else { 482 PetscCall(PetscSectionGetOffset(cellCoordSection, c, &off)); 483 if (dimEmbed != 3) { 484 for (i = 0; i < dof / dimEmbed; i++) { 485 y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(cx[off + i * dimEmbed + 0]); 486 y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(cx[off + i * dimEmbed + 1]) : 0.0); 487 y[cnt * 3 + 2] = (PetscVTUReal)0.0; 488 cnt++; 489 } 490 } else { 491 for (i = 0; i < dof; i++) y[cnt * 3 + i] = (PetscVTUReal)PetscRealPart(cx[off + i]); 492 cnt += dof / dimEmbed; 493 } 494 } 495 } 496 PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match"); 497 } else { 498 for (i = 0; i < piece.nvertices; i++) { 499 y[i * 3 + 0] = (PetscVTUReal)PetscRealPart(x[i * dimEmbed + 0]); 500 y[i * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[i * dimEmbed + 1]) : 0.); 501 y[i * 3 + 2] = (PetscVTUReal)((dimEmbed > 2) ? PetscRealPart(x[i * dimEmbed + 2]) : 0.); 502 } 503 } 504 } 505 nsend = piece.nvertices * 3; 506 PetscCall(TransferWrite(comm, viewer, fp, r, 0, copy ? (const void *)y : (const void *)x, buffer, nsend, MPIU_VTUREAL, tag)); 507 PetscCall(PetscFree(y)); 508 PetscCall(VecRestoreArrayRead(coords, &x)); 509 if (cellCoords) PetscCall(VecRestoreArrayRead(cellCoords, &cx)); 510 } 511 { /* Connectivity, offsets, types */ 512 PetscVTKInt *connectivity = NULL, *offsets = NULL; 513 PetscVTKType *types = NULL; 514 PetscCall(DMPlexGetVTKConnectivity(dm, localized, &piece, &connectivity, &offsets, &types)); 515 PetscCall(TransferWrite(comm, viewer, fp, r, 0, connectivity, buffer, piece.nconn, MPI_INT, tag)); 516 PetscCall(TransferWrite(comm, viewer, fp, r, 0, offsets, buffer, piece.ncells, MPI_INT, tag)); 517 PetscCall(TransferWrite(comm, viewer, fp, r, 0, types, buffer, piece.ncells, MPI_CHAR, tag)); 518 PetscCall(PetscFree3(connectivity, offsets, types)); 519 } 520 { /* Owners (cell data) */ 521 PetscVTKInt *owners; 522 PetscMPIInt orank; 523 524 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &orank)); 525 PetscCall(PetscMalloc1(piece.ncells, &owners)); 526 for (i = 0; i < piece.ncells; i++) owners[i] = orank; 527 PetscCall(TransferWrite(comm, viewer, fp, r, 0, owners, buffer, piece.ncells, MPI_INT, tag)); 528 PetscCall(PetscFree(owners)); 529 } 530 /* Cell data */ 531 for (link = vtk->link; link; link = link->next) { 532 Vec X = (Vec)link->vec; 533 DM dmX; 534 const PetscScalar *x; 535 PetscVTUReal *y; 536 PetscInt bs = 1, nfields, field; 537 PetscSection section = NULL; 538 539 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 540 PetscCall(VecGetDM(X, &dmX)); 541 if (!dmX) dmX = dm; 542 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 543 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 544 if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs)); 545 PetscCall(PetscSectionGetNumFields(section, &nfields)); 546 field = 0; 547 if (link->field >= 0) { 548 field = link->field; 549 nfields = field + 1; 550 } 551 PetscCall(VecGetArrayRead(X, &x)); 552 PetscCall(PetscMalloc1(piece.ncells * 3, &y)); 553 for (i = 0; field < (nfields ? nfields : 1); field++) { 554 PetscInt fbs, j; 555 PetscFV fv = NULL; 556 PetscObject f; 557 PetscClassId fClass; 558 PetscBool vector; 559 if (nfields && cEnd > cStart) { /* We have user-defined fields/components */ 560 PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs)); 561 } else fbs = bs; /* Say we have one field with 'bs' components */ 562 PetscCall(DMGetField(dmX, field, NULL, &f)); 563 PetscCall(PetscObjectGetClassId(f, &fClass)); 564 if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f; 565 vector = PETSC_FALSE; 566 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 567 vector = PETSC_TRUE; 568 for (j = 0; j < fbs; j++) { 569 const char *compName = NULL; 570 if (fv) { 571 PetscCall(PetscFVGetComponentName(fv, j, &compName)); 572 if (compName) break; 573 } 574 } 575 if (j < fbs) vector = PETSC_FALSE; 576 } 577 if (vector) { 578 PetscInt cnt, l; 579 for (l = 0; l < loops_per_scalar; l++) { 580 for (c = cStart, cnt = 0; c < cEnd; c++) { 581 const PetscScalar *xpoint; 582 PetscInt off, j; 583 584 if (hasLabel) { /* Ignore some cells */ 585 PetscInt value; 586 PetscCall(DMGetLabelValue(dmX, "vtk", c, &value)); 587 if (value != 1) continue; 588 } 589 if (nfields) { 590 PetscCall(PetscSectionGetFieldOffset(section, c, field, &off)); 591 } else { 592 PetscCall(PetscSectionGetOffset(section, c, &off)); 593 } 594 xpoint = &x[off]; 595 for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j])); 596 for (; j < 3; j++) y[cnt++] = 0.; 597 } 598 PetscCheck(cnt == piece.ncells * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match"); 599 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells * 3, MPIU_VTUREAL, tag)); 600 } 601 } else { 602 for (i = 0; i < fbs; i++) { 603 PetscInt cnt, l; 604 for (l = 0; l < loops_per_scalar; l++) { 605 for (c = cStart, cnt = 0; c < cEnd; c++) { 606 const PetscScalar *xpoint; 607 PetscInt off; 608 609 if (hasLabel) { /* Ignore some cells */ 610 PetscInt value; 611 PetscCall(DMGetLabelValue(dmX, "vtk", c, &value)); 612 if (value != 1) continue; 613 } 614 if (nfields) { 615 PetscCall(PetscSectionGetFieldOffset(section, c, field, &off)); 616 } else { 617 PetscCall(PetscSectionGetOffset(section, c, &off)); 618 } 619 xpoint = &x[off]; 620 y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i])); 621 } 622 PetscCheck(cnt == piece.ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match"); 623 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells, MPIU_VTUREAL, tag)); 624 } 625 } 626 } 627 } 628 PetscCall(PetscFree(y)); 629 PetscCall(VecRestoreArrayRead(X, &x)); 630 } 631 /* point data */ 632 for (link = vtk->link; link; link = link->next) { 633 Vec X = (Vec)link->vec; 634 DM dmX; 635 const PetscScalar *x; 636 PetscVTUReal *y; 637 PetscInt bs = 1, nfields, field; 638 PetscSection section = NULL; 639 640 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 641 PetscCall(VecGetDM(X, &dmX)); 642 if (!dmX) dmX = dm; 643 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 644 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 645 if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs)); 646 PetscCall(PetscSectionGetNumFields(section, &nfields)); 647 field = 0; 648 if (link->field >= 0) { 649 field = link->field; 650 nfields = field + 1; 651 } 652 PetscCall(VecGetArrayRead(X, &x)); 653 PetscCall(PetscMalloc1(piece.nvertices * 3, &y)); 654 for (i = 0; field < (nfields ? nfields : 1); field++) { 655 PetscInt fbs, j; 656 if (nfields && vEnd > vStart) { /* We have user-defined fields/components */ 657 PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs)); 658 } else fbs = bs; /* Say we have one field with 'bs' components */ 659 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 660 PetscInt cnt, l; 661 for (l = 0; l < loops_per_scalar; l++) { 662 if (!localized) { 663 for (v = vStart, cnt = 0; v < vEnd; v++) { 664 PetscInt off; 665 const PetscScalar *xpoint; 666 667 if (nfields) { 668 PetscCall(PetscSectionGetFieldOffset(section, v, field, &off)); 669 } else { 670 PetscCall(PetscSectionGetOffset(section, v, &off)); 671 } 672 xpoint = &x[off]; 673 for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j])); 674 for (; j < 3; j++) y[cnt++] = 0.; 675 } 676 } else { 677 for (c = cStart, cnt = 0; c < cEnd; c++) { 678 PetscInt *closure = NULL; 679 PetscInt closureSize, off; 680 681 PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure)); 682 for (v = 0, off = 0; v < closureSize * 2; v += 2) { 683 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 684 PetscInt voff; 685 const PetscScalar *xpoint; 686 687 if (nfields) { 688 PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff)); 689 } else { 690 PetscCall(PetscSectionGetOffset(section, closure[v], &voff)); 691 } 692 xpoint = &x[voff]; 693 for (j = 0; j < fbs; j++) y[cnt + off++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j])); 694 for (; j < 3; j++) y[cnt + off++] = 0.; 695 } 696 } 697 cnt += off; 698 PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure)); 699 } 700 } 701 PetscCheck(cnt == piece.nvertices * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match"); 702 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices * 3, MPIU_VTUREAL, tag)); 703 } 704 } else { 705 for (i = 0; i < fbs; i++) { 706 PetscInt cnt, l; 707 for (l = 0; l < loops_per_scalar; l++) { 708 if (!localized) { 709 for (v = vStart, cnt = 0; v < vEnd; v++) { 710 PetscInt off; 711 const PetscScalar *xpoint; 712 713 if (nfields) { 714 PetscCall(PetscSectionGetFieldOffset(section, v, field, &off)); 715 } else { 716 PetscCall(PetscSectionGetOffset(section, v, &off)); 717 } 718 xpoint = &x[off]; 719 y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i])); 720 } 721 } else { 722 for (c = cStart, cnt = 0; c < cEnd; c++) { 723 PetscInt *closure = NULL; 724 PetscInt closureSize, off; 725 726 PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure)); 727 for (v = 0, off = 0; v < closureSize * 2; v += 2) { 728 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 729 PetscInt voff; 730 const PetscScalar *xpoint; 731 732 if (nfields) { 733 PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff)); 734 } else { 735 PetscCall(PetscSectionGetOffset(section, closure[v], &voff)); 736 } 737 xpoint = &x[voff]; 738 y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i])); 739 } 740 } 741 cnt += off; 742 PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure)); 743 } 744 } 745 PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match"); 746 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices, MPIU_VTUREAL, tag)); 747 } 748 } 749 } 750 } 751 PetscCall(PetscFree(y)); 752 PetscCall(VecRestoreArrayRead(X, &x)); 753 } 754 } else if (rank == 0) { 755 PetscInt l; 756 757 PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag)); /* positions */ 758 PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nconn, MPI_INT, tag)); /* connectivity */ 759 PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag)); /* offsets */ 760 PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_CHAR, tag)); /* types */ 761 PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag)); /* owner rank (cells) */ 762 /* all cell data */ 763 for (link = vtk->link; link; link = link->next) { 764 Vec X = (Vec)link->vec; 765 PetscInt bs = 1, nfields, field; 766 DM dmX; 767 PetscSection section = NULL; 768 769 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 770 PetscCall(VecGetDM(X, &dmX)); 771 if (!dmX) dmX = dm; 772 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 773 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 774 if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs)); 775 PetscCall(PetscSectionGetNumFields(section, &nfields)); 776 field = 0; 777 if (link->field >= 0) { 778 field = link->field; 779 nfields = field + 1; 780 } 781 for (i = 0; field < (nfields ? nfields : 1); field++) { 782 PetscInt fbs, j; 783 PetscFV fv = NULL; 784 PetscObject f; 785 PetscClassId fClass; 786 PetscBool vector; 787 if (nfields && cEnd > cStart) { /* We have user-defined fields/components */ 788 PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs)); 789 } else fbs = bs; /* Say we have one field with 'bs' components */ 790 PetscCall(DMGetField(dmX, field, NULL, &f)); 791 PetscCall(PetscObjectGetClassId(f, &fClass)); 792 if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f; 793 vector = PETSC_FALSE; 794 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 795 vector = PETSC_TRUE; 796 for (j = 0; j < fbs; j++) { 797 const char *compName = NULL; 798 if (fv) { 799 PetscCall(PetscFVGetComponentName(fv, j, &compName)); 800 if (compName) break; 801 } 802 } 803 if (j < fbs) vector = PETSC_FALSE; 804 } 805 if (vector) { 806 for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells * 3, MPIU_VTUREAL, tag)); 807 } else { 808 for (i = 0; i < fbs; i++) { 809 for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPIU_VTUREAL, tag)); 810 } 811 } 812 } 813 } 814 /* all point data */ 815 for (link = vtk->link; link; link = link->next) { 816 Vec X = (Vec)link->vec; 817 DM dmX; 818 PetscInt bs = 1, nfields, field; 819 PetscSection section = NULL; 820 821 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 822 PetscCall(VecGetDM(X, &dmX)); 823 if (!dmX) dmX = dm; 824 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 825 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 826 if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs)); 827 PetscCall(PetscSectionGetNumFields(section, &nfields)); 828 field = 0; 829 if (link->field >= 0) { 830 field = link->field; 831 nfields = field + 1; 832 } 833 for (i = 0; field < (nfields ? nfields : 1); field++) { 834 PetscInt fbs; 835 if (nfields && vEnd > vStart) { /* We have user-defined fields/components */ 836 PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs)); 837 } else fbs = bs; /* Say we have one field with 'bs' components */ 838 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 839 for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag)); 840 } else { 841 for (i = 0; i < fbs; i++) { 842 for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices, MPIU_VTUREAL, tag)); 843 } 844 } 845 } 846 } 847 } 848 } 849 PetscCall(PetscFree(gpiece)); 850 PetscCall(PetscFree(buffer)); 851 PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n")); 852 PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n")); 853 PetscCall(PetscFClose(comm, fp)); 854 finalize: 855 /* this code sends to rank 0 that writes. 856 It may lead to very unbalanced log_view timings 857 of the next PETSc function logged. 858 Since this call is not performance critical, we 859 issue a barrier here to synchronize the processes */ 860 PetscCall(PetscBarrier((PetscObject)viewer)); 861 PetscFunctionReturn(PETSC_SUCCESS); 862 } 863