Actual source code: swarmpic_view.c
petsc-3.14.6 2021-03-30
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: }