Actual source code: plexvtu.c

petsc-main 2021-04-20
Report Typos and Errors
  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;

 32:   MPI_Comm_rank(comm,&rank);
 33:   if (rank == srank && rank != root) {
 34:     MPI_Send((void*)send,count,mpidatatype,root,tag,comm);
 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:       MPI_Recv(recv,count,mpidatatype,srank,tag,comm,&status);
 43:       MPI_Get_count(&status,mpidatatype,&nrecv);
 44:       if (count != nrecv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
 45:       buffer = recv;
 46:     }
 47:     PetscViewerVTKFWrite(viewer,fp,buffer,count,mpidatatype);
 48:   }
 49:   return(0);
 50: }

 52: static PetscErrorCode DMPlexGetVTKConnectivity(DM dm, PetscBool localized, PieceInfo *piece,PetscVTKInt **oconn,PetscVTKInt **ooffsets,PetscVTKType **otypes)
 53: {
 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;

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

 71:   countcell = 0;
 72:   countconn = 0;
 73:   for (c = cStart; c < cEnd; ++c) {
 74:     PetscInt nverts,dof = 0,celltype,startoffset,nC=0;

 76:     if (hasLabel) {
 77:       PetscInt value;

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

 90:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
 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:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
 99:     } else {
100:       for (nC = 0; nC < dof/dim; nC++) conn[countconn++] = startoffset + nC;
101:     }

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:       DMPlexReorderCell(dm, c, cone);
107:       for (i = 0; i < n; ++i) conn[s+i] = (int)cone[i];
108:     }

110:     offsets[countcell] = countconn;

112:     nverts = countconn - startoffset;
113:     DMPlexVTKGetCellType_Internal(dm,dim,nverts,&celltype);

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:   return(0);
124: }

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;

147:   PetscObjectGetComm((PetscObject)dm,&comm);
148: #if defined(PETSC_USE_COMPLEX)
149:   loops_per_scalar = 2;
150: #else
151:   loops_per_scalar = 1;
152: #endif
153:   MPI_Comm_size(comm,&size);
154:   MPI_Comm_rank(comm,&rank);
155:   PetscCommGetNewTag(comm,&tag);

157:   PetscFOpen(comm,vtk->filename,"wb",&fp);
158:   PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
159:   PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order);
160:   PetscFPrintf(comm,fp,"  <UnstructuredGrid>\n");

162:   DMGetCoordinateDim(dm, &dimEmbed);
163:   DMPlexGetVTKCellHeight(dm, &cellHeight);
164:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
165:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
166:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
167:   DMGetCoordinatesLocalized(dm, &localized);
168:   DMGetCoordinateSection(dm, &coordSection);

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;

180:       DMGetLabelValue(dm, "vtk", c, &value);
181:       if (value != 1) continue;
182:     }
183:     if (localized) {
184:       PetscSectionGetDof(coordSection, c, &dof);
185:     }
186:     if (!dof) {
187:       PetscInt *closure = NULL;
188:       PetscInt closureSize;

190:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
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:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
198:     } else {
199:       piece.nvertices += dof/dimEmbed;
200:       piece.nconn     += dof/dimEmbed;
201:     }
202:     piece.ncells++;
203:   }
204:   if (!rank) {PetscMalloc1(size,&gpiece);}
205:   MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);

