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: static PetscErrorCode private_PetscViewerCreate_XDMF(MPI_Comm comm, const char filename[], PetscViewer *v)
8: {
9: long int *bytes;
10: PetscContainer container;
11: PetscViewer viewer;
13: PetscFunctionBegin;
14: PetscCall(PetscViewerCreate(comm, &viewer));
15: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII));
16: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
17: PetscCall(PetscViewerFileSetName(viewer, filename));
19: PetscCall(PetscMalloc1(1, &bytes));
20: bytes[0] = 0;
21: PetscCall(PetscContainerCreate(comm, &container));
22: PetscCall(PetscContainerSetPointer(container, (void *)bytes));
23: PetscCall(PetscObjectCompose((PetscObject)viewer, "XDMFViewerContext", (PetscObject)container));
25: /* write xdmf header */
26: PetscCall(PetscViewerASCIIPrintf(viewer, "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n"));
27: PetscCall(PetscViewerASCIIPrintf(viewer, "<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude/\" Version=\"2.99\">\n"));
28: /* write xdmf domain */
29: PetscCall(PetscViewerASCIIPushTab(viewer));
30: PetscCall(PetscViewerASCIIPrintf(viewer, "<Domain>\n"));
31: *v = viewer;
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: static PetscErrorCode private_PetscViewerDestroy_XDMF(PetscViewer *v)
36: {
37: PetscViewer viewer;
38: DM dm = NULL;
39: long int *bytes;
40: PetscContainer container = NULL;
42: PetscFunctionBegin;
43: if (!v) PetscFunctionReturn(PETSC_SUCCESS);
44: viewer = *v;
46: PetscCall(PetscObjectQuery((PetscObject)viewer, "DMSwarm", (PetscObject *)&dm));
47: if (dm) {
48: PetscCall(PetscViewerASCIIPrintf(viewer, "</Grid>\n"));
49: PetscCall(PetscViewerASCIIPopTab(viewer));
50: }
52: /* close xdmf header */
53: PetscCall(PetscViewerASCIIPrintf(viewer, "</Domain>\n"));
54: PetscCall(PetscViewerASCIIPopTab(viewer));
55: PetscCall(PetscViewerASCIIPrintf(viewer, "</Xdmf>\n"));
57: PetscCall(PetscObjectQuery((PetscObject)viewer, "XDMFViewerContext", (PetscObject *)&container));
58: if (container) {
59: PetscCall(PetscContainerGetPointer(container, (void **)&bytes));
60: PetscCall(PetscFree(bytes));
61: PetscCall(PetscContainerDestroy(&container));
62: }
63: PetscCall(PetscViewerDestroy(&viewer));
64: *v = NULL;
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: static PetscErrorCode private_CreateDataFileNameXDMF(const char filename[], char dfilename[])
69: {
70: const char dot_xmf[] = ".xmf";
71: size_t len;
72: char viewername_minus_ext[PETSC_MAX_PATH_LEN];
73: PetscBool flg;
75: PetscFunctionBegin;
76: PetscCall(PetscStrendswith(filename, dot_xmf, &flg));
77: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "File extension must be %s", dot_xmf);
78: PetscCall(PetscStrncpy(viewername_minus_ext, filename, sizeof(viewername_minus_ext)));
79: PetscCall(PetscStrlen(filename, &len));
80: len -= sizeof(dot_xmf) - 1;
81: if (sizeof(viewername_minus_ext) > len) viewername_minus_ext[len] = '\0';
82: PetscCall(PetscSNPrintf(dfilename, PETSC_MAX_PATH_LEN - 1, "%s_swarm_fields.pbin", viewername_minus_ext));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode private_DMSwarmView_XDMF(DM dm, PetscViewer viewer)
87: {
88: PetscBool isswarm = PETSC_FALSE;
89: const char *viewername;
90: char datafile[PETSC_MAX_PATH_LEN];
91: char *datafilename;
92: PetscViewer fviewer;
93: PetscInt k, ng, dim;
94: Vec dvec;
95: long int *bytes = NULL;
96: PetscContainer container = NULL;
97: const char *dmname;
99: PetscFunctionBegin;
100: PetscCall(PetscObjectQuery((PetscObject)viewer, "XDMFViewerContext", (PetscObject *)&container));
101: if (container) {
102: PetscCall(PetscContainerGetPointer(container, (void **)&bytes));
103: } else SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Valid to find attached data XDMFViewerContext");
105: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMSWARM, &isswarm));
106: PetscCheck(isswarm, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Only valid for DMSwarm");
108: PetscCall(PetscObjectCompose((PetscObject)viewer, "DMSwarm", (PetscObject)dm));
110: PetscCall(PetscViewerASCIIPushTab(viewer));
111: PetscCall(PetscObjectGetName((PetscObject)dm, &dmname));
112: if (!dmname) PetscCall(DMGetOptionsPrefix(dm, &dmname));
113: if (!dmname) {
114: PetscCall(PetscViewerASCIIPrintf(viewer, "<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n"));
115: } else {
116: PetscCall(PetscViewerASCIIPrintf(viewer, "<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n", dmname));
117: }
119: /* create a sub-viewer for topology, geometry and all data fields */
120: /* name is viewer.name + "_swarm_fields.pbin" */
121: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)viewer), &fviewer));
122: PetscCall(PetscViewerSetType(fviewer, PETSCVIEWERBINARY));
123: PetscCall(PetscViewerBinarySetSkipHeader(fviewer, PETSC_TRUE));
124: PetscCall(PetscViewerBinarySetSkipInfo(fviewer, PETSC_TRUE));
125: PetscCall(PetscViewerFileSetMode(fviewer, FILE_MODE_WRITE));
127: PetscCall(PetscViewerFileGetName(viewer, &viewername));
128: PetscCall(private_CreateDataFileNameXDMF(viewername, datafile));
129: PetscCall(PetscViewerFileSetName(fviewer, datafile));
130: PetscCall(PetscStrrchr(datafile, '/', &datafilename));
132: PetscCall(DMSwarmGetSize(dm, &ng));
134: /* write topology header */
135: PetscCall(PetscViewerASCIIPushTab(viewer));
136: PetscCall(PetscViewerASCIIPrintf(viewer, "<Topology Dimensions=\"%" PetscInt_FMT "\" TopologyType=\"Mixed\">\n", ng));
137: PetscCall(PetscViewerASCIIPushTab(viewer));
138: PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Dimensions=\"%" PetscInt_FMT "\" Seek=\"%ld\">\n", ng * 3, bytes[0]));
139: PetscCall(PetscViewerASCIIPushTab(viewer));
140: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", datafilename));
141: PetscCall(PetscViewerASCIIPopTab(viewer));
142: PetscCall(PetscViewerASCIIPrintf(viewer, "</DataItem>\n"));
143: PetscCall(PetscViewerASCIIPopTab(viewer));
144: PetscCall(PetscViewerASCIIPrintf(viewer, "</Topology>\n"));
145: PetscCall(PetscViewerASCIIPopTab(viewer));
147: /* write topology data */
148: for (k = 0; k < ng; k++) {
149: PetscInt pvertex[3];
151: pvertex[0] = 1;
152: pvertex[1] = 1;
153: pvertex[2] = k;
154: PetscCall(PetscViewerBinaryWrite(fviewer, pvertex, 3, PETSC_INT));
155: }
156: bytes[0] += sizeof(PetscInt) * ng * 3;
158: /* write geometry header */
159: PetscCall(PetscViewerASCIIPushTab(viewer));
160: PetscCall(DMGetDimension(dm, &dim));
161: switch (dim) {
162: case 1:
163: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for 1D");
164: case 2:
165: PetscCall(PetscViewerASCIIPrintf(viewer, "<Geometry Type=\"XY\">\n"));
166: break;
167: case 3:
168: PetscCall(PetscViewerASCIIPrintf(viewer, "<Geometry Type=\"XYZ\">\n"));
169: break;
170: }
171: PetscCall(PetscViewerASCIIPushTab(viewer));
172: PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%" PetscInt_FMT " %" PetscInt_FMT "\" Seek=\"%ld\">\n", ng, dim, bytes[0]));
173: PetscCall(PetscViewerASCIIPushTab(viewer));
174: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", datafilename));
175: PetscCall(PetscViewerASCIIPopTab(viewer));
176: PetscCall(PetscViewerASCIIPrintf(viewer, "</DataItem>\n"));
177: PetscCall(PetscViewerASCIIPopTab(viewer));
178: PetscCall(PetscViewerASCIIPrintf(viewer, "</Geometry>\n"));
179: PetscCall(PetscViewerASCIIPopTab(viewer));
181: /* write geometry data */
182: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &dvec));
183: PetscCall(VecView(dvec, fviewer));
184: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &dvec));
185: bytes[0] += sizeof(PetscReal) * ng * dim;
187: PetscCall(PetscViewerDestroy(&fviewer));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: static PetscErrorCode private_VecView_Swarm_XDMF(Vec x, PetscViewer viewer)
192: {
193: long int *bytes = NULL;
194: PetscContainer container = NULL;
195: const char *viewername;
196: char datafile[PETSC_MAX_PATH_LEN];
197: char *datafilename;
198: PetscViewer fviewer;
199: PetscInt N, bs;
200: const char *vecname;
201: char fieldname[PETSC_MAX_PATH_LEN];
203: PetscFunctionBegin;
204: PetscCall(PetscObjectQuery((PetscObject)viewer, "XDMFViewerContext", (PetscObject *)&container));
205: PetscCheck(container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Unable to find attached data XDMFViewerContext");
206: PetscCall(PetscContainerGetPointer(container, (void **)&bytes));
207: PetscCall(PetscViewerFileGetName(viewer, &viewername));
208: PetscCall(private_CreateDataFileNameXDMF(viewername, datafile));
210: /* re-open a sub-viewer for all data fields */
211: /* name is viewer.name + "_swarm_fields.pbin" */
212: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)viewer), &fviewer));
213: PetscCall(PetscViewerSetType(fviewer, PETSCVIEWERBINARY));
214: PetscCall(PetscViewerBinarySetSkipHeader(fviewer, PETSC_TRUE));
215: PetscCall(PetscViewerBinarySetSkipInfo(fviewer, PETSC_TRUE));
216: PetscCall(PetscViewerFileSetMode(fviewer, FILE_MODE_APPEND));
217: PetscCall(PetscViewerFileSetName(fviewer, datafile));
218: PetscCall(PetscStrrchr(datafile, '/', &datafilename));
220: PetscCall(VecGetSize(x, &N));
221: PetscCall(VecGetBlockSize(x, &bs));
222: N = N / bs;
223: PetscCall(PetscObjectGetName((PetscObject)x, &vecname));
224: if (!vecname) {
225: PetscCall(PetscSNPrintf(fieldname, PETSC_MAX_PATH_LEN - 1, "swarmfield_%d", ((PetscObject)x)->tag));
226: } else {
227: PetscCall(PetscSNPrintf(fieldname, PETSC_MAX_PATH_LEN - 1, "%s", vecname));
228: }
230: /* write data header */
231: PetscCall(PetscViewerASCIIPushTab(viewer));
232: PetscCall(PetscViewerASCIIPrintf(viewer, "<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n", fieldname));
233: PetscCall(PetscViewerASCIIPushTab(viewer));
234: if (bs == 1) {
235: PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%" PetscInt_FMT "\" Seek=\"%ld\">\n", N, bytes[0]));
236: } else {
237: PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%" PetscInt_FMT " %" PetscInt_FMT "\" Seek=\"%ld\">\n", N, bs, bytes[0]));
238: }
239: PetscCall(PetscViewerASCIIPushTab(viewer));
240: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", datafilename));
241: PetscCall(PetscViewerASCIIPopTab(viewer));
242: PetscCall(PetscViewerASCIIPrintf(viewer, "</DataItem>\n"));
243: PetscCall(PetscViewerASCIIPopTab(viewer));
244: PetscCall(PetscViewerASCIIPrintf(viewer, "</Attribute>\n"));
245: PetscCall(PetscViewerASCIIPopTab(viewer));
247: /* write data */
248: PetscCall(VecView(x, fviewer));
249: bytes[0] += sizeof(PetscReal) * N * bs;
251: PetscCall(PetscViewerDestroy(&fviewer));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode private_ISView_Swarm_XDMF(IS is, PetscViewer viewer)
256: {
257: long int *bytes = NULL;
258: PetscContainer container = NULL;
259: const char *viewername;
260: char datafile[PETSC_MAX_PATH_LEN];
261: char *datafilename;
262: PetscViewer fviewer;
263: PetscInt N, bs;
264: const char *vecname;
265: char fieldname[PETSC_MAX_PATH_LEN];
267: PetscFunctionBegin;
268: PetscCall(PetscObjectQuery((PetscObject)viewer, "XDMFViewerContext", (PetscObject *)&container));
269: PetscCheck(container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Unable to find attached data XDMFViewerContext");
270: PetscCall(PetscContainerGetPointer(container, (void **)&bytes));
271: PetscCall(PetscViewerFileGetName(viewer, &viewername));
272: PetscCall(private_CreateDataFileNameXDMF(viewername, datafile));
274: /* re-open a sub-viewer for all data fields */
275: /* name is viewer.name + "_swarm_fields.pbin" */
276: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)viewer), &fviewer));
277: PetscCall(PetscViewerSetType(fviewer, PETSCVIEWERBINARY));
278: PetscCall(PetscViewerBinarySetSkipHeader(fviewer, PETSC_TRUE));
279: PetscCall(PetscViewerBinarySetSkipInfo(fviewer, PETSC_TRUE));
280: PetscCall(PetscViewerFileSetMode(fviewer, FILE_MODE_APPEND));
281: PetscCall(PetscViewerFileSetName(fviewer, datafile));
282: PetscCall(PetscStrrchr(datafile, '/', &datafilename));
284: PetscCall(ISGetSize(is, &N));
285: PetscCall(ISGetBlockSize(is, &bs));
286: N = N / bs;
287: PetscCall(PetscObjectGetName((PetscObject)is, &vecname));
288: if (!vecname) {
289: PetscCall(PetscSNPrintf(fieldname, PETSC_MAX_PATH_LEN - 1, "swarmfield_%d", ((PetscObject)is)->tag));
290: } else {
291: PetscCall(PetscSNPrintf(fieldname, PETSC_MAX_PATH_LEN - 1, "%s", vecname));
292: }
294: /* write data header */
295: PetscCall(PetscViewerASCIIPushTab(viewer));
296: PetscCall(PetscViewerASCIIPrintf(viewer, "<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n", fieldname));
297: PetscCall(PetscViewerASCIIPushTab(viewer));
298: if (bs == 1) {
299: PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%" PetscInt_FMT "\" Seek=\"%ld\">\n", N, bytes[0]));
300: } else {
301: PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%" PetscInt_FMT " %" PetscInt_FMT "\" Seek=\"%ld\">\n", N, bs, bytes[0]));
302: }
303: PetscCall(PetscViewerASCIIPushTab(viewer));
304: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", datafilename));
305: PetscCall(PetscViewerASCIIPopTab(viewer));
306: PetscCall(PetscViewerASCIIPrintf(viewer, "</DataItem>\n"));
307: PetscCall(PetscViewerASCIIPopTab(viewer));
308: PetscCall(PetscViewerASCIIPrintf(viewer, "</Attribute>\n"));
309: PetscCall(PetscViewerASCIIPopTab(viewer));
311: /* write data */
312: PetscCall(ISView(is, fviewer));
313: bytes[0] += sizeof(PetscInt) * N * bs;
315: PetscCall(PetscViewerDestroy(&fviewer));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@C
320: DMSwarmViewFieldsXDMF - Write a selection of DMSwarm fields to an XDMF3 file
322: Collective
324: Input Parameters:
325: + dm - the `DMSWARM`
326: . filename - the file name of the XDMF file (must have the extension .xmf)
327: . nfields - the number of fields to write into the XDMF file
328: - field_name_list - array of length nfields containing the textual name of fields to write
330: Level: beginner
332: Note:
333: Only fields registered with data type `PETSC_DOUBLE` or `PETSC_INT` can be written into the file
335: .seealso: `DM`, `DMSWARM`, `DMSwarmViewXDMF()`
336: @*/
337: PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm, const char filename[], PetscInt nfields, const char *field_name_list[])
338: {
339: Vec dvec;
340: PetscInt f, N;
341: PetscViewer viewer;
343: PetscFunctionBegin;
344: PetscCall(private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm), filename, &viewer));
345: PetscCall(private_DMSwarmView_XDMF(dm, viewer));
346: PetscCall(DMSwarmGetLocalSize(dm, &N));
347: for (f = 0; f < nfields; f++) {
348: void *data;
349: PetscDataType type;
351: PetscCall(DMSwarmGetField(dm, field_name_list[f], NULL, &type, &data));
352: PetscCall(DMSwarmRestoreField(dm, field_name_list[f], NULL, &type, &data));
353: if (type == PETSC_DOUBLE) {
354: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, field_name_list[f], &dvec));
355: PetscCall(PetscObjectSetName((PetscObject)dvec, field_name_list[f]));
356: PetscCall(private_VecView_Swarm_XDMF(dvec, viewer));
357: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, field_name_list[f], &dvec));
358: } else if (type == PETSC_INT) {
359: IS is;
360: const PetscInt *idx;
362: PetscCall(DMSwarmGetField(dm, field_name_list[f], NULL, &type, &data));
363: idx = (const PetscInt *)data;
365: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), N, idx, PETSC_USE_POINTER, &is));
366: PetscCall(PetscObjectSetName((PetscObject)is, field_name_list[f]));
367: PetscCall(private_ISView_Swarm_XDMF(is, viewer));
368: PetscCall(ISDestroy(&is));
369: PetscCall(DMSwarmRestoreField(dm, field_name_list[f], NULL, &type, &data));
370: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Can only write PETSC_INT and PETSC_DOUBLE");
371: }
372: PetscCall(private_PetscViewerDestroy_XDMF(&viewer));
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: /*@C
377: DMSwarmViewXDMF - Write `DMSWARM` fields to an XDMF3 file
379: Collective
381: Input Parameters:
382: + dm - the `DMSWARM`
383: - filename - the file name of the XDMF file (must have the extension .xmf)
385: Level: beginner
387: Note:
388: Only fields user registered with data type `PETSC_DOUBLE` or `PETSC_INT` will be written into the file
390: Developer Notes:
391: This should be removed and replaced with the standard use of `PetscViewer`
393: .seealso: `DM`, `DMSWARM`, `DMSwarmViewFieldsXDMF()`
394: @*/
395: PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm, const char filename[])
396: {
397: DM_Swarm *swarm = (DM_Swarm *)dm->data;
398: Vec dvec;
399: PetscInt f;
400: PetscViewer viewer;
402: PetscFunctionBegin;
403: PetscCall(private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm), filename, &viewer));
404: PetscCall(private_DMSwarmView_XDMF(dm, viewer));
405: for (f = 4; f < swarm->db->nfields; f++) { /* only examine user defined fields - the first 4 are internally created by DMSwarmPIC */
406: DMSwarmDataField field;
408: /* query field type - accept all those of type PETSC_DOUBLE */
409: field = swarm->db->field[f];
410: if (field->petsc_type == PETSC_DOUBLE) {
411: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, field->name, &dvec));
412: PetscCall(PetscObjectSetName((PetscObject)dvec, field->name));
413: PetscCall(private_VecView_Swarm_XDMF(dvec, viewer));
414: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, field->name, &dvec));
415: } else if (field->petsc_type == PETSC_INT) {
416: IS is;
417: PetscInt N;
418: const PetscInt *idx;
419: void *data;
421: PetscCall(DMSwarmGetLocalSize(dm, &N));
422: PetscCall(DMSwarmGetField(dm, field->name, NULL, NULL, &data));
423: idx = (const PetscInt *)data;
425: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), N, idx, PETSC_USE_POINTER, &is));
426: PetscCall(PetscObjectSetName((PetscObject)is, field->name));
427: PetscCall(private_ISView_Swarm_XDMF(is, viewer));
428: PetscCall(ISDestroy(&is));
429: PetscCall(DMSwarmRestoreField(dm, field->name, NULL, NULL, &data));
430: }
431: }
432: PetscCall(private_PetscViewerDestroy_XDMF(&viewer));
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }