Actual source code: plexvtu.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>

  4: typedef struct {
  5:   PetscInt nvertices;
  6:   PetscInt ncells;
  7:   PetscInt nconn;               /* number of entries in cell->vertex connectivity array */
  8: } PieceInfo;

 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

 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;

 30:   MPI_Comm_rank(comm,&rank);
 31:   if (rank == srank && rank != root) {
 32:     MPI_Send((void*)send,count,mpidatatype,root,tag,comm);
 33:   } else if (rank == root) {
 34:     const void *buffer;
 35:     if (root == srank) {        /* self */
 36:       buffer = send;
 37:     } else {
 38:       MPI_Status  status;
 39:       PetscMPIInt nrecv;
 40:       MPI_Recv(recv,count,mpidatatype,srank,tag,comm,&status);
 41:       MPI_Get_count(&status,mpidatatype,&nrecv);
 43:       buffer = recv;
 44:     }
 45:     PetscViewerVTKFWrite(viewer,fp,buffer,count,mpidatatype);
 46:   }
 47:   return 0;
 48: }

 50: static PetscErrorCode DMPlexGetVTKConnectivity(DM dm, PetscBool localized, PieceInfo *piece,PetscVTKInt **oconn,PetscVTKInt **ooffsets,PetscVTKType **otypes)
 51: {
 52:   PetscSection   coordSection;
 53:   PetscVTKInt    *conn,*offsets;
 54:   PetscVTKType   *types;
 55:   PetscInt       dim,vStart,vEnd,cStart,cEnd,pStart,pEnd,cellHeight,numLabelCells,hasLabel,c,v,countcell,countconn;

 57:   PetscMalloc3(piece->nconn,&conn,piece->ncells,&offsets,piece->ncells,&types);
 58:   DMGetCoordinateSection(dm, &coordSection);
 59:   DMGetDimension(dm,&dim);
 60:   DMPlexGetChart(dm,&pStart,&pEnd);
 61:   DMPlexGetVTKCellHeight(dm, &cellHeight);
 62:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
 63:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 64:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
 65:   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;

 67:   countcell = 0;
 68:   countconn = 0;
 69:   for (c = cStart; c < cEnd; ++c) {
 70:     PetscInt nverts,dof = 0,celltype,startoffset,nC=0;

 72:     if (hasLabel) {
 73:       PetscInt value;

 75:       DMGetLabelValue(dm, "vtk", c, &value);
 76:       if (value != 1) continue;
 77:     }
 78:     startoffset = countconn;
 79:     if (localized) {
 80:       PetscSectionGetDof(coordSection, c, &dof);
 81:     }
 82:     if (!dof) {
 83:       PetscInt *closure = NULL;
 84:       PetscInt  closureSize;

 86:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
 87:       for (v = 0; v < closureSize*2; v += 2) {
 88:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
 89:           if (!localized) conn[countconn++] = closure[v] - vStart;
 90:           else conn[countconn++] = startoffset + nC;
 91:           ++nC;
 92:         }
 93:       }
 94:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
 95:     } else {
 96:       for (nC = 0; nC < dof/dim; nC++) conn[countconn++] = startoffset + nC;
 97:     }

 99:     {
100:       PetscInt n = PetscMin(nC, 8), s = countconn - nC, i, cone[8];
101:       for (i = 0; i < n; ++i) cone[i] = conn[s+i];
102:       DMPlexReorderCell(dm, c, cone);
103:       for (i = 0; i < n; ++i) conn[s+i] = (int)cone[i];
104:     }

106:     offsets[countcell] = countconn;

108:     nverts = countconn - startoffset;
109:     DMPlexVTKGetCellType_Internal(dm,dim,nverts,&celltype);

111:     types[countcell] = celltype;
112:     countcell++;
113:   }
116:   *oconn    = conn;
117:   *ooffsets = offsets;
118:   *otypes   = types;
119:   return 0;
120: }

