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;

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

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

 24:   /* write xdmf header */
 25:   PetscViewerASCIIPrintf(viewer,"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n");
 26:   PetscViewerASCIIPrintf(viewer,"<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude/\" Version=\"2.99\">\n");
 27:   /* write xdmf domain */
 28:   PetscViewerASCIIPushTab(viewer);
 29:   PetscViewerASCIIPrintf(viewer,"<Domain>\n");
 30:   *v = viewer;
 31:   return 0;
 32: }

 34: PetscErrorCode private_PetscViewerDestroy_XDMF(PetscViewer *v)
 35: {
 36:   PetscViewer    viewer;
 37:   DM             dm = NULL;
 38:   long int       *bytes;
 39:   PetscContainer container = NULL;

 41:   if (!v) return 0;
 42:   viewer = *v;

 44:   PetscObjectQuery((PetscObject)viewer,"DMSwarm",(PetscObject*)&dm);
 45:   if (dm) {
 46:     PetscViewerASCIIPrintf(viewer,"</Grid>\n");
 47:     PetscViewerASCIIPopTab(viewer);
 48:   }

 50:   /* close xdmf header */
 51:   PetscViewerASCIIPrintf(viewer,"</Domain>\n");
 52:   PetscViewerASCIIPopTab(viewer);
 53:   PetscViewerASCIIPrintf(viewer,"</Xdmf>\n");

 55:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
 56:   if (container) {
 57:     PetscContainerGetPointer(container,(void**)&bytes);
 58:     PetscFree(bytes);
 59:     PetscContainerDestroy(&container);
 60:   }
 61:   PetscViewerDestroy(&viewer);
 62:   *v = NULL;
 63:   return 0;
 64: }

 66: PetscErrorCode private_CreateDataFileNameXDMF(const char filename[],char dfilename[])
 67: {
 68:   char           *ext;
 69:   PetscBool      flg;

 71:   PetscStrrchr(filename,'.',&ext);
 72:   PetscStrcmp("xmf",ext,&flg);
 73:   if (flg) {
 74:     size_t len;
 75:     char    viewername_minus_ext[PETSC_MAX_PATH_LEN];

 77:     PetscStrlen(filename,&len);
 78:     PetscStrncpy(viewername_minus_ext,filename,len-2);
 79:     PetscSNPrintf(dfilename,PETSC_MAX_PATH_LEN-1,"%s_swarm_fields.pbin",viewername_minus_ext);
 80:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"File extension must by .xmf");
 81:   return 0;
 82: }

 84: PetscErrorCode private_DMSwarmView_XDMF(DM dm,PetscViewer viewer)
 85: {
 86:   PetscBool      isswarm = PETSC_FALSE;
 87:   const char     *viewername;
 88:   char           datafile[PETSC_MAX_PATH_LEN];
 89:   char           *datafilename;
 90:   PetscViewer    fviewer;
 91:   PetscInt       k,ng,dim;
 92:   Vec            dvec;
 93:   long int       *bytes = NULL;
 94:   PetscContainer container = NULL;
 95:   const char     *dmname;

 97:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
 98:   if (container) {
 99:     PetscContainerGetPointer(container,(void**)&bytes);
100:   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");

102:   PetscObjectTypeCompare((PetscObject)dm,DMSWARM,&isswarm);

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

107:   PetscViewerASCIIPushTab(viewer);
108:   PetscObjectGetName((PetscObject)dm,&dmname);
109:   if (!dmname) {
110:     DMGetOptionsPrefix(dm,&dmname);
111:   }
112:   if (!dmname) {
113:     PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n");
114:   } else {
115:     PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n",dmname);
116:   }

118:   /* create a sub-viewer for topology, geometry and all data fields */
119:   /* name is viewer.name + "_swarm_fields.pbin" */
120:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
121:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
122:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
123:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
124:   PetscViewerFileSetMode(fviewer,FILE_MODE_WRITE);

126:   PetscViewerFileGetName(viewer,&viewername);
127:   private_CreateDataFileNameXDMF(viewername,datafile);
128:   PetscViewerFileSetName(fviewer,datafile);
129:   PetscStrrchr(datafile,'/',&datafilename);

131:   DMSwarmGetSize(dm,&ng);

133:   /* write topology header */
134:   PetscViewerASCIIPushTab(viewer);
135:   PetscViewerASCIIPrintf(viewer,"<Topology Dimensions=\"%D\" TopologyType=\"Mixed\">\n",ng);
136:   PetscViewerASCIIPushTab(viewer);
137:   PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Dimensions=\"%D\" Seek=\"%D\">\n",ng*3,bytes[0]);
138:   PetscViewerASCIIPushTab(viewer);
139:   PetscViewerASCIIPrintf(viewer,"%s\n",datafilename);
140:   PetscViewerASCIIPopTab(viewer);
141:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
142:   PetscViewerASCIIPopTab(viewer);
143:   PetscViewerASCIIPrintf(viewer,"</Topology>\n");
144:   PetscViewerASCIIPopTab(viewer);

146:   /* write topology data */
147:   for (k=0; k<ng; k++) {
148:     PetscInt pvertex[3];

150:     pvertex[0] = 1;
151:     pvertex[1] = 1;
152:     pvertex[2] = k;
153:     PetscViewerBinaryWrite(fviewer,pvertex,3,PETSC_INT);
154:   }
155:   bytes[0] += sizeof(PetscInt) * ng * 3;

157:   /* write geometry header */
158:   PetscViewerASCIIPushTab(viewer);
159:   DMGetDimension(dm,&dim);
160:   switch (dim) {
161:     case 1:
162:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for 1D");
163:     case 2:
164:       PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XY\">\n");
165:       break;
166:     case 3:
167:       PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XYZ\">\n");
168:       break;
169:   }
170:   PetscViewerASCIIPushTab(viewer);
171:   PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",ng,dim,bytes[0]);
172:   PetscViewerASCIIPushTab(viewer);
173:   PetscViewerASCIIPrintf(viewer,"%s\n",datafilename);
174:   PetscViewerASCIIPopTab(viewer);
175:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
176:   PetscViewerASCIIPopTab(viewer);
177:   PetscViewerASCIIPrintf(viewer,"</Geometry>\n");
178:   PetscViewerASCIIPopTab(viewer);

