Actual source code: swarmpic_view.c

  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;

 15:   PetscViewerCreate(comm,&viewer);
 16:   PetscViewerSetType(viewer,PETSCVIEWERASCII);
 17:   PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
 18:   PetscViewerFileSetName(viewer,filename);

 20:   PetscMalloc1(1,&bytes);
 21:   bytes[0] = 0;
 22:   PetscContainerCreate(comm,&container);
 23:   PetscContainerSetPointer(container,(void*)bytes);
 24:   PetscObjectCompose((PetscObject)viewer,"XDMFViewerContext",(PetscObject)container);

 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;

 45:   if (!v) return(0);
 46:   viewer = *v;

 48:   PetscObjectQuery((PetscObject)viewer,"DMSwarm",(PetscObject*)&dm);
 49:   if (dm) {
 50:     PetscViewerASCIIPrintf(viewer,"</Grid>\n");
 51:     PetscViewerASCIIPopTab(viewer);
 52:   }

 54:   /* close xdmf header */
 55:   PetscViewerASCIIPrintf(viewer,"</Domain>\n");
 56:   PetscViewerASCIIPopTab(viewer);
 57:   PetscViewerASCIIPrintf(viewer,"</Xdmf>\n");

 59:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
 60:   if (container) {
 61:     PetscContainerGetPointer(container,(void**)&bytes);
 62:     PetscFree(bytes);
 63:     PetscContainerDestroy(&container);
 64:   }

 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;

 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];

 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:   char           *datafilename;
 98:   PetscViewer    fviewer;
 99:   PetscInt       k,ng,dim;
100:   Vec            dvec;
101:   long int       *bytes = NULL;
102:   PetscContainer container = NULL;
103:   const char     *dmname;

107:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
108:   if (container) {
109:     PetscContainerGetPointer(container,(void**)&bytes);
110:   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");

112:   PetscObjectTypeCompare((PetscObject)dm,DMSWARM,&isswarm);
113:   if (!isswarm) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Only valid for DMSwarm");

115:   PetscObjectCompose((PetscObject)viewer,"DMSwarm",(PetscObject)dm);

117:   PetscViewerASCIIPushTab(viewer);
118:   PetscObjectGetName((PetscObject)dm,&dmname);
119:   if (!dmname) {
120:     DMGetOptionsPrefix(dm,&dmname);
121:   }
122:   if (!dmname) {
123:     PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n");
124:   } else {
125:     PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n",dmname);
126:   }

128:   /* create a sub-viewer for topology, geometry and all data fields */
129:   /* name is viewer.name + "_swarm_fields.pbin" */
130:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
131:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
132:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
133:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
134:   PetscViewerFileSetMode(fviewer,FILE_MODE_WRITE);

136:   PetscViewerFileGetName(viewer,&viewername);
137:   private_CreateDataFileNameXDMF(viewername,datafile);
138:   PetscViewerFileSetName(fviewer,datafile);
139:   PetscStrrchr(datafile,'/',&datafilename);

141:   DMSwarmGetSize(dm,&ng);

143:   /* write topology header */
144:   PetscViewerASCIIPushTab(viewer);
145:   PetscViewerASCIIPrintf(viewer,"<Topology Dimensions=\"%D\" TopologyType=\"Mixed\">\n",ng);
146:   PetscViewerASCIIPushTab(viewer);
147:   PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Dimensions=\"%D\" Seek=\"%D\">\n",ng*3,bytes[0]);
148:   PetscViewerASCIIPushTab(viewer);
149:   PetscViewerASCIIPrintf(viewer,"%s\n",datafilename);
150:   PetscViewerASCIIPopTab(viewer);
151:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
152:   PetscViewerASCIIPopTab(viewer);
153:   PetscViewerASCIIPrintf(viewer,"</Topology>\n");
154:   PetscViewerASCIIPopTab(viewer);

156:   /* write topology data */
157:   for (k=0; k<ng; k++) {
158:     PetscInt pvertex[3];

160:     pvertex[0] = 1;
161:     pvertex[1] = 1;
162:     pvertex[2] = k;
163:     PetscViewerBinaryWrite(fviewer,pvertex,3,PETSC_INT);
164:   }
165:   bytes[0] += sizeof(PetscInt) * ng * 3;

167:   /* write geometry header */
168:   PetscViewerASCIIPushTab(viewer);
169:   DMGetDimension(dm,&dim);
170:   switch (dim) {
171:     case 1:
172:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for 1D");
173:     case 2:
174:       PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XY\">\n");
175:       break;
176:     case 3:
177:       PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XYZ\">\n");
178:       break;
179:   }
180:   PetscViewerASCIIPushTab(viewer);
181:   PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",ng,dim,bytes[0]);
182:   PetscViewerASCIIPushTab(viewer);
183:   PetscViewerASCIIPrintf(viewer,"%s\n",datafilename);
184:   PetscViewerASCIIPopTab(viewer);
185:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
186:   PetscViewerASCIIPopTab(viewer);
187:   PetscViewerASCIIPrintf(viewer,"</Geometry>\n");
188:   PetscViewerASCIIPopTab(viewer);

