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) DMGetOptionsPrefix(dm, &dmname);
110:   if (!dmname) {
111:     PetscViewerASCIIPrintf(viewer, "<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n");
112:   } else {
113:     PetscViewerASCIIPrintf(viewer, "<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n", dmname);
114:   }

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

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

129:   DMSwarmGetSize(dm, &ng);

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

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

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

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

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

184:   PetscViewerDestroy(&fviewer);
185:   return 0;
186: }

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

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

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

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

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

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

247:   PetscViewerDestroy(&fviewer);
248:   return 0;
249: }

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

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

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

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

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

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

310:   PetscViewerDestroy(&fviewer);
311:   return 0;
312: }

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

317:    Collective on dm

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

325:    Level: beginner

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

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

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

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

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

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

370: /*@C
371:    DMSwarmViewXDMF - Write DMSwarm fields to an XDMF3 file

373:    Collective on dm

375:    Input parameters:
376: +  dm - the DMSwarm
377: -  filename - the file name of the XDMF file (must have the extension .xmf)

379:    Level: beginner

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

384:    Developer Notes:
385:      This should be removed and replaced with the standard use of PetscViewer

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

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

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

414:       DMSwarmGetLocalSize(dm, &N);
415:       DMSwarmGetField(dm, field->name, NULL, NULL, &data);
416:       idx = (const PetscInt *)data;

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