xref: /petsc/src/dm/impls/plex/tests/ex18.c (revision 2a27bf02e16e95ccba291fb88507e515be69a77f) 
1c4762a1bSJed Brown static char help[] = "Tests for parallel mesh loading and parallel topological interpolation\n\n";
2c4762a1bSJed Brown 
3c4762a1bSJed Brown #include <petsc/private/dmpleximpl.h>
4c4762a1bSJed Brown /* List of test meshes
5c4762a1bSJed Brown 
6c4762a1bSJed Brown Network
7c4762a1bSJed Brown -------
8c4762a1bSJed Brown Test 0 (2 ranks):
9c4762a1bSJed Brown 
10c4762a1bSJed Brown network=0:
11c4762a1bSJed Brown ---------
12c4762a1bSJed Brown   cell 0   cell 1   cell 2          nCells-1       (edge)
13c4762a1bSJed Brown 0 ------ 1 ------ 2 ------ 3 -- -- v --  -- nCells (vertex)
14c4762a1bSJed Brown 
15c4762a1bSJed Brown   vertex distribution:
16c4762a1bSJed Brown     rank 0: 0 1
17c4762a1bSJed Brown     rank 1: 2 3 ... nCells
18c4762a1bSJed Brown   cell(edge) distribution:
19c4762a1bSJed Brown     rank 0: 0 1
20c4762a1bSJed Brown     rank 1: 2 ... nCells-1
21c4762a1bSJed Brown 
22c4762a1bSJed Brown network=1:
23c4762a1bSJed Brown ---------
24c4762a1bSJed Brown                v2
25c4762a1bSJed Brown                 ^
26c4762a1bSJed Brown                 |
27c4762a1bSJed Brown                cell 2
28c4762a1bSJed Brown                 |
29c4762a1bSJed Brown  v0 --cell 0--> v3--cell 1--> v1
30c4762a1bSJed Brown 
31c4762a1bSJed Brown   vertex distribution:
32c4762a1bSJed Brown     rank 0: 0 1 3
33c4762a1bSJed Brown     rank 1: 2
34c4762a1bSJed Brown   cell(edge) distribution:
35c4762a1bSJed Brown     rank 0: 0 1
36c4762a1bSJed Brown     rank 1: 2
37c4762a1bSJed Brown 
38c4762a1bSJed Brown   example:
39c4762a1bSJed Brown     mpiexec -n 2 ./ex18 -distribute 1 -dim 1 -orig_dm_view -dist_dm_view -dist_dm_view -petscpartitioner_type parmetis -ncells 50
40c4762a1bSJed Brown 
41c4762a1bSJed Brown Triangle
42c4762a1bSJed Brown --------
43c4762a1bSJed Brown Test 0 (2 ranks):
44c4762a1bSJed Brown Two triangles sharing a face
45c4762a1bSJed Brown 
46c4762a1bSJed Brown         2
47c4762a1bSJed Brown       / | \
48c4762a1bSJed Brown      /  |  \
49c4762a1bSJed Brown     /   |   \
50c4762a1bSJed Brown    0  0 | 1  3
51c4762a1bSJed Brown     \   |   /
52c4762a1bSJed Brown      \  |  /
53c4762a1bSJed Brown       \ | /
54c4762a1bSJed Brown         1
55c4762a1bSJed Brown 
56c4762a1bSJed Brown   vertex distribution:
57c4762a1bSJed Brown     rank 0: 0 1
58c4762a1bSJed Brown     rank 1: 2 3
59c4762a1bSJed Brown   cell distribution:
60c4762a1bSJed Brown     rank 0: 0
61c4762a1bSJed Brown     rank 1: 1
62c4762a1bSJed Brown 
63c4762a1bSJed Brown Test 1 (3 ranks):
64c4762a1bSJed Brown Four triangles partitioned across 3 ranks
65c4762a1bSJed Brown 
66c4762a1bSJed Brown    0 _______ 3
67c4762a1bSJed Brown    | \     / |
68c4762a1bSJed Brown    |  \ 1 /  |
69c4762a1bSJed Brown    |   \ /   |
70c4762a1bSJed Brown    | 0  2  2 |
71c4762a1bSJed Brown    |   / \   |
72c4762a1bSJed Brown    |  / 3 \  |
73c4762a1bSJed Brown    | /     \ |
74c4762a1bSJed Brown    1 ------- 4
75c4762a1bSJed Brown 
76c4762a1bSJed Brown   vertex distribution:
77c4762a1bSJed Brown     rank 0: 0 1
78c4762a1bSJed Brown     rank 1: 2 3
79c4762a1bSJed Brown     rank 2: 4
80c4762a1bSJed Brown   cell distribution:
81c4762a1bSJed Brown     rank 0: 0
82c4762a1bSJed Brown     rank 1: 1
83c4762a1bSJed Brown     rank 2: 2 3
84c4762a1bSJed Brown 
85c4762a1bSJed Brown Test 2 (3 ranks):
86c4762a1bSJed Brown Four triangles partitioned across 3 ranks
87c4762a1bSJed Brown 
88c4762a1bSJed Brown    1 _______ 3
89c4762a1bSJed Brown    | \     / |
90c4762a1bSJed Brown    |  \ 1 /  |
91c4762a1bSJed Brown    |   \ /   |
92c4762a1bSJed Brown    | 0  0  2 |
93c4762a1bSJed Brown    |   / \   |
94c4762a1bSJed Brown    |  / 3 \  |
95c4762a1bSJed Brown    | /     \ |
96c4762a1bSJed Brown    2 ------- 4
97c4762a1bSJed Brown 
98c4762a1bSJed Brown   vertex distribution:
99c4762a1bSJed Brown     rank 0: 0 1
100c4762a1bSJed Brown     rank 1: 2 3
101c4762a1bSJed Brown     rank 2: 4
102c4762a1bSJed Brown   cell distribution:
103c4762a1bSJed Brown     rank 0: 0
104c4762a1bSJed Brown     rank 1: 1
105c4762a1bSJed Brown     rank 2: 2 3
106c4762a1bSJed Brown 
107c4762a1bSJed Brown Tetrahedron
108c4762a1bSJed Brown -----------
109c4762a1bSJed Brown Test 0:
110c4762a1bSJed Brown Two tets sharing a face
111c4762a1bSJed Brown 
112c4762a1bSJed Brown  cell   3 _______    cell
113c4762a1bSJed Brown  0    / | \      \   1
114c4762a1bSJed Brown      /  |  \      \
115c4762a1bSJed Brown     /   |   \      \
116c4762a1bSJed Brown    0----|----4-----2
117c4762a1bSJed Brown     \   |   /      /
118c4762a1bSJed Brown      \  |  /      /
119c4762a1bSJed Brown       \ | /      /
120c4762a1bSJed Brown         1-------
121c4762a1bSJed Brown    y
122c4762a1bSJed Brown    | x
123c4762a1bSJed Brown    |/
124c4762a1bSJed Brown    *----z
125c4762a1bSJed Brown 
126c4762a1bSJed Brown   vertex distribution:
127c4762a1bSJed Brown     rank 0: 0 1
128c4762a1bSJed Brown     rank 1: 2 3 4
129c4762a1bSJed Brown   cell distribution:
130c4762a1bSJed Brown     rank 0: 0
131c4762a1bSJed Brown     rank 1: 1
132c4762a1bSJed Brown 
133c4762a1bSJed Brown Quadrilateral
134c4762a1bSJed Brown -------------
135c4762a1bSJed Brown Test 0 (2 ranks):
136c4762a1bSJed Brown Two quads sharing a face
137c4762a1bSJed Brown 
138c4762a1bSJed Brown    3-------2-------5
139c4762a1bSJed Brown    |       |       |
140c4762a1bSJed Brown    |   0   |   1   |
141c4762a1bSJed Brown    |       |       |
142c4762a1bSJed Brown    0-------1-------4
143c4762a1bSJed Brown 
144c4762a1bSJed Brown   vertex distribution:
145c4762a1bSJed Brown     rank 0: 0 1 2
146c4762a1bSJed Brown     rank 1: 3 4 5
147c4762a1bSJed Brown   cell distribution:
148c4762a1bSJed Brown     rank 0: 0
149c4762a1bSJed Brown     rank 1: 1
150c4762a1bSJed Brown 
151c4762a1bSJed Brown TODO Test 1:
152c4762a1bSJed Brown A quad and a triangle sharing a face
153c4762a1bSJed Brown 
154c4762a1bSJed Brown    5-------4
155c4762a1bSJed Brown    |       | \
156c4762a1bSJed Brown    |   0   |  \
157c4762a1bSJed Brown    |       | 1 \
158c4762a1bSJed Brown    2-------3----6
159c4762a1bSJed Brown 
160c4762a1bSJed Brown Hexahedron
161c4762a1bSJed Brown ----------
162c4762a1bSJed Brown Test 0 (2 ranks):
163c4762a1bSJed Brown Two hexes sharing a face
164c4762a1bSJed Brown 
165c4762a1bSJed Brown cell   7-------------6-------------11 cell
166c4762a1bSJed Brown 0     /|            /|            /|     1
167c4762a1bSJed Brown      / |   F1      / |   F7      / |
168c4762a1bSJed Brown     /  |          /  |          /  |
169c4762a1bSJed Brown    4-------------5-------------10  |
170c4762a1bSJed Brown    |   |     F4  |   |     F10 |   |
171c4762a1bSJed Brown    |   |         |   |         |   |
172c4762a1bSJed Brown    |F5 |         |F3 |         |F9 |
173c4762a1bSJed Brown    |   |  F2     |   |   F8    |   |
174c4762a1bSJed Brown    |   3---------|---2---------|---9
175c4762a1bSJed Brown    |  /          |  /          |  /
176c4762a1bSJed Brown    | /   F0      | /    F6     | /
177c4762a1bSJed Brown    |/            |/            |/
178c4762a1bSJed Brown    0-------------1-------------8
179c4762a1bSJed Brown 
180c4762a1bSJed Brown   vertex distribution:
181c4762a1bSJed Brown     rank 0: 0 1 2 3 4 5
182c4762a1bSJed Brown     rank 1: 6 7 8 9 10 11
183c4762a1bSJed Brown   cell distribution:
184c4762a1bSJed Brown     rank 0: 0
185c4762a1bSJed Brown     rank 1: 1
186c4762a1bSJed Brown 
187c4762a1bSJed Brown */
188c4762a1bSJed Brown 
189c4762a1bSJed Brown typedef enum {NONE, CREATE, AFTER_CREATE, AFTER_DISTRIBUTE} InterpType;
190c4762a1bSJed Brown 
191c4762a1bSJed Brown typedef struct {
192c4762a1bSJed Brown   PetscInt   debug;                        /* The debugging level */
193c4762a1bSJed Brown   PetscInt   testNum;                      /* Indicates the mesh to create */
194c4762a1bSJed Brown   PetscInt   dim;                          /* The topological mesh dimension */
195c4762a1bSJed Brown   PetscBool  cellSimplex;                  /* Use simplices or hexes */
196c4762a1bSJed Brown   PetscBool  distribute;                   /* Distribute the mesh */
197c4762a1bSJed Brown   InterpType interpolate;                  /* Interpolate the mesh before or after DMPlexDistribute() */
198c4762a1bSJed Brown   PetscBool  useGenerator;                 /* Construct mesh with a mesh generator */
199c4762a1bSJed Brown   PetscBool  testOrientIF;                 /* Test for different original interface orientations */
200c4762a1bSJed Brown   PetscBool  testHeavy;                    /* Run the heavy PointSF test */
201c4762a1bSJed Brown   PetscBool  customView;                   /* Show results of DMPlexIsInterpolated() etc. */
202c4762a1bSJed Brown   PetscInt   ornt[2];                      /* Orientation of interface on rank 0 and rank 1 */
203c4762a1bSJed Brown   PetscInt   faces[3];                     /* Number of faces per dimension for generator */
204c4762a1bSJed Brown   PetscReal  coords[128];
205c4762a1bSJed Brown   PetscReal  coordsTol;
206c4762a1bSJed Brown   PetscInt   ncoords;
207c4762a1bSJed Brown   PetscInt   pointsToExpand[128];
208c4762a1bSJed Brown   PetscInt   nPointsToExpand;
209c4762a1bSJed Brown   PetscBool  testExpandPointsEmpty;
210c4762a1bSJed Brown   char       filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
211c4762a1bSJed Brown } AppCtx;
212c4762a1bSJed Brown 
213c4762a1bSJed Brown struct _n_PortableBoundary {
214c4762a1bSJed Brown   Vec coordinates;
215c4762a1bSJed Brown   PetscInt depth;
216c4762a1bSJed Brown   PetscSection *sections;
217c4762a1bSJed Brown };
218c4762a1bSJed Brown typedef struct _n_PortableBoundary * PortableBoundary;
219c4762a1bSJed Brown 
220956f8c0dSBarry Smith #if defined(PETSC_USE_LOG)
221c4762a1bSJed Brown static PetscLogStage  stage[3];
222956f8c0dSBarry Smith #endif
223c4762a1bSJed Brown 
224c4762a1bSJed Brown static PetscErrorCode DMPlexCheckPointSFHeavy(DM, PortableBoundary);
225c4762a1bSJed Brown static PetscErrorCode DMPlexSetOrientInterface_Private(DM,PetscBool);
226c4762a1bSJed Brown static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM, PortableBoundary *);
227c4762a1bSJed Brown static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM, IS, PetscSection, IS *);
228c4762a1bSJed Brown 
229c4762a1bSJed Brown static PetscErrorCode PortableBoundaryDestroy(PortableBoundary *bnd)
230c4762a1bSJed Brown {
231c4762a1bSJed Brown   PetscInt       d;
232c4762a1bSJed Brown   PetscErrorCode ierr;
233c4762a1bSJed Brown 
234c4762a1bSJed Brown   PetscFunctionBegin;
235c4762a1bSJed Brown   if (!*bnd) PetscFunctionReturn(0);
236c4762a1bSJed Brown   ierr = VecDestroy(&(*bnd)->coordinates);CHKERRQ(ierr);
237c4762a1bSJed Brown   for (d=0; d < (*bnd)->depth; d++) {
238c4762a1bSJed Brown     ierr = PetscSectionDestroy(&(*bnd)->sections[d]);CHKERRQ(ierr);
239c4762a1bSJed Brown   }
240c4762a1bSJed Brown   ierr = PetscFree((*bnd)->sections);CHKERRQ(ierr);
241c4762a1bSJed Brown   ierr = PetscFree(*bnd);CHKERRQ(ierr);
242c4762a1bSJed Brown   PetscFunctionReturn(0);
243c4762a1bSJed Brown }
244c4762a1bSJed Brown 
245c4762a1bSJed Brown static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
246c4762a1bSJed Brown {
247c4762a1bSJed Brown   const char    *interpTypes[4]  = {"none", "create", "after_create", "after_distribute"};
248c4762a1bSJed Brown   PetscInt       interp=NONE, dim;
249c4762a1bSJed Brown   PetscBool      flg1, flg2;
250c4762a1bSJed Brown   PetscErrorCode ierr;
251c4762a1bSJed Brown 
252c4762a1bSJed Brown   PetscFunctionBegin;
253c4762a1bSJed Brown   options->debug        = 0;
254c4762a1bSJed Brown   options->testNum      = 0;
255c4762a1bSJed Brown   options->dim          = 2;
256c4762a1bSJed Brown   options->cellSimplex  = PETSC_TRUE;
257c4762a1bSJed Brown   options->distribute   = PETSC_FALSE;
258c4762a1bSJed Brown   options->interpolate  = NONE;
259c4762a1bSJed Brown   options->useGenerator = PETSC_FALSE;
260c4762a1bSJed Brown   options->testOrientIF = PETSC_FALSE;
261c4762a1bSJed Brown   options->testHeavy    = PETSC_TRUE;
262c4762a1bSJed Brown   options->customView   = PETSC_FALSE;
263c4762a1bSJed Brown   options->testExpandPointsEmpty = PETSC_FALSE;
264c4762a1bSJed Brown   options->ornt[0]      = 0;
265c4762a1bSJed Brown   options->ornt[1]      = 0;
266c4762a1bSJed Brown   options->faces[0]     = 2;
267c4762a1bSJed Brown   options->faces[1]     = 2;
268c4762a1bSJed Brown   options->faces[2]     = 2;
269c4762a1bSJed Brown   options->filename[0]  = '\0';
270c4762a1bSJed Brown   options->coordsTol    = PETSC_DEFAULT;
271c4762a1bSJed Brown 
272c4762a1bSJed Brown   ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr);
273c4762a1bSJed Brown   ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex18.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr);
274c4762a1bSJed Brown   ierr = PetscOptionsBoundedInt("-testnum", "The mesh to create", "ex18.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr);
275c4762a1bSJed Brown   ierr = PetscOptionsBool("-cell_simplex", "Generate simplices if true, otherwise hexes", "ex18.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr);
276c4762a1bSJed Brown   ierr = PetscOptionsBool("-distribute", "Distribute the mesh", "ex18.c", options->distribute, &options->distribute, NULL);CHKERRQ(ierr);
277c4762a1bSJed Brown   ierr = PetscOptionsEList("-interpolate", "Type of mesh interpolation (none, create, after_create, after_distribute)", "ex18.c", interpTypes, 4, interpTypes[options->interpolate], &interp, NULL);CHKERRQ(ierr);
278c4762a1bSJed Brown   options->interpolate = (InterpType) interp;
279c4762a1bSJed Brown   if (!options->distribute && options->interpolate == AFTER_DISTRIBUTE) SETERRQ(comm, PETSC_ERR_SUP, "-interpolate after_distribute  needs  -distribute 1");
280c4762a1bSJed Brown   ierr = PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex18.c", options->useGenerator, &options->useGenerator, NULL);CHKERRQ(ierr);
281c4762a1bSJed Brown   options->ncoords = 128;
282c4762a1bSJed Brown   ierr = PetscOptionsRealArray("-view_vertices_from_coords", "Print DAG points corresponding to vertices with given coordinates", "ex18.c", options->coords, &options->ncoords, NULL);CHKERRQ(ierr);
283c4762a1bSJed Brown   ierr = PetscOptionsReal("-view_vertices_from_coords_tol", "Tolerance for -view_vertices_from_coords", "ex18.c", options->coordsTol, &options->coordsTol, NULL);CHKERRQ(ierr);
284c4762a1bSJed Brown   options->nPointsToExpand = 128;
285c4762a1bSJed Brown   ierr = PetscOptionsIntArray("-test_expand_points", "Expand given array of DAG point using DMPlexGetConeRecursive() and print results", "ex18.c", options->pointsToExpand, &options->nPointsToExpand, NULL);CHKERRQ(ierr);
286c4762a1bSJed Brown   if (options->nPointsToExpand) {
287c4762a1bSJed Brown     ierr = PetscOptionsBool("-test_expand_points_empty", "For -test_expand_points, rank 0 will have empty input array", "ex18.c", options->testExpandPointsEmpty, &options->testExpandPointsEmpty, NULL);CHKERRQ(ierr);
288c4762a1bSJed Brown   }
289c4762a1bSJed Brown   ierr = PetscOptionsBool("-test_heavy", "Run the heavy PointSF test", "ex18.c", options->testHeavy, &options->testHeavy, NULL);CHKERRQ(ierr);
290c4762a1bSJed Brown   ierr = PetscOptionsBool("-custom_view", "Custom DMPlex view", "ex18.c", options->customView, &options->customView, NULL);CHKERRQ(ierr);
291c4762a1bSJed Brown   ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex18.c", options->dim, &options->dim, &flg1,1,3);CHKERRQ(ierr);
292c4762a1bSJed Brown   dim = 3;
293c4762a1bSJed Brown   ierr = PetscOptionsIntArray("-faces", "Number of faces per dimension", "ex18.c", options->faces, &dim, &flg2);CHKERRQ(ierr);
294c4762a1bSJed Brown   if (flg2) {
295c4762a1bSJed Brown     if (flg1 && dim != options->dim) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "specified -dim %D is not equal to length %D of -faces (note that -dim can be omitted)", options->dim, dim);
296c4762a1bSJed Brown     options->dim = dim;
297c4762a1bSJed Brown   }
298589a23caSBarry Smith   ierr = PetscOptionsString("-filename", "The mesh file", "ex18.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr);
299c4762a1bSJed Brown   ierr = PetscOptionsBoundedInt("-rotate_interface_0", "Rotation (relative orientation) of interface on rank 0; implies -interpolate create -distribute 0", "ex18.c", options->ornt[0], &options->ornt[0], &options->testOrientIF,0);CHKERRQ(ierr);
300c4762a1bSJed Brown   ierr = PetscOptionsBoundedInt("-rotate_interface_1", "Rotation (relative orientation) of interface on rank 1; implies -interpolate create -distribute 0", "ex18.c", options->ornt[1], &options->ornt[1], &flg2,0);CHKERRQ(ierr);
301c4762a1bSJed Brown   if (flg2 != options->testOrientIF) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "neither or both -rotate_interface_0 -rotate_interface_1 must be set");
302c4762a1bSJed Brown   if (options->testOrientIF) {
303c4762a1bSJed Brown     PetscInt i;
304c4762a1bSJed Brown     for (i=0; i<2; i++) {
305c4762a1bSJed Brown       if (options->ornt[i] >= 10) options->ornt[i] = -(options->ornt[i]-10);  /* 11 12 13 become -1 -2 -3 */
306c4762a1bSJed Brown     }
307c4762a1bSJed Brown     options->filename[0]  = 0;
308c4762a1bSJed Brown     options->useGenerator = PETSC_FALSE;
309c4762a1bSJed Brown     options->dim          = 3;
310c4762a1bSJed Brown     options->cellSimplex  = PETSC_TRUE;
311c4762a1bSJed Brown     options->interpolate  = CREATE;
312c4762a1bSJed Brown     options->distribute   = PETSC_FALSE;
313c4762a1bSJed Brown   }
314c4762a1bSJed Brown   ierr = PetscOptionsEnd();
315c4762a1bSJed Brown   PetscFunctionReturn(0);
316c4762a1bSJed Brown }
317c4762a1bSJed Brown 
318c4762a1bSJed Brown static PetscErrorCode CreateMesh_1D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
319c4762a1bSJed Brown {
320c4762a1bSJed Brown   PetscInt       testNum = user->testNum;
321c4762a1bSJed Brown   PetscMPIInt    rank,size;
322c4762a1bSJed Brown   PetscErrorCode ierr;
323c4762a1bSJed Brown   PetscInt       numCorners=2,i;
324c4762a1bSJed Brown   PetscInt       numCells,numVertices,network;
325a4a685f2SJacob Faibussowitsch   PetscInt       *cells;
326c4762a1bSJed Brown   PetscReal      *coords;
327c4762a1bSJed Brown 
328c4762a1bSJed Brown   PetscFunctionBegin;
329ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
330ffc4695bSBarry Smith   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
331c4762a1bSJed Brown   if (size > 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for <=2 processes",testNum);
332c4762a1bSJed Brown 
333c4762a1bSJed Brown   numCells = 3;
334c4762a1bSJed Brown   ierr = PetscOptionsGetInt(NULL, NULL, "-ncells", &numCells, NULL);CHKERRQ(ierr);
335c4762a1bSJed Brown   if (numCells < 3) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test ncells must >=3",numCells);
336c4762a1bSJed Brown 
337c4762a1bSJed Brown   if (size == 1) {
338a4a685f2SJacob Faibussowitsch     PetscReal *dcoords;
339c4762a1bSJed Brown     numVertices = numCells + 1;
340c4762a1bSJed Brown     ierr = PetscMalloc2(2*numCells,&cells,2*numVertices,&dcoords);CHKERRQ(ierr);
341c4762a1bSJed Brown     for (i=0; i<numCells; i++) {
342c4762a1bSJed Brown       cells[2*i] = i; cells[2*i+1] = i + 1;
343c4762a1bSJed Brown       dcoords[2*i] = i; dcoords[2*i+1] = i + 1;
344c4762a1bSJed Brown     }
345c4762a1bSJed Brown 
346a4a685f2SJacob Faibussowitsch     ierr = DMPlexCreateFromCellListPetsc(comm, user->dim, numCells, numVertices, numCorners, PETSC_FALSE, cells, user->dim, dcoords, dm);CHKERRQ(ierr);
347c4762a1bSJed Brown     ierr = PetscFree2(cells,coords);CHKERRQ(ierr);
348c4762a1bSJed Brown     PetscFunctionReturn(0);
349c4762a1bSJed Brown   }
350c4762a1bSJed Brown 
351c4762a1bSJed Brown   network = 0;
352c4762a1bSJed Brown   ierr = PetscOptionsGetInt(NULL, NULL, "-network_case", &network, NULL);CHKERRQ(ierr);
353c4762a1bSJed Brown   if (network == 0) {
354c4762a1bSJed Brown     switch (rank) {
355c4762a1bSJed Brown     case 0:
356c4762a1bSJed Brown     {
357c4762a1bSJed Brown       numCells    = 2;
358c4762a1bSJed Brown       numVertices = numCells;
359c4762a1bSJed Brown       ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr);
360c4762a1bSJed Brown       cells[0] = 0; cells[1] = 1;
361c4762a1bSJed Brown       cells[2] = 1; cells[3] = 2;
362c4762a1bSJed Brown       coords[0] = 0.; coords[1] = 1.;
363c4762a1bSJed Brown       coords[2] = 1.; coords[3] = 2.;
364c4762a1bSJed Brown     }
365c4762a1bSJed Brown     break;
366c4762a1bSJed Brown     case 1:
367c4762a1bSJed Brown     {
368c4762a1bSJed Brown       numCells    -= 2;
369c4762a1bSJed Brown       numVertices = numCells + 1;
370c4762a1bSJed Brown       ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr);
371c4762a1bSJed Brown       for (i=0; i<numCells; i++) {
372c4762a1bSJed Brown         cells[2*i] = 2+i; cells[2*i+1] = 2 + i + 1;
373c4762a1bSJed Brown         coords[2*i] = 2+i; coords[2*i+1] = 2 + i + 1;
374c4762a1bSJed Brown       }
375c4762a1bSJed Brown     }
376c4762a1bSJed Brown     break;
377c4762a1bSJed Brown     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
378c4762a1bSJed Brown     }
379c4762a1bSJed Brown   } else { /* network_case = 1 */
380c4762a1bSJed Brown     /* ----------------------- */
381c4762a1bSJed Brown     switch (rank) {
382c4762a1bSJed Brown     case 0:
383c4762a1bSJed Brown     {
384c4762a1bSJed Brown       numCells    = 2;
385c4762a1bSJed Brown       numVertices = 3;
386c4762a1bSJed Brown       ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr);
387c4762a1bSJed Brown       cells[0] = 0; cells[1] = 3;
388c4762a1bSJed Brown       cells[2] = 3; cells[3] = 1;
389c4762a1bSJed Brown     }
390c4762a1bSJed Brown     break;
391c4762a1bSJed Brown     case 1:
392c4762a1bSJed Brown     {
393c4762a1bSJed Brown       numCells    = 1;
394c4762a1bSJed Brown       numVertices = 1;
395c4762a1bSJed Brown       ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr);
396c4762a1bSJed Brown       cells[0] = 3; cells[1] = 2;
397c4762a1bSJed Brown     }
398c4762a1bSJed Brown     break;
399c4762a1bSJed Brown     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
400c4762a1bSJed Brown     }
401c4762a1bSJed Brown   }
40225b6865aSVaclav Hapla   ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, PETSC_FALSE, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
403c4762a1bSJed Brown   ierr = PetscFree2(cells,coords);CHKERRQ(ierr);
404c4762a1bSJed Brown   PetscFunctionReturn(0);
405c4762a1bSJed Brown }
406c4762a1bSJed Brown 
407c4762a1bSJed Brown static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
408c4762a1bSJed Brown {
409c4762a1bSJed Brown   PetscInt       testNum = user->testNum, p;
410c4762a1bSJed Brown   PetscMPIInt    rank, size;
411c4762a1bSJed Brown   PetscErrorCode ierr;
412c4762a1bSJed Brown 
413c4762a1bSJed Brown   PetscFunctionBegin;
414ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
415ffc4695bSBarry Smith   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
416c4762a1bSJed Brown   switch (testNum) {
417c4762a1bSJed Brown   case 0:
418c4762a1bSJed Brown     if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum);
419c4762a1bSJed Brown     switch (rank) {
420c4762a1bSJed Brown       case 0:
421c4762a1bSJed Brown       {
422c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
423a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {0, 1, 2};
424c4762a1bSJed Brown         PetscReal      coords[4] = {-0.5, 0.5, 0.0, 0.0};
425c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
426c4762a1bSJed Brown 
42725b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
428c4762a1bSJed Brown         for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
429c4762a1bSJed Brown       }
430c4762a1bSJed Brown       break;
431c4762a1bSJed Brown       case 1:
432c4762a1bSJed Brown       {
433c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
434a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {1, 3, 2};
435c4762a1bSJed Brown         PetscReal      coords[4] = {0.0, 1.0, 0.5, 0.5};
436c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
437c4762a1bSJed Brown 
43825b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
439c4762a1bSJed Brown         for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
440c4762a1bSJed Brown       }
441c4762a1bSJed Brown       break;
442c4762a1bSJed Brown       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
443c4762a1bSJed Brown     }
444c4762a1bSJed Brown     break;
445c4762a1bSJed Brown   case 1:
446c4762a1bSJed Brown     if (size != 3) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum);
447c4762a1bSJed Brown     switch (rank) {
448c4762a1bSJed Brown       case 0:
449c4762a1bSJed Brown       {
450c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
451a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {0, 1, 2};
452c4762a1bSJed Brown         PetscReal      coords[4] = {0.0, 1.0, 0.0, 0.0};
453c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
454c4762a1bSJed Brown 
45525b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
456c4762a1bSJed Brown         for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
457c4762a1bSJed Brown       }
458c4762a1bSJed Brown       break;
459c4762a1bSJed Brown       case 1:
460c4762a1bSJed Brown       {
461c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
462a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {0, 2, 3};
463c4762a1bSJed Brown         PetscReal      coords[4] = {0.5, 0.5, 1.0, 1.0};
464c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
465c4762a1bSJed Brown 
46625b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
467c4762a1bSJed Brown         for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
468c4762a1bSJed Brown       }
469c4762a1bSJed Brown       break;
470c4762a1bSJed Brown       case 2:
471c4762a1bSJed Brown       {
472c4762a1bSJed Brown         const PetscInt numCells  = 2, numVertices = 1, numCorners = 3;
473a4a685f2SJacob Faibussowitsch         const PetscInt cells[6]  = {2, 4, 3, 2, 1, 4};
474c4762a1bSJed Brown         PetscReal      coords[2] = {1.0, 0.0};
475c4762a1bSJed Brown         PetscInt       markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
476c4762a1bSJed Brown 
47725b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
478c4762a1bSJed Brown         for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
479c4762a1bSJed Brown       }
480c4762a1bSJed Brown       break;
481c4762a1bSJed Brown       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
482c4762a1bSJed Brown     }
483c4762a1bSJed Brown     break;
484c4762a1bSJed Brown   case 2:
485c4762a1bSJed Brown     if (size != 3) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum);
486c4762a1bSJed Brown     switch (rank) {
487c4762a1bSJed Brown       case 0:
488c4762a1bSJed Brown       {
489c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
490a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {1, 2, 0};
491c4762a1bSJed Brown         PetscReal      coords[4] = {0.5, 0.5, 0.0, 1.0};
492c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
493c4762a1bSJed Brown 
49425b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
495c4762a1bSJed Brown         for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
496c4762a1bSJed Brown       }
497c4762a1bSJed Brown       break;
498c4762a1bSJed Brown       case 1:
499c4762a1bSJed Brown       {
500c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
501a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {1, 0, 3};
502c4762a1bSJed Brown         PetscReal      coords[4] = {0.0, 0.0, 1.0, 1.0};
503c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
504c4762a1bSJed Brown 
50525b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
506c4762a1bSJed Brown         for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
507c4762a1bSJed Brown       }
508c4762a1bSJed Brown       break;
509c4762a1bSJed Brown       case 2:
510c4762a1bSJed Brown       {
511c4762a1bSJed Brown         const PetscInt numCells  = 2, numVertices = 1, numCorners = 3;
512a4a685f2SJacob Faibussowitsch         const PetscInt cells[6]  = {0, 4, 3, 0, 2, 4};
513c4762a1bSJed Brown         PetscReal      coords[2] = {1.0, 0.0};
514c4762a1bSJed Brown         PetscInt       markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
515c4762a1bSJed Brown 
51625b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
517c4762a1bSJed Brown         for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
518c4762a1bSJed Brown       }
519c4762a1bSJed Brown       break;
520c4762a1bSJed Brown       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
521c4762a1bSJed Brown     }
522c4762a1bSJed Brown     break;
523c4762a1bSJed Brown   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum);
524c4762a1bSJed Brown   }
525c4762a1bSJed Brown   PetscFunctionReturn(0);
526c4762a1bSJed Brown }
527c4762a1bSJed Brown 
528c4762a1bSJed Brown static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
529c4762a1bSJed Brown {
530c4762a1bSJed Brown   PetscInt       testNum = user->testNum, p;
531c4762a1bSJed Brown   PetscMPIInt    rank, size;
532c4762a1bSJed Brown   PetscErrorCode ierr;
533c4762a1bSJed Brown 
534c4762a1bSJed Brown   PetscFunctionBegin;
535ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
536ffc4695bSBarry Smith   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
537c4762a1bSJed Brown   switch (testNum) {
538c4762a1bSJed Brown   case 0:
539c4762a1bSJed Brown     if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum);
540c4762a1bSJed Brown     switch (rank) {
541c4762a1bSJed Brown       case 0:
542c4762a1bSJed Brown       {
543c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 4;
544a4a685f2SJacob Faibussowitsch         const PetscInt cells[4]  = {0, 2, 1, 3};
545c4762a1bSJed Brown         PetscReal      coords[6] = {0.0, 0.0, -0.5,  0.0, -0.5, 0.0};
546c4762a1bSJed Brown         PetscInt       markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1};
547c4762a1bSJed Brown 
54825b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
549c4762a1bSJed Brown         for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
550c4762a1bSJed Brown       }
551c4762a1bSJed Brown       break;
552c4762a1bSJed Brown       case 1:
553c4762a1bSJed Brown       {
554c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 3, numCorners = 4;
555a4a685f2SJacob Faibussowitsch         const PetscInt cells[4]  = {1, 2, 4, 3};
556c4762a1bSJed Brown         PetscReal      coords[9] = {1.0, 0.0, 0.0,  0.0, 0.5, 0.0,  0.0, 0.0, 0.5};
557c4762a1bSJed Brown         PetscInt       markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1};
558c4762a1bSJed Brown 
55925b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
560c4762a1bSJed Brown         for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
561c4762a1bSJed Brown       }
562c4762a1bSJed Brown       break;
563c4762a1bSJed Brown       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
564c4762a1bSJed Brown     }
565c4762a1bSJed Brown     break;
566c4762a1bSJed Brown   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum);
567c4762a1bSJed Brown   }
568c4762a1bSJed Brown   if (user->testOrientIF) {
569c4762a1bSJed Brown     PetscInt start;
570c4762a1bSJed Brown     PetscBool reverse;
571c4762a1bSJed Brown     PetscInt ifp[] = {8, 6};
572c4762a1bSJed Brown 
573c4762a1bSJed Brown     ierr = PetscObjectSetName((PetscObject) *dm, "Mesh before orientation");CHKERRQ(ierr);
574c4762a1bSJed Brown     ierr = DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view");CHKERRQ(ierr);
575c4762a1bSJed Brown     /* rotate interface face ifp[rank] by given orientation ornt[rank] */
576c4762a1bSJed Brown     ierr = DMPlexFixFaceOrientations_Translate_Private(user->ornt[rank], &start, &reverse);CHKERRQ(ierr);
577c4762a1bSJed Brown     ierr = DMPlexOrientCell_Internal(*dm, ifp[rank], start, reverse);CHKERRQ(ierr);
578c4762a1bSJed Brown     ierr = DMPlexCheckFaces(*dm, 0);CHKERRQ(ierr);
579c4762a1bSJed Brown     ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr);
580c4762a1bSJed Brown     ierr = PetscPrintf(comm, "Orientation test PASSED\n");CHKERRQ(ierr);
581c4762a1bSJed Brown   }
582c4762a1bSJed Brown   PetscFunctionReturn(0);
583c4762a1bSJed Brown }
584c4762a1bSJed Brown 
585c4762a1bSJed Brown static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
586c4762a1bSJed Brown {
587c4762a1bSJed Brown   PetscInt       testNum = user->testNum, p;
588c4762a1bSJed Brown   PetscMPIInt    rank, size;
589c4762a1bSJed Brown   PetscErrorCode ierr;
590c4762a1bSJed Brown 
591c4762a1bSJed Brown   PetscFunctionBegin;
592ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
593ffc4695bSBarry Smith   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
594c4762a1bSJed Brown   switch (testNum) {
595c4762a1bSJed Brown   case 0:
596c4762a1bSJed Brown     if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum);
597c4762a1bSJed Brown     switch (rank) {
598c4762a1bSJed Brown       case 0:
599c4762a1bSJed Brown       {
600c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 3, numCorners = 4;
601a4a685f2SJacob Faibussowitsch         const PetscInt cells[4]  = {0, 1, 2, 3};
602c4762a1bSJed Brown         PetscReal      coords[6] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0};
603c4762a1bSJed Brown         PetscInt       markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1};
604c4762a1bSJed Brown 
60525b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
606c4762a1bSJed Brown         for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
607c4762a1bSJed Brown       }
608c4762a1bSJed Brown       break;
609c4762a1bSJed Brown       case 1:
610c4762a1bSJed Brown       {
611c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 3, numCorners = 4;
612a4a685f2SJacob Faibussowitsch         const PetscInt cells[4]  = {1, 4, 5, 2};
613c4762a1bSJed Brown         PetscReal      coords[6] = {-0.5, 1.0, 0.5, 0.0, 0.5, 1.0};
614c4762a1bSJed Brown         PetscInt       markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1};
615c4762a1bSJed Brown 
61625b6865aSVaclav Hapla         ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
617c4762a1bSJed Brown         for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
618c4762a1bSJed Brown       }
619c4762a1bSJed Brown       break;
620c4762a1bSJed Brown       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
621c4762a1bSJed Brown     }
622c4762a1bSJed Brown     break;
623c4762a1bSJed Brown   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum);
624c4762a1bSJed Brown   }
625c4762a1bSJed Brown   PetscFunctionReturn(0);
626c4762a1bSJed Brown }
627c4762a1bSJed Brown 
628c4762a1bSJed Brown static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
629c4762a1bSJed Brown {
630c4762a1bSJed Brown   PetscInt       testNum = user->testNum, p;
631c4762a1bSJed Brown   PetscMPIInt    rank, size;
632c4762a1bSJed Brown   PetscErrorCode ierr;
633c4762a1bSJed Brown 
634c4762a1bSJed Brown   PetscFunctionBegin;
635ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
636ffc4695bSBarry Smith   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
637c4762a1bSJed Brown   switch (testNum) {
638c4762a1bSJed Brown   case 0:
639c4762a1bSJed Brown     if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum);
640c4762a1bSJed Brown     switch (rank) {
641c4762a1bSJed Brown     case 0:
642c4762a1bSJed Brown     {
643c4762a1bSJed Brown       const PetscInt numCells  = 1, numVertices = 6, numCorners = 8;
644a4a685f2SJacob Faibussowitsch       const PetscInt cells[8]  = {0, 3, 2, 1, 4, 5, 6, 7};
645c4762a1bSJed Brown       PetscReal      coords[6*3] = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0, -0.5,0.0,1.0, 0.0,0.0,1.0};
646c4762a1bSJed Brown       PetscInt       markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1};
647c4762a1bSJed Brown 
64825b6865aSVaclav Hapla       ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
649c4762a1bSJed Brown       for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
650c4762a1bSJed Brown     }
651c4762a1bSJed Brown     break;
652c4762a1bSJed Brown     case 1:
653c4762a1bSJed Brown     {
654c4762a1bSJed Brown       const PetscInt numCells  = 1, numVertices = 6, numCorners = 8;
655a4a685f2SJacob Faibussowitsch       const PetscInt cells[8]  = {1, 2, 9, 8, 5, 10, 11, 6};
656c4762a1bSJed Brown       PetscReal      coords[6*3] = {0.0,1.0,1.0, -0.5,1.0,1.0, 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0,  0.5,1.0,1.0};
657c4762a1bSJed Brown       PetscInt       markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1};
658c4762a1bSJed Brown 
65925b6865aSVaclav Hapla       ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr);
660c4762a1bSJed Brown       for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
661c4762a1bSJed Brown     }
662c4762a1bSJed Brown     break;
663c4762a1bSJed Brown     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
664c4762a1bSJed Brown     }
665c4762a1bSJed Brown   break;
666c4762a1bSJed Brown   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum);
667c4762a1bSJed Brown   }
668c4762a1bSJed Brown   PetscFunctionReturn(0);
669c4762a1bSJed Brown }
670c4762a1bSJed Brown 
671c4762a1bSJed Brown static PetscErrorCode CustomView(DM dm, PetscViewer v)
672c4762a1bSJed Brown {
673c4762a1bSJed Brown   DMPlexInterpolatedFlag interpolated;
674c4762a1bSJed Brown   PetscBool              distributed;
675c4762a1bSJed Brown   PetscErrorCode         ierr;
676c4762a1bSJed Brown 
677c4762a1bSJed Brown   PetscFunctionBegin;
678c4762a1bSJed Brown   ierr = DMPlexIsDistributed(dm, &distributed);CHKERRQ(ierr);
679c4762a1bSJed Brown   ierr = DMPlexIsInterpolatedCollective(dm, &interpolated);CHKERRQ(ierr);
680c4762a1bSJed Brown   ierr = PetscViewerASCIIPrintf(v, "DMPlexIsDistributed: %s\n", PetscBools[distributed]);CHKERRQ(ierr);
681c4762a1bSJed Brown   ierr = PetscViewerASCIIPrintf(v, "DMPlexIsInterpolatedCollective: %s\n", DMPlexInterpolatedFlags[interpolated]);CHKERRQ(ierr);
682c4762a1bSJed Brown   PetscFunctionReturn(0);
683c4762a1bSJed Brown }
684c4762a1bSJed Brown 
685c4762a1bSJed Brown static PetscErrorCode CreateMeshFromFile(MPI_Comm comm, AppCtx *user, DM *dm, DM *serialDM)
686c4762a1bSJed Brown {
687c4762a1bSJed Brown   const char    *filename       = user->filename;
688c4762a1bSJed Brown   PetscBool      testHeavy      = user->testHeavy;
689c4762a1bSJed Brown   PetscBool      interpCreate   = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE;
690c4762a1bSJed Brown   PetscBool      distributed    = PETSC_FALSE;
691c4762a1bSJed Brown   PetscErrorCode ierr;
692c4762a1bSJed Brown 
693c4762a1bSJed Brown   PetscFunctionBegin;
694c4762a1bSJed Brown   *serialDM = NULL;
695c4762a1bSJed Brown   if (testHeavy && interpCreate) {ierr = DMPlexSetOrientInterface_Private(NULL, PETSC_FALSE);CHKERRQ(ierr);}
696c4762a1bSJed Brown   ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr);
697c4762a1bSJed Brown   ierr = DMPlexCreateFromFile(comm, filename, interpCreate, dm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */
698c4762a1bSJed Brown   ierr = PetscLogStagePop();CHKERRQ(ierr);
699c4762a1bSJed Brown   if (testHeavy && interpCreate) {ierr = DMPlexSetOrientInterface_Private(NULL, PETSC_TRUE);CHKERRQ(ierr);}
700c4762a1bSJed Brown   ierr = DMPlexIsDistributed(*dm, &distributed);CHKERRQ(ierr);
701c4762a1bSJed Brown   ierr = PetscPrintf(comm, "DMPlexCreateFromFile produced %s mesh.\n", distributed ? "distributed" : "serial");CHKERRQ(ierr);
702c4762a1bSJed Brown   if (testHeavy && distributed) {
703c4762a1bSJed Brown     ierr = PetscOptionsSetValue(NULL, "-dm_plex_hdf5_force_sequential", NULL);CHKERRQ(ierr);
704c4762a1bSJed Brown     ierr = DMPlexCreateFromFile(comm, filename, interpCreate, serialDM);CHKERRQ(ierr);
705c4762a1bSJed Brown     ierr = DMPlexIsDistributed(*serialDM, &distributed);CHKERRQ(ierr);
706c4762a1bSJed Brown     if (distributed) SETERRQ(comm, PETSC_ERR_PLIB, "unable to create a serial DM from file");
707c4762a1bSJed Brown   }
708c4762a1bSJed Brown   ierr = DMGetDimension(*dm, &user->dim);CHKERRQ(ierr);
709c4762a1bSJed Brown   PetscFunctionReturn(0);
710c4762a1bSJed Brown }
711c4762a1bSJed Brown 
712c4762a1bSJed Brown static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
713c4762a1bSJed Brown {
714c4762a1bSJed Brown   PetscPartitioner part;
715c4762a1bSJed Brown   PortableBoundary boundary     = NULL;
716c4762a1bSJed Brown   DM             serialDM       = NULL;
717c4762a1bSJed Brown   PetscBool      cellSimplex    = user->cellSimplex;
718c4762a1bSJed Brown   PetscBool      useGenerator   = user->useGenerator;
719c4762a1bSJed Brown   PetscBool      interpCreate   = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE;
720c4762a1bSJed Brown   PetscBool      interpSerial   = user->interpolate == AFTER_CREATE ? PETSC_TRUE : PETSC_FALSE;
721c4762a1bSJed Brown   PetscBool      interpParallel = user->interpolate == AFTER_DISTRIBUTE ? PETSC_TRUE : PETSC_FALSE;
722c4762a1bSJed Brown   PetscBool      testHeavy      = user->testHeavy;
723c4762a1bSJed Brown   PetscMPIInt    rank;
724c4762a1bSJed Brown   PetscErrorCode ierr;
725c4762a1bSJed Brown 
726c4762a1bSJed Brown   PetscFunctionBegin;
727ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
728c4762a1bSJed Brown   if (user->filename[0]) {
729c4762a1bSJed Brown     ierr = CreateMeshFromFile(comm, user, dm, &serialDM);CHKERRQ(ierr);
730c4762a1bSJed Brown   } else if (useGenerator) {
731c4762a1bSJed Brown     ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr);
732c4762a1bSJed Brown     ierr = DMPlexCreateBoxMesh(comm, user->dim, cellSimplex, user->faces, NULL, NULL, NULL, interpCreate, dm);CHKERRQ(ierr);
733c4762a1bSJed Brown     ierr = PetscLogStagePop();CHKERRQ(ierr);
734c4762a1bSJed Brown   } else {
735c4762a1bSJed Brown     ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr);
736c4762a1bSJed Brown     switch (user->dim) {
737c4762a1bSJed Brown     case 1:
738c4762a1bSJed Brown       ierr = CreateMesh_1D(comm, interpCreate, user, dm);CHKERRQ(ierr);
739c4762a1bSJed Brown       break;
740c4762a1bSJed Brown     case 2:
741c4762a1bSJed Brown       if (cellSimplex) {
742c4762a1bSJed Brown         ierr = CreateSimplex_2D(comm, interpCreate, user, dm);CHKERRQ(ierr);
743c4762a1bSJed Brown       } else {
744c4762a1bSJed Brown         ierr = CreateQuad_2D(comm, interpCreate, user, dm);CHKERRQ(ierr);
745c4762a1bSJed Brown       }
746c4762a1bSJed Brown       break;
747c4762a1bSJed Brown     case 3:
748c4762a1bSJed Brown       if (cellSimplex) {
749c4762a1bSJed Brown         ierr = CreateSimplex_3D(comm, interpCreate, user, dm);CHKERRQ(ierr);
750c4762a1bSJed Brown       } else {
751c4762a1bSJed Brown         ierr = CreateHex_3D(comm, interpCreate, user, dm);CHKERRQ(ierr);
752c4762a1bSJed Brown       }
753c4762a1bSJed Brown       break;
754c4762a1bSJed Brown     default:
755c4762a1bSJed Brown       SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", user->dim);
756c4762a1bSJed Brown     }
757c4762a1bSJed Brown     ierr = PetscLogStagePop();CHKERRQ(ierr);
758c4762a1bSJed Brown   }
759c4762a1bSJed Brown   if (user->ncoords % user->dim) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "length of coordinates array %D must be divisable by spatial dimension %D", user->ncoords, user->dim);
760c4762a1bSJed Brown   ierr = PetscObjectSetName((PetscObject) *dm, "Original Mesh");CHKERRQ(ierr);
761c4762a1bSJed Brown   ierr = DMViewFromOptions(*dm, NULL, "-orig_dm_view");CHKERRQ(ierr);
762c4762a1bSJed Brown 
763c4762a1bSJed Brown   if (interpSerial) {
764c4762a1bSJed Brown     DM idm;
765c4762a1bSJed Brown 
766c4762a1bSJed Brown     if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE);CHKERRQ(ierr);}
767c4762a1bSJed Brown     ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr);
768c4762a1bSJed Brown     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */
769c4762a1bSJed Brown     ierr = PetscLogStagePop();CHKERRQ(ierr);
770c4762a1bSJed Brown     if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE);CHKERRQ(ierr);}
771c4762a1bSJed Brown     ierr = DMDestroy(dm);CHKERRQ(ierr);
772c4762a1bSJed Brown     *dm = idm;
773c4762a1bSJed Brown     ierr = PetscObjectSetName((PetscObject) *dm, "Interpolated Mesh");CHKERRQ(ierr);
774c4762a1bSJed Brown     ierr = DMViewFromOptions(*dm, NULL, "-intp_dm_view");CHKERRQ(ierr);
775c4762a1bSJed Brown   }
776c4762a1bSJed Brown 
777c4762a1bSJed Brown   /* Set partitioner options */
778c4762a1bSJed Brown   ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr);
779c4762a1bSJed Brown   if (part) {
780c4762a1bSJed Brown     ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE);CHKERRQ(ierr);
781c4762a1bSJed Brown     ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
782c4762a1bSJed Brown   }
783c4762a1bSJed Brown 
784c4762a1bSJed Brown   if (user->customView) {ierr = CustomView(*dm, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);}
785c4762a1bSJed Brown   if (testHeavy) {
786c4762a1bSJed Brown     PetscBool distributed;
787c4762a1bSJed Brown 
788c4762a1bSJed Brown     ierr = DMPlexIsDistributed(*dm, &distributed);CHKERRQ(ierr);
789c4762a1bSJed Brown     if (!serialDM && !distributed) {
790c4762a1bSJed Brown       serialDM = *dm;
791c4762a1bSJed Brown       ierr = PetscObjectReference((PetscObject)*dm);CHKERRQ(ierr);
792c4762a1bSJed Brown     }
793c4762a1bSJed Brown     if (serialDM) {
794c4762a1bSJed Brown       ierr = DMPlexGetExpandedBoundary_Private(serialDM, &boundary);CHKERRQ(ierr);
795c4762a1bSJed Brown     }
796c4762a1bSJed Brown     if (boundary) {
797c4762a1bSJed Brown       /* check DM which has been created in parallel and already interpolated */
798c4762a1bSJed Brown       ierr = DMPlexCheckPointSFHeavy(*dm, boundary);CHKERRQ(ierr);
799c4762a1bSJed Brown     }
800c4762a1bSJed Brown     /* Orient interface because it could be deliberately skipped above. It is idempotent. */
801c4762a1bSJed Brown     ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr);
802c4762a1bSJed Brown   }
803c4762a1bSJed Brown   if (user->distribute) {
804c4762a1bSJed Brown     DM               pdm = NULL;
805c4762a1bSJed Brown 
806c4762a1bSJed Brown     /* Redistribute mesh over processes using that partitioner */
807c4762a1bSJed Brown     ierr = PetscLogStagePush(stage[1]);CHKERRQ(ierr);
808c4762a1bSJed Brown     ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr);
809c4762a1bSJed Brown     ierr = PetscLogStagePop();CHKERRQ(ierr);
810c4762a1bSJed Brown     if (pdm) {
811c4762a1bSJed Brown       ierr = DMDestroy(dm);CHKERRQ(ierr);
812c4762a1bSJed Brown       *dm  = pdm;
813c4762a1bSJed Brown       ierr = PetscObjectSetName((PetscObject) *dm, "Redistributed Mesh");CHKERRQ(ierr);
814c4762a1bSJed Brown       ierr = DMViewFromOptions(*dm, NULL, "-dist_dm_view");CHKERRQ(ierr);
815c4762a1bSJed Brown     }
816c4762a1bSJed Brown 
817c4762a1bSJed Brown     if (interpParallel) {
818c4762a1bSJed Brown       DM idm;
819c4762a1bSJed Brown 
820c4762a1bSJed Brown       if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE);CHKERRQ(ierr);}
821c4762a1bSJed Brown       ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr);
822c4762a1bSJed Brown       ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */
823c4762a1bSJed Brown       ierr = PetscLogStagePop();CHKERRQ(ierr);
824c4762a1bSJed Brown       if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE);CHKERRQ(ierr);}
825c4762a1bSJed Brown       ierr = DMDestroy(dm);CHKERRQ(ierr);
826c4762a1bSJed Brown       *dm = idm;
827c4762a1bSJed Brown       ierr = PetscObjectSetName((PetscObject) *dm, "Interpolated Redistributed Mesh");CHKERRQ(ierr);
828c4762a1bSJed Brown       ierr = DMViewFromOptions(*dm, NULL, "-intp_dm_view");CHKERRQ(ierr);
829c4762a1bSJed Brown     }
830c4762a1bSJed Brown   }
831c4762a1bSJed Brown   if (testHeavy) {
832c4762a1bSJed Brown     if (boundary) {
833c4762a1bSJed Brown       ierr = DMPlexCheckPointSFHeavy(*dm, boundary);CHKERRQ(ierr);
834c4762a1bSJed Brown     }
835c4762a1bSJed Brown     /* Orient interface because it could be deliberately skipped above. It is idempotent. */
836c4762a1bSJed Brown     ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr);
837c4762a1bSJed Brown   }
838c4762a1bSJed Brown 
839c4762a1bSJed Brown   ierr = PetscObjectSetName((PetscObject) *dm, "Parallel Mesh");CHKERRQ(ierr);
840c4762a1bSJed Brown   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
841c4762a1bSJed Brown   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
842c4762a1bSJed Brown 
843c4762a1bSJed Brown   if (user->customView) {ierr = CustomView(*dm, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);}
844c4762a1bSJed Brown   ierr = DMDestroy(&serialDM);CHKERRQ(ierr);
845c4762a1bSJed Brown   ierr = PortableBoundaryDestroy(&boundary);CHKERRQ(ierr);
846c4762a1bSJed Brown   PetscFunctionReturn(0);
847c4762a1bSJed Brown }
848c4762a1bSJed Brown 
849c4762a1bSJed Brown PETSC_STATIC_INLINE PetscErrorCode coord2str(char buf[], size_t len, PetscInt dim, PetscReal *coords, PetscReal tol)
850c4762a1bSJed Brown {
851c4762a1bSJed Brown   PetscErrorCode ierr;
852c4762a1bSJed Brown 
853c4762a1bSJed Brown   PetscFunctionBegin;
854c4762a1bSJed Brown   if (dim > 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "dim must be less than or equal 3");
855c4762a1bSJed Brown   if (tol >= 1e-3) {
856c4762a1bSJed Brown     switch (dim) {
857c4762a1bSJed Brown       case 1: ierr = PetscSNPrintf(buf,len,"(%12.3f)",(double)coords[0]);CHKERRQ(ierr);
858c4762a1bSJed Brown       case 2: ierr = PetscSNPrintf(buf,len,"(%12.3f, %12.3f)",(double)coords[0],(double)coords[1]);CHKERRQ(ierr);
859c4762a1bSJed Brown       case 3: ierr = PetscSNPrintf(buf,len,"(%12.3f, %12.3f, %12.3f)",(double)coords[0],(double)coords[1],(double)coords[2]);CHKERRQ(ierr);
860c4762a1bSJed Brown     }
861c4762a1bSJed Brown   } else {
862c4762a1bSJed Brown     switch (dim) {
863c4762a1bSJed Brown       case 1: ierr = PetscSNPrintf(buf,len,"(%12.6f)",(double)coords[0]);CHKERRQ(ierr);
864c4762a1bSJed Brown       case 2: ierr = PetscSNPrintf(buf,len,"(%12.6f, %12.6f)",(double)coords[0],(double)coords[1]);CHKERRQ(ierr);
865c4762a1bSJed Brown       case 3: ierr = PetscSNPrintf(buf,len,"(%12.6f, %12.6f, %12.6f)",(double)coords[0],(double)coords[1],(double)coords[2]);CHKERRQ(ierr);
866c4762a1bSJed Brown     }
867c4762a1bSJed Brown   }
868c4762a1bSJed Brown   PetscFunctionReturn(0);
869c4762a1bSJed Brown }
870c4762a1bSJed Brown 
871c4762a1bSJed Brown static PetscErrorCode ViewVerticesFromCoords(DM dm, PetscInt npoints, PetscReal coords[], PetscReal tol, PetscViewer viewer)
872c4762a1bSJed Brown {
873c4762a1bSJed Brown   PetscInt       dim, i;
874c4762a1bSJed Brown   PetscInt       *points;
875c4762a1bSJed Brown   char           coordstr[128];
876c4762a1bSJed Brown   MPI_Comm       comm;
877c4762a1bSJed Brown   PetscMPIInt    rank;
878c4762a1bSJed Brown   PetscErrorCode ierr;
879c4762a1bSJed Brown 
880c4762a1bSJed Brown   PetscFunctionBegin;
881c4762a1bSJed Brown   ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
882ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
883c4762a1bSJed Brown   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
884c4762a1bSJed Brown   ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
885c4762a1bSJed Brown   ierr = PetscMalloc1(npoints, &points);CHKERRQ(ierr);
886c4762a1bSJed Brown   ierr = DMPlexFindVertices(dm, npoints, coords, tol, points);CHKERRQ(ierr);
887c4762a1bSJed Brown   for (i=0; i < npoints; i++) {
888c4762a1bSJed Brown     ierr = coord2str(coordstr, sizeof(coordstr), dim, &coords[i*dim], tol);CHKERRQ(ierr);
889c4762a1bSJed Brown     if (!rank && i) {ierr = PetscViewerASCIISynchronizedPrintf(viewer, "-----\n");CHKERRQ(ierr);}
890c4762a1bSJed Brown     ierr = PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %s --> points[%D] = %D\n", rank, coordstr, i, points[i]);CHKERRQ(ierr);
891c4762a1bSJed Brown     ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
892c4762a1bSJed Brown   }
893c4762a1bSJed Brown   ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
894c4762a1bSJed Brown   ierr = PetscFree(points);CHKERRQ(ierr);
895c4762a1bSJed Brown   PetscFunctionReturn(0);
896c4762a1bSJed Brown }
897c4762a1bSJed Brown 
898c4762a1bSJed Brown static PetscErrorCode TestExpandPoints(DM dm, AppCtx *user)
899c4762a1bSJed Brown {
900c4762a1bSJed Brown   IS                is;
901c4762a1bSJed Brown   PetscSection      *sects;
902c4762a1bSJed Brown   IS                *iss;
903c4762a1bSJed Brown   PetscInt          d,depth;
904c4762a1bSJed Brown   PetscMPIInt       rank;
905c4762a1bSJed Brown   PetscErrorCode    ierr;
906c4762a1bSJed Brown   PetscViewer       viewer=PETSC_VIEWER_STDOUT_WORLD, sviewer;
907c4762a1bSJed Brown 
908c4762a1bSJed Brown   PetscFunctionBegin;
909ffc4695bSBarry Smith   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRMPI(ierr);
910c4762a1bSJed Brown   if (user->testExpandPointsEmpty && !rank) {
911c4762a1bSJed Brown     ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &is);CHKERRQ(ierr);
912c4762a1bSJed Brown   } else {
913c4762a1bSJed Brown     ierr = ISCreateGeneral(PETSC_COMM_SELF, user->nPointsToExpand, user->pointsToExpand, PETSC_USE_POINTER, &is);CHKERRQ(ierr);
914c4762a1bSJed Brown   }
915c4762a1bSJed Brown   ierr = DMPlexGetConeRecursive(dm, is, &depth, &iss, &sects);CHKERRQ(ierr);
916c4762a1bSJed Brown   ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr);
917c4762a1bSJed Brown   ierr = PetscViewerASCIIPrintf(sviewer, "[%d] ==========================\n",rank);CHKERRQ(ierr);
918c4762a1bSJed Brown   for (d=depth-1; d>=0; d--) {
919c4762a1bSJed Brown     IS          checkIS;
920c4762a1bSJed Brown     PetscBool   flg;
921c4762a1bSJed Brown 
922c4762a1bSJed Brown     ierr = PetscViewerASCIIPrintf(sviewer, "depth %D ---------------\n",d);CHKERRQ(ierr);
923c4762a1bSJed Brown     ierr = PetscSectionView(sects[d], sviewer);CHKERRQ(ierr);
924c4762a1bSJed Brown     ierr = ISView(iss[d], sviewer);CHKERRQ(ierr);
925c4762a1bSJed Brown     /* check reverse operation */
926c4762a1bSJed Brown     if (d < depth-1) {
927c4762a1bSJed Brown       ierr = DMPlexExpandedConesToFaces_Private(dm, iss[d], sects[d], &checkIS);CHKERRQ(ierr);
928c4762a1bSJed Brown       ierr = ISEqualUnsorted(checkIS, iss[d+1], &flg);CHKERRQ(ierr);
929c4762a1bSJed Brown       if (!flg) SETERRQ(PetscObjectComm((PetscObject) checkIS), PETSC_ERR_PLIB, "DMPlexExpandedConesToFaces_Private produced wrong IS");
930c4762a1bSJed Brown       ierr = ISDestroy(&checkIS);CHKERRQ(ierr);
931c4762a1bSJed Brown     }
932c4762a1bSJed Brown   }
933c4762a1bSJed Brown   ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr);
934c4762a1bSJed Brown   ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
935c4762a1bSJed Brown   ierr = DMPlexRestoreConeRecursive(dm, is, &depth, &iss, &sects);CHKERRQ(ierr);
936c4762a1bSJed Brown   ierr = ISDestroy(&is);CHKERRQ(ierr);
937c4762a1bSJed Brown   PetscFunctionReturn(0);
938c4762a1bSJed Brown }
939c4762a1bSJed Brown 
940c4762a1bSJed Brown static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM dm, IS is, PetscSection section, IS *newis)
941c4762a1bSJed Brown {
942c4762a1bSJed Brown   PetscInt          n,n1,ncone,numCoveredPoints,o,p,q,start,end;
943c4762a1bSJed Brown   const PetscInt    *coveredPoints;
944c4762a1bSJed Brown   const PetscInt    *arr, *cone;
945c4762a1bSJed Brown   PetscInt          *newarr;
946c4762a1bSJed Brown   PetscErrorCode    ierr;
947c4762a1bSJed Brown 
948c4762a1bSJed Brown   PetscFunctionBegin;
949c4762a1bSJed Brown   ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr);
950c4762a1bSJed Brown   ierr = PetscSectionGetStorageSize(section, &n1);CHKERRQ(ierr);
951c4762a1bSJed Brown   ierr = PetscSectionGetChart(section, &start, &end);CHKERRQ(ierr);
952c4762a1bSJed Brown   if (n != n1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "IS size = %D != %D = section storage size\n", n, n1);
953c4762a1bSJed Brown   ierr = ISGetIndices(is, &arr);CHKERRQ(ierr);
954c4762a1bSJed Brown   ierr = PetscMalloc1(end-start, &newarr);CHKERRQ(ierr);
955c4762a1bSJed Brown   for (q=start; q<end; q++) {
956c4762a1bSJed Brown     ierr = PetscSectionGetDof(section, q, &ncone);CHKERRQ(ierr);
957c4762a1bSJed Brown     ierr = PetscSectionGetOffset(section, q, &o);CHKERRQ(ierr);
958c4762a1bSJed Brown     cone = &arr[o];
959c4762a1bSJed Brown     if (ncone == 1) {
960c4762a1bSJed Brown       numCoveredPoints = 1;
961c4762a1bSJed Brown       p = cone[0];
962c4762a1bSJed Brown     } else {
963c4762a1bSJed Brown       PetscInt i;
964c4762a1bSJed Brown       p = PETSC_MAX_INT;
965c4762a1bSJed Brown       for (i=0; i<ncone; i++) if (cone[i] < 0) {p = -1; break;}
966c4762a1bSJed Brown       if (p >= 0) {
967c4762a1bSJed Brown         ierr = DMPlexGetJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr);
968c4762a1bSJed Brown         if (numCoveredPoints > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "more than one covered points for section point %D",q);
969c4762a1bSJed Brown         if (numCoveredPoints) p = coveredPoints[0];
970c4762a1bSJed Brown         else                  p = -2;
971c4762a1bSJed Brown         ierr = DMPlexRestoreJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr);
972c4762a1bSJed Brown       }
973c4762a1bSJed Brown     }
974c4762a1bSJed Brown     newarr[q-start] = p;
975c4762a1bSJed Brown   }
976c4762a1bSJed Brown   ierr = ISRestoreIndices(is, &arr);CHKERRQ(ierr);
977c4762a1bSJed Brown   ierr = ISCreateGeneral(PETSC_COMM_SELF, end-start, newarr, PETSC_OWN_POINTER, newis);CHKERRQ(ierr);
978c4762a1bSJed Brown   PetscFunctionReturn(0);
979c4762a1bSJed Brown }
980c4762a1bSJed Brown 
981c4762a1bSJed Brown static PetscErrorCode DMPlexExpandedVerticesToFaces_Private(DM dm, IS boundary_expanded_is, PetscInt depth, PetscSection sections[], IS *boundary_is)
982c4762a1bSJed Brown {
983c4762a1bSJed Brown   PetscInt          d;
984c4762a1bSJed Brown   IS                is,newis;
985c4762a1bSJed Brown   PetscErrorCode    ierr;
986c4762a1bSJed Brown 
987c4762a1bSJed Brown   PetscFunctionBegin;
988c4762a1bSJed Brown   is = boundary_expanded_is;
989c4762a1bSJed Brown   ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr);
990c4762a1bSJed Brown   for (d = 0; d < depth-1; ++d) {
991c4762a1bSJed Brown     ierr = DMPlexExpandedConesToFaces_Private(dm, is, sections[d], &newis);CHKERRQ(ierr);
992c4762a1bSJed Brown     ierr = ISDestroy(&is);CHKERRQ(ierr);
993c4762a1bSJed Brown     is = newis;
994c4762a1bSJed Brown   }
995c4762a1bSJed Brown   *boundary_is = is;
996c4762a1bSJed Brown   PetscFunctionReturn(0);
997c4762a1bSJed Brown }
998c4762a1bSJed Brown 
999c4762a1bSJed Brown #define CHKERRQI(incall,ierr) if (ierr) {incall = PETSC_FALSE; }
1000c4762a1bSJed Brown 
1001c4762a1bSJed Brown static PetscErrorCode DMLabelViewFromOptionsOnComm_Private(DMLabel label, const char optionname[], MPI_Comm comm)
1002c4762a1bSJed Brown {
1003c4762a1bSJed Brown   PetscErrorCode    ierr;
1004c4762a1bSJed Brown   PetscViewer       viewer;
1005c4762a1bSJed Brown   PetscBool         flg;
1006c4762a1bSJed Brown   static PetscBool  incall = PETSC_FALSE;
1007c4762a1bSJed Brown   PetscViewerFormat format;
1008c4762a1bSJed Brown 
1009c4762a1bSJed Brown   PetscFunctionBegin;
1010c4762a1bSJed Brown   if (incall) PetscFunctionReturn(0);
1011c4762a1bSJed Brown   incall = PETSC_TRUE;
1012c4762a1bSJed Brown   ierr   = PetscOptionsGetViewer(comm,((PetscObject)label)->options,((PetscObject)label)->prefix,optionname,&viewer,&format,&flg);CHKERRQI(incall,ierr);
1013c4762a1bSJed Brown   if (flg) {
1014c4762a1bSJed Brown     ierr = PetscViewerPushFormat(viewer,format);CHKERRQI(incall,ierr);
1015c4762a1bSJed Brown     ierr = DMLabelView(label, viewer);CHKERRQI(incall,ierr);
1016c4762a1bSJed Brown     ierr = PetscViewerPopFormat(viewer);CHKERRQI(incall,ierr);
1017c4762a1bSJed Brown     ierr = PetscViewerDestroy(&viewer);CHKERRQI(incall,ierr);
1018c4762a1bSJed Brown   }
1019c4762a1bSJed Brown   incall = PETSC_FALSE;
1020c4762a1bSJed Brown   PetscFunctionReturn(0);
1021c4762a1bSJed Brown }
1022c4762a1bSJed Brown 
1023c4762a1bSJed Brown /* TODO: this is hotfixing DMLabelGetStratumIS() - it should be fixed systematically instead */
1024c4762a1bSJed Brown PETSC_STATIC_INLINE PetscErrorCode DMLabelGetStratumISOnComm_Private(DMLabel label, PetscInt value, MPI_Comm comm, IS *is)
1025c4762a1bSJed Brown {
1026c4762a1bSJed Brown   IS                tmpis;
1027c4762a1bSJed Brown   PetscErrorCode    ierr;
1028c4762a1bSJed Brown 
1029c4762a1bSJed Brown   PetscFunctionBegin;
1030c4762a1bSJed Brown   ierr = DMLabelGetStratumIS(label, value, &tmpis);CHKERRQ(ierr);
1031c4762a1bSJed Brown   if (!tmpis) {ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &tmpis);CHKERRQ(ierr);}
1032c4762a1bSJed Brown   ierr = ISOnComm(tmpis, comm, PETSC_COPY_VALUES, is);CHKERRQ(ierr);
1033c4762a1bSJed Brown   ierr = ISDestroy(&tmpis);CHKERRQ(ierr);
1034c4762a1bSJed Brown   PetscFunctionReturn(0);
1035c4762a1bSJed Brown }
1036c4762a1bSJed Brown 
1037c4762a1bSJed Brown /* currently only for simple PetscSection without fields or constraints */
1038c4762a1bSJed Brown static PetscErrorCode PetscSectionReplicate_Private(MPI_Comm comm, PetscMPIInt rootrank, PetscSection sec0, PetscSection *secout)
1039c4762a1bSJed Brown {
1040c4762a1bSJed Brown   PetscSection      sec;
1041c4762a1bSJed Brown   PetscInt          chart[2], p;
1042c4762a1bSJed Brown   PetscInt          *dofarr;
1043c4762a1bSJed Brown   PetscMPIInt       rank;
1044c4762a1bSJed Brown   PetscErrorCode    ierr;
1045c4762a1bSJed Brown 
1046c4762a1bSJed Brown   PetscFunctionBegin;
1047ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
1048c4762a1bSJed Brown   if (rank == rootrank) {
1049c4762a1bSJed Brown     ierr = PetscSectionGetChart(sec0, &chart[0], &chart[1]);CHKERRQ(ierr);
1050c4762a1bSJed Brown   }
1051ffc4695bSBarry Smith   ierr = MPI_Bcast(chart, 2, MPIU_INT, rootrank, comm);CHKERRMPI(ierr);
1052c4762a1bSJed Brown   ierr = PetscMalloc1(chart[1]-chart[0], &dofarr);CHKERRQ(ierr);
1053c4762a1bSJed Brown   if (rank == rootrank) {
1054c4762a1bSJed Brown     for (p = chart[0]; p < chart[1]; p++) {
1055c4762a1bSJed Brown       ierr = PetscSectionGetDof(sec0, p, &dofarr[p-chart[0]]);CHKERRQ(ierr);
1056c4762a1bSJed Brown     }
1057c4762a1bSJed Brown   }
1058ffc4695bSBarry Smith   ierr = MPI_Bcast(dofarr, chart[1]-chart[0], MPIU_INT, rootrank, comm);CHKERRMPI(ierr);
1059c4762a1bSJed Brown   ierr = PetscSectionCreate(comm, &sec);CHKERRQ(ierr);
1060c4762a1bSJed Brown   ierr = PetscSectionSetChart(sec, chart[0], chart[1]);CHKERRQ(ierr);
1061c4762a1bSJed Brown   for (p = chart[0]; p < chart[1]; p++) {
1062c4762a1bSJed Brown     ierr = PetscSectionSetDof(sec, p, dofarr[p-chart[0]]);CHKERRQ(ierr);
1063c4762a1bSJed Brown   }
1064c4762a1bSJed Brown   ierr = PetscSectionSetUp(sec);CHKERRQ(ierr);
1065c4762a1bSJed Brown   ierr = PetscFree(dofarr);CHKERRQ(ierr);
1066c4762a1bSJed Brown   *secout = sec;
1067c4762a1bSJed Brown   PetscFunctionReturn(0);
1068c4762a1bSJed Brown }
1069c4762a1bSJed Brown 
1070c4762a1bSJed Brown static PetscErrorCode VecToPetscReal_Private(Vec vec, PetscReal *rvals[])
1071c4762a1bSJed Brown {
1072c4762a1bSJed Brown   PetscInt          n;
1073c4762a1bSJed Brown   const PetscScalar *svals;
1074c4762a1bSJed Brown   PetscErrorCode    ierr;
1075c4762a1bSJed Brown 
1076c4762a1bSJed Brown   PetscFunctionBegin;
1077c4762a1bSJed Brown   ierr = VecGetLocalSize(vec, &n);CHKERRQ(ierr);
1078c4762a1bSJed Brown   ierr = VecGetArrayRead(vec, &svals);CHKERRQ(ierr);
1079c4762a1bSJed Brown   ierr = PetscMalloc1(n, rvals);CHKERRQ(ierr);
1080c4762a1bSJed Brown #if defined(PETSC_USE_COMPLEX)
1081c4762a1bSJed Brown   {
1082c4762a1bSJed Brown     PetscInt i;
1083c4762a1bSJed Brown     for (i=0; i<n; i++) (*rvals)[i] = PetscRealPart(svals[i]);
1084c4762a1bSJed Brown   }
1085c4762a1bSJed Brown #else
1086c4762a1bSJed Brown   ierr = PetscMemcpy(*rvals, svals, n*sizeof(PetscReal));
1087c4762a1bSJed Brown #endif
1088c4762a1bSJed Brown   ierr = VecRestoreArrayRead(vec, &svals);CHKERRQ(ierr);
1089c4762a1bSJed Brown   PetscFunctionReturn(0);
1090c4762a1bSJed Brown }
1091c4762a1bSJed Brown 
1092c4762a1bSJed Brown static PetscErrorCode DMPlexExpandedVerticesCoordinatesToFaces_Private(DM ipdm, PortableBoundary bnd, IS *face_is)
1093c4762a1bSJed Brown {
1094c4762a1bSJed Brown   PetscInt            dim, ncoords, npoints;
1095c4762a1bSJed Brown   PetscReal           *rcoords;
1096c4762a1bSJed Brown   PetscInt            *points;
1097c4762a1bSJed Brown   IS                  faces_expanded_is;
1098c4762a1bSJed Brown   PetscErrorCode      ierr;
1099c4762a1bSJed Brown 
1100c4762a1bSJed Brown   PetscFunctionBegin;
1101c4762a1bSJed Brown   ierr = DMGetDimension(ipdm, &dim);CHKERRQ(ierr);
1102c4762a1bSJed Brown   ierr = VecGetLocalSize(bnd->coordinates, &ncoords);CHKERRQ(ierr);
1103c4762a1bSJed Brown   ierr = VecToPetscReal_Private(bnd->coordinates, &rcoords);CHKERRQ(ierr);
1104c4762a1bSJed Brown   npoints = ncoords / dim;
1105c4762a1bSJed Brown   ierr = PetscMalloc1(npoints, &points);CHKERRQ(ierr);
1106c4762a1bSJed Brown   ierr = DMPlexFindVertices(ipdm, npoints, rcoords, 0.0, points);CHKERRQ(ierr);
1107c4762a1bSJed Brown   ierr = ISCreateGeneral(PETSC_COMM_SELF, npoints, points, PETSC_OWN_POINTER, &faces_expanded_is);CHKERRQ(ierr);
1108c4762a1bSJed Brown   ierr = DMPlexExpandedVerticesToFaces_Private(ipdm, faces_expanded_is, bnd->depth, bnd->sections, face_is);CHKERRQ(ierr);
1109c4762a1bSJed Brown   ierr = PetscFree(rcoords);CHKERRQ(ierr);
1110c4762a1bSJed Brown   ierr = ISDestroy(&faces_expanded_is);CHKERRQ(ierr);
1111c4762a1bSJed Brown   PetscFunctionReturn(0);
1112c4762a1bSJed Brown }
1113c4762a1bSJed Brown 
1114c4762a1bSJed Brown /* hack disabling DMPlexOrientInterface() call in DMPlexInterpolate() via -dm_plex_interpolate_orient_interfaces option */
1115c4762a1bSJed Brown static PetscErrorCode DMPlexSetOrientInterface_Private(DM dm, PetscBool enable)
1116c4762a1bSJed Brown {
1117c4762a1bSJed Brown   PetscOptions      options = NULL;
1118c4762a1bSJed Brown   const char        *prefix = NULL;
1119c4762a1bSJed Brown   const char        opt[] = "-dm_plex_interpolate_orient_interfaces";
1120c4762a1bSJed Brown   char              prefix_opt[512];
1121c4762a1bSJed Brown   PetscBool         flg, set;
1122c4762a1bSJed Brown   static PetscBool  wasSetTrue = PETSC_FALSE;
1123c4762a1bSJed Brown   PetscErrorCode    ierr;
1124c4762a1bSJed Brown 
1125c4762a1bSJed Brown   PetscFunctionBegin;
1126c4762a1bSJed Brown   if (dm) {
1127c4762a1bSJed Brown     ierr = PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);CHKERRQ(ierr);
1128c4762a1bSJed Brown     options = ((PetscObject)dm)->options;
1129c4762a1bSJed Brown   }
1130c4762a1bSJed Brown   ierr = PetscStrcpy(prefix_opt, "-");CHKERRQ(ierr);
1131c4762a1bSJed Brown   ierr = PetscStrlcat(prefix_opt, prefix, sizeof(prefix_opt));CHKERRQ(ierr);
1132c4762a1bSJed Brown   ierr = PetscStrlcat(prefix_opt, &opt[1], sizeof(prefix_opt));CHKERRQ(ierr);
1133c4762a1bSJed Brown   ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr);
1134c4762a1bSJed Brown   if (!enable) {
1135c4762a1bSJed Brown     if (set && flg) wasSetTrue = PETSC_TRUE;
1136c4762a1bSJed Brown     ierr = PetscOptionsSetValue(options, prefix_opt, "0");CHKERRQ(ierr);
1137c4762a1bSJed Brown   } else if (set && !flg) {
1138c4762a1bSJed Brown     if (wasSetTrue) {
1139c4762a1bSJed Brown       ierr = PetscOptionsSetValue(options, prefix_opt, "1");CHKERRQ(ierr);
1140c4762a1bSJed Brown     } else {
1141c4762a1bSJed Brown       /* default is PETSC_TRUE */
1142c4762a1bSJed Brown       ierr = PetscOptionsClearValue(options, prefix_opt);CHKERRQ(ierr);
1143c4762a1bSJed Brown     }
1144c4762a1bSJed Brown     wasSetTrue = PETSC_FALSE;
1145c4762a1bSJed Brown   }
114676bd3646SJed Brown   if (PetscDefined(USE_DEBUG)) {
1147c4762a1bSJed Brown     ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr);
1148c4762a1bSJed Brown     if (PetscUnlikely(set && flg != enable)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "PetscOptionsSetValue did not have the desired effect");
1149c4762a1bSJed Brown   }
1150c4762a1bSJed Brown   PetscFunctionReturn(0);
1151c4762a1bSJed Brown }
1152c4762a1bSJed Brown 
1153c4762a1bSJed Brown /* get coordinate description of the whole-domain boundary */
1154c4762a1bSJed Brown static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM dm, PortableBoundary *boundary)
1155c4762a1bSJed Brown {
1156c4762a1bSJed Brown   PortableBoundary       bnd0, bnd;
1157c4762a1bSJed Brown   MPI_Comm               comm;
1158c4762a1bSJed Brown   DM                     idm;
1159c4762a1bSJed Brown   DMLabel                label;
1160c4762a1bSJed Brown   PetscInt               d;
1161c4762a1bSJed Brown   const char             boundaryName[] = "DMPlexDistributeInterpolateMarkInterface_boundary";
1162c4762a1bSJed Brown   IS                     boundary_is;
1163c4762a1bSJed Brown   IS                     *boundary_expanded_iss;
1164c4762a1bSJed Brown   PetscMPIInt            rootrank = 0;
1165c4762a1bSJed Brown   PetscMPIInt            rank, size;
1166c4762a1bSJed Brown   PetscInt               value = 1;
1167c4762a1bSJed Brown   DMPlexInterpolatedFlag intp;
1168c4762a1bSJed Brown   PetscBool              flg;
1169c4762a1bSJed Brown   PetscErrorCode         ierr;
1170c4762a1bSJed Brown 
1171c4762a1bSJed Brown   PetscFunctionBegin;
1172c4762a1bSJed Brown   ierr = PetscNew(&bnd);CHKERRQ(ierr);
1173c4762a1bSJed Brown   ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
1174ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
1175ffc4695bSBarry Smith   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
1176c4762a1bSJed Brown   ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr);
1177c4762a1bSJed Brown   if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "serial DM (all points on one rank) needed");
1178c4762a1bSJed Brown 
1179c4762a1bSJed Brown   /* interpolate serial DM if not yet interpolated */
1180c4762a1bSJed Brown   ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr);
1181c4762a1bSJed Brown   if (intp == DMPLEX_INTERPOLATED_FULL) {
1182c4762a1bSJed Brown     idm = dm;
1183c4762a1bSJed Brown     ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);
1184c4762a1bSJed Brown   } else {
1185c4762a1bSJed Brown     ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
1186c4762a1bSJed Brown     ierr = DMViewFromOptions(idm, NULL, "-idm_view");CHKERRQ(ierr);
1187c4762a1bSJed Brown   }
1188c4762a1bSJed Brown 
1189c4762a1bSJed Brown   /* mark whole-domain boundary of the serial DM */
1190c4762a1bSJed Brown   ierr = DMLabelCreate(PETSC_COMM_SELF, boundaryName, &label);CHKERRQ(ierr);
1191c4762a1bSJed Brown   ierr = DMAddLabel(idm, label);CHKERRQ(ierr);
1192c4762a1bSJed Brown   ierr = DMPlexMarkBoundaryFaces(idm, value, label);CHKERRQ(ierr);
1193c4762a1bSJed Brown   ierr = DMLabelViewFromOptionsOnComm_Private(label, "-idm_boundary_view", comm);CHKERRQ(ierr);
1194c4762a1bSJed Brown   ierr = DMLabelGetStratumIS(label, value, &boundary_is);CHKERRQ(ierr);
1195c4762a1bSJed Brown 
1196c4762a1bSJed Brown   /* translate to coordinates */
1197c4762a1bSJed Brown   ierr = PetscNew(&bnd0);CHKERRQ(ierr);
1198c4762a1bSJed Brown   ierr = DMGetCoordinatesLocalSetUp(idm);CHKERRQ(ierr);
1199c4762a1bSJed Brown   if (rank == rootrank) {
1200c4762a1bSJed Brown     ierr = DMPlexGetConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr);
1201c4762a1bSJed Brown     ierr = DMGetCoordinatesLocalTuple(dm, boundary_expanded_iss[0], NULL, &bnd0->coordinates);CHKERRQ(ierr);
1202c4762a1bSJed Brown     /* self-check */
1203c4762a1bSJed Brown     {
1204c4762a1bSJed Brown       IS is0;
1205c4762a1bSJed Brown       ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(idm, bnd0, &is0);CHKERRQ(ierr);
1206c4762a1bSJed Brown       ierr = ISEqual(is0, boundary_is, &flg);CHKERRQ(ierr);
1207c4762a1bSJed Brown       if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMPlexExpandedVerticesCoordinatesToFaces_Private produced a wrong IS");
1208c4762a1bSJed Brown       ierr = ISDestroy(&is0);CHKERRQ(ierr);
1209c4762a1bSJed Brown     }
1210c4762a1bSJed Brown   } else {
1211c4762a1bSJed Brown     ierr = VecCreateSeq(PETSC_COMM_SELF, 0, &bnd0->coordinates);CHKERRQ(ierr);
1212c4762a1bSJed Brown   }
1213c4762a1bSJed Brown 
1214c4762a1bSJed Brown   {
1215c4762a1bSJed Brown     Vec         tmp;
1216c4762a1bSJed Brown     VecScatter  sc;
1217c4762a1bSJed Brown     IS          xis;
1218c4762a1bSJed Brown     PetscInt    n;
1219c4762a1bSJed Brown 
1220c4762a1bSJed Brown     /* just convert seq vectors to mpi vector */
1221c4762a1bSJed Brown     ierr = VecGetLocalSize(bnd0->coordinates, &n);CHKERRQ(ierr);
1222ffc4695bSBarry Smith     ierr = MPI_Bcast(&n, 1, MPIU_INT, rootrank, comm);CHKERRMPI(ierr);
1223c4762a1bSJed Brown     if (rank == rootrank) {
1224c4762a1bSJed Brown       ierr = VecCreateMPI(comm, n, n, &tmp);CHKERRQ(ierr);
1225c4762a1bSJed Brown     } else {
1226c4762a1bSJed Brown       ierr = VecCreateMPI(comm, 0, n, &tmp);CHKERRQ(ierr);
1227c4762a1bSJed Brown     }
1228c4762a1bSJed Brown     ierr = VecCopy(bnd0->coordinates, tmp);CHKERRQ(ierr);
1229c4762a1bSJed Brown     ierr = VecDestroy(&bnd0->coordinates);CHKERRQ(ierr);
1230c4762a1bSJed Brown     bnd0->coordinates = tmp;
1231c4762a1bSJed Brown 
1232c4762a1bSJed Brown     /* replicate coordinates from root rank to all ranks */
1233c4762a1bSJed Brown     ierr = VecCreateMPI(comm, n, n*size, &bnd->coordinates);CHKERRQ(ierr);
1234c4762a1bSJed Brown     ierr = ISCreateStride(comm, n, 0, 1, &xis);CHKERRQ(ierr);
1235c4762a1bSJed Brown     ierr = VecScatterCreate(bnd0->coordinates, xis, bnd->coordinates, NULL, &sc);CHKERRQ(ierr);
1236c4762a1bSJed Brown     ierr = VecScatterBegin(sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
1237c4762a1bSJed Brown     ierr = VecScatterEnd(  sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
1238c4762a1bSJed Brown     ierr = VecScatterDestroy(&sc);CHKERRQ(ierr);
1239c4762a1bSJed Brown     ierr = ISDestroy(&xis);CHKERRQ(ierr);
1240c4762a1bSJed Brown   }
1241c4762a1bSJed Brown   bnd->depth = bnd0->depth;
1242ffc4695bSBarry Smith   ierr = MPI_Bcast(&bnd->depth, 1, MPIU_INT, rootrank, comm);CHKERRMPI(ierr);
1243c4762a1bSJed Brown   ierr = PetscMalloc1(bnd->depth, &bnd->sections);CHKERRQ(ierr);
1244c4762a1bSJed Brown   for (d=0; d<bnd->depth; d++) {
1245c4762a1bSJed Brown     ierr = PetscSectionReplicate_Private(comm, rootrank, (rank == rootrank) ? bnd0->sections[d] : NULL, &bnd->sections[d]);CHKERRQ(ierr);
1246c4762a1bSJed Brown   }
1247c4762a1bSJed Brown 
1248c4762a1bSJed Brown   if (rank == rootrank) {
1249c4762a1bSJed Brown     ierr = DMPlexRestoreConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr);
1250c4762a1bSJed Brown   }
1251c4762a1bSJed Brown   ierr = PortableBoundaryDestroy(&bnd0);CHKERRQ(ierr);
1252c4762a1bSJed Brown   ierr = DMRemoveLabelBySelf(idm, &label, PETSC_TRUE);CHKERRQ(ierr);
1253c4762a1bSJed Brown   ierr = DMLabelDestroy(&label);CHKERRQ(ierr);
1254c4762a1bSJed Brown   ierr = ISDestroy(&boundary_is);CHKERRQ(ierr);
1255c4762a1bSJed Brown   ierr = DMDestroy(&idm);CHKERRQ(ierr);
1256c4762a1bSJed Brown   *boundary = bnd;
1257c4762a1bSJed Brown   PetscFunctionReturn(0);
1258c4762a1bSJed Brown }
1259c4762a1bSJed Brown 
1260c4762a1bSJed Brown /* get faces of inter-partition interface */
1261c4762a1bSJed Brown static PetscErrorCode DMPlexGetInterfaceFaces_Private(DM ipdm, IS boundary_faces_is, IS *interface_faces_is)
1262c4762a1bSJed Brown {
1263c4762a1bSJed Brown   MPI_Comm               comm;
1264c4762a1bSJed Brown   DMLabel                label;
1265c4762a1bSJed Brown   IS                     part_boundary_faces_is;
1266c4762a1bSJed Brown   const char             partBoundaryName[] = "DMPlexDistributeInterpolateMarkInterface_partBoundary";
1267c4762a1bSJed Brown   PetscInt               value = 1;
1268c4762a1bSJed Brown   DMPlexInterpolatedFlag intp;
1269c4762a1bSJed Brown   PetscErrorCode         ierr;
1270c4762a1bSJed Brown 
1271c4762a1bSJed Brown   PetscFunctionBegin;
1272c4762a1bSJed Brown   ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr);
1273c4762a1bSJed Brown   ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr);
1274c4762a1bSJed Brown   if (intp != DMPLEX_INTERPOLATED_FULL) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex");
1275c4762a1bSJed Brown 
1276c4762a1bSJed Brown   /* get ipdm partition boundary (partBoundary) */
1277c4762a1bSJed Brown   ierr = DMLabelCreate(PETSC_COMM_SELF, partBoundaryName, &label);CHKERRQ(ierr);
1278c4762a1bSJed Brown   ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr);
1279c4762a1bSJed Brown   ierr = DMPlexMarkBoundaryFaces(ipdm, value, label);CHKERRQ(ierr);
1280c4762a1bSJed Brown   ierr = DMLabelViewFromOptionsOnComm_Private(label, "-ipdm_part_boundary_view", comm);CHKERRQ(ierr);
1281c4762a1bSJed Brown   ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, &part_boundary_faces_is);CHKERRQ(ierr);
1282c4762a1bSJed Brown   ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr);
1283c4762a1bSJed Brown   ierr = DMLabelDestroy(&label);CHKERRQ(ierr);
1284c4762a1bSJed Brown 
1285c4762a1bSJed Brown   /* remove ipdm whole-domain boundary (boundary_faces_is) from ipdm partition boundary (part_boundary_faces_is), resulting just in inter-partition interface */
1286c4762a1bSJed Brown   ierr = ISDifference(part_boundary_faces_is,boundary_faces_is,interface_faces_is);CHKERRQ(ierr);
1287c4762a1bSJed Brown   ierr = ISDestroy(&part_boundary_faces_is);CHKERRQ(ierr);
1288c4762a1bSJed Brown   PetscFunctionReturn(0);
1289c4762a1bSJed Brown }
1290c4762a1bSJed Brown 
1291c4762a1bSJed Brown /* compute inter-partition interface including edges and vertices */
1292c4762a1bSJed Brown static PetscErrorCode DMPlexComputeCompleteInterface_Private(DM ipdm, IS interface_faces_is, IS *interface_is)
1293c4762a1bSJed Brown {
1294c4762a1bSJed Brown   DMLabel                label;
1295c4762a1bSJed Brown   PetscInt               value = 1;
1296c4762a1bSJed Brown   const char             interfaceName[] = "DMPlexDistributeInterpolateMarkInterface_interface";
1297c4762a1bSJed Brown   DMPlexInterpolatedFlag intp;
1298c4762a1bSJed Brown   MPI_Comm               comm;
1299c4762a1bSJed Brown   PetscErrorCode         ierr;
1300c4762a1bSJed Brown 
1301c4762a1bSJed Brown   PetscFunctionBegin;
1302c4762a1bSJed Brown   ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr);
1303c4762a1bSJed Brown   ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr);
1304c4762a1bSJed Brown   if (intp != DMPLEX_INTERPOLATED_FULL) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex");
1305c4762a1bSJed Brown 
1306c4762a1bSJed Brown   ierr = DMLabelCreate(PETSC_COMM_SELF, interfaceName, &label);CHKERRQ(ierr);
1307c4762a1bSJed Brown   ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr);
1308c4762a1bSJed Brown   ierr = DMLabelSetStratumIS(label, value, interface_faces_is);CHKERRQ(ierr);
1309c4762a1bSJed Brown   ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_faces_view", comm);CHKERRQ(ierr);
1310c4762a1bSJed Brown   ierr = DMPlexLabelComplete(ipdm, label);CHKERRQ(ierr);
1311c4762a1bSJed Brown   ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_view", comm);CHKERRQ(ierr);
1312c4762a1bSJed Brown   ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, interface_is);CHKERRQ(ierr);
1313c4762a1bSJed Brown   ierr = PetscObjectSetName((PetscObject)*interface_is, "interface_is");CHKERRQ(ierr);
1314c4762a1bSJed Brown   ierr = ISViewFromOptions(*interface_is, NULL, "-interface_is_view");CHKERRQ(ierr);
1315c4762a1bSJed Brown   ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr);
1316c4762a1bSJed Brown   ierr = DMLabelDestroy(&label);CHKERRQ(ierr);
1317c4762a1bSJed Brown   PetscFunctionReturn(0);
1318c4762a1bSJed Brown }
1319c4762a1bSJed Brown 
1320c4762a1bSJed Brown static PetscErrorCode PointSFGetOutwardInterfacePoints(PetscSF sf, IS *is)
1321c4762a1bSJed Brown {
1322c4762a1bSJed Brown   PetscInt        n;
1323c4762a1bSJed Brown   const PetscInt  *arr;
1324c4762a1bSJed Brown   PetscErrorCode  ierr;
1325c4762a1bSJed Brown 
1326c4762a1bSJed Brown   PetscFunctionBegin;
1327c4762a1bSJed Brown   ierr = PetscSFGetGraph(sf, NULL, &n, &arr, NULL);CHKERRQ(ierr);
1328c4762a1bSJed Brown   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_USE_POINTER, is);CHKERRQ(ierr);
1329c4762a1bSJed Brown   PetscFunctionReturn(0);
1330c4762a1bSJed Brown }
1331c4762a1bSJed Brown 
1332c4762a1bSJed Brown static PetscErrorCode PointSFGetInwardInterfacePoints(PetscSF sf, IS *is)
1333c4762a1bSJed Brown {
1334c4762a1bSJed Brown   PetscInt        n;
1335c4762a1bSJed Brown   const PetscInt  *rootdegree;
1336c4762a1bSJed Brown   PetscInt        *arr;
1337c4762a1bSJed Brown   PetscErrorCode  ierr;
1338c4762a1bSJed Brown 
1339c4762a1bSJed Brown   PetscFunctionBegin;
1340c4762a1bSJed Brown   ierr = PetscSFSetUp(sf);CHKERRQ(ierr);
1341c4762a1bSJed Brown   ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr);
1342c4762a1bSJed Brown   ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr);
1343c4762a1bSJed Brown   ierr = PetscSFComputeMultiRootOriginalNumbering(sf, rootdegree, &n, &arr);CHKERRQ(ierr);
1344c4762a1bSJed Brown   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_OWN_POINTER, is);CHKERRQ(ierr);
1345c4762a1bSJed Brown   PetscFunctionReturn(0);
1346c4762a1bSJed Brown }
1347c4762a1bSJed Brown 
1348c4762a1bSJed Brown static PetscErrorCode PointSFGetInterfacePoints_Private(PetscSF pointSF, IS *is)
1349c4762a1bSJed Brown {
1350c4762a1bSJed Brown   IS pointSF_out_is, pointSF_in_is;
1351c4762a1bSJed Brown   PetscErrorCode  ierr;
1352c4762a1bSJed Brown 
1353c4762a1bSJed Brown   PetscFunctionBegin;
1354c4762a1bSJed Brown   ierr = PointSFGetOutwardInterfacePoints(pointSF, &pointSF_out_is);CHKERRQ(ierr);
1355c4762a1bSJed Brown   ierr = PointSFGetInwardInterfacePoints(pointSF, &pointSF_in_is);CHKERRQ(ierr);
1356c4762a1bSJed Brown   ierr = ISExpand(pointSF_out_is, pointSF_in_is, is);CHKERRQ(ierr);
1357c4762a1bSJed Brown   ierr = ISDestroy(&pointSF_out_is);CHKERRQ(ierr);
1358c4762a1bSJed Brown   ierr = ISDestroy(&pointSF_in_is);CHKERRQ(ierr);
1359c4762a1bSJed Brown   PetscFunctionReturn(0);
1360c4762a1bSJed Brown }
1361c4762a1bSJed Brown 
1362c4762a1bSJed Brown #define MYCHKERRQ(ierr) do {if (PetscUnlikely(ierr)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PointSF is wrong. Unable to show details!");} while (0)
1363c4762a1bSJed Brown 
1364c4762a1bSJed Brown static PetscErrorCode ViewPointsWithType_Internal(DM dm, IS pointsIS, PetscViewer v)
1365c4762a1bSJed Brown {
1366c4762a1bSJed Brown   DMLabel         label;
1367c4762a1bSJed Brown   PetscSection    coordsSection;
1368c4762a1bSJed Brown   Vec             coordsVec;
1369c4762a1bSJed Brown   PetscScalar     *coordsScalar;
1370c4762a1bSJed Brown   PetscInt        coneSize, depth, dim, i, p, npoints;
1371c4762a1bSJed Brown   const PetscInt  *points;
1372c4762a1bSJed Brown   PetscErrorCode  ierr;
1373c4762a1bSJed Brown 
1374c4762a1bSJed Brown   PetscFunctionBegin;
1375c4762a1bSJed Brown   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1376c4762a1bSJed Brown   ierr = DMGetCoordinateSection(dm, &coordsSection);CHKERRQ(ierr);
1377c4762a1bSJed Brown   ierr = DMGetCoordinatesLocal(dm, &coordsVec);CHKERRQ(ierr);
1378c4762a1bSJed Brown   ierr = VecGetArray(coordsVec, &coordsScalar);CHKERRQ(ierr);
1379c4762a1bSJed Brown   ierr = ISGetLocalSize(pointsIS, &npoints);CHKERRQ(ierr);
1380c4762a1bSJed Brown   ierr = ISGetIndices(pointsIS, &points);CHKERRQ(ierr);
1381c4762a1bSJed Brown   ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr);
1382c4762a1bSJed Brown   ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr);
1383c4762a1bSJed Brown   for (i=0; i<npoints; i++) {
1384c4762a1bSJed Brown     p = points[i];
1385c4762a1bSJed Brown     ierr = DMLabelGetValue(label, p, &depth);CHKERRQ(ierr);
1386c4762a1bSJed Brown     if (!depth) {
1387c4762a1bSJed Brown       PetscInt        c, n, o;
1388c4762a1bSJed Brown       PetscReal       coords[3];
1389c4762a1bSJed Brown       char            coordstr[128];
1390c4762a1bSJed Brown 
1391c4762a1bSJed Brown       ierr = PetscSectionGetDof(coordsSection, p, &n);CHKERRQ(ierr);
1392c4762a1bSJed Brown       ierr = PetscSectionGetOffset(coordsSection, p, &o);CHKERRQ(ierr);
1393c4762a1bSJed Brown       for (c=0; c<n; c++) coords[c] = PetscRealPart(coordsScalar[o+c]);
1394c4762a1bSJed Brown       ierr = coord2str(coordstr, sizeof(coordstr), n, coords, 1.0);CHKERRQ(ierr);
1395c4762a1bSJed Brown       ierr = PetscViewerASCIISynchronizedPrintf(v, "vertex %D w/ coordinates %s\n", p, coordstr);CHKERRQ(ierr);
1396c4762a1bSJed Brown     } else {
1397c4762a1bSJed Brown       char            entityType[16];
1398c4762a1bSJed Brown 
1399c4762a1bSJed Brown       switch (depth) {
1400c4762a1bSJed Brown         case 1: ierr = PetscStrcpy(entityType, "edge");CHKERRQ(ierr); break;
1401c4762a1bSJed Brown         case 2: ierr = PetscStrcpy(entityType, "face");CHKERRQ(ierr); break;
1402c4762a1bSJed Brown         case 3: ierr = PetscStrcpy(entityType, "cell");CHKERRQ(ierr); break;
1403*2a27bf02SStefano Zampini         default: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for depth <= 3");
1404c4762a1bSJed Brown       }
1405c4762a1bSJed Brown       if (depth == dim && dim < 3) {
1406c4762a1bSJed Brown         ierr = PetscStrlcat(entityType, " (cell)", sizeof(entityType));CHKERRQ(ierr);
1407c4762a1bSJed Brown       }
1408c4762a1bSJed Brown       ierr = PetscViewerASCIISynchronizedPrintf(v, "%s %D\n", entityType, p);CHKERRQ(ierr);
1409c4762a1bSJed Brown     }
1410c4762a1bSJed Brown     ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
1411c4762a1bSJed Brown     if (coneSize) {
1412c4762a1bSJed Brown       const PetscInt *cone;
1413c4762a1bSJed Brown       IS             coneIS;
1414c4762a1bSJed Brown 
1415c4762a1bSJed Brown       ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
1416c4762a1bSJed Brown       ierr = ISCreateGeneral(PETSC_COMM_SELF, coneSize, cone, PETSC_USE_POINTER, &coneIS);CHKERRQ(ierr);
1417c4762a1bSJed Brown       ierr = ViewPointsWithType_Internal(dm, coneIS, v);CHKERRQ(ierr);
1418c4762a1bSJed Brown       ierr = ISDestroy(&coneIS);CHKERRQ(ierr);
1419c4762a1bSJed Brown     }
1420c4762a1bSJed Brown   }
1421c4762a1bSJed Brown   ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr);
1422c4762a1bSJed Brown   ierr = VecRestoreArray(coordsVec, &coordsScalar);CHKERRQ(ierr);
1423c4762a1bSJed Brown   ierr = ISRestoreIndices(pointsIS, &points);CHKERRQ(ierr);
1424c4762a1bSJed Brown   PetscFunctionReturn(0);
1425c4762a1bSJed Brown }
1426c4762a1bSJed Brown 
1427c4762a1bSJed Brown static PetscErrorCode ViewPointsWithType(DM dm, IS points, PetscViewer v)
1428c4762a1bSJed Brown {
1429c4762a1bSJed Brown   PetscBool       flg;
1430c4762a1bSJed Brown   PetscInt        npoints;
1431c4762a1bSJed Brown   PetscMPIInt     rank;
1432c4762a1bSJed Brown   PetscErrorCode  ierr;
1433c4762a1bSJed Brown 
1434c4762a1bSJed Brown   PetscFunctionBegin;
1435c4762a1bSJed Brown   ierr = PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &flg);CHKERRQ(ierr);
1436c4762a1bSJed Brown   if (!flg) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for ASCII viewer");CHKERRQ(ierr);
1437ffc4695bSBarry Smith   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);CHKERRMPI(ierr);
1438c4762a1bSJed Brown   ierr = PetscViewerASCIIPushSynchronized(v);CHKERRQ(ierr);
1439c4762a1bSJed Brown   ierr = ISGetLocalSize(points, &npoints);CHKERRQ(ierr);
1440c4762a1bSJed Brown   if (npoints) {
1441c4762a1bSJed Brown     ierr = PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", rank);CHKERRQ(ierr);
1442c4762a1bSJed Brown     ierr = ViewPointsWithType_Internal(dm, points, v);CHKERRQ(ierr);
1443c4762a1bSJed Brown   }
1444c4762a1bSJed Brown   ierr = PetscViewerFlush(v);CHKERRQ(ierr);
1445c4762a1bSJed Brown   ierr = PetscViewerASCIIPopSynchronized(v);CHKERRQ(ierr);
1446c4762a1bSJed Brown   PetscFunctionReturn(0);
1447c4762a1bSJed Brown }
1448c4762a1bSJed Brown 
1449c4762a1bSJed Brown static PetscErrorCode DMPlexComparePointSFWithInterface_Private(DM ipdm, IS interface_is)
1450c4762a1bSJed Brown {
1451c4762a1bSJed Brown   PetscSF         pointsf;
1452c4762a1bSJed Brown   IS              pointsf_is;
1453c4762a1bSJed Brown   PetscBool       flg;
1454c4762a1bSJed Brown   MPI_Comm        comm;
1455c4762a1bSJed Brown   PetscMPIInt     size;
1456c4762a1bSJed Brown   PetscErrorCode  ierr;
1457c4762a1bSJed Brown 
1458c4762a1bSJed Brown   PetscFunctionBegin;
1459c4762a1bSJed Brown   ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr);
1460ffc4695bSBarry Smith   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
1461c4762a1bSJed Brown   ierr = DMGetPointSF(ipdm, &pointsf);CHKERRQ(ierr);
1462c4762a1bSJed Brown   if (pointsf) {
1463c4762a1bSJed Brown     PetscInt nroots;
1464c4762a1bSJed Brown     ierr = PetscSFGetGraph(pointsf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
1465c4762a1bSJed Brown     if (nroots < 0) pointsf = NULL; /* uninitialized SF */
1466c4762a1bSJed Brown   }
1467c4762a1bSJed Brown   if (!pointsf) {
1468c4762a1bSJed Brown     PetscInt N=0;
1469c4762a1bSJed Brown     if (interface_is) {ierr = ISGetSize(interface_is, &N);CHKERRQ(ierr);}
1470c4762a1bSJed Brown     if (N) SETERRQ(comm, PETSC_ERR_PLIB, "interface_is should be NULL or empty for PointSF being NULL");
1471c4762a1bSJed Brown     PetscFunctionReturn(0);
1472c4762a1bSJed Brown   }
1473c4762a1bSJed Brown 
1474c4762a1bSJed Brown   /* get PointSF points as IS pointsf_is */
1475c4762a1bSJed Brown   ierr = PointSFGetInterfacePoints_Private(pointsf, &pointsf_is);CHKERRQ(ierr);
1476c4762a1bSJed Brown 
1477c4762a1bSJed Brown   /* compare pointsf_is with interface_is */
1478c4762a1bSJed Brown   ierr = ISEqual(interface_is, pointsf_is, &flg);CHKERRQ(ierr);
1479ffc4695bSBarry Smith   ierr = MPI_Allreduce(MPI_IN_PLACE,&flg,1,MPIU_BOOL,MPI_LAND,comm);CHKERRMPI(ierr);
1480c4762a1bSJed Brown   if (!flg) {
1481c4762a1bSJed Brown     IS pointsf_extra_is, pointsf_missing_is;
1482c4762a1bSJed Brown     PetscViewer errv = PETSC_VIEWER_STDERR_(comm);
1483c4762a1bSJed Brown     ierr = ISDifference(interface_is, pointsf_is, &pointsf_missing_is);MYCHKERRQ(ierr);
1484c4762a1bSJed Brown     ierr = ISDifference(pointsf_is, interface_is, &pointsf_extra_is);MYCHKERRQ(ierr);
1485c4762a1bSJed Brown     ierr = PetscViewerASCIIPrintf(errv, "Points missing in PointSF:\n");MYCHKERRQ(ierr);
1486c4762a1bSJed Brown     ierr = ViewPointsWithType(ipdm, pointsf_missing_is, errv);MYCHKERRQ(ierr);
1487c4762a1bSJed Brown     ierr = PetscViewerASCIIPrintf(errv, "Extra points in PointSF:\n");MYCHKERRQ(ierr);
1488c4762a1bSJed Brown     ierr = ViewPointsWithType(ipdm, pointsf_extra_is, errv);MYCHKERRQ(ierr);
1489c4762a1bSJed Brown     ierr = ISDestroy(&pointsf_extra_is);MYCHKERRQ(ierr);
1490c4762a1bSJed Brown     ierr = ISDestroy(&pointsf_missing_is);MYCHKERRQ(ierr);
1491c4762a1bSJed Brown     SETERRQ(comm, PETSC_ERR_PLIB, "PointSF is wrong! See details above.");
1492c4762a1bSJed Brown   }
1493c4762a1bSJed Brown   ierr = ISDestroy(&pointsf_is);CHKERRQ(ierr);
1494c4762a1bSJed Brown   PetscFunctionReturn(0);
1495c4762a1bSJed Brown }
1496c4762a1bSJed Brown 
1497c4762a1bSJed Brown /* remove faces & edges from label, leave just vertices */
1498c4762a1bSJed Brown static PetscErrorCode DMPlexISFilterVertices_Private(DM dm, IS points)
1499c4762a1bSJed Brown {
1500c4762a1bSJed Brown   PetscInt        vStart, vEnd;
1501c4762a1bSJed Brown   MPI_Comm        comm;
1502c4762a1bSJed Brown   PetscErrorCode  ierr;
1503c4762a1bSJed Brown 
1504c4762a1bSJed Brown   PetscFunctionBegin;
1505c4762a1bSJed Brown   ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
1506c4762a1bSJed Brown   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
1507c4762a1bSJed Brown   ierr = ISGeneralFilter(points, vStart, vEnd);CHKERRQ(ierr);
1508c4762a1bSJed Brown   PetscFunctionReturn(0);
1509c4762a1bSJed Brown }
1510c4762a1bSJed Brown 
1511c4762a1bSJed Brown /*
1512c4762a1bSJed Brown   DMPlexCheckPointSFHeavy - Thorough test that the PointSF after parallel DMPlexInterpolate() includes exactly all interface points.
1513c4762a1bSJed Brown 
1514c4762a1bSJed Brown   Collective
1515c4762a1bSJed Brown 
1516c4762a1bSJed Brown   Input Parameters:
1517c4762a1bSJed Brown . dm - The DMPlex object
1518c4762a1bSJed Brown 
1519c4762a1bSJed Brown   Notes:
1520c4762a1bSJed Brown   The input DMPlex must be serial (one partition has all points, the other partitions have no points).
1521c4762a1bSJed Brown   This is a heavy test which involves DMPlexInterpolate() if the input DM is not interpolated yet, and depends on having a representation of the whole-domain boundary (PortableBoundary), which can be obtained only by DMPlexGetExpandedBoundary_Private() (which involves DMPlexInterpolate() of a sequential DM).
1522c4762a1bSJed Brown   This is mainly intended for debugging/testing purposes.
1523c4762a1bSJed Brown 
1524c4762a1bSJed Brown   Algorithm:
1525c4762a1bSJed Brown   1. boundary faces of the serial version of the whole mesh are found using DMPlexMarkBoundaryFaces()
1526c4762a1bSJed Brown   2. boundary faces are translated into vertices using DMPlexGetConeRecursive() and these are translated into coordinates - this description (aka PortableBoundary) is completely independent of partitioning and point numbering
1527c4762a1bSJed Brown   3. the mesh is distributed or loaded in parallel
1528c4762a1bSJed Brown   4. boundary faces of the distributed mesh are reconstructed from PortableBoundary using DMPlexFindVertices()
1529c4762a1bSJed Brown   5. partition boundary faces of the parallel mesh are found using DMPlexMarkBoundaryFaces()
1530c4762a1bSJed Brown   6. partition interfaces are computed as set difference of partition boundary faces minus the reconstructed boundary
1531c4762a1bSJed Brown   7. check that interface covered by PointSF (union of inward and outward points) is equal to the partition interface for each rank, otherwise print the difference and throw an error
1532c4762a1bSJed Brown 
1533c4762a1bSJed Brown   Level: developer
1534c4762a1bSJed Brown 
1535c4762a1bSJed Brown .seealso: DMGetPointSF(), DMPlexCheckSymmetry(), DMPlexCheckSkeleton(), DMPlexCheckFaces()
1536c4762a1bSJed Brown */
1537c4762a1bSJed Brown static PetscErrorCode DMPlexCheckPointSFHeavy(DM dm, PortableBoundary bnd)
1538c4762a1bSJed Brown {
1539c4762a1bSJed Brown   DM                     ipdm=NULL;
1540c4762a1bSJed Brown   IS                     boundary_faces_is, interface_faces_is, interface_is;
1541c4762a1bSJed Brown   DMPlexInterpolatedFlag intp;
1542c4762a1bSJed Brown   MPI_Comm               comm;
1543c4762a1bSJed Brown   PetscErrorCode         ierr;
1544c4762a1bSJed Brown 
1545c4762a1bSJed Brown   PetscFunctionBegin;
1546c4762a1bSJed Brown   ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
1547c4762a1bSJed Brown 
1548c4762a1bSJed Brown   ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr);
1549c4762a1bSJed Brown   if (intp == DMPLEX_INTERPOLATED_FULL) {
1550c4762a1bSJed Brown     ipdm = dm;
1551c4762a1bSJed Brown   } else {
1552c4762a1bSJed Brown     /* create temporary interpolated DM if input DM is not interpolated */
1553c4762a1bSJed Brown     ierr = DMPlexSetOrientInterface_Private(dm, PETSC_FALSE);CHKERRQ(ierr);
1554c4762a1bSJed Brown     ierr = DMPlexInterpolate(dm, &ipdm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexComparePointSFWithInterface_Private() below */
1555c4762a1bSJed Brown     ierr = DMPlexSetOrientInterface_Private(dm, PETSC_TRUE);CHKERRQ(ierr);
1556c4762a1bSJed Brown   }
1557c4762a1bSJed Brown   ierr = DMViewFromOptions(ipdm, NULL, "-ipdm_view");CHKERRQ(ierr);
1558c4762a1bSJed Brown 
1559c4762a1bSJed Brown   /* recover ipdm whole-domain boundary faces from the expanded vertices coordinates */
1560c4762a1bSJed Brown   ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(ipdm, bnd, &boundary_faces_is);CHKERRQ(ierr);
1561c4762a1bSJed Brown   /* get inter-partition interface faces (interface_faces_is)*/
1562c4762a1bSJed Brown   ierr = DMPlexGetInterfaceFaces_Private(ipdm, boundary_faces_is, &interface_faces_is);CHKERRQ(ierr);
1563c4762a1bSJed Brown   /* compute inter-partition interface including edges and vertices (interface_is) */
1564c4762a1bSJed Brown   ierr = DMPlexComputeCompleteInterface_Private(ipdm, interface_faces_is, &interface_is);CHKERRQ(ierr);
1565c4762a1bSJed Brown   /* destroy immediate ISs */
1566c4762a1bSJed Brown   ierr = ISDestroy(&boundary_faces_is);CHKERRQ(ierr);
1567c4762a1bSJed Brown   ierr = ISDestroy(&interface_faces_is);CHKERRQ(ierr);
1568c4762a1bSJed Brown 
1569c4762a1bSJed Brown   /* for uninterpolated case, keep just vertices in interface */
1570c4762a1bSJed Brown   if (!intp) {
1571c4762a1bSJed Brown     ierr = DMPlexISFilterVertices_Private(ipdm, interface_is);CHKERRQ(ierr);
1572c4762a1bSJed Brown     ierr = DMDestroy(&ipdm);CHKERRQ(ierr);
1573c4762a1bSJed Brown   }
1574c4762a1bSJed Brown 
1575c4762a1bSJed Brown   /* compare PointSF with the boundary reconstructed from coordinates */
1576c4762a1bSJed Brown   ierr = DMPlexComparePointSFWithInterface_Private(dm, interface_is);CHKERRQ(ierr);
1577c4762a1bSJed Brown   ierr = PetscPrintf(comm, "DMPlexCheckPointSFHeavy PASSED\n");CHKERRQ(ierr);
1578c4762a1bSJed Brown   ierr = ISDestroy(&interface_is);CHKERRQ(ierr);
1579c4762a1bSJed Brown   PetscFunctionReturn(0);
1580c4762a1bSJed Brown }
1581c4762a1bSJed Brown 
1582c4762a1bSJed Brown int main(int argc, char **argv)
1583c4762a1bSJed Brown {
1584c4762a1bSJed Brown   DM             dm;
1585c4762a1bSJed Brown   AppCtx         user;
1586c4762a1bSJed Brown   PetscErrorCode ierr;
1587c4762a1bSJed Brown 
1588c4762a1bSJed Brown   ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
1589c4762a1bSJed Brown   ierr = PetscLogStageRegister("create",&stage[0]);CHKERRQ(ierr);
1590c4762a1bSJed Brown   ierr = PetscLogStageRegister("distribute",&stage[1]);CHKERRQ(ierr);
1591c4762a1bSJed Brown   ierr = PetscLogStageRegister("interpolate",&stage[2]);CHKERRQ(ierr);
1592c4762a1bSJed Brown   ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
1593c4762a1bSJed Brown   ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
1594c4762a1bSJed Brown   if (user.nPointsToExpand) {
1595c4762a1bSJed Brown     ierr = TestExpandPoints(dm, &user);CHKERRQ(ierr);
1596c4762a1bSJed Brown   }
1597c4762a1bSJed Brown   if (user.ncoords) {
1598c4762a1bSJed Brown     ierr = ViewVerticesFromCoords(dm, user.ncoords/user.dim, user.coords, user.coordsTol, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1599c4762a1bSJed Brown   }
1600c4762a1bSJed Brown   ierr = DMDestroy(&dm);CHKERRQ(ierr);
1601c4762a1bSJed Brown   ierr = PetscFinalize();
1602c4762a1bSJed Brown   return ierr;
1603c4762a1bSJed Brown }
1604c4762a1bSJed Brown 
1605c4762a1bSJed Brown /*TEST
1606c4762a1bSJed Brown 
1607c4762a1bSJed Brown   testset:
1608c4762a1bSJed Brown     nsize: 2
1609c4762a1bSJed Brown     args: -dm_view ascii::ascii_info_detail
1610c4762a1bSJed Brown     args: -dm_plex_check_all
1611c4762a1bSJed Brown     test:
1612c4762a1bSJed Brown       suffix: 1_tri_dist0
1613c4762a1bSJed Brown       args: -distribute 0 -interpolate {{none create}separate output}
1614c4762a1bSJed Brown     test:
1615c4762a1bSJed Brown       suffix: 1_tri_dist1
1616c4762a1bSJed Brown       args: -distribute 1 -interpolate {{none create after_distribute}separate output}
1617c4762a1bSJed Brown     test:
1618c4762a1bSJed Brown       suffix: 1_quad_dist0
1619c4762a1bSJed Brown       args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output}
1620c4762a1bSJed Brown     test:
1621c4762a1bSJed Brown       suffix: 1_quad_dist1
1622c4762a1bSJed Brown       args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output}
1623c4762a1bSJed Brown     test:
1624c4762a1bSJed Brown       suffix: 1_1d_dist1
1625c4762a1bSJed Brown       args: -dim 1 -distribute 1
1626c4762a1bSJed Brown 
1627c4762a1bSJed Brown   testset:
1628c4762a1bSJed Brown     nsize: 3
1629c4762a1bSJed Brown     args: -testnum 1 -interpolate create
1630c4762a1bSJed Brown     args: -dm_plex_check_all
1631c4762a1bSJed Brown     test:
1632c4762a1bSJed Brown       suffix: 2
1633c4762a1bSJed Brown       args: -dm_view ascii::ascii_info_detail
1634c4762a1bSJed Brown     test:
1635c4762a1bSJed Brown       suffix: 2a
1636c4762a1bSJed Brown       args: -dm_plex_check_cones_conform_on_interfaces_verbose
1637c4762a1bSJed Brown     test:
1638c4762a1bSJed Brown       suffix: 2b
1639c4762a1bSJed Brown       args: -test_expand_points 0,1,2,5,6
1640c4762a1bSJed Brown     test:
1641c4762a1bSJed Brown       suffix: 2c
1642c4762a1bSJed Brown       args: -test_expand_points 0,1,2,5,6 -test_expand_points_empty
1643c4762a1bSJed Brown 
1644c4762a1bSJed Brown   testset:
1645c4762a1bSJed Brown     # the same as 1% for 3D
1646c4762a1bSJed Brown     nsize: 2
1647c4762a1bSJed Brown     args: -dim 3 -dm_view ascii::ascii_info_detail
1648c4762a1bSJed Brown     args: -dm_plex_check_all
1649c4762a1bSJed Brown     test:
1650c4762a1bSJed Brown       suffix: 4_tet_dist0
1651c4762a1bSJed Brown       args: -distribute 0 -interpolate {{none create}separate output}
1652c4762a1bSJed Brown     test:
1653c4762a1bSJed Brown       suffix: 4_tet_dist1
1654c4762a1bSJed Brown       args: -distribute 1 -interpolate {{none create after_distribute}separate output}
1655c4762a1bSJed Brown     test:
1656c4762a1bSJed Brown       suffix: 4_hex_dist0
1657c4762a1bSJed Brown       args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output}
1658c4762a1bSJed Brown     test:
1659c4762a1bSJed Brown       suffix: 4_hex_dist1
1660c4762a1bSJed Brown       args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output}
1661c4762a1bSJed Brown 
1662c4762a1bSJed Brown   test:
1663c4762a1bSJed Brown     # the same as 4_tet_dist0 but test different initial orientations
1664c4762a1bSJed Brown     suffix: 4_tet_test_orient
1665c4762a1bSJed Brown     nsize: 2
1666c4762a1bSJed Brown     args: -dim 3 -distribute 0
1667c4762a1bSJed Brown     args: -dm_plex_check_all
1668c4762a1bSJed Brown     args: -rotate_interface_0 {{0 1 2 11 12 13}}
1669c4762a1bSJed Brown     args: -rotate_interface_1 {{0 1 2 11 12 13}}
1670c4762a1bSJed Brown 
1671c4762a1bSJed Brown   testset:
1672c4762a1bSJed Brown     requires: exodusii
1673c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/TwoQuads.exo
1674c4762a1bSJed Brown     args: -dm_view ascii::ascii_info_detail
1675c4762a1bSJed Brown     args: -dm_plex_check_all
1676c4762a1bSJed Brown     args: -custom_view
1677c4762a1bSJed Brown     test:
1678c4762a1bSJed Brown       suffix: 5_seq
1679c4762a1bSJed Brown       nsize: 1
1680c4762a1bSJed Brown       args: -distribute 0 -interpolate {{none create}separate output}
1681c4762a1bSJed Brown     test:
168230602db0SMatthew G. Knepley       # Detail viewing in a non-distributed mesh is broken because the DMLabelView() is collective, but the label is not shared
1683c4762a1bSJed Brown       suffix: 5_dist0
1684c4762a1bSJed Brown       nsize: 2
168530602db0SMatthew G. Knepley       args: -distribute 0 -interpolate {{none create}separate output} -dm_view
1686c4762a1bSJed Brown     test:
1687c4762a1bSJed Brown       suffix: 5_dist1
1688c4762a1bSJed Brown       nsize: 2
1689c4762a1bSJed Brown       args: -distribute 1 -interpolate {{none create after_distribute}separate output}
1690c4762a1bSJed Brown 
1691c4762a1bSJed Brown   testset:
1692c4762a1bSJed Brown     nsize: {{1 2 4}}
1693c4762a1bSJed Brown     args: -use_generator
1694c4762a1bSJed Brown     args: -dm_plex_check_all
1695c4762a1bSJed Brown     args: -distribute -interpolate none
1696c4762a1bSJed Brown     test:
1697c4762a1bSJed Brown       suffix: 6_tri
1698c4762a1bSJed Brown       requires: triangle
1699c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 1 -dm_plex_generator triangle
1700c4762a1bSJed Brown     test:
1701c4762a1bSJed Brown       suffix: 6_quad
1702c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 0
1703c4762a1bSJed Brown     test:
1704c4762a1bSJed Brown       suffix: 6_tet
1705c4762a1bSJed Brown       requires: ctetgen
1706c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen
1707c4762a1bSJed Brown     test:
1708c4762a1bSJed Brown       suffix: 6_hex
1709c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 0
1710c4762a1bSJed Brown   testset:
1711c4762a1bSJed Brown     nsize: {{1 2 4}}
1712c4762a1bSJed Brown     args: -use_generator
1713c4762a1bSJed Brown     args: -dm_plex_check_all
1714c4762a1bSJed Brown     args: -distribute -interpolate create
1715c4762a1bSJed Brown     test:
1716c4762a1bSJed Brown       suffix: 6_int_tri
1717c4762a1bSJed Brown       requires: triangle
1718c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 1 -dm_plex_generator triangle
1719c4762a1bSJed Brown     test:
1720c4762a1bSJed Brown       suffix: 6_int_quad
1721c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 0
1722c4762a1bSJed Brown     test:
1723c4762a1bSJed Brown       suffix: 6_int_tet
1724c4762a1bSJed Brown       requires: ctetgen
1725c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen
1726c4762a1bSJed Brown     test:
1727c4762a1bSJed Brown       suffix: 6_int_hex
1728c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 0
1729c4762a1bSJed Brown   testset:
1730c4762a1bSJed Brown     nsize: {{2 4}}
1731c4762a1bSJed Brown     args: -use_generator
1732c4762a1bSJed Brown     args: -dm_plex_check_all
1733c4762a1bSJed Brown     args: -distribute -interpolate after_distribute
1734c4762a1bSJed Brown     test:
1735c4762a1bSJed Brown       suffix: 6_parint_tri
1736c4762a1bSJed Brown       requires: triangle
1737c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 1 -dm_plex_generator triangle
1738c4762a1bSJed Brown     test:
1739c4762a1bSJed Brown       suffix: 6_parint_quad
1740c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 0
1741c4762a1bSJed Brown     test:
1742c4762a1bSJed Brown       suffix: 6_parint_tet
1743c4762a1bSJed Brown       requires: ctetgen
1744c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen
1745c4762a1bSJed Brown     test:
1746c4762a1bSJed Brown       suffix: 6_parint_hex
1747c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 0
1748c4762a1bSJed Brown 
1749c4762a1bSJed Brown   testset: # 7 EXODUS
1750c4762a1bSJed Brown     requires: exodusii
1751c4762a1bSJed Brown     args: -dm_plex_check_all
1752c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo
1753c4762a1bSJed Brown     args: -distribute
1754c4762a1bSJed Brown     test: # seq load, simple partitioner
1755c4762a1bSJed Brown       suffix: 7_exo
1756c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1757c4762a1bSJed Brown       args: -interpolate none
1758c4762a1bSJed Brown     test: # seq load, seq interpolation, simple partitioner
1759c4762a1bSJed Brown       suffix: 7_exo_int_simple
1760c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1761c4762a1bSJed Brown       args: -interpolate create
1762c4762a1bSJed Brown     test: # seq load, seq interpolation, metis partitioner
1763c4762a1bSJed Brown       suffix: 7_exo_int_metis
1764c4762a1bSJed Brown       requires: parmetis
1765c4762a1bSJed Brown       nsize: {{2 4 5}}
1766c4762a1bSJed Brown       args: -interpolate create
1767c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1768c4762a1bSJed Brown     test: # seq load, simple partitioner, par interpolation
1769c4762a1bSJed Brown       suffix: 7_exo_simple_int
1770c4762a1bSJed Brown       nsize: {{2 4 5}}
1771c4762a1bSJed Brown       args: -interpolate after_distribute
1772c4762a1bSJed Brown     test: # seq load, metis partitioner, par interpolation
1773c4762a1bSJed Brown       suffix: 7_exo_metis_int
1774c4762a1bSJed Brown       requires: parmetis
1775c4762a1bSJed Brown       nsize: {{2 4 5}}
1776c4762a1bSJed Brown       args: -interpolate after_distribute
1777c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1778c4762a1bSJed Brown 
1779c4762a1bSJed Brown   testset: # 7 HDF5 SEQUANTIAL LOAD
1780c4762a1bSJed Brown     requires: hdf5 !complex
1781c4762a1bSJed Brown     args: -dm_plex_check_all
1782c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf
1783c4762a1bSJed Brown     args: -dm_plex_hdf5_force_sequential
1784c4762a1bSJed Brown     args: -distribute
1785c4762a1bSJed Brown     test: # seq load, simple partitioner
1786c4762a1bSJed Brown       suffix: 7_seq_hdf5_simple
1787c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1788c4762a1bSJed Brown       args: -interpolate none
1789c4762a1bSJed Brown     test: # seq load, seq interpolation, simple partitioner
1790c4762a1bSJed Brown       suffix: 7_seq_hdf5_int_simple
1791c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1792c4762a1bSJed Brown       args: -interpolate after_create
1793c4762a1bSJed Brown     test: # seq load, seq interpolation, metis partitioner
1794c4762a1bSJed Brown       nsize: {{2 4 5}}
1795c4762a1bSJed Brown       suffix: 7_seq_hdf5_int_metis
1796c4762a1bSJed Brown       requires: parmetis
1797c4762a1bSJed Brown       args: -interpolate after_create
1798c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1799c4762a1bSJed Brown     test: # seq load, simple partitioner, par interpolation
1800c4762a1bSJed Brown       suffix: 7_seq_hdf5_simple_int
1801c4762a1bSJed Brown       nsize: {{2 4 5}}
1802c4762a1bSJed Brown       args: -interpolate after_distribute
1803c4762a1bSJed Brown     test: # seq load, metis partitioner, par interpolation
1804c4762a1bSJed Brown       nsize: {{2 4 5}}
1805c4762a1bSJed Brown       suffix: 7_seq_hdf5_metis_int
1806c4762a1bSJed Brown       requires: parmetis
1807c4762a1bSJed Brown       args: -interpolate after_distribute
1808c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1809c4762a1bSJed Brown 
1810c4762a1bSJed Brown   testset: # 7 HDF5 PARALLEL LOAD
1811c4762a1bSJed Brown     requires: hdf5 !complex
1812c4762a1bSJed Brown     nsize: {{2 4 5}}
1813c4762a1bSJed Brown     args: -dm_plex_check_all
1814c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf
1815c4762a1bSJed Brown     test: # par load
1816c4762a1bSJed Brown       suffix: 7_par_hdf5
1817c4762a1bSJed Brown       args: -interpolate none
1818c4762a1bSJed Brown     test: # par load, par interpolation
1819c4762a1bSJed Brown       suffix: 7_par_hdf5_int
1820c4762a1bSJed Brown       args: -interpolate after_create
1821c4762a1bSJed Brown     test: # par load, parmetis repartitioner
1822c4762a1bSJed Brown       TODO: Parallel partitioning of uninterpolated meshes not supported
1823c4762a1bSJed Brown       suffix: 7_par_hdf5_parmetis
1824c4762a1bSJed Brown       requires: parmetis
1825c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1826c4762a1bSJed Brown       args: -interpolate none
1827c4762a1bSJed Brown     test: # par load, par interpolation, parmetis repartitioner
1828c4762a1bSJed Brown       suffix: 7_par_hdf5_int_parmetis
1829c4762a1bSJed Brown       requires: parmetis
1830c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1831c4762a1bSJed Brown       args: -interpolate after_create
1832c4762a1bSJed Brown     test: # par load, parmetis partitioner, par interpolation
1833c4762a1bSJed Brown       TODO: Parallel partitioning of uninterpolated meshes not supported
1834c4762a1bSJed Brown       suffix: 7_par_hdf5_parmetis_int
1835c4762a1bSJed Brown       requires: parmetis
1836c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1837c4762a1bSJed Brown       args: -interpolate after_distribute
1838c4762a1bSJed Brown 
1839c4762a1bSJed Brown     test:
1840c4762a1bSJed Brown       suffix: 7_hdf5_hierarch
1841c4762a1bSJed Brown       requires: hdf5 ptscotch !complex
1842c4762a1bSJed Brown       nsize: {{2 3 4}separate output}
1843c4762a1bSJed Brown       args: -distribute
1844c4762a1bSJed Brown       args: -interpolate after_create
1845c4762a1bSJed Brown       args: -petscpartitioner_type matpartitioning -petscpartitioner_view ::ascii_info
1846c4762a1bSJed Brown       args: -mat_partitioning_type hierarch -mat_partitioning_hierarchical_nfineparts 2
1847c4762a1bSJed Brown       args: -mat_partitioning_hierarchical_coarseparttype ptscotch -mat_partitioning_hierarchical_fineparttype ptscotch
1848c4762a1bSJed Brown 
1849c4762a1bSJed Brown   test:
1850c4762a1bSJed Brown     suffix: 8
1851c4762a1bSJed Brown     requires: hdf5 !complex
1852c4762a1bSJed Brown     nsize: 4
1853c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf
1854c4762a1bSJed Brown     args: -distribute 0 -interpolate after_create
1855c4762a1bSJed Brown     args: -view_vertices_from_coords 0.,1.,0.,-0.5,1.,0.,0.583,-0.644,0.,-2.,-2.,-2. -view_vertices_from_coords_tol 1e-3
1856c4762a1bSJed Brown     args: -dm_plex_check_all
1857c4762a1bSJed Brown     args: -custom_view
1858c4762a1bSJed Brown 
1859c4762a1bSJed Brown   testset: # 9 HDF5 SEQUANTIAL LOAD
1860c4762a1bSJed Brown     requires: hdf5 !complex datafilespath
1861c4762a1bSJed Brown     args: -dm_plex_check_all
1862c4762a1bSJed Brown     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
1863c4762a1bSJed Brown     args: -dm_plex_hdf5_force_sequential
1864c4762a1bSJed Brown     args: -distribute
1865c4762a1bSJed Brown     test: # seq load, simple partitioner
1866c4762a1bSJed Brown       suffix: 9_seq_hdf5_simple
1867c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1868c4762a1bSJed Brown       args: -interpolate none
1869c4762a1bSJed Brown     test: # seq load, seq interpolation, simple partitioner
1870c4762a1bSJed Brown       suffix: 9_seq_hdf5_int_simple
1871c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1872c4762a1bSJed Brown       args: -interpolate after_create
1873c4762a1bSJed Brown     test: # seq load, seq interpolation, metis partitioner
1874c4762a1bSJed Brown       nsize: {{2 4 5}}
1875c4762a1bSJed Brown       suffix: 9_seq_hdf5_int_metis
1876c4762a1bSJed Brown       requires: parmetis
1877c4762a1bSJed Brown       args: -interpolate after_create
1878c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1879c4762a1bSJed Brown     test: # seq load, simple partitioner, par interpolation
1880c4762a1bSJed Brown       suffix: 9_seq_hdf5_simple_int
1881c4762a1bSJed Brown       nsize: {{2 4 5}}
1882c4762a1bSJed Brown       args: -interpolate after_distribute
1883c4762a1bSJed Brown     test: # seq load, simple partitioner, par interpolation
1884c4762a1bSJed Brown       # This is like 9_seq_hdf5_simple_int but testing error output of DMPlexCheckPointSFHeavy().
1885c4762a1bSJed Brown       # Once 9_seq_hdf5_simple_int gets fixed, this one gets broken.
1886c4762a1bSJed Brown       # We can then provide an intentionally broken mesh instead.
1887c4762a1bSJed Brown       TODO: This test is broken because PointSF is fixed.
1888c4762a1bSJed Brown       suffix: 9_seq_hdf5_simple_int_err
1889c4762a1bSJed Brown       nsize: 4
1890c4762a1bSJed Brown       args: -interpolate after_distribute
1891c4762a1bSJed Brown       filter: sed -e "/PETSC ERROR/,$$d"
1892c4762a1bSJed Brown     test: # seq load, metis partitioner, par interpolation
1893c4762a1bSJed Brown       nsize: {{2 4 5}}
1894c4762a1bSJed Brown       suffix: 9_seq_hdf5_metis_int
1895c4762a1bSJed Brown       requires: parmetis
1896c4762a1bSJed Brown       args: -interpolate after_distribute
1897c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1898c4762a1bSJed Brown 
1899c4762a1bSJed Brown   testset: # 9 HDF5 PARALLEL LOAD
1900c4762a1bSJed Brown     requires: hdf5 !complex datafilespath
1901c4762a1bSJed Brown     nsize: {{2 4 5}}
1902c4762a1bSJed Brown     args: -dm_plex_check_all
1903c4762a1bSJed Brown     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
1904c4762a1bSJed Brown     test: # par load
1905c4762a1bSJed Brown       suffix: 9_par_hdf5
1906c4762a1bSJed Brown       args: -interpolate none
1907c4762a1bSJed Brown     test: # par load, par interpolation
1908c4762a1bSJed Brown       suffix: 9_par_hdf5_int
1909c4762a1bSJed Brown       args: -interpolate after_create
1910c4762a1bSJed Brown     test: # par load, parmetis repartitioner
1911c4762a1bSJed Brown       TODO: Parallel partitioning of uninterpolated meshes not supported
1912c4762a1bSJed Brown       suffix: 9_par_hdf5_parmetis
1913c4762a1bSJed Brown       requires: parmetis
1914c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1915c4762a1bSJed Brown       args: -interpolate none
1916c4762a1bSJed Brown     test: # par load, par interpolation, parmetis repartitioner
1917c4762a1bSJed Brown       suffix: 9_par_hdf5_int_parmetis
1918c4762a1bSJed Brown       requires: parmetis
1919c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1920c4762a1bSJed Brown       args: -interpolate after_create
1921c4762a1bSJed Brown     test: # par load, parmetis partitioner, par interpolation
1922c4762a1bSJed Brown       TODO: Parallel partitioning of uninterpolated meshes not supported
1923c4762a1bSJed Brown       suffix: 9_par_hdf5_parmetis_int
1924c4762a1bSJed Brown       requires: parmetis
1925c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1926c4762a1bSJed Brown       args: -interpolate after_distribute
1927c4762a1bSJed Brown 
1928c4762a1bSJed Brown   testset: # 10 HDF5 PARALLEL LOAD
1929c4762a1bSJed Brown     requires: hdf5 !complex datafilespath
1930c4762a1bSJed Brown     nsize: {{2 4 7}}
1931c4762a1bSJed Brown     args: -dm_plex_check_all
1932c4762a1bSJed Brown     args: -filename ${DATAFILESPATH}/meshes/cube-hexahedra-refined2.h5 -dm_plex_create_from_hdf5_xdmf -dm_plex_hdf5_topology_path /topo -dm_plex_hdf5_geometry_path /geom
1933c4762a1bSJed Brown     test: # par load, par interpolation
1934c4762a1bSJed Brown       suffix: 10_par_hdf5_int
1935c4762a1bSJed Brown       args: -interpolate after_create
1936c4762a1bSJed Brown TEST*/
1937