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*) §ion);
244: if (!section) DMGetLocalSection(dmX, §ion);
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*) §ion);
346: if (!section) DMGetLocalSection(dmX, §ion);
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*) §ion);
519: if (!section) DMGetLocalSection(dmX, §ion);
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*) §ion);
624: if (!section) DMGetLocalSection(dmX, §ion);
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*) §ion);
757: if (!section) DMGetLocalSection(dmX, §ion);
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*) §ion);
815: if (!section) DMGetLocalSection(dmX, §ion);
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: }