190:   /* write geometry data */
191:   DMSwarmCreateGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);
192:   VecView(dvec,fviewer);
193:   DMSwarmDestroyGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);
194:   bytes[0] += sizeof(PetscReal) * ng * dim;

196:   PetscViewerDestroy(&fviewer);

198:   return(0);
199: }

201: PetscErrorCode private_VecView_Swarm_XDMF(Vec x,PetscViewer viewer)
202: {
203:   long int       *bytes = NULL;
204:   PetscContainer container = NULL;
205:   const char     *viewername;
206:   char           datafile[PETSC_MAX_PATH_LEN];
207:   char           *datafilename;
208:   PetscViewer    fviewer;
209:   PetscInt       N,bs;
210:   const char     *vecname;
211:   char           fieldname[PETSC_MAX_PATH_LEN];

215:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
216:   if (container) {
217:     PetscContainerGetPointer(container,(void**)&bytes);
218:   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");

220:   PetscViewerFileGetName(viewer,&viewername);
221:   private_CreateDataFileNameXDMF(viewername,datafile);

223:   /* re-open a sub-viewer for all data fields */
224:   /* name is viewer.name + "_swarm_fields.pbin" */
225:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
226:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
227:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
228:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
229:   PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);
230:   PetscViewerFileSetName(fviewer,datafile);
231:   PetscStrrchr(datafile,'/',&datafilename);

233:   VecGetSize(x,&N);
234:   VecGetBlockSize(x,&bs);
235:   N = N/bs;
236:   PetscObjectGetName((PetscObject)x,&vecname);
237:   if (!vecname) {
238:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)x)->tag);
239:   } else {
240:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);
241:   }

243:   /* write data header */
244:   PetscViewerASCIIPushTab(viewer);
245:   PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);
246:   PetscViewerASCIIPushTab(viewer);
247:   if (bs == 1) {
248:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);
249:   } else {
250:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);
251:   }
252:   PetscViewerASCIIPushTab(viewer);
253:   PetscViewerASCIIPrintf(viewer,"%s\n",datafilename);
254:   PetscViewerASCIIPopTab(viewer);
255:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
256:   PetscViewerASCIIPopTab(viewer);
257:   PetscViewerASCIIPrintf(viewer,"</Attribute>\n");
258:   PetscViewerASCIIPopTab(viewer);

260:   /* write data */
261:   VecView(x,fviewer);
262:   bytes[0] += sizeof(PetscReal) * N * bs;

264:   PetscViewerDestroy(&fviewer);

266:   return(0);
267: }

269: PetscErrorCode private_ISView_Swarm_XDMF(IS is,PetscViewer viewer)
270: {
271:   long int       *bytes = NULL;
272:   PetscContainer container = NULL;
273:   const char     *viewername;
274:   char           datafile[PETSC_MAX_PATH_LEN];
275:   char           *datafilename;
276:   PetscViewer    fviewer;
277:   PetscInt       N,bs;
278:   const char     *vecname;
279:   char           fieldname[PETSC_MAX_PATH_LEN];

283:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
284:   if (container) {
285:     PetscContainerGetPointer(container,(void**)&bytes);
286:   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");

288:   PetscViewerFileGetName(viewer,&viewername);
289:   private_CreateDataFileNameXDMF(viewername,datafile);

291:   /* re-open a sub-viewer for all data fields */
292:   /* name is viewer.name + "_swarm_fields.pbin" */
293:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
294:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
295:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
296:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
297:   PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);
298:   PetscViewerFileSetName(fviewer,datafile);
299:   PetscStrrchr(datafile,'/',&datafilename);

301:   ISGetSize(is,&N);
302:   ISGetBlockSize(is,&bs);
303:   N = N/bs;
304:   PetscObjectGetName((PetscObject)is,&vecname);
305:   if (!vecname) {
306:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)is)->tag);
307:   } else {
308:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);
309:   }

311:   /* write data header */
312:   PetscViewerASCIIPushTab(viewer);
313:   PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);
314:   PetscViewerASCIIPushTab(viewer);
315:   if (bs == 1) {
316:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);
317:   } else {
318:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);
319:   }
320:   PetscViewerASCIIPushTab(viewer);
321:   PetscViewerASCIIPrintf(viewer,"%s\n",datafilename);
322:   PetscViewerASCIIPopTab(viewer);
323:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
324:   PetscViewerASCIIPopTab(viewer);
325:   PetscViewerASCIIPrintf(viewer,"</Attribute>\n");
326:   PetscViewerASCIIPopTab(viewer);