207:   /*
208:    * Write file header
209:    */
210:   if (!rank) {
211:     PetscInt boffset = 0;

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

230:       /*
231:        * Cell Data headers
232:        */
233:       PetscFPrintf(PETSC_COMM_SELF,fp,"      <CellData>\n");
234:       PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
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:           PetscObjectGetName((PetscObject)X,&vecname);
246:         }
247:         VecGetDM(X, &dmX);
248:         if (!dmX) dmX = dm;
249:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
250:         if (!section) {DMGetLocalSection(dmX, &section);}
251:         if (cEnd > cStart) { PetscSectionGetDof(section,cStart,&bs); }
252:         PetscSectionGetNumFields(section,&nfields);
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:             PetscSectionGetFieldDof(section,cStart,field,&fbs);
268:             PetscSectionGetFieldName(section,field,&fieldname);
269:           } else fbs = bs;      /* Say we have one field with 'bs' components */
270:           DMGetField(dmX,field,NULL,&f);
271:           PetscObjectGetClassId(f,&fClass);
272:           if (fClass == PETSCFV_CLASSID) {
273:             fv = (PetscFV) f;
274:           }
275:           if (nfields && !fieldname) {
276:             PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);
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:                 PetscFVGetComponentName(fv,j,&compName);
287:                 if (compName) break;
288:               }
289:             }
290:             if (j < fbs) vector = PETSC_FALSE;
291:           }
292:           if (vector) {
293: #if defined(PETSC_USE_COMPLEX)
294:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
295:             boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
296:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
297:             boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
298: #else
299:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
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:                 PetscFVGetComponentName(fv,j,&compName);
309:               }
310:               if (compName) {
311:                 PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName);
312:               }
313:               else if (fbs > 1) {
314:                 PetscSNPrintf(finalname,255,"%s%s.%D",vecname,fieldname,j);
315:               } else {
316:                 PetscSNPrintf(finalname,255,"%s%s",vecname,fieldname);
317:               }
318: #if defined(PETSC_USE_COMPLEX)
319:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
320:               boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
321:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
322:               boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
323: #else
324:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
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:       PetscFPrintf(PETSC_COMM_SELF,fp,"      </CellData>\n");

335:       /*
336:        * Point Data headers
337:        */
338:       PetscFPrintf(PETSC_COMM_SELF,fp,"      <PointData>\n");
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:           PetscObjectGetName((PetscObject)X,&vecname);
348:         }
349:         VecGetDM(X, &dmX);
350:         if (!dmX) dmX = dm;
351:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
352:         if (!section) {DMGetLocalSection(dmX, &section);}
353:         if (vEnd > vStart) { PetscSectionGetDof(section,vStart,&bs); }
354:         PetscSectionGetNumFields(section,&nfields);
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:             PetscSectionGetFieldDof(section,vStart,field,&fbs);
366:             PetscSectionGetFieldName(section,field,&fieldname);
367:           } else fbs = bs;      /* Say we have one field with 'bs' components */
368:           if (nfields && !fieldname) {
369:             PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);
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:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
376:             boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
377:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
378:             boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
379: #else
380:             PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
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:               PetscSectionGetComponentName(section,field,j,&compName);
388:               PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName);
389: #if defined(PETSC_USE_COMPLEX)
390:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
391:               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
392:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
393:               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
394: #else
395:               PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
396:               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
397: #endif
398:             }
399:           }
400:         }
401:       }
402:       PetscFPrintf(PETSC_COMM_SELF,fp,"      </PointData>\n");
403:       PetscFPrintf(PETSC_COMM_SELF,fp,"    </Piece>\n");
404:     }
405:   }

407:   PetscFPrintf(comm,fp,"  </UnstructuredGrid>\n");
408:   PetscFPrintf(comm,fp,"  <AppendedData encoding=\"raw\">\n");
409:   PetscFPrintf(comm,fp,"_");

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:     PetscMalloc(maxsize,&buffer);
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;