180:   /* write geometry data */
181:   DMSwarmCreateGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);
182:   VecView(dvec,fviewer);
183:   DMSwarmDestroyGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);
184:   bytes[0] += sizeof(PetscReal) * ng * dim;

186:   PetscViewerDestroy(&fviewer);
187:   return 0;
188: }

190: PetscErrorCode private_VecView_Swarm_XDMF(Vec x,PetscViewer viewer)
191: {
192:   long int       *bytes = NULL;
193:   PetscContainer container = NULL;
194:   const char     *viewername;
195:   char           datafile[PETSC_MAX_PATH_LEN];
196:   char           *datafilename;
197:   PetscViewer    fviewer;
198:   PetscInt       N,bs;
199:   const char     *vecname;
200:   char           fieldname[PETSC_MAX_PATH_LEN];

202:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
204:   PetscContainerGetPointer(container,(void**)&bytes);
205:   PetscViewerFileGetName(viewer,&viewername);
206:   private_CreateDataFileNameXDMF(viewername,datafile);

208:   /* re-open a sub-viewer for all data fields */
209:   /* name is viewer.name + "_swarm_fields.pbin" */
210:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
211:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
212:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
213:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
214:   PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);
215:   PetscViewerFileSetName(fviewer,datafile);
216:   PetscStrrchr(datafile,'/',&datafilename);

218:   VecGetSize(x,&N);
219:   VecGetBlockSize(x,&bs);
220:   N = N/bs;
221:   PetscObjectGetName((PetscObject)x,&vecname);
222:   if (!vecname) {
223:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)x)->tag);
224:   } else {
225:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);
226:   }

228:   /* write data header */
229:   PetscViewerASCIIPushTab(viewer);
230:   PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);
231:   PetscViewerASCIIPushTab(viewer);
232:   if (bs == 1) {
233:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);
234:   } else {
235:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);
236:   }
237:   PetscViewerASCIIPushTab(viewer);
238:   PetscViewerASCIIPrintf(viewer,"%s\n",datafilename);
239:   PetscViewerASCIIPopTab(viewer);
240:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
241:   PetscViewerASCIIPopTab(viewer);
242:   PetscViewerASCIIPrintf(viewer,"</Attribute>\n");
243:   PetscViewerASCIIPopTab(viewer);

245:   /* write data */
246:   VecView(x,fviewer);
247:   bytes[0] += sizeof(PetscReal) * N * bs;

249:   PetscViewerDestroy(&fviewer);
250:   return 0;
251: }

253: PetscErrorCode private_ISView_Swarm_XDMF(IS is,PetscViewer viewer)
254: {
255:   long int       *bytes = NULL;
256:   PetscContainer container = NULL;
257:   const char     *viewername;
258:   char           datafile[PETSC_MAX_PATH_LEN];
259:   char           *datafilename;
260:   PetscViewer    fviewer;
261:   PetscInt       N,bs;
262:   const char     *vecname;
263:   char           fieldname[PETSC_MAX_PATH_LEN];

265:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
267:   PetscContainerGetPointer(container,(void**)&bytes);
268:   PetscViewerFileGetName(viewer,&viewername);
269:   private_CreateDataFileNameXDMF(viewername,datafile);

271:   /* re-open a sub-viewer for all data fields */
272:   /* name is viewer.name + "_swarm_fields.pbin" */
273:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
274:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
275:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
276:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
277:   PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);
278:   PetscViewerFileSetName(fviewer,datafile);
279:   PetscStrrchr(datafile,'/',&datafilename);