328:   /* write data */
329:   ISView(is,fviewer);
330:   bytes[0] += sizeof(PetscInt) * N * bs;

332:   PetscViewerDestroy(&fviewer);

334:   return(0);
335: }

337: /*@C
338:    DMSwarmViewFieldsXDMF - Write a selection of DMSwarm fields to an XDMF3 file

340:    Collective on dm

342:    Input parameters:
343: +  dm - the DMSwarm
344: .  filename - the file name of the XDMF file (must have the extension .xmf)
345: .  nfields - the number of fields to write into the XDMF file
346: -  field_name_list - array of length nfields containing the textual name of fields to write

348:    Level: beginner

350:    Notes:
351:    Only fields registered with data type PETSC_DOUBLE or PETSC_INT can be written into the file

353: .seealso: DMSwarmViewXDMF()
354: @*/
355: PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm,const char filename[],PetscInt nfields,const char *field_name_list[])
356: {
358:   Vec            dvec;
359:   PetscInt       f,N;
360:   PetscViewer    viewer;

363:   private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);
364:   private_DMSwarmView_XDMF(dm,viewer);
365:   DMSwarmGetLocalSize(dm,&N);
366:   for (f=0; f<nfields; f++) {
367:     void          *data;
368:     PetscDataType type;

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

373:     if (type == PETSC_DOUBLE) {
374:       DMSwarmCreateGlobalVectorFromField(dm,field_name_list[f],&dvec);
375:       PetscObjectSetName((PetscObject)dvec,field_name_list[f]);
376:       private_VecView_Swarm_XDMF(dvec,viewer);
377:       DMSwarmDestroyGlobalVectorFromField(dm,field_name_list[f],&dvec);
378:     } else if (type == PETSC_INT) {
379:       IS is;
380:       const PetscInt *idx;

382:       DMSwarmGetField(dm,field_name_list[f],NULL,&type,&data);
383:       idx = (const PetscInt*)data;

385:       ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);
386:       PetscObjectSetName((PetscObject)is,field_name_list[f]);
387:       private_ISView_Swarm_XDMF(is,viewer);
388:       ISDestroy(&is);
389:       DMSwarmRestoreField(dm,field_name_list[f],NULL,&type,&data);
390:     } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Can only write PETSC_INT and PETSC_DOUBLE");

392:   }
393:   private_PetscViewerDestroy_XDMF(&viewer);
394:   return(0);
395: }

397: /*@C
398:    DMSwarmViewXDMF - Write DMSwarm fields to an XDMF3 file

400:    Collective on dm

402:    Input parameters:
403: +  dm - the DMSwarm
404: -  filename - the file name of the XDMF file (must have the extension .xmf)

406:    Level: beginner

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

411: .seealso: DMSwarmViewFieldsXDMF()
412: @*/
413: PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm,const char filename[])
414: {
415:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
417:   Vec            dvec;
418:   PetscInt       f;
419:   PetscViewer    viewer;

422:   private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);
423:   private_DMSwarmView_XDMF(dm,viewer);
424:   for (f=4; f<swarm->db->nfields; f++) { /* only examine user defined fields - the first 4 are internally created by DMSwarmPIC */
425:     DMSwarmDataField field;

427:     /* query field type - accept all those of type PETSC_DOUBLE */
428:     field = swarm->db->field[f];
429:     if (field->petsc_type == PETSC_DOUBLE) {
430:       DMSwarmCreateGlobalVectorFromField(dm,field->name,&dvec);
431:       PetscObjectSetName((PetscObject)dvec,field->name);
432:       private_VecView_Swarm_XDMF(dvec,viewer);
433:       DMSwarmDestroyGlobalVectorFromField(dm,field->name,&dvec);
434:     } else if (field->petsc_type == PETSC_INT) {
435:       IS             is;
436:       PetscInt       N;
437:       const PetscInt *idx;
438:       void           *data;

440:       DMSwarmGetLocalSize(dm,&N);
441:       DMSwarmGetField(dm,field->name,NULL,NULL,&data);
442:       idx = (const PetscInt*)data;

444:       ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);
445:       PetscObjectSetName((PetscObject)is,field->name);
446:       private_ISView_Swarm_XDMF(is,viewer);
447:       ISDestroy(&is);
448:       DMSwarmRestoreField(dm,field->name,NULL,NULL,&data);
449:     }
450:   }
451:   private_PetscViewerDestroy_XDMF(&viewer);
452:   return(0);
453: }