429:         DMGetCoordinatesLocal(dm,&coords);
430:         VecGetArrayRead(coords,&x);
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:           PetscMalloc1(piece.nvertices*3,&y);
438:           if (localized) {
439:             PetscInt cnt;
440:             for (c=cStart,cnt=0; c<cEnd; c++) {
441:               PetscInt off, dof;

443:               PetscSectionGetDof(coordSection, c, &dof);
444:               if (!dof) {
445:                 PetscInt *closure = NULL;
446:                 PetscInt closureSize;

448:                 DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
449:                 for (v = 0; v < closureSize*2; v += 2) {
450:                   if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
451:                     PetscSectionGetOffset(coordSection, closure[v], &off);
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:                 DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
465:               } else {
466:                 PetscSectionGetOffset(coordSection, c, &off);
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:         TransferWrite(comm,viewer,fp,r,0,copy ? (const void *) y : (const void *) x,buffer,nsend,MPIU_VTUREAL,tag);
493:         PetscFree(y);
494:         VecRestoreArrayRead(coords,&x);
495:       }
496:       {                           /* Connectivity, offsets, types */
497:         PetscVTKInt  *connectivity = NULL, *offsets = NULL;
498:         PetscVTKType *types = NULL;
499:         DMPlexGetVTKConnectivity(dm,localized,&piece,&connectivity,&offsets,&types);
500:         TransferWrite(comm,viewer,fp,r,0,connectivity,buffer,piece.nconn,MPI_INT,tag);
501:         TransferWrite(comm,viewer,fp,r,0,offsets,buffer,piece.ncells,MPI_INT,tag);
502:         TransferWrite(comm,viewer,fp,r,0,types,buffer,piece.ncells,MPI_CHAR,tag);
503:         PetscFree3(connectivity,offsets,types);
504:       }
505:       {                         /* Owners (cell data) */
506:         PetscVTKInt *owners;
507:         PetscMalloc1(piece.ncells,&owners);
508:         for (i=0; i<piece.ncells; i++) owners[i] = rank;
509:         TransferWrite(comm,viewer,fp,r,0,owners,buffer,piece.ncells,MPI_INT,tag);
510:         PetscFree(owners);
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;

521:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
522:         VecGetDM(X, &dmX);
523:         if (!dmX) dmX = dm;
524:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
525:         if (!section) {DMGetLocalSection(dmX, &section);}
526:         if (cEnd > cStart) { PetscSectionGetDof(section,cStart,&bs); }
527:         PetscSectionGetNumFields(section,&nfields);
528:         field = 0;
529:         if (link->field >= 0) {
530:           field = link->field;
531:           nfields = field + 1;
532:         }
533:         VecGetArrayRead(X,&x);
534:         PetscMalloc1(piece.ncells*3,&y);
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:             PetscSectionGetFieldDof(section,cStart,field,&fbs);
543:           } else fbs = bs;      /* Say we have one field with 'bs' components */
544:           DMGetField(dmX,field,NULL,&f);
545:           PetscObjectGetClassId(f,&fClass);
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:                 PetscFVGetComponentName(fv,j,&compName);
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;

568:                 if (hasLabel) {     /* Ignore some cells */
569:                   PetscInt value;
570:                   DMGetLabelValue(dmX, "vtk", c, &value);
571:                   if (value != 1) continue;
572:                 }
573:                 if (nfields) {
574:                   PetscSectionGetFieldOffset(section, c, field, &off);
575:                 } else {
576:                   PetscSectionGetOffset(section, c, &off);
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:               TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.ncells*3,MPIU_VTUREAL,tag);
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;

595:                   if (hasLabel) {     /* Ignore some cells */
596:                     PetscInt value;
597:                     DMGetLabelValue(dmX, "vtk", c, &value);
598:                     if (value != 1) continue;
599:                   }
600:                   if (nfields) {
601:                     PetscSectionGetFieldOffset(section, c, field, &off);
602:                   } else {
603:                     PetscSectionGetOffset(section, c, &off);
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:                 TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.ncells,MPIU_VTUREAL,tag);
610:               }
611:             }
612:           }
613:         }
614:         PetscFree(y);
615:         VecRestoreArrayRead(X,&x);
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;

626:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
627:         VecGetDM(X, &dmX);
628:         if (!dmX) dmX = dm;
629:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
630:         if (!section) {DMGetLocalSection(dmX, &section);}
631:         if (vEnd > vStart) { PetscSectionGetDof(section,vStart,&bs); }
632:         PetscSectionGetNumFields(section,&nfields);
633:         field = 0;
634:         if (link->field >= 0) {
635:           field = link->field;
636:           nfields = field + 1;
637:         }
638:         VecGetArrayRead(X,&x);
639:         PetscMalloc1(piece.nvertices*3,&y);
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:             PetscSectionGetFieldDof(section,vStart,field,&fbs);
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;