281:   ISGetSize(is,&N);
282:   ISGetBlockSize(is,&bs);
283:   N = N/bs;
284:   PetscObjectGetName((PetscObject)is,&vecname);
285:   if (!vecname) {
286:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)is)->tag);
287:   } else {
288:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);
289:   }

291:   /* write data header */
292:   PetscViewerASCIIPushTab(viewer);
293:   PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);
294:   PetscViewerASCIIPushTab(viewer);
295:   if (bs == 1) {
296:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);
297:   } else {
298:     PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);
299:   }
300:   PetscViewerASCIIPushTab(viewer);
301:   PetscViewerASCIIPrintf(viewer,"%s\n",datafilename);
302:   PetscViewerASCIIPopTab(viewer);
303:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
304:   PetscViewerASCIIPopTab(viewer);
305:   PetscViewerASCIIPrintf(viewer,"</Attribute>\n");
306:   PetscViewerASCIIPopTab(viewer);

308:   /* write data */
309:   ISView(is,fviewer);
310:   bytes[0] += sizeof(PetscInt) * N * bs;

312:   PetscViewerDestroy(&fviewer);
313:   return 0;
314: }

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

319:    Collective on dm

321:    Input parameters:
322: +  dm - the DMSwarm
323: .  filename - the file name of the XDMF file (must have the extension .xmf)
324: .  nfields - the number of fields to write into the XDMF file
325: -  field_name_list - array of length nfields containing the textual name of fields to write

327:    Level: beginner

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

332: .seealso: DMSwarmViewXDMF()
333: @*/
334: PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm,const char filename[],PetscInt nfields,const char *field_name_list[])
335: {
336:   Vec            dvec;
337:   PetscInt       f,N;
338:   PetscViewer    viewer;

340:   private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);
341:   private_DMSwarmView_XDMF(dm,viewer);
342:   DMSwarmGetLocalSize(dm,&N);
343:   for (f=0; f<nfields; f++) {
344:     void          *data;
345:     PetscDataType type;

347:     DMSwarmGetField(dm,field_name_list[f],NULL,&type,&data);
348:     DMSwarmRestoreField(dm,field_name_list[f],NULL,&type,&data);
349:     if (type == PETSC_DOUBLE) {
350:       DMSwarmCreateGlobalVectorFromField(dm,field_name_list[f],&dvec);
351:       PetscObjectSetName((PetscObject)dvec,field_name_list[f]);
352:       private_VecView_Swarm_XDMF(dvec,viewer);
353:       DMSwarmDestroyGlobalVectorFromField(dm,field_name_list[f],&dvec);
354:     } else if (type == PETSC_INT) {
355:       IS is;
356:       const PetscInt *idx;

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

361:       ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);
362:       PetscObjectSetName((PetscObject)is,field_name_list[f]);
363:       private_ISView_Swarm_XDMF(is,viewer);
364:       ISDestroy(&is);
365:       DMSwarmRestoreField(dm,field_name_list[f],NULL,&type,&data);
366:     } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Can only write PETSC_INT and PETSC_DOUBLE");

368:   }
369:   private_PetscViewerDestroy_XDMF(&viewer);
370:   return 0;
371: }

373: /*@C
374:    DMSwarmViewXDMF - Write DMSwarm fields to an XDMF3 file

376:    Collective on dm

378:    Input parameters:
379: +  dm - the DMSwarm
380: -  filename - the file name of the XDMF file (must have the extension .xmf)

382:    Level: beginner

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

387:    Developer Notes:
388:      This should be removed and replaced with the standard use of PetscViewer

390: .seealso: DMSwarmViewFieldsXDMF()
391: @*/
392: PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm,const char filename[])
393: {
394:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
395:   Vec            dvec;
396:   PetscInt       f;
397:   PetscViewer    viewer;

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

404:     /* query field type - accept all those of type PETSC_DOUBLE */
405:     field = swarm->db->field[f];
406:     if (field->petsc_type == PETSC_DOUBLE) {
407:       DMSwarmCreateGlobalVectorFromField(dm,field->name,&dvec);
408:       PetscObjectSetName((PetscObject)dvec,field->name);
409:       private_VecView_Swarm_XDMF(dvec,viewer);
410:       DMSwarmDestroyGlobalVectorFromField(dm,field->name,&dvec);
411:     } else if (field->petsc_type == PETSC_INT) {
412:       IS             is;
413:       PetscInt       N;
414:       const PetscInt *idx;
415:       void           *data;

417:       DMSwarmGetLocalSize(dm,&N);
418:       DMSwarmGetField(dm,field->name,NULL,NULL,&data);
419:       idx = (const PetscInt*)data;

421:       ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);
422:       PetscObjectSetName((PetscObject)is,field->name);
423:       private_ISView_Swarm_XDMF(is,viewer);
424:       ISDestroy(&is);
425:       DMSwarmRestoreField(dm,field->name,NULL,NULL,&data);
426:     }
427:   }
428:   private_PetscViewerDestroy_XDMF(&viewer);
429:   return 0;
430: }