Actual source code: swarmpic_view.c

petsc-3.8.4 2018-03-24
Report Typos and Errors
  1:  #include <petscdm.h>
  2:  #include <petscdmda.h>
  3:  #include <petscdmswarm.h>
  4:  #include <petsc/private/dmswarmimpl.h>
  5: #include "data_bucket.h"

  7: PetscErrorCode private_PetscViewerCreate_XDMF(MPI_Comm comm,const char filename[],PetscViewer *v)
  8: {
  9:   long int *bytes;
 10:   PetscContainer container;
 11:   PetscViewer viewer;
 13: 
 15:   PetscViewerCreate(comm,&viewer);
 16:   PetscViewerSetType(viewer,PETSCVIEWERASCII);
 17:   PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
 18:   PetscViewerFileSetName(viewer,filename);
 19: 
 20:   PetscMalloc1(1,&bytes);
 21:   bytes[0] = 0;
 22:   PetscContainerCreate(comm,&container);
 23:   PetscContainerSetPointer(container,(void*)bytes);
 24:   PetscObjectCompose((PetscObject)viewer,"XDMFViewerContext",(PetscObject)container);
 25: 
 26:   /* write xdmf header */
 27:   PetscViewerASCIIPrintf(viewer,"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n");
 28:   PetscViewerASCIIPrintf(viewer,"<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\" Version=\"2.99\">\n");
 29:   /* write xdmf domain */
 30:   PetscViewerASCIIPushTab(viewer);
 31:   PetscViewerASCIIPrintf(viewer,"<Domain>\n");
 32:   *v = viewer;
 33:   return(0);
 34: }

 36: PetscErrorCode private_PetscViewerDestroy_XDMF(PetscViewer *v)
 37: {
 38:   PetscViewer viewer;
 39:   DM dm = NULL;
 40:   long int *bytes;
 41:   PetscContainer container = NULL;
 43: 
 45:   if (!v) return(0);
 46:   viewer = *v;
 47: 
 48:   PetscObjectQuery((PetscObject)viewer,"DMSwarm",(PetscObject*)&dm);
 49:   if (dm) {
 50:     PetscViewerASCIIPrintf(viewer,"</Grid>\n");
 51:     PetscViewerASCIIPopTab(viewer);
 52:   }
 53: 
 54:   /* close xdmf header */
 55:   PetscViewerASCIIPrintf(viewer,"</Domain>\n");
 56:   PetscViewerASCIIPopTab(viewer);
 57:   PetscViewerASCIIPrintf(viewer,"</Xdmf>\n");
 58: 
 59:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
 60:   if (container) {
 61:     PetscContainerGetPointer(container,(void**)&bytes);
 62:     PetscFree(bytes);
 63:     PetscContainerDestroy(&container);
 64:   }
 65: 
 66:   PetscViewerDestroy(&viewer);
 67:   *v = NULL;
 68:   return(0);
 69: }

 71: PetscErrorCode private_CreateDataFileNameXDMF(const char filename[],char dfilename[])
 72: {
 73:   char *ext;
 74:   PetscBool flg;
 76: 
 78:   PetscStrrchr(filename,'.',&ext);
 79:   PetscStrcmp("xmf",ext,&flg);
 80:   if (flg) {
 81:     size_t len;
 82:     char viewername_minus_ext[PETSC_MAX_PATH_LEN];
 83: 
 84:     PetscStrlen(filename,&len);
 85:     PetscSNPrintf(viewername_minus_ext,len-3,"%s",filename);
 86:     PetscSNPrintf(dfilename,PETSC_MAX_PATH_LEN-1,"%s_swarm_fields.pbin",viewername_minus_ext);
 87:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"File extension must by .xmf");

 89:   return(0);
 90: }

 92: PetscErrorCode private_DMSwarmView_XDMF(DM dm,PetscViewer viewer)
 93: {
 94:   PetscBool isswarm = PETSC_FALSE;
 95:   const char *viewername;
 96:   char datafile[PETSC_MAX_PATH_LEN];
 97:   PetscViewer fviewer;
 98:   PetscInt k,ng,dim;
 99:   Vec dvec;
100:   long int *bytes = NULL;
101:   PetscContainer container = NULL;
102:   const char *dmname;
104: 
106:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
107:   if (container) {
108:     PetscContainerGetPointer(container,(void**)&bytes);
109:   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");
110: 
111:   PetscObjectTypeCompare((PetscObject)dm,DMSWARM,&isswarm);
112:   if (!isswarm) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Only valid for DMSwarm");
113: 
114:   PetscObjectCompose((PetscObject)viewer,"DMSwarm",(PetscObject)dm);
115: 
116:   PetscViewerASCIIPushTab(viewer);
117:   PetscObjectGetName((PetscObject)dm,&dmname);
118:   if (!dmname) {
119:     DMGetOptionsPrefix(dm,&dmname);
120:   }
121:   if (!dmname) {
122:     PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n");
123:   } else {
124:     PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n",dmname);
125:   }
126: 
127:   /* create a sub-viewer for topology, geometry and all data fields */
128:   /* name is viewer.name + "_swarm_fields.pbin" */
129:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
130:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
131:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
132:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
133:   PetscViewerFileSetMode(fviewer,FILE_MODE_WRITE);
134: 
135:   PetscViewerFileGetName(viewer,&viewername);
136:   private_CreateDataFileNameXDMF(viewername,datafile);
137:   PetscViewerFileSetName(fviewer,datafile);
138: 
139:   DMSwarmGetSize(dm,&ng);
140: 
141:   /* write topology header */
142:   PetscViewerASCIIPushTab(viewer);
143:   PetscViewerASCIIPrintf(viewer,"<Topology Dimensions=\"%D\" TopologyType=\"Mixed\">\n",ng);
144:   PetscViewerASCIIPushTab(viewer);
145:   PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Dimensions=\"%D\" Seek=\"%D\">\n",ng*3,bytes[0]);
146:   PetscViewerASCIIPushTab(viewer);
147:   PetscViewerASCIIPrintf(viewer,"%s\n",datafile);
148:   PetscViewerASCIIPopTab(viewer);
149:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
150:   PetscViewerASCIIPopTab(viewer);
151:   PetscViewerASCIIPrintf(viewer,"</Topology>\n");
152:   PetscViewerASCIIPopTab(viewer);
153: 
154:   /* write topology data */
155:   for (k=0; k<ng; k++) {
156:     PetscInt pvertex[3];
157: 
158:     pvertex[0] = 1;
159:     pvertex[1] = 1;
160:     pvertex[2] = k;
161:     PetscViewerBinaryWrite(fviewer,pvertex,3,PETSC_INT,PETSC_FALSE);
162:   }
163:   bytes[0] += sizeof(PetscInt) * ng * 3;
164: 
165:   /* write geometry header */
166:   PetscViewerASCIIPushTab(viewer);
167:   DMGetDimension(dm,&dim);
168:   switch (dim) {
169:     case 1:
170:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for 1D");
171:       break;
172:     case 2:
173:       PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XY\">\n");
174:       break;
175:     case 3:
176:       PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XYZ\">\n");
177:       break;
178:   }
179:   PetscViewerASCIIPushTab(viewer);
180:   PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",ng,dim,bytes[0]);
181:   PetscViewerASCIIPushTab(viewer);
182:   PetscViewerASCIIPrintf(viewer,"%s\n",datafile);
183:   PetscViewerASCIIPopTab(viewer);
184:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
185:   PetscViewerASCIIPopTab(viewer);
186:   PetscViewerASCIIPrintf(viewer,"</Geometry>\n");
187:   PetscViewerASCIIPopTab(viewer);
188: 
189:   /* write geometry data */
190:   DMSwarmCreateGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);
191:   VecView(dvec,fviewer);
192:   DMSwarmDestroyGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);
193:   bytes[0] += sizeof(PetscReal) * ng * dim;
194: 
195:   PetscViewerDestroy(&fviewer);
196: 
197:   return(0);
198: }

