xref: /petsc/src/dm/impls/plex/tests/ex55.c (revision 5cca0b87407614a21a26cb11f48462082081e6c3)
1*5cca0b87SVaclav Hapla static char help[] = "Load and save the mesh and fields to HDF5 and ExodusII\n\n";
2*5cca0b87SVaclav Hapla 
3*5cca0b87SVaclav Hapla #include <petscdmplex.h>
4*5cca0b87SVaclav Hapla #include <petscviewerhdf5.h>
5*5cca0b87SVaclav Hapla #include <petscsf.h>
6*5cca0b87SVaclav Hapla 
7*5cca0b87SVaclav Hapla typedef struct {
8*5cca0b87SVaclav Hapla   PetscBool compare;                      /* Compare the meshes using DMPlexEqual() */
9*5cca0b87SVaclav Hapla   PetscBool distribute;                   /* Distribute the mesh */
10*5cca0b87SVaclav Hapla   PetscBool interpolate;                  /* Generate intermediate mesh elements */
11*5cca0b87SVaclav Hapla   char      filename[PETSC_MAX_PATH_LEN]; /* Mesh filename */
12*5cca0b87SVaclav Hapla   PetscViewerFormat format;               /* Format to write and read */
13*5cca0b87SVaclav Hapla   PetscBool second_write_read;            /* Write and read for the 2nd time */
14*5cca0b87SVaclav Hapla } AppCtx;
15*5cca0b87SVaclav Hapla 
16*5cca0b87SVaclav Hapla static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
17*5cca0b87SVaclav Hapla {
18*5cca0b87SVaclav Hapla   PetscErrorCode ierr;
19*5cca0b87SVaclav Hapla 
20*5cca0b87SVaclav Hapla   PetscFunctionBeginUser;
21*5cca0b87SVaclav Hapla   options->compare = PETSC_FALSE;
22*5cca0b87SVaclav Hapla   options->distribute = PETSC_TRUE;
23*5cca0b87SVaclav Hapla   options->interpolate = PETSC_FALSE;
24*5cca0b87SVaclav Hapla   options->filename[0] = '\0';
25*5cca0b87SVaclav Hapla   options->format = PETSC_VIEWER_DEFAULT;
26*5cca0b87SVaclav Hapla   options->second_write_read = PETSC_FALSE;
27*5cca0b87SVaclav Hapla 
28*5cca0b87SVaclav Hapla   ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr);
29*5cca0b87SVaclav Hapla   ierr = PetscOptionsBool("-compare", "Compare the meshes using DMPlexEqual()", "ex5.c", options->compare, &options->compare, NULL);CHKERRQ(ierr);
30*5cca0b87SVaclav Hapla   ierr = PetscOptionsBool("-distribute", "Distribute the mesh", "ex5.c", options->distribute, &options->distribute, NULL);CHKERRQ(ierr);
31*5cca0b87SVaclav Hapla   ierr = PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex5.c", options->interpolate, &options->interpolate, NULL);CHKERRQ(ierr);
32*5cca0b87SVaclav Hapla   ierr = PetscOptionsString("-filename", "The mesh file", "ex5.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr);
33*5cca0b87SVaclav Hapla   ierr = PetscOptionsEnum("-format", "Format to write and read", "ex5.c", PetscViewerFormats, (PetscEnum)options->format, (PetscEnum*)&options->format, NULL);CHKERRQ(ierr);
34*5cca0b87SVaclav Hapla   ierr = PetscOptionsBool("-second_write_read", "Write and read for the 2nd time", "ex5.c", options->second_write_read, &options->second_write_read, NULL);CHKERRQ(ierr);
35*5cca0b87SVaclav Hapla   ierr = PetscOptionsEnd();
36*5cca0b87SVaclav Hapla   PetscFunctionReturn(0);
37*5cca0b87SVaclav Hapla };
38*5cca0b87SVaclav Hapla 
39*5cca0b87SVaclav Hapla static PetscErrorCode DMPlexWriteAndReadHDF5(DM dm, const char filename[], PetscViewerFormat format, const char prefix[], DM *dm_new)
40*5cca0b87SVaclav Hapla {
41*5cca0b87SVaclav Hapla   DM             dmnew;
42*5cca0b87SVaclav Hapla   PetscViewer    v;
43*5cca0b87SVaclav Hapla   PetscErrorCode ierr;
44*5cca0b87SVaclav Hapla 
45*5cca0b87SVaclav Hapla   PetscFunctionBeginUser;
46*5cca0b87SVaclav Hapla   ierr = PetscViewerHDF5Open(PetscObjectComm((PetscObject) dm), filename, FILE_MODE_WRITE, &v);CHKERRQ(ierr);
47*5cca0b87SVaclav Hapla   ierr = PetscViewerPushFormat(v, format);CHKERRQ(ierr);
48*5cca0b87SVaclav Hapla   ierr = DMView(dm, v);CHKERRQ(ierr);
49*5cca0b87SVaclav Hapla 
50*5cca0b87SVaclav Hapla   ierr = PetscViewerFileSetMode(v, FILE_MODE_READ);CHKERRQ(ierr);
51*5cca0b87SVaclav Hapla   ierr = DMCreate(PETSC_COMM_WORLD, &dmnew);CHKERRQ(ierr);
52*5cca0b87SVaclav Hapla   ierr = DMSetType(dmnew, DMPLEX);CHKERRQ(ierr);
53*5cca0b87SVaclav Hapla   ierr = DMSetOptionsPrefix(dmnew, prefix);CHKERRQ(ierr);
54*5cca0b87SVaclav Hapla   ierr = DMLoad(dmnew, v);CHKERRQ(ierr);
55*5cca0b87SVaclav Hapla   ierr = PetscObjectSetName((PetscObject)dmnew,"Mesh_new");CHKERRQ(ierr);
56*5cca0b87SVaclav Hapla 
57*5cca0b87SVaclav Hapla   ierr = PetscViewerPopFormat(v);CHKERRQ(ierr);
58*5cca0b87SVaclav Hapla   ierr = PetscViewerDestroy(&v);CHKERRQ(ierr);
59*5cca0b87SVaclav Hapla   *dm_new = dmnew;
60*5cca0b87SVaclav Hapla   PetscFunctionReturn(0);
61*5cca0b87SVaclav Hapla }
62*5cca0b87SVaclav Hapla 
63*5cca0b87SVaclav Hapla int main(int argc, char **argv)
64*5cca0b87SVaclav Hapla {
65*5cca0b87SVaclav Hapla   DM             dm, dmnew;
66*5cca0b87SVaclav Hapla   PetscPartitioner part;
67*5cca0b87SVaclav Hapla   AppCtx         user;
68*5cca0b87SVaclav Hapla   PetscBool      flg;
69*5cca0b87SVaclav Hapla   PetscErrorCode ierr;
70*5cca0b87SVaclav Hapla 
71*5cca0b87SVaclav Hapla   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
72*5cca0b87SVaclav Hapla   ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
73*5cca0b87SVaclav Hapla   ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, user.filename, user.interpolate, &dm);CHKERRQ(ierr);
74*5cca0b87SVaclav Hapla   ierr = DMSetOptionsPrefix(dm,"orig_");CHKERRQ(ierr);
75*5cca0b87SVaclav Hapla   ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
76*5cca0b87SVaclav Hapla 
77*5cca0b87SVaclav Hapla   if (user.distribute) {
78*5cca0b87SVaclav Hapla     DM dmdist;
79*5cca0b87SVaclav Hapla 
80*5cca0b87SVaclav Hapla     ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr);
81*5cca0b87SVaclav Hapla     ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
82*5cca0b87SVaclav Hapla     ierr = DMPlexDistribute(dm, 0, NULL, &dmdist);CHKERRQ(ierr);
83*5cca0b87SVaclav Hapla     if (dmdist) {
84*5cca0b87SVaclav Hapla       ierr = DMDestroy(&dm);CHKERRQ(ierr);
85*5cca0b87SVaclav Hapla       dm   = dmdist;
86*5cca0b87SVaclav Hapla     }
87*5cca0b87SVaclav Hapla   }
88*5cca0b87SVaclav Hapla 
89*5cca0b87SVaclav Hapla   ierr = DMSetOptionsPrefix(dm,NULL);CHKERRQ(ierr);
90*5cca0b87SVaclav Hapla   ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
91*5cca0b87SVaclav Hapla   ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
92*5cca0b87SVaclav Hapla 
93*5cca0b87SVaclav Hapla   ierr = DMPlexWriteAndReadHDF5(dm, "dmdist.h5", user.format, "new_", &dmnew);CHKERRQ(ierr);
94*5cca0b87SVaclav Hapla 
95*5cca0b87SVaclav Hapla   if (user.second_write_read) {
96*5cca0b87SVaclav Hapla     ierr = DMDestroy(&dm);CHKERRQ(ierr);
97*5cca0b87SVaclav Hapla     dm = dmnew;
98*5cca0b87SVaclav Hapla     ierr = DMPlexWriteAndReadHDF5(dm, "dmdist.h5", user.format, "new_", &dmnew);CHKERRQ(ierr);
99*5cca0b87SVaclav Hapla   }
100*5cca0b87SVaclav Hapla 
101*5cca0b87SVaclav Hapla   ierr = DMViewFromOptions(dmnew, NULL, "-dm_view");CHKERRQ(ierr);
102*5cca0b87SVaclav Hapla   /* TODO: Is it still true? */
103*5cca0b87SVaclav Hapla   /* The NATIVE format for coordiante viewing is killing parallel output, since we have a local vector. Map it to global, and it will work. */
104*5cca0b87SVaclav Hapla 
105*5cca0b87SVaclav Hapla   /* This currently makes sense only for sequential meshes. */
106*5cca0b87SVaclav Hapla   if (user.compare) {
107*5cca0b87SVaclav Hapla     ierr = DMPlexEqual(dmnew, dm, &flg);CHKERRQ(ierr);
108*5cca0b87SVaclav Hapla     if (flg) {ierr = PetscPrintf(PETSC_COMM_WORLD,"DMs equal\n");CHKERRQ(ierr);}
109*5cca0b87SVaclav Hapla     else     {ierr = PetscPrintf(PETSC_COMM_WORLD,"DMs are not equal\n");CHKERRQ(ierr);}
110*5cca0b87SVaclav Hapla   }
111*5cca0b87SVaclav Hapla 
112*5cca0b87SVaclav Hapla   ierr = DMDestroy(&dm);CHKERRQ(ierr);
113*5cca0b87SVaclav Hapla   ierr = DMDestroy(&dmnew);CHKERRQ(ierr);
114*5cca0b87SVaclav Hapla   ierr = PetscFinalize();
115*5cca0b87SVaclav Hapla   return ierr;
116*5cca0b87SVaclav Hapla }
117*5cca0b87SVaclav Hapla 
118*5cca0b87SVaclav Hapla /*TEST
119*5cca0b87SVaclav Hapla   build:
120*5cca0b87SVaclav Hapla     requires: hdf5
121*5cca0b87SVaclav Hapla   # Idempotence of saving/loading
122*5cca0b87SVaclav Hapla   #   Have to replace Exodus file, which is creating uninterpolated edges
123*5cca0b87SVaclav Hapla   test:
124*5cca0b87SVaclav Hapla     suffix: 0
125*5cca0b87SVaclav Hapla     requires: exodusii broken
126*5cca0b87SVaclav Hapla     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/Rect-tri3.exo -dm_view ascii::ascii_info_detail
127*5cca0b87SVaclav Hapla     args: -format hdf5_petsc -compare
128*5cca0b87SVaclav Hapla   test:
129*5cca0b87SVaclav Hapla     suffix: 1
130*5cca0b87SVaclav Hapla     requires: exodusii parmetis !define(PETSC_USE_64BIT_INDICES) broken
131*5cca0b87SVaclav Hapla     nsize: 2
132*5cca0b87SVaclav Hapla     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/Rect-tri3.exo -dm_view ascii::ascii_info_detail
133*5cca0b87SVaclav Hapla     args: -petscpartitioner_type parmetis
134*5cca0b87SVaclav Hapla     args: -format hdf5_petsc -new_dm_view ascii::ascii_info_detail
135*5cca0b87SVaclav Hapla   testset:
136*5cca0b87SVaclav Hapla     requires: exodusii
137*5cca0b87SVaclav Hapla     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo
138*5cca0b87SVaclav Hapla     args: -petscpartitioner_type simple
139*5cca0b87SVaclav Hapla     args: -dm_view ascii::ascii_info_detail
140*5cca0b87SVaclav Hapla     args: -new_dm_view ascii::ascii_info_detail
141*5cca0b87SVaclav Hapla     test:
142*5cca0b87SVaclav Hapla       suffix: 2
143*5cca0b87SVaclav Hapla       nsize: {{1 2 4 8}separate output}
144*5cca0b87SVaclav Hapla       args: -format {{default hdf5_petsc}separate output}
145*5cca0b87SVaclav Hapla       args: -interpolate {{0 1}separate output}
146*5cca0b87SVaclav Hapla     test:
147*5cca0b87SVaclav Hapla       suffix: 2a
148*5cca0b87SVaclav Hapla       nsize: {{1 2 4 8}separate output}
149*5cca0b87SVaclav Hapla       args: -format {{hdf5_xdmf hdf5_viz}separate output}
150*5cca0b87SVaclav Hapla   test:
151*5cca0b87SVaclav Hapla     suffix: 3
152*5cca0b87SVaclav Hapla     requires: exodusii
153*5cca0b87SVaclav Hapla     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -compare
154*5cca0b87SVaclav Hapla 
155*5cca0b87SVaclav Hapla   # Load HDF5 file in XDMF format in parallel, write, read dm1, write, read dm2, and compare dm1 and dm2
156*5cca0b87SVaclav Hapla   test:
157*5cca0b87SVaclav Hapla     suffix: 4
158*5cca0b87SVaclav Hapla     requires: !complex
159*5cca0b87SVaclav Hapla     nsize: {{1 2 3 4 8}}
160*5cca0b87SVaclav Hapla     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5
161*5cca0b87SVaclav Hapla     args: -dm_plex_create_from_hdf5_xdmf -distribute 0 -format hdf5_xdmf -second_write_read -compare
162*5cca0b87SVaclav Hapla 
163*5cca0b87SVaclav Hapla   testset:
164*5cca0b87SVaclav Hapla     # the same data and settings as dm_impls_plex_tests-ex18_9%
165*5cca0b87SVaclav Hapla     requires: hdf5 !complex datafilespath
166*5cca0b87SVaclav Hapla     #TODO DMPlexCheckPointSF() fails for nsize 4
167*5cca0b87SVaclav Hapla     nsize: {{1 2}}
168*5cca0b87SVaclav Hapla     args: -dm_plex_check_symmetry -dm_plex_check_skeleton -dm_plex_check_geometry
169*5cca0b87SVaclav Hapla     args: -filename ${DATAFILESPATH}/meshes/cube-hexahedra-refined.h5 -dm_plex_create_from_hdf5_xdmf -dm_plex_hdf5_topology_path /cells -dm_plex_hdf5_geometry_path /coordinates
170*5cca0b87SVaclav Hapla     args: -format hdf5_xdmf -second_write_read -compare
171*5cca0b87SVaclav Hapla     test:
172*5cca0b87SVaclav Hapla       suffix: 9_hdf5_seqload
173*5cca0b87SVaclav Hapla       args: -distribute -petscpartitioner_type simple
174*5cca0b87SVaclav Hapla       args: -interpolate {{0 1}}
175*5cca0b87SVaclav Hapla       args: -dm_plex_hdf5_force_sequential
176*5cca0b87SVaclav Hapla     test:
177*5cca0b87SVaclav Hapla       suffix: 9_hdf5_seqload_metis
178*5cca0b87SVaclav Hapla       requires: parmetis
179*5cca0b87SVaclav Hapla       args: -distribute -petscpartitioner_type parmetis
180*5cca0b87SVaclav Hapla       args: -interpolate 1
181*5cca0b87SVaclav Hapla       args: -dm_plex_hdf5_force_sequential
182*5cca0b87SVaclav Hapla     test:
183*5cca0b87SVaclav Hapla       suffix: 9_hdf5
184*5cca0b87SVaclav Hapla       args: -interpolate 1 -petscpartitioner_type simple
185*5cca0b87SVaclav Hapla     test:
186*5cca0b87SVaclav Hapla       suffix: 9_hdf5_repart
187*5cca0b87SVaclav Hapla       requires: parmetis
188*5cca0b87SVaclav Hapla       args: -distribute -petscpartitioner_type parmetis
189*5cca0b87SVaclav Hapla       args: -interpolate 1
190*5cca0b87SVaclav Hapla     test:
191*5cca0b87SVaclav Hapla       TODO: Parallel partitioning of uninterpolated meshes not supported
192*5cca0b87SVaclav Hapla       suffix: 9_hdf5_repart_ppu
193*5cca0b87SVaclav Hapla       requires: parmetis
194*5cca0b87SVaclav Hapla       args: -distribute -petscpartitioner_type parmetis
195*5cca0b87SVaclav Hapla       args: -interpolate 0
196*5cca0b87SVaclav Hapla 
197*5cca0b87SVaclav Hapla   # reproduce PetscSFView() crash - fixed, left as regression test
198*5cca0b87SVaclav Hapla   test:
199*5cca0b87SVaclav Hapla     suffix: new_dm_view
200*5cca0b87SVaclav Hapla     requires: exodusii
201*5cca0b87SVaclav Hapla     nsize: 2
202*5cca0b87SVaclav Hapla     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/TwoQuads.exo -new_dm_view ascii:ex5_new_dm_view.log:ascii_info_detail
203*5cca0b87SVaclav Hapla TEST*/
204