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, §ion);
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(§ion);
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, §ion);
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(§ion);
416: DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf);
417: /* Load as local vector */
418: DMGetLocalSection(sdm, §ion);
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*/