200: PetscErrorCode private_VecView_Swarm_XDMF(Vec x,PetscViewer viewer)
201: {
202:   long int *bytes = NULL;
203:   PetscContainer container = NULL;
204:   const char *viewername;
205:   char datafile[PETSC_MAX_PATH_LEN];
206:   PetscViewer fviewer;
207:   PetscInt N,bs;
208:   const char *vecname;
209:   char fieldname[PETSC_MAX_PATH_LEN];
211: 
213:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
214:   if (container) {
215:     PetscContainerGetPointer(container,(void**)&bytes);
216:   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");
217: 
218:   PetscViewerFileGetName(viewer,&viewername);
219:   private_CreateDataFileNameXDMF(viewername,datafile);
220: 
221:   /* re-open a sub-viewer for all data fields */
222:   /* name is viewer.name + "_swarm_fields.pbin" */
223:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
224:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
225:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
226:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
227:   PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);
228:   PetscViewerFileSetName(fviewer,datafile);
229: 
230:   VecGetSize(x,&N);
231:   VecGetBlockSize(x,&bs);
232:   N = N/bs;
233:   PetscObjectGetName((PetscObject)x,&vecname);
234:   if (!vecname) {
235:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)x)->tag);
236:   } else {
237:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);
238:   }
239: 
240:   /* write data header */
241:   PetscViewerASCIIPushTab(viewer);
242:   PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);
243:   PetscViewerASCIIPushTab(viewer);
244:   if (bs == 1) {
245:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);
246:   } else {
247:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);
248:   }
249:   PetscViewerASCIIPushTab(viewer);
250:   PetscViewerASCIIPrintf(viewer,"%s\n",datafile);
251:   PetscViewerASCIIPopTab(viewer);
252:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
253:   PetscViewerASCIIPopTab(viewer);
254:   PetscViewerASCIIPrintf(viewer,"</Attribute>\n");
255:   PetscViewerASCIIPopTab(viewer);
256: 
257:   /* write data */
258:   VecView(x,fviewer);
259:   bytes[0] += sizeof(PetscReal) * N * bs;
260: 
261:   PetscViewerDestroy(&fviewer);
262: 
263:   return(0);
264: }