122: /*
123:   Write all fields that have been provided to the viewer
124:   Multi-block XML format with binary appended data.
125: */
126: PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer)
127: {
128:   MPI_Comm                 comm;
129:   PetscSection             coordSection;
130:   PetscViewer_VTK          *vtk = (PetscViewer_VTK*)viewer->data;
131:   PetscViewerVTKObjectLink link;
132:   FILE                     *fp;
133:   PetscMPIInt              rank,size,tag;
134:   PetscInt                 dimEmbed,cellHeight,cStart,cEnd,vStart,vEnd,numLabelCells,hasLabel,c,v,r,i;
135:   PetscBool                localized;
136:   PieceInfo                piece,*gpiece = NULL;
137:   void                     *buffer = NULL;
138:   const char               *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
139:   PetscInt                 loops_per_scalar;

141:   PetscObjectGetComm((PetscObject)dm,&comm);
142: #if defined(PETSC_USE_COMPLEX)
143:   loops_per_scalar = 2;
144: #else
145:   loops_per_scalar = 1;
146: #endif
147:   MPI_Comm_size(comm,&size);
148:   MPI_Comm_rank(comm,&rank);
149:   PetscCommGetNewTag(comm,&tag);

151:   PetscFOpen(comm,vtk->filename,"wb",&fp);
152:   PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
153:   PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order);
154:   PetscFPrintf(comm,fp,"  <UnstructuredGrid>\n");

156:   DMGetCoordinateDim(dm, &dimEmbed);
157:   DMPlexGetVTKCellHeight(dm, &cellHeight);
158:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
159:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
160:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
161:   DMGetCoordinatesLocalized(dm, &localized);
162:   DMGetCoordinateSection(dm, &coordSection);

164:   hasLabel        = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
165:   piece.nvertices = 0;
166:   piece.ncells    = 0;
167:   piece.nconn     = 0;
168:   if (!localized) piece.nvertices = vEnd - vStart;
169:   for (c = cStart; c < cEnd; ++c) {
170:     PetscInt dof = 0;
171:     if (hasLabel) {
172:       PetscInt value;

174:       DMGetLabelValue(dm, "vtk", c, &value);
175:       if (value != 1) continue;
176:     }
177:     if (localized) {
178:       PetscSectionGetDof(coordSection, c, &dof);
179:     }
180:     if (!dof) {
181:       PetscInt *closure = NULL;
182:       PetscInt closureSize;

184:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
185:       for (v = 0; v < closureSize*2; v += 2) {
186:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
187:           piece.nconn++;
188:           if (localized) piece.nvertices++;
189:         }
190:       }
191:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
192:     } else {
193:       piece.nvertices += dof/dimEmbed;
194:       piece.nconn     += dof/dimEmbed;
195:     }
196:     piece.ncells++;
197:   }
198:   if (rank == 0) PetscMalloc1(size,&gpiece);
199:   MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);

201:   /*
202:    * Write file header
203:    */
204:   if (rank == 0) {
205:     PetscInt boffset = 0;

207:     for (r=0; r<size; r++) {
208:       PetscFPrintf(PETSC_COMM_SELF,fp,"    <Piece NumberOfPoints=\"%D\" NumberOfCells=\"%D\">\n",gpiece[r].nvertices,gpiece[r].ncells);
209:       /* Coordinate positions */
210:       PetscFPrintf(PETSC_COMM_SELF,fp,"      <Points>\n");
211:       PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
212:       boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
213:       PetscFPrintf(PETSC_COMM_SELF,fp,"      </Points>\n");
214:       /* Cell connectivity */
215:       PetscFPrintf(PETSC_COMM_SELF,fp,"      <Cells>\n");
216:       PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
217:       boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int);
218:       PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"offsets\"      NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
219:       boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int);
220:       PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"UInt8\" Name=\"types\"        NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
221:       boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int);
222:       PetscFPrintf(PETSC_COMM_SELF,fp,"      </Cells>\n");

