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"
6 
private_PetscViewerCreate_XDMF(MPI_Comm comm,const char filename[],PetscViewer * v)7 PetscErrorCode private_PetscViewerCreate_XDMF(MPI_Comm comm,const char filename[],PetscViewer *v)
8 {
9   long int       *bytes;
10   PetscContainer container;
11   PetscViewer    viewer;
12   PetscErrorCode ierr;
13 
14   PetscFunctionBegin;
15   ierr = PetscViewerCreate(comm,&viewer);CHKERRQ(ierr);
16   ierr = PetscViewerSetType(viewer,PETSCVIEWERASCII);CHKERRQ(ierr);
17   ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
18   ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr);
19 
20   ierr = PetscMalloc1(1,&bytes);CHKERRQ(ierr);
21   bytes[0] = 0;
22   ierr = PetscContainerCreate(comm,&container);CHKERRQ(ierr);
23   ierr = PetscContainerSetPointer(container,(void*)bytes);CHKERRQ(ierr);
24   ierr = PetscObjectCompose((PetscObject)viewer,"XDMFViewerContext",(PetscObject)container);CHKERRQ(ierr);
25 
26   /* write xdmf header */
27   ierr = PetscViewerASCIIPrintf(viewer,"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n");CHKERRQ(ierr);
28   ierr = PetscViewerASCIIPrintf(viewer,"<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude/\" Version=\"2.99\">\n");CHKERRQ(ierr);
29   /* write xdmf domain */
30   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
31   ierr = PetscViewerASCIIPrintf(viewer,"<Domain>\n");CHKERRQ(ierr);
32   *v = viewer;
33   PetscFunctionReturn(0);
34 }
35 
private_PetscViewerDestroy_XDMF(PetscViewer * v)36 PetscErrorCode private_PetscViewerDestroy_XDMF(PetscViewer *v)
37 {
38   PetscViewer    viewer;
39   DM             dm = NULL;
40   long int       *bytes;
41   PetscContainer container = NULL;
42   PetscErrorCode ierr;
43 
44   PetscFunctionBegin;
45   if (!v) PetscFunctionReturn(0);
46   viewer = *v;
47 
48   ierr = PetscObjectQuery((PetscObject)viewer,"DMSwarm",(PetscObject*)&dm);CHKERRQ(ierr);
49   if (dm) {
50     ierr = PetscViewerASCIIPrintf(viewer,"</Grid>\n");CHKERRQ(ierr);
51     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
52   }
53 
54   /* close xdmf header */
55   ierr = PetscViewerASCIIPrintf(viewer,"</Domain>\n");CHKERRQ(ierr);
56   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
57   ierr = PetscViewerASCIIPrintf(viewer,"</Xdmf>\n");CHKERRQ(ierr);
58 
59   ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr);
60   if (container) {
61     ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr);
62     ierr = PetscFree(bytes);CHKERRQ(ierr);
63     ierr = PetscContainerDestroy(&container);CHKERRQ(ierr);
64   }
65 
66   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
67   *v = NULL;
68   PetscFunctionReturn(0);
69 }
70 
private_CreateDataFileNameXDMF(const char filename[],char dfilename[])71 PetscErrorCode private_CreateDataFileNameXDMF(const char filename[],char dfilename[])
72 {
73   char           *ext;
74   PetscBool      flg;
75   PetscErrorCode ierr;
76 
77   PetscFunctionBegin;
78   ierr = PetscStrrchr(filename,'.',&ext);CHKERRQ(ierr);
79   ierr = PetscStrcmp("xmf",ext,&flg);CHKERRQ(ierr);
80   if (flg) {
81     size_t len;
82     char    viewername_minus_ext[PETSC_MAX_PATH_LEN];
83 
84     ierr = PetscStrlen(filename,&len);CHKERRQ(ierr);
85     ierr = PetscStrncpy(viewername_minus_ext,filename,len-2);CHKERRQ(ierr);
86     ierr = PetscSNPrintf(dfilename,PETSC_MAX_PATH_LEN-1,"%s_swarm_fields.pbin",viewername_minus_ext);CHKERRQ(ierr);
87   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"File extension must by .xmf");
88 
89   PetscFunctionReturn(0);
90 }
91 
private_DMSwarmView_XDMF(DM dm,PetscViewer viewer)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;
103   PetscErrorCode ierr;
104 
105   PetscFunctionBegin;
106   ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr);
107   if (container) {
108     ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr);
109   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");
110 
111   ierr = PetscObjectTypeCompare((PetscObject)dm,DMSWARM,&isswarm);CHKERRQ(ierr);
112   if (!isswarm) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Only valid for DMSwarm");
113 
114   ierr = PetscObjectCompose((PetscObject)viewer,"DMSwarm",(PetscObject)dm);CHKERRQ(ierr);
115 
116   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
117   ierr = PetscObjectGetName((PetscObject)dm,&dmname);CHKERRQ(ierr);
118   if (!dmname) {
119     ierr = DMGetOptionsPrefix(dm,&dmname);CHKERRQ(ierr);
120   }
121   if (!dmname) {
122     ierr = PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n");CHKERRQ(ierr);
123   } else {
124     ierr = PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n",dmname);CHKERRQ(ierr);
125   }
126 
127   /* create a sub-viewer for topology, geometry and all data fields */
128   /* name is viewer.name + "_swarm_fields.pbin" */
129   ierr = PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);CHKERRQ(ierr);
130   ierr = PetscViewerSetType(fviewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
131   ierr = PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);CHKERRQ(ierr);
132   ierr = PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);CHKERRQ(ierr);
133   ierr = PetscViewerFileSetMode(fviewer,FILE_MODE_WRITE);CHKERRQ(ierr);
134 
135   ierr = PetscViewerFileGetName(viewer,&viewername);CHKERRQ(ierr);
136   ierr = private_CreateDataFileNameXDMF(viewername,datafile);CHKERRQ(ierr);
137   ierr = PetscViewerFileSetName(fviewer,datafile);CHKERRQ(ierr);
138 
139   ierr = DMSwarmGetSize(dm,&ng);CHKERRQ(ierr);
140 
141   /* write topology header */
142   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
143   ierr = PetscViewerASCIIPrintf(viewer,"<Topology Dimensions=\"%D\" TopologyType=\"Mixed\">\n",ng);CHKERRQ(ierr);
144   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
145   ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Dimensions=\"%D\" Seek=\"%D\">\n",ng*3,bytes[0]);CHKERRQ(ierr);
146   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
147   ierr = PetscViewerASCIIPrintf(viewer,"%s\n",datafile);CHKERRQ(ierr);
148   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
149   ierr = PetscViewerASCIIPrintf(viewer,"</DataItem>\n");CHKERRQ(ierr);
150   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
151   ierr = PetscViewerASCIIPrintf(viewer,"</Topology>\n");CHKERRQ(ierr);
152   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
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     ierr = PetscViewerBinaryWrite(fviewer,pvertex,3,PETSC_INT);CHKERRQ(ierr);
162   }
163   bytes[0] += sizeof(PetscInt) * ng * 3;
164 
165   /* write geometry header */
166   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
167   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
168   switch (dim) {
169     case 1:
170       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for 1D");
171       break;
172     case 2:
173       ierr = PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XY\">\n");CHKERRQ(ierr);
174       break;
175     case 3:
176       ierr = PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XYZ\">\n");CHKERRQ(ierr);
177       break;
178   }
179   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
180   ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",ng,dim,bytes[0]);CHKERRQ(ierr);
181   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
182   ierr = PetscViewerASCIIPrintf(viewer,"%s\n",datafile);CHKERRQ(ierr);
183   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
184   ierr = PetscViewerASCIIPrintf(viewer,"</DataItem>\n");CHKERRQ(ierr);
185   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
186   ierr = PetscViewerASCIIPrintf(viewer,"</Geometry>\n");CHKERRQ(ierr);
187   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
188 
189   /* write geometry data */
190   ierr = DMSwarmCreateGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);CHKERRQ(ierr);
191   ierr = VecView(dvec,fviewer);CHKERRQ(ierr);
192   ierr = DMSwarmDestroyGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);CHKERRQ(ierr);
193   bytes[0] += sizeof(PetscReal) * ng * dim;
194 
195   ierr = PetscViewerDestroy(&fviewer);CHKERRQ(ierr);
196 
197   PetscFunctionReturn(0);
198 }
199 
private_VecView_Swarm_XDMF(Vec x,PetscViewer viewer)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];
210   PetscErrorCode ierr;
211 
212   PetscFunctionBegin;
213   ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr);
214   if (container) {
215     ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr);
216   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");
217 
218   ierr = PetscViewerFileGetName(viewer,&viewername);CHKERRQ(ierr);
219   ierr = private_CreateDataFileNameXDMF(viewername,datafile);CHKERRQ(ierr);
220 
221   /* re-open a sub-viewer for all data fields */
222   /* name is viewer.name + "_swarm_fields.pbin" */
223   ierr = PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);CHKERRQ(ierr);
224   ierr = PetscViewerSetType(fviewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
225   ierr = PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);CHKERRQ(ierr);
226   ierr = PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);CHKERRQ(ierr);
227   ierr = PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);CHKERRQ(ierr);
228   ierr = PetscViewerFileSetName(fviewer,datafile);CHKERRQ(ierr);
229 
230   ierr = VecGetSize(x,&N);CHKERRQ(ierr);
231   ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr);
232   N = N/bs;
233   ierr = PetscObjectGetName((PetscObject)x,&vecname);CHKERRQ(ierr);
234   if (!vecname) {
235     ierr = PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)x)->tag);CHKERRQ(ierr);
236   } else {
237     ierr = PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);CHKERRQ(ierr);
238   }
239 
240   /* write data header */
241   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
242   ierr = PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);CHKERRQ(ierr);
243   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
244   if (bs == 1) {
245     ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);CHKERRQ(ierr);
246   } else {
247     ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);CHKERRQ(ierr);
248   }
249   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
250   ierr = PetscViewerASCIIPrintf(viewer,"%s\n",datafile);CHKERRQ(ierr);
251   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
252   ierr = PetscViewerASCIIPrintf(viewer,"</DataItem>\n");CHKERRQ(ierr);
253   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
254   ierr = PetscViewerASCIIPrintf(viewer,"</Attribute>\n");CHKERRQ(ierr);
255   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
256 
257   /* write data */
258   ierr = VecView(x,fviewer);CHKERRQ(ierr);
259   bytes[0] += sizeof(PetscReal) * N * bs;
260 
261   ierr = PetscViewerDestroy(&fviewer);CHKERRQ(ierr);
262 
263   PetscFunctionReturn(0);
264 }
265 
private_ISView_Swarm_XDMF(IS is,PetscViewer viewer)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];
276   PetscErrorCode ierr;
277 
278   PetscFunctionBegin;
279   ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr);
280   if (container) {
281     ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr);
282   } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");
283 
284   ierr = PetscViewerFileGetName(viewer,&viewername);CHKERRQ(ierr);
285   ierr = private_CreateDataFileNameXDMF(viewername,datafile);CHKERRQ(ierr);
286 
287   /* re-open a sub-viewer for all data fields */
288   /* name is viewer.name + "_swarm_fields.pbin" */
289   ierr = PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);CHKERRQ(ierr);
290   ierr = PetscViewerSetType(fviewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
291   ierr = PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);CHKERRQ(ierr);
292   ierr = PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);CHKERRQ(ierr);
293   ierr = PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);CHKERRQ(ierr);
294   ierr = PetscViewerFileSetName(fviewer,datafile);CHKERRQ(ierr);
295 
296   ierr = ISGetSize(is,&N);CHKERRQ(ierr);
297   ierr = ISGetBlockSize(is,&bs);CHKERRQ(ierr);
298   N = N/bs;
299   ierr = PetscObjectGetName((PetscObject)is,&vecname);CHKERRQ(ierr);
300   if (!vecname) {
301     ierr = PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)is)->tag);CHKERRQ(ierr);
302   } else {
303     ierr = PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);CHKERRQ(ierr);
304   }
305 
306   /* write data header */
307   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
308   ierr = PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);CHKERRQ(ierr);
309   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
310   if (bs == 1) {
311     ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);CHKERRQ(ierr);
312   } else {
313     ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);CHKERRQ(ierr);
314   }
315   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
316   ierr = PetscViewerASCIIPrintf(viewer,"%s\n",datafile);CHKERRQ(ierr);
317   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
318   ierr = PetscViewerASCIIPrintf(viewer,"</DataItem>\n");CHKERRQ(ierr);
319   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
320   ierr = PetscViewerASCIIPrintf(viewer,"</Attribute>\n");CHKERRQ(ierr);
321   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
322 
323   /* write data */
324   ierr = ISView(is,fviewer);CHKERRQ(ierr);
325   bytes[0] += sizeof(PetscInt) * N * bs;
326 
327   ierr = PetscViewerDestroy(&fviewer);CHKERRQ(ierr);
328 
329   PetscFunctionReturn(0);
330 }
331 
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
344 
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 @*/
DMSwarmViewFieldsXDMF(DM dm,const char filename[],PetscInt nfields,const char * field_name_list[])350 PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm,const char filename[],PetscInt nfields,const char *field_name_list[])
351 {
352   PetscErrorCode ierr;
353   Vec            dvec;
354   PetscInt       f,N;
355   PetscViewer    viewer;
356 
357   PetscFunctionBegin;
358   ierr = private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);CHKERRQ(ierr);
359   ierr = private_DMSwarmView_XDMF(dm,viewer);CHKERRQ(ierr);
360   ierr = DMSwarmGetLocalSize(dm,&N);CHKERRQ(ierr);
361   for (f=0; f<nfields; f++) {
362     void          *data;
363     PetscDataType type;
364 
365     ierr = DMSwarmGetField(dm,field_name_list[f],NULL,&type,&data);CHKERRQ(ierr);
366     ierr = DMSwarmRestoreField(dm,field_name_list[f],NULL,&type,&data);CHKERRQ(ierr);
367 
368     if (type == PETSC_DOUBLE) {
369       ierr = DMSwarmCreateGlobalVectorFromField(dm,field_name_list[f],&dvec);CHKERRQ(ierr);
370       ierr = PetscObjectSetName((PetscObject)dvec,field_name_list[f]);CHKERRQ(ierr);
371       ierr = private_VecView_Swarm_XDMF(dvec,viewer);CHKERRQ(ierr);
372       ierr = DMSwarmDestroyGlobalVectorFromField(dm,field_name_list[f],&dvec);CHKERRQ(ierr);
373     } else if (type == PETSC_INT) {
374       IS is;
375       const PetscInt *idx;
376 
377       ierr = DMSwarmGetField(dm,field_name_list[f],NULL,&type,&data);CHKERRQ(ierr);
378       idx = (const PetscInt*)data;
379 
380       ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);CHKERRQ(ierr);
381       ierr = PetscObjectSetName((PetscObject)is,field_name_list[f]);CHKERRQ(ierr);
382       ierr = private_ISView_Swarm_XDMF(is,viewer);CHKERRQ(ierr);
383       ierr = ISDestroy(&is);CHKERRQ(ierr);
384       ierr = DMSwarmRestoreField(dm,field_name_list[f],NULL,&type,&data);CHKERRQ(ierr);
385     } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Can only write PETSC_INT and PETSC_DOUBLE");
386 
387   }
388   ierr = private_PetscViewerDestroy_XDMF(&viewer);CHKERRQ(ierr);
389   PetscFunctionReturn(0);
390 }
391 
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
402 
403    Notes:
404    Only fields user registered with data type PETSC_DOUBLE or PETSC_INT will be written into the file
405 
406 .seealso: DMSwarmViewFieldsXDMF()
407 @*/
DMSwarmViewXDMF(DM dm,const char filename[])408 PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm,const char filename[])
409 {
410   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
411   PetscErrorCode ierr;
412   Vec            dvec;
413   PetscInt       f;
414   PetscViewer    viewer;
415 
416   PetscFunctionBegin;
417   ierr = private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm),filename,&viewer);CHKERRQ(ierr);
418   ierr = private_DMSwarmView_XDMF(dm,viewer);CHKERRQ(ierr);
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       ierr = DMSwarmCreateGlobalVectorFromField(dm,field->name,&dvec);CHKERRQ(ierr);
426       ierr = PetscObjectSetName((PetscObject)dvec,field->name);CHKERRQ(ierr);
427       ierr = private_VecView_Swarm_XDMF(dvec,viewer);CHKERRQ(ierr);
428       ierr = DMSwarmDestroyGlobalVectorFromField(dm,field->name,&dvec);CHKERRQ(ierr);
429     } else if (field->petsc_type == PETSC_INT) {
430       IS             is;
431       PetscInt       N;
432       const PetscInt *idx;
433       void           *data;
434 
435       ierr = DMSwarmGetLocalSize(dm,&N);CHKERRQ(ierr);
436       ierr = DMSwarmGetField(dm,field->name,NULL,NULL,&data);CHKERRQ(ierr);
437       idx = (const PetscInt*)data;
438 
439       ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),N,idx,PETSC_USE_POINTER,&is);CHKERRQ(ierr);
440       ierr = PetscObjectSetName((PetscObject)is,field->name);CHKERRQ(ierr);
441       ierr = private_ISView_Swarm_XDMF(is,viewer);CHKERRQ(ierr);
442       ierr = ISDestroy(&is);CHKERRQ(ierr);
443       ierr = DMSwarmRestoreField(dm,field->name,NULL,NULL,&data);CHKERRQ(ierr);
444     }
445   }
446   ierr = private_PetscViewerDestroy_XDMF(&viewer);CHKERRQ(ierr);
447   PetscFunctionReturn(0);
448 }
449