Actual source code: swarmpic_view.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  1:  #include <petscdm.h>
  2:  #include <petscdmda.h>
  3:  #include <petscdmswarm.h>
  4:  #include <petsc/private/dmswarmimpl.h>
  5: #include "../src/dm/impls/swarm/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:     PetscStrncpy(viewername_minus_ext,filename,len-2);
 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: PetscErrorCode private_ISView_Swarm_XDMF(IS is,PetscViewer viewer)
267: {
268:   long int       *bytes = NULL;
269:   PetscContainer container = NULL;
270:   const char     *viewername;
271:   char           datafile[PETSC_MAX_PATH_LEN];
272:   PetscViewer    fviewer;
273:   PetscInt       N,bs;
274:   const char     *vecname;
275:   char           fieldname[PETSC_MAX_PATH_LEN];
277: 
279:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
280:   if (container) {
281:     PetscContainerGetPointer(container,(void**)&bytes);
282:   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");
283: 
284:   PetscViewerFileGetName(viewer,&viewername);
285:   private_CreateDataFileNameXDMF(viewername,datafile);
286: 
287:   /* re-open a sub-viewer for all data fields */
288:   /* name is viewer.name + "_swarm_fields.pbin" */
289:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
290:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
291:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
292:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
293:   PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);
294:   PetscViewerFileSetName(fviewer,datafile);
295: 
296:   ISGetSize(is,&N);
297:   ISGetBlockSize(is,&bs);
298:   N = N/bs;
299:   PetscObjectGetName((PetscObject)is,&vecname);
300:   if (!vecname) {
301:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)is)->tag);
302:   } else {
303:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);
304:   }
305: 
306:   /* write data header */
307:   PetscViewerASCIIPushTab(viewer);
308:   PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);
309:   PetscViewerASCIIPushTab(viewer);
310:   if (bs == 1) {
311:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);
312:   } else {
313:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);
314:   }
315:   PetscViewerASCIIPushTab(viewer);
316:   PetscViewerASCIIPrintf(viewer,"%s\n",datafile);
317:   PetscViewerASCIIPopTab(viewer);
318:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
319:   PetscViewerASCIIPopTab(viewer);
320:   PetscViewerASCIIPrintf(viewer,"</Attribute>\n");
321:   PetscViewerASCIIPopTab(viewer);
322: 
323:   /* write data */
324:   ISView(is,fviewer);
325:   bytes[0] += sizeof(PetscInt) * N * bs;
326: 
327:   PetscViewerDestroy(&fviewer);
328: 
329:   return(0);
330: }

332: /*@C
333:    DMSwarmViewFieldsXDMF - Write a selection of DMSwarm fields to an XDMF3 file
334:  
335:    Collective on DM
336:  
337:    Input parameters:
338: +  dm - the DMSwarm
339: .  filename - the file name of the XDMF file (must have the extension .xmf)
340: .  nfields - the number of fields to write into the XDMF file
341: -  field_name_list - array of length nfields containing the textual name of fields to write
342:  
343:    Level: beginner

345:    Notes:
346:    Only fields registered with data type PETSC_DOUBLE or PETSC_INT can be written into the file
347:  
348: .seealso: DMSwarmViewXDMF()
349: @*/
350: PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm,const char filename[],PetscInt nfields,const char *field_name_list[])
351: {
353:   Vec            dvec;
354:   PetscInt       f,N;
355:   PetscViewer    viewer;
356: 
358:   private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);
359:   private_DMSwarmView_XDMF(dm,viewer);
360:   DMSwarmGetLocalSize(dm,&N);
361:   for (f=0; f<nfields; f++) {
362:     void          *data;
363:     PetscDataType type;

365:     DMSwarmGetField(dm,field_name_list[f],NULL,&type,&data);
366:     DMSwarmRestoreField(dm,field_name_list[f],NULL,&type,&data);

368:     if (type == PETSC_DOUBLE) {
369:       DMSwarmCreateGlobalVectorFromField(dm,field_name_list[f],&dvec);
370:       PetscObjectSetName((PetscObject)dvec,field_name_list[f]);
371:       private_VecView_Swarm_XDMF(dvec,viewer);
372:       DMSwarmDestroyGlobalVectorFromField(dm,field_name_list[f],&dvec);
373:     } else if (type == PETSC_INT) {
374:       IS is;
375:       const PetscInt *idx;

377:       DMSwarmGetField(dm,field_name_list[f],NULL,&type,&data);
378:       idx = (const PetscInt*)data;
379: 
380:       ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);
381:       PetscObjectSetName((PetscObject)is,field_name_list[f]);
382:       private_ISView_Swarm_XDMF(is,viewer);
383:       ISDestroy(&is);
384:       DMSwarmRestoreField(dm,field_name_list[f],NULL,&type,&data);
385:     } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Can only write PETSC_INT and PETSC_DOUBLE");

387:   }
388:   private_PetscViewerDestroy_XDMF(&viewer);
389:   return(0);
390: }

392: /*@C
393:    DMSwarmViewXDMF - Write DMSwarm fields to an XDMF3 file
394:  
395:    Collective on DM
396:  
397:    Input parameters:
398: +  dm - the DMSwarm
399: -  filename - the file name of the XDMF file (must have the extension .xmf)
400:  
401:    Level: beginner

403:    Notes:
404:    Only fields user registered with data type PETSC_DOUBLE or PETSC_INT will be written into the file

406: .seealso: DMSwarmViewFieldsXDMF()
407: @*/
408: PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm,const char filename[])
409: {
410:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
412:   Vec            dvec;
413:   PetscInt       f;
414:   PetscViewer    viewer;
415: 
417:   private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);
418:   private_DMSwarmView_XDMF(dm,viewer);
419:   for (f=4; f<swarm->db->nfields; f++) { /* only examine user defined fields - the first 4 are internally created by DMSwarmPIC */
420:     DMSwarmDataField field;
421: 
422:     /* query field type - accept all those of type PETSC_DOUBLE */
423:     field = swarm->db->field[f];
424:     if (field->petsc_type == PETSC_DOUBLE) {
425:       DMSwarmCreateGlobalVectorFromField(dm,field->name,&dvec);
426:       PetscObjectSetName((PetscObject)dvec,field->name);
427:       private_VecView_Swarm_XDMF(dvec,viewer);
428:       DMSwarmDestroyGlobalVectorFromField(dm,field->name,&dvec);
429:     } else if (field->petsc_type == PETSC_INT) {
430:       IS             is;
431:       PetscInt       N;
432:       const PetscInt *idx;
433:       void           *data;
434: 
435:       DMSwarmGetLocalSize(dm,&N);
436:       DMSwarmGetField(dm,field->name,NULL,NULL,&data);
437:       idx = (const PetscInt*)data;
438: 
439:       ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);
440:       PetscObjectSetName((PetscObject)is,field->name);
441:       private_ISView_Swarm_XDMF(is,viewer);
442:       ISDestroy(&is);
443:       DMSwarmRestoreField(dm,field->name,NULL,NULL,&data);
444:     }
445:   }
446:   private_PetscViewerDestroy_XDMF(&viewer);
447:   return(0);
448: }