Actual source code: ex21.c

  1: static char help[] = "Tests save/load of plex/section/vec in HDF5.\n\n";

  3: #include <petscdmshell.h>
  4: #include <petscdmplex.h>
  5: #include <petscsection.h>
  6: #include <petscsf.h>
  7: #include <petsclayouthdf5.h>

  9: /* A six-element mesh

 11: =====================
 12:  Save on 2 processes
 13: =====================

 15: exampleDMPlex: Local numbering:

 17:              7---17--8---18--9--19--(12)(24)(13)
 18:              |       |       |       |       |
 19:   rank 0:   20   0  21   1  22   2  (25) (3)(26)
 20:              |       |       |       |       |
 21:              4---14--5---15--6--16--(10)(23)(11)

 23:                            (13)(25)--8--17---9--18--10--19--11
 24:                              |       |       |       |       |
 25:   rank 1:                  (26) (3) 20   0   21  1  22   2  23
 26:                              |       |       |       |       |
 27:                            (12)(24)--4--14---5--15---6--16---7

 29: exampleDMPlex: globalPointNumbering:

 31:              9--23--10--24--11--25--16--32--17--33--18--34--19
 32:              |       |       |       |       |       |       |
 33:             26   0  27   1  28   2  35   3  36   4  37   5  38
 34:              |       |       |       |       |       |       |
 35:              6--20---7--21---8--22--12--29--13--30--14--31--15

 37: exampleSectionDM:
 38:   - includesConstraints = TRUE for local section (default)
 39:   - includesConstraints = FALSE for global section (default)

 41: exampleSectionDM: Dofs (Field 0):

 43:              0---0---0---0---0---0---2---0---0---0---0---0---0
 44:              |       |       |       |       |       |       |
 45:              0   0   0   0   0   0   0   2   0   0   0   0   0
 46:              |       |       |       |       |       |       |
 47:              0---0---0---0---0---0---0---0---0---0---0---0---0

 49: exampleSectionDM: Dofs (Field 1):      constrained
 50:                                       /
 51:              0---0---0---0---0---0---1---0---0---0---0---0---0
 52:              |       |       |       |       |       |       |
 53:              0   0   0   0   0   0   2   0   0   1   0   0   0
 54:              |       |       |       |       |       |       |
 55:              0---0---0---0---0---0---0---0---0---0---0---0---0

 57: exampleSectionDM: Offsets (total) in global section:

 59:              0---0---0---0---0---0---3---5---5---5---5---5---5
 60:              |       |       |       |       |       |       |
 61:              0   0   0   0   0   0   5   0   7   2   7   3   7
 62:              |       |       |       |       |       |       |
 63:              0---0---0---0---0---0---3---5---3---5---3---5---3

 65: exampleVec: Values (Field 0):          (1.3, 1.4)
 66:                                       /
 67:              +-------+-------+-------*-------+-------+-------+
 68:              |       |       |       |       |       |       |
 69:              |       |       |       |   * (1.0, 1.1)|       |
 70:              |       |       |       |       |       |       |
 71:              +-------+-------+-------+-------+-------+-------+

 73: exampleVec: Values (Field 1):          (1.5,) constrained
 74:                                       /
 75:              +-------+-------+-------*-------+-------+-------+
 76:              |       |       |       |       |       |       |
 77:              |       |    (1.6, 1.7) *       |   * (1.2,)    |
 78:              |       |       |       |       |       |       |
 79:              +-------+-------+-------+-------+-------+-------+

 81: exampleVec: as global vector

 83:   rank 0: []
 84:   rakn 1: [1.0, 1.1, 1.2, 1.3, 1.4, 1.6, 1.7]

 86: =====================
 87:  Load on 3 Processes
 88: =====================

 90: exampleDMPlex: Loaded/Distributed:

 92:              5--13---6--14--(8)(18)(10)
 93:              |       |       |       |
 94:   rank 0:   15   0   16  1  (19)(2)(20)
 95:              |       |       |       |
 96:              3--11---4--12--(7)(17)-(9)

 98:                     (9)(21)--5--15---7--18-(12)(24)(13)
 99:                      |       |       |       |       |
100:   rank 1:          (22) (2) 16   0  19   1 (25) (3)(26)
101:                      |       |       |       |       |
102:                     (8)(20)--4--14---6--17-(10)(23)(11)

104:                                +-> (10)(19)--6--13---7--14---8
105:                        permute |     |       |       |       |
106:   rank 2:                      +-> (20) (2) 15   0  16   1  17
107:                                      |       |       |       |
108:                                     (9)(18)--3--11---4--12---5

110: exampleSectionDM:
111:   - includesConstraints = TRUE for local section (default)
112:   - includesConstraints = FALSE for global section (default)

114: exampleVec: as local vector:

116:   rank 0: [1.3, 1.4, 1.5, 1.6, 1.7]
117:   rank 1: [1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7]
118:   rank 2: [1.2, 1.0, 1.1, 1.6, 1.7, 1.3, 1.4, 1.5]

120: exampleVec: as global vector:

122:   rank 0: []
123:   rank 1: [1.0, 1.1, 1.3, 1.4, 1.6, 1.7]
124:   rank 2: [1.2]

126: */

