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