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, PetscMPIInt 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(MPI_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(MPI_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) conn[countconn++] = closure[v] - vStart; 91 else conn[countconn++] = startoffset + nC; 92 ++nC; 93 } 94 } 95 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 96 } else { 97 for (nC = 0; nC < dof / dim; nC++) conn[countconn++] = startoffset + nC; 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) conn[s + i] = (int)cone[i]; 105 } 106 107 offsets[countcell] = countconn; 108 109 nverts = countconn - startoffset; 110 PetscCall(DMPlexVTKGetCellType_Internal(dm, dim, nverts, &celltype)); 111 112 types[countcell] = celltype; 113 countcell++; 114 } 115 PetscCheck(countcell == piece->ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent cell count"); 116 PetscCheck(countconn == piece->nconn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent connectivity count"); 117 *oconn = conn; 118 *ooffsets = offsets; 119 *otypes = types; 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 PETSC_INTERN PetscErrorCode DMPlexGetNonEmptyComm_Private(DM dm, MPI_Comm *comm) 124 { 125 DM_Plex *mesh = (DM_Plex *)dm->data; 126 127 PetscFunctionBegin; 128 if (mesh->nonempty_comm == MPI_COMM_SELF) { /* Not yet setup */ 129 PetscInt cStart, cEnd, cellHeight; 130 MPI_Comm dmcomm = PetscObjectComm((PetscObject)dm); 131 PetscMPIInt color, rank; 132 133 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 134 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 135 color = (cStart < cEnd) ? 0 : 1; 136 PetscCallMPI(MPI_Comm_rank(dmcomm, &rank)); 137 PetscCallMPI(MPI_Comm_split(dmcomm, color, rank, &mesh->nonempty_comm)); 138 if (color == 1) { 139 PetscCallMPI(MPI_Comm_free(&mesh->nonempty_comm)); 140 mesh->nonempty_comm = MPI_COMM_NULL; 141 } 142 } 143 *comm = mesh->nonempty_comm; 144 PetscFunctionReturn(PETSC_SUCCESS); 145 } 146 147 /* 148 Write all fields that have been provided to the viewer 149 Multi-block XML format with binary appended data. 150 */ 151 PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm, PetscViewer viewer) 152 { 153 MPI_Comm comm; 154 PetscSection coordSection, cellCoordSection; 155 PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data; 156 PetscViewerVTKObjectLink link; 157 FILE *fp; 158 PetscMPIInt rank, size, tag; 159 PetscInt dimEmbed, cellHeight, cStart, cEnd, vStart, vEnd, numLabelCells, hasLabel, c, v, r, i; 160 PetscBool localized; 161 PieceInfo piece, *gpiece = NULL; 162 void *buffer = NULL; 163 const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian"; 164 PetscInt loops_per_scalar; 165 166 PetscFunctionBegin; 167 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 168 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 169 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 170 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 171 PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells)); 172 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 173 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 174 PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection)); 175 PetscCall(PetscCommGetNewTag(PetscObjectComm((PetscObject)dm), &tag)); 176 177 PetscCall(DMPlexGetNonEmptyComm_Private(dm, &comm)); 178 #if defined(PETSC_USE_COMPLEX) 179 loops_per_scalar = 2; 180 #else 181 loops_per_scalar = 1; 182 #endif 183 if (comm == MPI_COMM_NULL) goto finalize; 184 PetscCallMPI(MPI_Comm_size(comm, &size)); 185 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 186 187 PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp)); 188 PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n")); 189 PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order)); 190 PetscCall(PetscFPrintf(comm, fp, " <UnstructuredGrid>\n")); 191 192 hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE; 193 piece.nvertices = 0; 194 piece.ncells = 0; 195 piece.nconn = 0; 196 if (!localized) piece.nvertices = vEnd - vStart; 197 for (c = cStart; c < cEnd; ++c) { 198 PetscInt dof = 0; 199 if (hasLabel) { 200 PetscInt value; 201 202 PetscCall(DMGetLabelValue(dm, "vtk", c, &value)); 203 if (value != 1) continue; 204 } 205 if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof)); 206 if (!dof) { 207 PetscInt *closure = NULL; 208 PetscInt closureSize; 209 210 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 211 for (v = 0; v < closureSize * 2; v += 2) { 212 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 213 piece.nconn++; 214 if (localized) piece.nvertices++; 215 } 216 } 217 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 218 } else { 219 piece.nvertices += dof / dimEmbed; 220 piece.nconn += dof / dimEmbed; 221 } 222 piece.ncells++; 223 } 224 if (rank == 0) PetscCall(PetscMalloc1(size, &gpiece)); 225 PetscCallMPI(MPI_Gather((PetscInt *)&piece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, (PetscInt *)gpiece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, 0, comm)); 226 227 /* 228 * Write file header 229 */ 230 if (rank == 0) { 231 PetscInt64 boffset = 0; 232 233 for (r = 0; r < size; r++) { 234 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Piece NumberOfPoints=\"%" PetscInt_FMT "\" NumberOfCells=\"%" PetscInt_FMT "\">\n", gpiece[r].nvertices, gpiece[r].ncells)); 235 /* Coordinate positions */ 236 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Points>\n")); 237 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset)); 238 boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 239 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Points>\n")); 240 /* Cell connectivity */ 241 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <Cells>\n")); 242 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset)); 243 boffset += gpiece[r].nconn * sizeof(PetscVTKInt) + (gpiece[r].nconn ? sizeof(PetscInt64) : 0); 244 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset)); 245 boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 246 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset)); 247 boffset += gpiece[r].ncells * sizeof(unsigned char) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 248 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Cells>\n")); 249 250 /* 251 * Cell Data headers 252 */ 253 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <CellData>\n")); 254 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset)); 255 boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 256 /* all the vectors */ 257 for (link = vtk->link; link; link = link->next) { 258 Vec X = (Vec)link->vec; 259 DM dmX = NULL; 260 PetscInt bs = 1, nfields, field; 261 const char *vecname = ""; 262 PetscSection section; 263 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 264 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. */ 265 PetscCall(PetscObjectGetName((PetscObject)X, &vecname)); 266 } 267 PetscCall(VecGetDM(X, &dmX)); 268 if (!dmX) dmX = dm; 269 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 270 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 271 if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs)); 272 PetscCall(PetscSectionGetNumFields(section, &nfields)); 273 field = 0; 274 if (link->field >= 0) { 275 field = link->field; 276 nfields = field + 1; 277 } 278 for (i = 0; field < (nfields ? nfields : 1); field++) { 279 PetscInt fbs, j; 280 PetscFV fv = NULL; 281 PetscObject f; 282 PetscClassId fClass; 283 const char *fieldname = NULL; 284 char buf[256]; 285 PetscBool vector; 286 if (nfields) { /* We have user-defined fields/components */ 287 PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs)); 288 PetscCall(PetscSectionGetFieldName(section, field, &fieldname)); 289 } else fbs = bs; /* Say we have one field with 'bs' components */ 290 PetscCall(DMGetField(dmX, field, NULL, &f)); 291 PetscCall(PetscObjectGetClassId(f, &fClass)); 292 if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f; 293 if (nfields && !fieldname) { 294 PetscCall(PetscSNPrintf(buf, sizeof(buf), "CellField%" PetscInt_FMT, field)); 295 fieldname = buf; 296 } 297 vector = PETSC_FALSE; 298 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 299 vector = PETSC_TRUE; 300 PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Cell vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs); 301 for (j = 0; j < fbs; j++) { 302 const char *compName = NULL; 303 if (fv) { 304 PetscCall(PetscFVGetComponentName(fv, j, &compName)); 305 if (compName) break; 306 } 307 } 308 if (j < fbs) vector = PETSC_FALSE; 309 } 310 if (vector) { 311 #if defined(PETSC_USE_COMPLEX) 312 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 313 boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 314 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 315 boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 316 #else 317 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 318 boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 319 #endif 320 i += fbs; 321 } else { 322 for (j = 0; j < fbs; j++) { 323 const char *compName = NULL; 324 char finalname[256]; 325 if (fv) PetscCall(PetscFVGetComponentName(fv, j, &compName)); 326 if (compName) { 327 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName)); 328 } else if (fbs > 1) { 329 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%" PetscInt_FMT, vecname, fieldname, j)); 330 } else { 331 PetscCall(PetscSNPrintf(finalname, 255, "%s%s", vecname, fieldname)); 332 } 333 #if defined(PETSC_USE_COMPLEX) 334 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 335 boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 336 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 337 boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 338 #else 339 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 340 boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0); 341 #endif 342 i++; 343 } 344 } 345 } 346 //PetscCheck(i == bs,comm,PETSC_ERR_PLIB,"Total number of field components %" PetscInt_FMT " != block size %" PetscInt_FMT,i,bs); 347 } 348 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </CellData>\n")); 349 350 /* 351 * Point Data headers 352 */ 353 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " <PointData>\n")); 354 for (link = vtk->link; link; link = link->next) { 355 Vec X = (Vec)link->vec; 356 DM dmX; 357 PetscInt bs = 1, nfields, field; 358 const char *vecname = ""; 359 PetscSection section; 360 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 361 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. */ 362 PetscCall(PetscObjectGetName((PetscObject)X, &vecname)); 363 } 364 PetscCall(VecGetDM(X, &dmX)); 365 if (!dmX) dmX = dm; 366 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 367 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 368 if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs)); 369 PetscCall(PetscSectionGetNumFields(section, &nfields)); 370 field = 0; 371 if (link->field >= 0) { 372 field = link->field; 373 nfields = field + 1; 374 } 375 for (i = 0; field < (nfields ? nfields : 1); field++) { 376 PetscInt fbs, j; 377 const char *fieldname = NULL; 378 char buf[256]; 379 if (nfields) { /* We have user-defined fields/components */ 380 PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs)); 381 PetscCall(PetscSectionGetFieldName(section, field, &fieldname)); 382 } else fbs = bs; /* Say we have one field with 'bs' components */ 383 if (nfields && !fieldname) { 384 PetscCall(PetscSNPrintf(buf, sizeof(buf), "PointField%" PetscInt_FMT, field)); 385 fieldname = buf; 386 } 387 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 388 PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Point vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs); 389 #if defined(PETSC_USE_COMPLEX) 390 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 391 boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 392 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 393 boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 394 #else 395 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 396 boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 397 #endif 398 } else { 399 for (j = 0; j < fbs; j++) { 400 const char *compName = NULL; 401 char finalname[256]; 402 PetscCall(PetscSectionGetComponentName(section, field, j, &compName)); 403 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName)); 404 #if defined(PETSC_USE_COMPLEX) 405 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 406 boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 407 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 408 boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 409 #else 410 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset)); 411 boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0); 412 #endif 413 } 414 } 415 } 416 } 417 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </PointData>\n")); 418 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, " </Piece>\n")); 419 } 420 } 421 422 PetscCall(PetscFPrintf(comm, fp, " </UnstructuredGrid>\n")); 423 PetscCall(PetscFPrintf(comm, fp, " <AppendedData encoding=\"raw\">\n")); 424 PetscCall(PetscFPrintf(comm, fp, "_")); 425 426 if (rank == 0) { 427 PetscInt maxsize = 0; 428 for (r = 0; r < size; r++) { 429 maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nvertices * 3 * sizeof(PetscVTUReal))); 430 maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].ncells * 3 * sizeof(PetscVTUReal))); 431 maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nconn * sizeof(PetscVTKInt))); 432 } 433 PetscCall(PetscMalloc(maxsize, &buffer)); 434 } 435 for (r = 0; r < size; r++) { 436 if (r == rank) { 437 PetscInt nsend; 438 { /* Position */ 439 const PetscScalar *x, *cx = NULL; 440 PetscVTUReal *y = NULL; 441 Vec coords, cellCoords; 442 PetscBool copy; 443 444 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 445 PetscCall(VecGetArrayRead(coords, &x)); 446 PetscCall(DMGetCellCoordinatesLocal(dm, &cellCoords)); 447 if (cellCoords) PetscCall(VecGetArrayRead(cellCoords, &cx)); 448 #if defined(PETSC_USE_COMPLEX) 449 copy = PETSC_TRUE; 450 #else 451 copy = (PetscBool)(dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal))); 452 #endif 453 if (copy) { 454 PetscCall(PetscMalloc1(piece.nvertices * 3, &y)); 455 if (localized) { 456 PetscInt cnt; 457 for (c = cStart, cnt = 0; c < cEnd; c++) { 458 PetscInt off, dof; 459 460 PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof)); 461 if (!dof) { 462 PetscInt *closure = NULL; 463 PetscInt closureSize; 464 465 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 466 for (v = 0; v < closureSize * 2; v += 2) { 467 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 468 PetscCall(PetscSectionGetOffset(coordSection, closure[v], &off)); 469 if (dimEmbed != 3) { 470 y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]); 471 y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[off + 1]) : 0.0); 472 y[cnt * 3 + 2] = (PetscVTUReal)0.0; 473 } else { 474 y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]); 475 y[cnt * 3 + 1] = (PetscVTUReal)PetscRealPart(x[off + 1]); 476 y[cnt * 3 + 2] = (PetscVTUReal)PetscRealPart(x[off + 2]); 477 } 478 cnt++; 479 } 480 } 481 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 482 } else { 483 PetscCall(PetscSectionGetOffset(cellCoordSection, c, &off)); 484 if (dimEmbed != 3) { 485 for (i = 0; i < dof / dimEmbed; i++) { 486 y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(cx[off + i * dimEmbed + 0]); 487 y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(cx[off + i * dimEmbed + 1]) : 0.0); 488 y[cnt * 3 + 2] = (PetscVTUReal)0.0; 489 cnt++; 490 } 491 } else { 492 for (i = 0; i < dof; i++) y[cnt * 3 + i] = (PetscVTUReal)PetscRealPart(cx[off + i]); 493 cnt += dof / dimEmbed; 494 } 495 } 496 } 497 PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match"); 498 } else { 499 for (i = 0; i < piece.nvertices; i++) { 500 y[i * 3 + 0] = (PetscVTUReal)PetscRealPart(x[i * dimEmbed + 0]); 501 y[i * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[i * dimEmbed + 1]) : 0.); 502 y[i * 3 + 2] = (PetscVTUReal)((dimEmbed > 2) ? PetscRealPart(x[i * dimEmbed + 2]) : 0.); 503 } 504 } 505 } 506 nsend = piece.nvertices * 3; 507 PetscCall(TransferWrite(comm, viewer, fp, r, 0, copy ? (const void *)y : (const void *)x, buffer, nsend, MPIU_VTUREAL, tag)); 508 PetscCall(PetscFree(y)); 509 PetscCall(VecRestoreArrayRead(coords, &x)); 510 if (cellCoords) PetscCall(VecRestoreArrayRead(cellCoords, &cx)); 511 } 512 { /* Connectivity, offsets, types */ 513 PetscVTKInt *connectivity = NULL, *offsets = NULL; 514 PetscVTKType *types = NULL; 515 PetscCall(DMPlexGetVTKConnectivity(dm, localized, &piece, &connectivity, &offsets, &types)); 516 PetscCall(TransferWrite(comm, viewer, fp, r, 0, connectivity, buffer, piece.nconn, MPI_INT, tag)); 517 PetscCall(TransferWrite(comm, viewer, fp, r, 0, offsets, buffer, piece.ncells, MPI_INT, tag)); 518 PetscCall(TransferWrite(comm, viewer, fp, r, 0, types, buffer, piece.ncells, MPI_CHAR, tag)); 519 PetscCall(PetscFree3(connectivity, offsets, types)); 520 } 521 { /* Owners (cell data) */ 522 PetscVTKInt *owners; 523 PetscMPIInt orank; 524 525 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &orank)); 526 PetscCall(PetscMalloc1(piece.ncells, &owners)); 527 for (i = 0; i < piece.ncells; i++) owners[i] = orank; 528 PetscCall(TransferWrite(comm, viewer, fp, r, 0, owners, buffer, piece.ncells, MPI_INT, tag)); 529 PetscCall(PetscFree(owners)); 530 } 531 /* Cell data */ 532 for (link = vtk->link; link; link = link->next) { 533 Vec X = (Vec)link->vec; 534 DM dmX; 535 const PetscScalar *x; 536 PetscVTUReal *y; 537 PetscInt bs = 1, nfields, field; 538 PetscSection section = NULL; 539 540 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 541 PetscCall(VecGetDM(X, &dmX)); 542 if (!dmX) dmX = dm; 543 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 544 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 545 if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs)); 546 PetscCall(PetscSectionGetNumFields(section, &nfields)); 547 field = 0; 548 if (link->field >= 0) { 549 field = link->field; 550 nfields = field + 1; 551 } 552 PetscCall(VecGetArrayRead(X, &x)); 553 PetscCall(PetscMalloc1(piece.ncells * 3, &y)); 554 for (i = 0; field < (nfields ? nfields : 1); field++) { 555 PetscInt fbs, j; 556 PetscFV fv = NULL; 557 PetscObject f; 558 PetscClassId fClass; 559 PetscBool vector; 560 if (nfields && cEnd > cStart) { /* We have user-defined fields/components */ 561 PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs)); 562 } else fbs = bs; /* Say we have one field with 'bs' components */ 563 PetscCall(DMGetField(dmX, field, NULL, &f)); 564 PetscCall(PetscObjectGetClassId(f, &fClass)); 565 if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f; 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 (i = 0; 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 (i = 0; field < (nfields ? nfields : 1); field++) { 783 PetscInt fbs, j; 784 PetscFV fv = NULL; 785 PetscObject f; 786 PetscClassId fClass; 787 PetscBool vector; 788 if (nfields && cEnd > cStart) { /* We have user-defined fields/components */ 789 PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs)); 790 } else fbs = bs; /* Say we have one field with 'bs' components */ 791 PetscCall(DMGetField(dmX, field, NULL, &f)); 792 PetscCall(PetscObjectGetClassId(f, &fClass)); 793 if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f; 794 vector = PETSC_FALSE; 795 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 796 vector = PETSC_TRUE; 797 for (j = 0; j < fbs; j++) { 798 const char *compName = NULL; 799 if (fv) { 800 PetscCall(PetscFVGetComponentName(fv, j, &compName)); 801 if (compName) break; 802 } 803 } 804 if (j < fbs) vector = PETSC_FALSE; 805 } 806 if (vector) { 807 for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells * 3, MPIU_VTUREAL, tag)); 808 } else { 809 for (i = 0; i < fbs; i++) { 810 for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPIU_VTUREAL, tag)); 811 } 812 } 813 } 814 } 815 /* all point data */ 816 for (link = vtk->link; link; link = link->next) { 817 Vec X = (Vec)link->vec; 818 DM dmX; 819 PetscInt bs = 1, nfields, field; 820 PetscSection section = NULL; 821 822 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 823 PetscCall(VecGetDM(X, &dmX)); 824 if (!dmX) dmX = dm; 825 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion)); 826 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 827 if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs)); 828 PetscCall(PetscSectionGetNumFields(section, &nfields)); 829 field = 0; 830 if (link->field >= 0) { 831 field = link->field; 832 nfields = field + 1; 833 } 834 for (i = 0; field < (nfields ? nfields : 1); field++) { 835 PetscInt fbs; 836 if (nfields && vEnd > vStart) { /* We have user-defined fields/components */ 837 PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs)); 838 } else fbs = bs; /* Say we have one field with 'bs' components */ 839 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 840 for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag)); 841 } else { 842 for (i = 0; i < fbs; i++) { 843 for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices, MPIU_VTUREAL, tag)); 844 } 845 } 846 } 847 } 848 } 849 } 850 PetscCall(PetscFree(gpiece)); 851 PetscCall(PetscFree(buffer)); 852 PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n")); 853 PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n")); 854 PetscCall(PetscFClose(comm, fp)); 855 finalize: 856 /* this code sends to rank 0 that writes. 857 It may lead to very unbalanced log_view timings 858 of the next PETSc function logged. 859 Since this call is not performance critical, we 860 issue a barrier here to synchronize the processes */ 861 PetscCall(PetscBarrier((PetscObject)viewer)); 862 PetscFunctionReturn(PETSC_SUCCESS); 863 } 864