Actual source code: plexvtu.c
petsc-3.14.6 2021-03-30
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*) §ion);
250: if (!section) {DMGetLocalSection(dmX, §ion);}
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*) §ion);
352: if (!section) {DMGetLocalSection(dmX, §ion);}
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*) §ion);
525: if (!section) {DMGetLocalSection(dmX, §ion);}
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*) §ion);
630: if (!section) {DMGetLocalSection(dmX, §ion);}
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*) §ion);
763: if (!section) {DMGetLocalSection(dmX, §ion);}
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*) §ion);
821: if (!section) {DMGetLocalSection(dmX, §ion);}
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: }