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