224:       /*
225:        * Cell Data headers
226:        */
227:       PetscFPrintf(PETSC_COMM_SELF,fp,"      <CellData>\n");
228:       PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
229:       boffset += gpiece[r].ncells*sizeof(int) + sizeof(int);
230:       /* all the vectors */
231:       for (link=vtk->link; link; link=link->next) {
232:         Vec          X = (Vec)link->vec;
233:         DM           dmX = NULL;
234:         PetscInt     bs = 1,nfields,field;
235:         const char   *vecname = "";
236:         PetscSection section;
237:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
238:         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. */
239:           PetscObjectGetName((PetscObject)X,&vecname);
240:         }
241:         VecGetDM(X, &dmX);
242:         if (!dmX) dmX = dm;
243:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
244:         if (!section) DMGetLocalSection(dmX, &section);
245:         if (cEnd > cStart) PetscSectionGetDof(section,cStart,&bs);
246:         PetscSectionGetNumFields(section,&nfields);
247:         field = 0;
248:         if (link->field >= 0) {
249:           field = link->field;
250:           nfields = field + 1;
251:         }
252:         for (i=0; field<(nfields?nfields:1); field++) {
253:           PetscInt     fbs,j;
254:           PetscFV      fv = NULL;
255:           PetscObject  f;
256:           PetscClassId fClass;
257:           const char *fieldname = NULL;
258:           char       buf[256];
259:           PetscBool    vector;
260:           if (nfields) {        /* We have user-defined fields/components */
261:             PetscSectionGetFieldDof(section,cStart,field,&fbs);
262:             PetscSectionGetFieldName(section,field,&fieldname);
263:           } else fbs = bs;      /* Say we have one field with 'bs' components */
264:           DMGetField(dmX,field,NULL,&f);
265:           PetscObjectGetClassId(f,&fClass);
266:           if (fClass == PETSCFV_CLASSID) {
267:             fv = (PetscFV) f;
268:           }
269:           if (nfields && !fieldname) {
270:             PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);
271:             fieldname = buf;
272:           }
273:           vector = PETSC_FALSE;
274:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
275:             vector = PETSC_TRUE;
277:             for (j = 0; j < fbs; j++) {
278:               const char *compName = NULL;
279:               if (fv) {
280:                 PetscFVGetComponentName(fv,j,&compName);
281:                 if (compName) break;
282:               }
283:             }
284:             if (j < fbs) vector = PETSC_FALSE;
285:           }
286:           if (vector) {
287: #if defined(PETSC_USE_COMPLEX)
288:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
289:             boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
290:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
291:             boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
292: #else
293:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
294:             boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
295: #endif
296:             i+=fbs;
297:           } else {
298:             for (j=0; j<fbs; j++) {
299:               const char *compName = NULL;
300:               char finalname[256];
301:               if (fv) {
302:                 PetscFVGetComponentName(fv,j,&compName);
303:               }
304:               if (compName) {
305:                 PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName);
306:               }
307:               else if (fbs > 1) {
308:                 PetscSNPrintf(finalname,255,"%s%s.%D",vecname,fieldname,j);
309:               } else {
310:                 PetscSNPrintf(finalname,255,"%s%s",vecname,fieldname);
311:               }
312: #if defined(PETSC_USE_COMPLEX)
313:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
314:               boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
315:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
316:               boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
317: #else
318:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
319:               boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
320: #endif
321:               i++;
322:             }
323:           }
324:         }
326:       }
327:       PetscFPrintf(PETSC_COMM_SELF,fp,"      </CellData>\n");

