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