128: typedef struct {
129:   char      fname[PETSC_MAX_PATH_LEN]; /* Output mesh filename */
130:   PetscBool shell;                     /* Use DMShell to wrap sections */
131: } AppCtx;

133: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
134: {
135:   PetscBool flg;

137:   options->fname[0] = '\0';
138:   PetscOptionsBegin(comm, "", "DMPlex View/Load Test Options", "DMPLEX");
139:   PetscOptionsString("-fname", "The output mesh file", "ex12.c", options->fname, options->fname, sizeof(options->fname), &flg);
140:   PetscOptionsBool("-shell", "Use DMShell to wrap sections", "ex12.c", options->shell, &options->shell, NULL);
141:   PetscOptionsEnd();
142:   return 0;
143: }

145: int main(int argc, char **argv)
146: {
147:   MPI_Comm          comm;
148:   PetscMPIInt       size, rank, mycolor;
149:   const char        exampleDMPlexName[]    = "exampleDMPlex";
150:   const char        exampleSectionDMName[] = "exampleSectionDM";
151:   const char        exampleVecName[]       = "exampleVec";
152:   PetscScalar       constraintValue        = 1.5;
153:   PetscViewerFormat format                 = PETSC_VIEWER_HDF5_PETSC;
154:   AppCtx            user;

157:   PetscInitialize(&argc, &argv, NULL, help);
158:   ProcessOptions(PETSC_COMM_WORLD, &user);
159:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
160:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

163:   /* Save */
164:   mycolor = (PetscMPIInt)(rank >= 2);
165:   MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm);
166:   if (mycolor == 0) {
167:     DM          dm;
168:     PetscViewer viewer;

170:     PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer);
171:     /* Save exampleDMPlex */
172:     {
173:       DM             pdm;
174:       const PetscInt faces[2] = {6, 1};
175:       PetscSF        sf;
176:       PetscInt       overlap = 1;

178:       DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, &dm);
179:       DMPlexDistribute(dm, overlap, &sf, &pdm);
180:       if (pdm) {
181:         DMDestroy(&dm);
182:         dm = pdm;
183:       }
184:       PetscSFDestroy(&sf);
185:       PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
186:       PetscViewerPushFormat(viewer, format);
187:       DMPlexTopologyView(dm, viewer);
188:       DMPlexLabelsView(dm, viewer);
189:       PetscViewerPopFormat(viewer);
190:     }
191:     /* Save coordinates */
192:     PetscViewerPushFormat(viewer, format);
193:     DMPlexCoordinatesView(dm, viewer);
194:     PetscViewerPopFormat(viewer);
195:     /* Save exampleVec */
196:     {
197:       PetscInt     pStart = -1, pEnd = -1;
198:       DM           sdm;
199:       PetscSection section, gsection;
200:       PetscBool    includesConstraints = PETSC_FALSE;
201:       Vec          vec;
202:       PetscScalar *array = NULL;

204:       /* Create section */
205:       PetscSectionCreate(comm, &section);
206:       PetscSectionSetNumFields(section, 2);
207:       DMPlexGetChart(dm, &pStart, &pEnd);
208:       PetscSectionSetChart(section, pStart, pEnd);
209:       switch (rank) {
210:       case 0:
211:         PetscSectionSetDof(section, 3, 2);
212:         PetscSectionSetDof(section, 12, 3);
213:         PetscSectionSetDof(section, 25, 2);
214:         PetscSectionSetConstraintDof(section, 12, 1);
215:         PetscSectionSetFieldDof(section, 3, 0, 2);
216:         PetscSectionSetFieldDof(section, 12, 0, 2);
217:         PetscSectionSetFieldDof(section, 12, 1, 1);
218:         PetscSectionSetFieldDof(section, 25, 1, 2);
219:         PetscSectionSetFieldConstraintDof(section, 12, 1, 1);
220:         break;
221:       case 1:
222:         PetscSectionSetDof(section, 0, 2);
223:         PetscSectionSetDof(section, 1, 1);
224:         PetscSectionSetDof(section, 8, 3);
225:         PetscSectionSetDof(section, 20, 2);
226:         PetscSectionSetConstraintDof(section, 8, 1);
227:         PetscSectionSetFieldDof(section, 0, 0, 2);
228:         PetscSectionSetFieldDof(section, 8, 0, 2);
229:         PetscSectionSetFieldDof(section, 1, 1, 1);
230:         PetscSectionSetFieldDof(section, 8, 1, 1);
231:         PetscSectionSetFieldDof(section, 20, 1, 2);
232:         PetscSectionSetFieldConstraintDof(section, 8, 1, 1);
233:         break;
234:       }
235:       PetscSectionSetUp(section);
236:       {
237:         const PetscInt indices[]  = {2};
238:         const PetscInt indices1[] = {0};

240:         switch (rank) {
241:         case 0:
242:           PetscSectionSetConstraintIndices(section, 12, indices);
243:           PetscSectionSetFieldConstraintIndices(section, 12, 1, indices1);
244:           break;
245:         case 1:
246:           PetscSectionSetConstraintIndices(section, 8, indices);
247:           PetscSectionSetFieldConstraintIndices(section, 8, 1, indices1);
248:           break;
249:         }
250:       }
251:       if (user.shell) {
252:         PetscSF sf;

254:         DMShellCreate(comm, &sdm);
255:         DMGetPointSF(dm, &sf);
256:         DMSetPointSF(sdm, sf);
257:       } else {
258:         DMClone(dm, &sdm);
259:       }
260:       PetscObjectSetName((PetscObject)sdm, exampleSectionDMName);
261:       DMSetLocalSection(sdm, section);
262:       PetscSectionDestroy(&section);
263:       DMPlexSectionView(dm, viewer, sdm);
264:       /* Create global vector */
265:       DMGetGlobalSection(sdm, &gsection);
266:       PetscSectionGetIncludesConstraints(gsection, &includesConstraints);
267:       if (user.shell) {
268:         PetscInt n = -1;

270:         VecCreate(comm, &vec);
271:         if (includesConstraints) PetscSectionGetStorageSize(gsection, &n);
272:         else PetscSectionGetConstrainedStorageSize(gsection, &n);
273:         VecSetSizes(vec, n, PETSC_DECIDE);
274:         VecSetUp(vec);
275:       } else {
276:         DMGetGlobalVector(sdm, &vec);
277:       }
278:       PetscObjectSetName((PetscObject)vec, exampleVecName);
279:       VecGetArrayWrite(vec, &array);
280:       if (includesConstraints) {
281:         switch (rank) {
282:         case 0:
283:           break;
284:         case 1:
285:           array[0] = 1.0;
286:           array[1] = 1.1;
287:           array[2] = 1.2;
288:           array[3] = 1.3;
289:           array[4] = 1.4;
290:           array[5] = 1.5;
291:           array[6] = 1.6;
292:           array[7] = 1.7;
293:           break;
294:         }
295:       } else {
296:         switch (rank) {
297:         case 0:
298:           break;
299:         case 1:
300:           array[0] = 1.0;
301:           array[1] = 1.1;
302:           array[2] = 1.2;
303:           array[3] = 1.3;
304:           array[4] = 1.4;
305:           array[5] = 1.6;
306:           array[6] = 1.7;
307:           break;
308:         }
309:       }
310:       VecRestoreArrayWrite(vec, &array);
311:       DMPlexGlobalVectorView(dm, viewer, sdm, vec);
312:       if (user.shell) {
313:         VecDestroy(&vec);
314:       } else {
315:         DMRestoreGlobalVector(sdm, &vec);
316:       }
317:       DMDestroy(&sdm);
318:     }
319:     PetscViewerDestroy(&viewer);
320:     DMDestroy(&dm);
321:   }
322:   MPI_Comm_free(&comm);
323:   /* Load */
324:   mycolor = (PetscMPIInt)(rank >= 3);
325:   MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm);
326:   if (mycolor == 0) {
327:     DM          dm;
328:     PetscSF     sfXC;
329:     PetscViewer viewer;

331:     PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer);
332:     /* Load exampleDMPlex */
333:     {
334:       PetscSF sfXB, sfBC;

336:       DMCreate(comm, &dm);
337:       DMSetType(dm, DMPLEX);
338:       PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
339:       /* sfXB: X -> B                         */
340:       /* X: set of globalPointNumbers, [0, N) */
341:       /* B: loaded naive in-memory plex       */
342:       PetscViewerPushFormat(viewer, format);
343:       DMPlexTopologyLoad(dm, viewer, &sfXB);
344:       PetscViewerPopFormat(viewer);
345:       PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
346:       {
347:         DM               distributedDM;
348:         PetscInt         overlap = 1;
349:         PetscPartitioner part;

351:         DMPlexGetPartitioner(dm, &part);
352:         PetscPartitionerSetFromOptions(part);
353:         /* sfBC: B -> C                    */
354:         /* B: loaded naive in-memory plex  */
355:         /* C: redistributed good in-memory */
356:         DMPlexDistribute(dm, overlap, &sfBC, &distributedDM);
357:         if (distributedDM) {
358:           DMDestroy(&dm);
359:           dm = distributedDM;
360:         }
361:         PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
362:       }
363:       /* sfXC: X -> C */
364:       PetscSFCompose(sfXB, sfBC, &sfXC);
365:       PetscSFDestroy(&sfXB);
366:       PetscSFDestroy(&sfBC);
367:     }
368:     /* Load labels */
369:     PetscViewerPushFormat(viewer, format);
370:     DMPlexLabelsLoad(dm, viewer, sfXC);
371:     PetscViewerPopFormat(viewer);
372:     /* Load coordinates */
373:     PetscViewerPushFormat(viewer, format);
374:     DMPlexCoordinatesLoad(dm, viewer, sfXC);
375:     PetscViewerPopFormat(viewer);
376:     PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)");
377:     DMViewFromOptions(dm, NULL, "-dm_view");
378:     PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
379:     /* Load exampleVec */
380:     {
381:       DM           sdm;
382:       PetscSection section, gsection;
383:       IS           perm;
384:       PetscBool    includesConstraints = PETSC_FALSE;
385:       Vec          vec;
386:       PetscSF      lsf, gsf;

388:       if (user.shell) {
389:         PetscSF sf;

391:         DMShellCreate(comm, &sdm);
392:         DMGetPointSF(dm, &sf);
393:         DMSetPointSF(sdm, sf);
394:       } else {
395:         DMClone(dm, &sdm);
396:       }
397:       PetscObjectSetName((PetscObject)sdm, exampleSectionDMName);
398:       PetscSectionCreate(comm, &section);
399:       {
400:         PetscInt  pStart = -1, pEnd = -1, p = -1;
401:         PetscInt *pinds = NULL;

403:         DMPlexGetChart(dm, &pStart, &pEnd);
404:         PetscMalloc1(pEnd - pStart, &pinds);
405:         for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p;
406:         if (rank == 2) {
407:           pinds[10] = 20;
408:           pinds[20] = 10;
409:         }
410:         ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm);
411:       }
412:       PetscSectionSetPermutation(section, perm);
413:       ISDestroy(&perm);
414:       DMSetLocalSection(sdm, section);
415:       PetscSectionDestroy(&section);
416:       DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf);
417:       /* Load as local vector */
418:       DMGetLocalSection(sdm, &section);
419:       PetscObjectSetName((PetscObject)section, "Load: local section");
420:       PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm));
421:       PetscSectionGetIncludesConstraints(section, &includesConstraints);
422:       if (user.shell) {
423:         PetscInt m = -1;

425:         VecCreate(comm, &vec);
426:         if (includesConstraints) PetscSectionGetStorageSize(section, &m);
427:         else PetscSectionGetConstrainedStorageSize(section, &m);
428:         VecSetSizes(vec, m, PETSC_DECIDE);
429:         VecSetUp(vec);
430:       } else {
431:         DMGetLocalVector(sdm, &vec);
432:       }
433:       PetscObjectSetName((PetscObject)vec, exampleVecName);
434:       VecSet(vec, constraintValue);
435:       DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec);
436:       PetscSFDestroy(&lsf);
437:       if (user.shell) {
438:         VecView(vec, PETSC_VIEWER_STDOUT_(comm));
439:         VecDestroy(&vec);
440:       } else {
441:         DMRestoreLocalVector(sdm, &vec);
442:       }
443:       /* Load as global vector */
444:       DMGetGlobalSection(sdm, &gsection);
445:       PetscObjectSetName((PetscObject)gsection, "Load: global section");
446:       PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm));
447:       PetscSectionGetIncludesConstraints(gsection, &includesConstraints);
448:       if (user.shell) {
449:         PetscInt m = -1;

451:         VecCreate(comm, &vec);
452:         if (includesConstraints) PetscSectionGetStorageSize(gsection, &m);
453:         else PetscSectionGetConstrainedStorageSize(gsection, &m);
454:         VecSetSizes(vec, m, PETSC_DECIDE);
455:         VecSetUp(vec);
456:       } else {
457:         DMGetGlobalVector(sdm, &vec);
458:       }
459:       PetscObjectSetName((PetscObject)vec, exampleVecName);
460:       DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec);
461:       PetscSFDestroy(&gsf);
462:       VecView(vec, PETSC_VIEWER_STDOUT_(comm));
463:       if (user.shell) {
464:         VecDestroy(&vec);
465:       } else {
466:         DMRestoreGlobalVector(sdm, &vec);
467:       }
468:       DMDestroy(&sdm);
469:     }
470:     PetscViewerDestroy(&viewer);
471:     PetscSFDestroy(&sfXC);
472:     DMDestroy(&dm);
473:   }
474:   MPI_Comm_free(&comm);

476:   /* Finalize */
477:   PetscFinalize();
478:   return 0;
479: }

481: /*TEST

483:   build:
484:     requires: hdf5
485:   testset:
486:     suffix: 0
487:     requires: !complex
488:     nsize: 4
489:     args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail
490:     args: -dm_plex_view_hdf5_storage_version 2.0.0
491:     test:
492:       suffix: parmetis
493:       requires: parmetis
494:       args: -petscpartitioner_type parmetis
495:     test:
496:       suffix: ptscotch
497:       requires: ptscotch
498:       args: -petscpartitioner_type ptscotch

500: TEST*/