329:       /*
330:        * Point Data headers
331:        */
332:       PetscFPrintf(PETSC_COMM_SELF,fp,"      <PointData>\n");
333:       for (link=vtk->link; link; link=link->next) {
334:         Vec          X = (Vec)link->vec;
335:         DM           dmX;
336:         PetscInt     bs = 1,nfields,field;
337:         const char   *vecname = "";
338:         PetscSection section;
339:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
340:         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. */
341:           PetscObjectGetName((PetscObject)X,&vecname);
342:         }
343:         VecGetDM(X, &dmX);
344:         if (!dmX) dmX = dm;
345:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
346:         if (!section) DMGetLocalSection(dmX, &section);
347:         if (vEnd > vStart) PetscSectionGetDof(section,vStart,&bs);
348:         PetscSectionGetNumFields(section,&nfields);
349:         field = 0;
350:         if (link->field >= 0) {
351:           field = link->field;
352:           nfields = field + 1;
353:         }
354:         for (i=0; field<(nfields?nfields:1); field++) {
355:           PetscInt   fbs,j;
356:           const char *fieldname = NULL;
357:           char       buf[256];
358:           if (nfields) {        /* We have user-defined fields/components */
359:             PetscSectionGetFieldDof(section,vStart,field,&fbs);
360:             PetscSectionGetFieldName(section,field,&fieldname);
361:           } else fbs = bs;      /* Say we have one field with 'bs' components */
362:           if (nfields && !fieldname) {
363:             PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);
364:             fieldname = buf;
365:           }
366:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
368: #if defined(PETSC_USE_COMPLEX)
369:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
370:             boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
371:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
372:             boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
373: #else
374:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
375:             boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
376: #endif
377:           } else {
378:             for (j=0; j<fbs; j++) {
379:               const char *compName = NULL;
380:               char finalname[256];
381:               PetscSectionGetComponentName(section,field,j,&compName);
382:               PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName);
383: #if defined(PETSC_USE_COMPLEX)
384:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
385:               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
386:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
387:               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
388: #else
389:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
390:               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
391: #endif
392:             }
393:           }
394:         }
395:       }
396:       PetscFPrintf(PETSC_COMM_SELF,fp,"      </PointData>\n");
397:       PetscFPrintf(PETSC_COMM_SELF,fp,"    </Piece>\n");
398:     }
399:   }

401:   PetscFPrintf(comm,fp,"  </UnstructuredGrid>\n");
402:   PetscFPrintf(comm,fp,"  <AppendedData encoding=\"raw\">\n");
403:   PetscFPrintf(comm,fp,"_");

