Actual source code: swarmpic_view.c

petsc-3.14.6 2021-03-30
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;

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

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");

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

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

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:   }

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

135:   PetscViewerFileGetName(viewer,&viewername);
136:   private_CreateDataFileNameXDMF(viewername,datafile);
137:   PetscViewerFileSetName(fviewer,datafile);

139:   DMSwarmGetSize(dm,&ng);

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

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

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

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:     case 2:
172:       PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XY\">\n");
173:       break;
174:     case 3:
175:       PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XYZ\">\n");
176:       break;
177:   }
178:   PetscViewerASCIIPushTab(viewer);
179:   PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",ng,dim,bytes[0]);
180:   PetscViewerASCIIPushTab(viewer);
181:   PetscViewerASCIIPrintf(viewer,"%s\n",datafile);
182:   PetscViewerASCIIPopTab(viewer);
183:   PetscViewerASCIIPrintf(viewer,"</DataItem>\n");
184:   PetscViewerASCIIPopTab(viewer);
185:   PetscViewerASCIIPrintf(viewer,"</Geometry>\n");
186:   PetscViewerASCIIPopTab(viewer);

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

194:   PetscViewerDestroy(&fviewer);

196:   return(0);
197: }

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

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

217:   PetscViewerFileGetName(viewer,&viewername);
218:   private_CreateDataFileNameXDMF(viewername,datafile);

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

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

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

256:   /* write data */
257:   VecView(x,fviewer);
258:   bytes[0] += sizeof(PetscReal) * N * bs;

260:   PetscViewerDestroy(&fviewer);

262:   return(0);
263: }

265: PetscErrorCode private_ISView_Swarm_XDMF(IS is,PetscViewer viewer)
266: {
267:   long int       *bytes = NULL;
268:   PetscContainer container = NULL;
269:   const char     *viewername;
270:   char           datafile[PETSC_MAX_PATH_LEN];
271:   PetscViewer    fviewer;
272:   PetscInt       N,bs;
273:   const char     *vecname;
274:   char           fieldname[PETSC_MAX_PATH_LEN];

278:   PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);
279:   if (container) {
280:     PetscContainerGetPointer(container,(void**)&bytes);
281:   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");

283:   PetscViewerFileGetName(viewer,&viewername);
284:   private_CreateDataFileNameXDMF(viewername,datafile);

286:   /* re-open a sub-viewer for all data fields */
287:   /* name is viewer.name + "_swarm_fields.pbin" */
288:   PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);
289:   PetscViewerSetType(fviewer,PETSCVIEWERBINARY);
290:   PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);
291:   PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);
292:   PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);
293:   PetscViewerFileSetName(fviewer,datafile);

295:   ISGetSize(is,&N);
296:   ISGetBlockSize(is,&bs);
297:   N = N/bs;
298:   PetscObjectGetName((PetscObject)is,&vecname);
299:   if (!vecname) {
300:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)is)->tag);
301:   } else {
302:     PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);
303:   }

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

322:   /* write data */
323:   ISView(is,fviewer);
324:   bytes[0] += sizeof(PetscInt) * N * bs;

326:   PetscViewerDestroy(&fviewer);

328:   return(0);
329: }

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

334:    Collective on dm

336:    Input parameters:
337: +  dm - the DMSwarm
338: .  filename - the file name of the XDMF file (must have the extension .xmf)
339: .  nfields - the number of fields to write into the XDMF file
340: -  field_name_list - array of length nfields containing the textual name of fields to write

342:    Level: beginner

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

347: .seealso: DMSwarmViewXDMF()
348: @*/
349: PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm,const char filename[],PetscInt nfields,const char *field_name_list[])
350: {
352:   Vec            dvec;
353:   PetscInt       f,N;
354:   PetscViewer    viewer;

357:   private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);
358:   private_DMSwarmView_XDMF(dm,viewer);
359:   DMSwarmGetLocalSize(dm,&N);
360:   for (f=0; f<nfields; f++) {
361:     void          *data;
362:     PetscDataType type;

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

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

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

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

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

391: /*@C
392:    DMSwarmViewXDMF - Write DMSwarm fields to an XDMF3 file

394:    Collective on dm

396:    Input parameters:
397: +  dm - the DMSwarm
398: -  filename - the file name of the XDMF file (must have the extension .xmf)

400:    Level: beginner

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

405: .seealso: DMSwarmViewFieldsXDMF()
406: @*/
407: PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm,const char filename[])
408: {
409:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
411:   Vec            dvec;
412:   PetscInt       f;
413:   PetscViewer    viewer;

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

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

434:       DMSwarmGetLocalSize(dm,&N);
435:       DMSwarmGetField(dm,field->name,NULL,NULL,&data);
436:       idx = (const PetscInt*)data;

438:       ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);
439:       PetscObjectSetName((PetscObject)is,field->name);
440:       private_ISView_Swarm_XDMF(is,viewer);
441:       ISDestroy(&is);
442:       DMSwarmRestoreField(dm,field->name,NULL,NULL,&data);
443:     }
444:   }
445:   private_PetscViewerDestroy_XDMF(&viewer);
446:   return(0);
447: }