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