653:                   if (nfields) {
654:                     PetscSectionGetFieldOffset(section,v,field,&off);
655:                   } else {
656:                     PetscSectionGetOffset(section,v,&off);
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;

669:                   DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
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;

675:                       if (nfields) {
676:                         PetscSectionGetFieldOffset(section,closure[v],field,&voff);
677:                       } else {
678:                         PetscSectionGetOffset(section,closure[v],&voff);
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:                   DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
689:                 }
690:               }
691:               if (cnt != piece.nvertices*3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
692:               TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.nvertices*3,MPIU_VTUREAL,tag);
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;

703:                     if (nfields) {
704:                       PetscSectionGetFieldOffset(section,v,field,&off);
705:                     } else {
706:                       PetscSectionGetOffset(section,v,&off);
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;

716:                     DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
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;

722:                         if (nfields) {
723:                           PetscSectionGetFieldOffset(section,closure[v],field,&voff);
724:                         } else {
725:                           PetscSectionGetOffset(section,closure[v],&voff);
726:                         }
727:                         xpoint         = &x[voff];
728:                         y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
729:                       }
730:                     }
731:                     cnt += off;
732:                     DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
733:                   }
734:                 }
735:                 if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
736:                 TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.nvertices,MPIU_VTUREAL,tag);
737:               }
738:             }
739:           }
740:         }
741:         PetscFree(y);
742:         VecRestoreArrayRead(X,&x);
743:       }
744:     } else if (!rank) {
745:       PetscInt l;

747:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag); /* positions */
748:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,MPI_INT,tag); /* connectivity */
749:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag); /* offsets */
750:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_CHAR,tag); /* types */
751:       TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag); /* 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;

759:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
760:         VecGetDM(X, &dmX);
761:         if (!dmX) dmX = dm;
762:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
763:         if (!section) {DMGetLocalSection(dmX, &section);}
764:         if (cEnd > cStart) { PetscSectionGetDof(section,cStart,&bs); }
765:         PetscSectionGetNumFields(section,&nfields);
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:             PetscSectionGetFieldDof(section,cStart,field,&fbs);
779:           } else fbs = bs;      /* Say we have one field with 'bs' components */
780:           DMGetField(dmX,field,NULL,&f);
781:           PetscObjectGetClassId(f,&fClass);
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:                 PetscFVGetComponentName(fv,j,&compName);
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:               TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells*3,MPIU_VTUREAL,tag);
800:             }
801:           } else {
802:             for (i=0; i<fbs; i++) {
803:               for (l = 0; l < loops_per_scalar; l++) {
804:                 TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPIU_VTUREAL,tag);
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;

817:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
818:         VecGetDM(X, &dmX);
819:         if (!dmX) dmX = dm;
820:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
821:         if (!section) {DMGetLocalSection(dmX, &section);}
822:         if (vEnd > vStart) { PetscSectionGetDof(section,vStart,&bs); }
823:         PetscSectionGetNumFields(section,&nfields);
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:             PetscSectionGetFieldDof(section,vStart,field,&fbs);
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:               TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag);
837:             }
838:           } else {
839:             for (i=0; i<fbs; i++) {
840:               for (l = 0; l < loops_per_scalar; l++) {
841:                 TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,MPIU_VTUREAL,tag);
842:               }
843:             }
844:           }
845:         }
846:       }
847:     }
848:   }
849:   PetscFree(gpiece);
850:   PetscFree(buffer);
851:   PetscFPrintf(comm,fp,"\n  </AppendedData>\n");
852:   PetscFPrintf(comm,fp,"</VTKFile>\n");
853:   PetscFClose(comm,fp);
854:   return(0);
855: }