405:   if (rank == 0) {
406:     PetscInt maxsize = 0;
407:     for (r=0; r<size; r++) {
408:       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscVTUReal)));
409:       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*3*sizeof(PetscVTUReal)));
410:       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt)));
411:     }
412:     PetscMalloc(maxsize,&buffer);
413:   }
414:   for (r=0; r<size; r++) {
415:     if (r == rank) {
416:       PetscInt nsend;
417:       {                         /* Position */
418:         const PetscScalar *x;
419:         PetscVTUReal      *y = NULL;
420:         Vec               coords;
421:         PetscBool         copy;

423:         DMGetCoordinatesLocal(dm,&coords);
424:         VecGetArrayRead(coords,&x);
425: #if defined(PETSC_USE_COMPLEX)
426:         copy = PETSC_TRUE;
427: #else
428:         copy = (PetscBool) (dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal)));
429: #endif
430:         if (copy) {
431:           PetscMalloc1(piece.nvertices*3,&y);
432:           if (localized) {
433:             PetscInt cnt;
434:             for (c=cStart,cnt=0; c<cEnd; c++) {
435:               PetscInt off, dof;

437:               PetscSectionGetDof(coordSection, c, &dof);
438:               if (!dof) {
439:                 PetscInt *closure = NULL;
440:                 PetscInt closureSize;

442:                 DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
443:                 for (v = 0; v < closureSize*2; v += 2) {
444:                   if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
445:                     PetscSectionGetOffset(coordSection, closure[v], &off);
446:                     if (dimEmbed != 3) {
447:                       y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off+0]);
448:                       y[cnt*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[off+1]) : 0.0);
449:                       y[cnt*3+2] = (PetscVTUReal) 0.0;
450:                     } else {
451:                       y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off+0]);
452:                       y[cnt*3+1] = (PetscVTUReal) PetscRealPart(x[off+1]);
453:                       y[cnt*3+2] = (PetscVTUReal) PetscRealPart(x[off+2]);
454:                     }
455:                     cnt++;
456:                   }
457:                 }
458:                 DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
459:               } else {
460:                 PetscSectionGetOffset(coordSection, c, &off);
461:                 if (dimEmbed != 3) {
462:                   for (i=0; i<dof/dimEmbed; i++) {
463:                     y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off + i*dimEmbed + 0]);
464:                     y[cnt*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[off + i*dimEmbed + 1]) : 0.0);
465:                     y[cnt*3+2] = (PetscVTUReal) 0.0;
466:                     cnt++;
467:                   }
468:                 } else {
469:                   for (i=0; i<dof; i ++) {
470:                     y[cnt*3+i] = (PetscVTUReal) PetscRealPart(x[off + i]);
471:                   }
472:                   cnt += dof/dimEmbed;
473:                 }
474:               }
475:             }
477:           } else {
478:             for (i=0; i<piece.nvertices; i++) {
479:               y[i*3+0] = (PetscVTUReal) PetscRealPart(x[i*dimEmbed+0]);
480:               y[i*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[i*dimEmbed+1]) : 0.);
481:               y[i*3+2] = (PetscVTUReal) ((dimEmbed > 2) ? PetscRealPart(x[i*dimEmbed+2]) : 0.);
482:             }
483:           }
484:         }
485:         nsend = piece.nvertices*3;
486:         TransferWrite(comm,viewer,fp,r,0,copy ? (const void *) y : (const void *) x,buffer,nsend,MPIU_VTUREAL,tag);
487:         PetscFree(y);
488:         VecRestoreArrayRead(coords,&x);
489:       }
490:       {                           /* Connectivity, offsets, types */
491:         PetscVTKInt  *connectivity = NULL, *offsets = NULL;
492:         PetscVTKType *types = NULL;
493:         DMPlexGetVTKConnectivity(dm,localized,&piece,&connectivity,&offsets,&types);
494:         TransferWrite(comm,viewer,fp,r,0,connectivity,buffer,piece.nconn,MPI_INT,tag);
495:         TransferWrite(comm,viewer,fp,r,0,offsets,buffer,piece.ncells,MPI_INT,tag);
496:         TransferWrite(comm,viewer,fp,r,0,types,buffer,piece.ncells,MPI_CHAR,tag);
497:         PetscFree3(connectivity,offsets,types);
498:       }
499:       {                         /* Owners (cell data) */
500:         PetscVTKInt *owners;
501:         PetscMalloc1(piece.ncells,&owners);
502:         for (i=0; i<piece.ncells; i++) owners[i] = rank;
503:         TransferWrite(comm,viewer,fp,r,0,owners,buffer,piece.ncells,MPI_INT,tag);
504:         PetscFree(owners);
505:       }
506:       /* Cell data */
507:       for (link=vtk->link; link; link=link->next) {
508:         Vec               X = (Vec)link->vec;
509:         DM                dmX;
510:         const PetscScalar *x;
511:         PetscVTUReal      *y;
512:         PetscInt          bs = 1, nfields, field;
513:         PetscSection      section = NULL;

515:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
516:         VecGetDM(X, &dmX);
517:         if (!dmX) dmX = dm;
518:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
519:         if (!section) DMGetLocalSection(dmX, &section);
520:         if (cEnd > cStart) PetscSectionGetDof(section,cStart,&bs);
521:         PetscSectionGetNumFields(section,&nfields);
522:         field = 0;
523:         if (link->field >= 0) {
524:           field = link->field;
525:           nfields = field + 1;
526:         }
527:         VecGetArrayRead(X,&x);
528:         PetscMalloc1(piece.ncells*3,&y);
529:         for (i=0; field<(nfields?nfields:1); field++) {
530:           PetscInt     fbs,j;
531:           PetscFV      fv = NULL;
532:           PetscObject  f;
533:           PetscClassId fClass;
534:           PetscBool    vector;
535:           if (nfields && cEnd > cStart) {        /* We have user-defined fields/components */
536:             PetscSectionGetFieldDof(section,cStart,field,&fbs);
537:           } else fbs = bs;      /* Say we have one field with 'bs' components */
538:           DMGetField(dmX,field,NULL,&f);
539:           PetscObjectGetClassId(f,&fClass);
540:           if (fClass == PETSCFV_CLASSID) {
541:             fv = (PetscFV) f;
542:           }
543:           vector = PETSC_FALSE;
544:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
545:             vector = PETSC_TRUE;
546:             for (j = 0; j < fbs; j++) {
547:               const char *compName = NULL;
548:               if (fv) {
549:                 PetscFVGetComponentName(fv,j,&compName);
550:                 if (compName) break;
551:               }
552:             }
553:             if (j < fbs) vector = PETSC_FALSE;
554:           }
555:           if (vector) {
556:             PetscInt cnt, l;
557:             for (l = 0; l < loops_per_scalar; l++) {
558:               for (c=cStart,cnt=0; c<cEnd; c++) {
559:                 const PetscScalar *xpoint;
560:                 PetscInt off, j;

562:                 if (hasLabel) {     /* Ignore some cells */
563:                   PetscInt value;
564:                   DMGetLabelValue(dmX, "vtk", c, &value);
565:                   if (value != 1) continue;
566:                 }
567:                 if (nfields) {
568:                   PetscSectionGetFieldOffset(section, c, field, &off);
569:                 } else {
570:                   PetscSectionGetOffset(section, c, &off);
571:                 }
572:                 xpoint = &x[off];
573:                 for (j = 0; j < fbs; j++) {
574:                   y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
575:                 }
576:                 for (; j < 3; j++) y[cnt++] = 0.;
577:               }
579:               TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.ncells*3,MPIU_VTUREAL,tag);
580:             }
581:           } else {
582:             for (i=0; i<fbs; i++) {
583:               PetscInt cnt, l;
584:               for (l = 0; l < loops_per_scalar; l++) {
585:                 for (c=cStart,cnt=0; c<cEnd; c++) {
586:                   const PetscScalar *xpoint;
587:                   PetscInt off;

589:                   if (hasLabel) {     /* Ignore some cells */
590:                     PetscInt value;
591:                     DMGetLabelValue(dmX, "vtk", c, &value);
592:                     if (value != 1) continue;
593:                   }
594:                   if (nfields) {
595:                     PetscSectionGetFieldOffset(section, c, field, &off);
596:                   } else {
597:                     PetscSectionGetOffset(section, c, &off);
598:                   }
599:                   xpoint   = &x[off];
600:                   y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
601:                 }
603:                 TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.ncells,MPIU_VTUREAL,tag);
604:               }
605:             }
606:           }
607:         }
608:         PetscFree(y);
609:         VecRestoreArrayRead(X,&x);
610:       }
611:       /* point data */
612:       for (link=vtk->link; link; link=link->next) {
613:         Vec               X = (Vec)link->vec;
614:         DM                dmX;
615:         const PetscScalar *x;
616:         PetscVTUReal      *y;
617:         PetscInt          bs = 1, nfields, field;
618:         PetscSection      section = NULL;

620:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
621:         VecGetDM(X, &dmX);
622:         if (!dmX) dmX = dm;
623:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
624:         if (!section) DMGetLocalSection(dmX, &section);
625:         if (vEnd > vStart) PetscSectionGetDof(section,vStart,&bs);
626:         PetscSectionGetNumFields(section,&nfields);
627:         field = 0;
628:         if (link->field >= 0) {
629:           field = link->field;
630:           nfields = field + 1;
631:         }
632:         VecGetArrayRead(X,&x);
633:         PetscMalloc1(piece.nvertices*3,&y);
634:         for (i=0; field<(nfields?nfields:1); field++) {
635:           PetscInt fbs,j;
636:           if (nfields && vEnd > vStart) {        /* We have user-defined fields/components */
637:             PetscSectionGetFieldDof(section,vStart,field,&fbs);
638:           } else fbs = bs;      /* Say we have one field with 'bs' components */
639:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
640:             PetscInt cnt, l;
641:             for (l = 0; l < loops_per_scalar; l++) {
642:               if (!localized) {
643:                 for (v=vStart,cnt=0; v<vEnd; v++) {
644:                   PetscInt          off;
645:                   const PetscScalar *xpoint;

647:                   if (nfields) {
648:                     PetscSectionGetFieldOffset(section,v,field,&off);
649:                   } else {
650:                     PetscSectionGetOffset(section,v,&off);
651:                   }
652:                   xpoint = &x[off];
653:                   for (j = 0; j < fbs; j++) {
654:                     y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
655:                   }
656:                   for (; j < 3; j++) y[cnt++] = 0.;
657:                 }
658:               } else {
659:                 for (c=cStart,cnt=0; c<cEnd; c++) {
660:                   PetscInt *closure = NULL;
661:                   PetscInt  closureSize, off;

663:                   DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
664:                   for (v = 0, off = 0; v < closureSize*2; v += 2) {
665:                     if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
666:                       PetscInt    voff;
667:                       const PetscScalar *xpoint;

669:                       if (nfields) {
670:                         PetscSectionGetFieldOffset(section,closure[v],field,&voff);
671:                       } else {
672:                         PetscSectionGetOffset(section,closure[v],&voff);
673:                       }
674:                       xpoint = &x[voff];
675:                       for (j = 0; j < fbs; j++) {
676:                         y[cnt + off++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
677:                       }
678:                       for (; j < 3; j++) y[cnt + off++] = 0.;
679:                     }
680:                   }
681:                   cnt += off;
682:                   DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
683:                 }
684:               }
686:               TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.nvertices*3,MPIU_VTUREAL,tag);
687:             }
688:           } else {
689:             for (i=0; i<fbs; i++) {
690:               PetscInt cnt, l;
691:               for (l = 0; l < loops_per_scalar; l++) {
692:                 if (!localized) {
693:                   for (v=vStart,cnt=0; v<vEnd; v++) {
694:                     PetscInt          off;
695:                     const PetscScalar *xpoint;

697:                     if (nfields) {
698:                       PetscSectionGetFieldOffset(section,v,field,&off);
699:                     } else {
700:                       PetscSectionGetOffset(section,v,&off);
701:                     }
702:                     xpoint   = &x[off];
703:                     y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
704:                   }
705:                 } else {
706:                   for (c=cStart,cnt=0; c<cEnd; c++) {
707:                     PetscInt *closure = NULL;
708:                     PetscInt  closureSize, off;

710:                     DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
711:                     for (v = 0, off = 0; v < closureSize*2; v += 2) {
712:                       if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
713:                         PetscInt    voff;
714:                         const PetscScalar *xpoint;

716:                         if (nfields) {
717:                           PetscSectionGetFieldOffset(section,closure[v],field,&voff);
718:                         } else {
719:                           PetscSectionGetOffset(section,closure[v],&voff);
720:                         }
721:                         xpoint         = &x[voff];
722:                         y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
723:                       }
724:                     }
725:                     cnt += off;
726:                     DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
727:                   }
728:                 }
730:                 TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.nvertices,MPIU_VTUREAL,tag);
731:               }
732:             }
733:           }
734:         }
735:         PetscFree(y);
736:         VecRestoreArrayRead(X,&x);
737:       }
738:     } else if (rank == 0) {
739:       PetscInt l;

741:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag); /* positions */
742:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,MPI_INT,tag); /* connectivity */
743:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag); /* offsets */
744:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_CHAR,tag); /* types */
745:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag); /* owner rank (cells) */
746:       /* all cell data */
747:       for (link=vtk->link; link; link=link->next) {
748:         Vec          X = (Vec)link->vec;
749:         PetscInt     bs = 1, nfields, field;
750:         DM           dmX;
751:         PetscSection section = NULL;

753:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
754:         VecGetDM(X, &dmX);
755:         if (!dmX) dmX = dm;
756:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
757:         if (!section) DMGetLocalSection(dmX, &section);
758:         if (cEnd > cStart) PetscSectionGetDof(section,cStart,&bs);
759:         PetscSectionGetNumFields(section,&nfields);
760:         field = 0;
761:         if (link->field >= 0) {
762:           field = link->field;
763:           nfields = field + 1;
764:         }
765:         for (i=0; field<(nfields?nfields:1); field++) {
766:           PetscInt     fbs,j;
767:           PetscFV      fv = NULL;
768:           PetscObject  f;
769:           PetscClassId fClass;
770:           PetscBool    vector;
771:           if (nfields && cEnd > cStart) {        /* We have user-defined fields/components */
772:             PetscSectionGetFieldDof(section,cStart,field,&fbs);
773:           } else fbs = bs;      /* Say we have one field with 'bs' components */
774:           DMGetField(dmX,field,NULL,&f);
775:           PetscObjectGetClassId(f,&fClass);
776:           if (fClass == PETSCFV_CLASSID) {
777:             fv = (PetscFV) f;
778:           }
779:           vector = PETSC_FALSE;
780:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
781:             vector = PETSC_TRUE;
782:             for (j = 0; j < fbs; j++) {
783:               const char *compName = NULL;
784:               if (fv) {
785:                 PetscFVGetComponentName(fv,j,&compName);
786:                 if (compName) break;
787:               }
788:             }
789:             if (j < fbs) vector = PETSC_FALSE;
790:           }
791:           if (vector) {
792:             for (l = 0; l < loops_per_scalar; l++) {
793:               TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells*3,MPIU_VTUREAL,tag);
794:             }
795:           } else {
796:             for (i=0; i<fbs; i++) {
797:               for (l = 0; l < loops_per_scalar; l++) {
798:                 TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPIU_VTUREAL,tag);
799:               }
800:             }
801:           }
802:         }
803:       }
804:       /* all point data */
805:       for (link=vtk->link; link; link=link->next) {
806:         Vec          X = (Vec)link->vec;
807:         DM           dmX;
808:         PetscInt     bs = 1, nfields, field;
809:         PetscSection section = NULL;

811:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
812:         VecGetDM(X, &dmX);
813:         if (!dmX) dmX = dm;
814:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
815:         if (!section) DMGetLocalSection(dmX, &section);
816:         if (vEnd > vStart) PetscSectionGetDof(section,vStart,&bs);
817:         PetscSectionGetNumFields(section,&nfields);
818:         field = 0;
819:         if (link->field >= 0) {
820:           field = link->field;
821:           nfields = field + 1;
822:         }
823:         for (i=0; field<(nfields?nfields:1); field++) {
824:           PetscInt   fbs;
825:           if (nfields && vEnd > vStart) {        /* We have user-defined fields/components */
826:             PetscSectionGetFieldDof(section,vStart,field,&fbs);
827:           } else fbs = bs;      /* Say we have one field with 'bs' components */
828:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
829:             for (l = 0; l < loops_per_scalar; l++) {
830:               TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag);
831:             }
832:           } else {
833:             for (i=0; i<fbs; i++) {
834:               for (l = 0; l < loops_per_scalar; l++) {
835:                 TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,MPIU_VTUREAL,tag);
836:               }
837:             }
838:           }
839:         }
840:       }
841:     }
842:   }
843:   PetscFree(gpiece);
844:   PetscFree(buffer);
845:   PetscFPrintf(comm,fp,"\n  </AppendedData>\n");
846:   PetscFPrintf(comm,fp,"</VTKFile>\n");
847:   PetscFClose(comm,fp);
848:   return 0;
849: }