266: /*@C
267:    DMSwarmViewFieldsXDMF - Write a selection of DMSwarm fields to an XDMF3 file
268:  
269:    Collective on DM
270:  
271:    Input parameters:
272: +  dm - the DMSwarm
273: .  filename - the file name of the XDMF file (must have the extension .xmf)
274: .  nfields - the number of fields to write into the XDMF file
275: -  field_name_list - array of length nfields containing the textual name of fields to write
276:  
277:    Level: beginner

279:    Notes:
280:    Only fields registered with data type PETSC_REAL can be written into the file
281:  
282: .seealso: DMSwarmViewXDMF()
283: @*/
284: PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm,const char filename[],PetscInt nfields,const char *field_name_list[])
285: {
287:   Vec dvec;
288:   PetscInt f;
289:   PetscViewer viewer;
290: 
292:   private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);
293:   private_DMSwarmView_XDMF(dm,viewer);
294:   for (f=0; f<nfields; f++) {
295:     DMSwarmCreateGlobalVectorFromField(dm,field_name_list[f],&dvec);
296:     PetscObjectSetName((PetscObject)dvec,field_name_list[f]);
297:     private_VecView_Swarm_XDMF(dvec,viewer);
298:     DMSwarmDestroyGlobalVectorFromField(dm,field_name_list[f],&dvec);
299:   }
300:   private_PetscViewerDestroy_XDMF(&viewer);
301:   return(0);
302: }

304: /*@C
305:    DMSwarmViewXDMF - Write DMSwarm fields to an XDMF3 file
306:  
307:    Collective on DM
308:  
309:    Input parameters:
310: +  dm - the DMSwarm
311: -  filename - the file name of the XDMF file (must have the extension .xmf)
312:  
313:    Level: beginner

315:    Notes:
316:    Only fields registered with data type PETSC_REAL will be written into the file

318: .seealso: DMSwarmViewFieldsXDMF()
319: @*/
320: PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm,const char filename[])
321: {
322:   DM_Swarm *swarm = (DM_Swarm*)dm->data;
324:   Vec dvec;
325:   PetscInt f;
326:   PetscViewer viewer;
327: 
329:   private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);
330:   private_DMSwarmView_XDMF(dm,viewer);
331:   for (f=4; f<swarm->db->nfields; f++) { /* only examine user defined fields - the first 4 are internally created by DMSwarmPIC */
332:     DataField field;
333: 
334:     /* query field type - accept all those of type PETSC_DOUBLE */
335:     field = swarm->db->field[f];
336:     if (field->petsc_type != PETSC_DOUBLE) continue;
337: 
338:     DMSwarmCreateGlobalVectorFromField(dm,field->name,&dvec);
339:     PetscObjectSetName((PetscObject)dvec,field->name);
340:     private_VecView_Swarm_XDMF(dvec,viewer);
341:     DMSwarmDestroyGlobalVectorFromField(dm,field->name,&dvec);
342:   }
343:   private_PetscViewerDestroy_XDMF(&viewer);
344:   return(0);
345: }