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