xref: /petsc/src/dm/impls/plex/tests/ex18.c (revision b122ec5aa1bd4469eb4e0673542fb7de3f411254)
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 */
204d3e1f4ccSVaclav Hapla   PetscScalar 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 
233c4762a1bSJed Brown   PetscFunctionBegin;
234c4762a1bSJed Brown   if (!*bnd) PetscFunctionReturn(0);
2355f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&(*bnd)->coordinates));
236c4762a1bSJed Brown   for (d=0; d < (*bnd)->depth; d++) {
2375f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscSectionDestroy(&(*bnd)->sections[d]));
238c4762a1bSJed Brown   }
2395f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscFree((*bnd)->sections));
2405f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscFree(*bnd));
241c4762a1bSJed Brown   PetscFunctionReturn(0);
242c4762a1bSJed Brown }
243c4762a1bSJed Brown 
244c4762a1bSJed Brown static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
245c4762a1bSJed Brown {
246c4762a1bSJed Brown   const char    *interpTypes[4]  = {"none", "create", "after_create", "after_distribute"};
247c4762a1bSJed Brown   PetscInt       interp=NONE, dim;
248c4762a1bSJed Brown   PetscBool      flg1, flg2;
249c4762a1bSJed Brown   PetscErrorCode ierr;
250c4762a1bSJed Brown 
251c4762a1bSJed Brown   PetscFunctionBegin;
252c4762a1bSJed Brown   options->debug        = 0;
253c4762a1bSJed Brown   options->testNum      = 0;
254c4762a1bSJed Brown   options->dim          = 2;
255c4762a1bSJed Brown   options->cellSimplex  = PETSC_TRUE;
256c4762a1bSJed Brown   options->distribute   = PETSC_FALSE;
257c4762a1bSJed Brown   options->interpolate  = NONE;
258c4762a1bSJed Brown   options->useGenerator = PETSC_FALSE;
259c4762a1bSJed Brown   options->testOrientIF = PETSC_FALSE;
260c4762a1bSJed Brown   options->testHeavy    = PETSC_TRUE;
261c4762a1bSJed Brown   options->customView   = PETSC_FALSE;
262c4762a1bSJed Brown   options->testExpandPointsEmpty = PETSC_FALSE;
263c4762a1bSJed Brown   options->ornt[0]      = 0;
264c4762a1bSJed Brown   options->ornt[1]      = 0;
265c4762a1bSJed Brown   options->faces[0]     = 2;
266c4762a1bSJed Brown   options->faces[1]     = 2;
267c4762a1bSJed Brown   options->faces[2]     = 2;
268c4762a1bSJed Brown   options->filename[0]  = '\0';
269c4762a1bSJed Brown   options->coordsTol    = PETSC_DEFAULT;
270c4762a1bSJed Brown 
271c4762a1bSJed Brown   ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr);
2725f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsBoundedInt("-debug", "The debugging level", "ex18.c", options->debug, &options->debug, NULL,0));
2735f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsBoundedInt("-testnum", "The mesh to create", "ex18.c", options->testNum, &options->testNum, NULL,0));
2745f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsBool("-cell_simplex", "Generate simplices if true, otherwise hexes", "ex18.c", options->cellSimplex, &options->cellSimplex, NULL));
2755f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsBool("-distribute", "Distribute the mesh", "ex18.c", options->distribute, &options->distribute, NULL));
2765f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsEList("-interpolate", "Type of mesh interpolation (none, create, after_create, after_distribute)", "ex18.c", interpTypes, 4, interpTypes[options->interpolate], &interp, NULL));
277c4762a1bSJed Brown   options->interpolate = (InterpType) interp;
2782c71b3e2SJacob Faibussowitsch   PetscCheckFalse(!options->distribute && options->interpolate == AFTER_DISTRIBUTE,comm, PETSC_ERR_SUP, "-interpolate after_distribute  needs  -distribute 1");
2795f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex18.c", options->useGenerator, &options->useGenerator, NULL));
280c4762a1bSJed Brown   options->ncoords = 128;
2815f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsScalarArray("-view_vertices_from_coords", "Print DAG points corresponding to vertices with given coordinates", "ex18.c", options->coords, &options->ncoords, NULL));
2825f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsReal("-view_vertices_from_coords_tol", "Tolerance for -view_vertices_from_coords", "ex18.c", options->coordsTol, &options->coordsTol, NULL));
283c4762a1bSJed Brown   options->nPointsToExpand = 128;
2845f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsIntArray("-test_expand_points", "Expand given array of DAG point using DMPlexGetConeRecursive() and print results", "ex18.c", options->pointsToExpand, &options->nPointsToExpand, NULL));
285c4762a1bSJed Brown   if (options->nPointsToExpand) {
2865f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscOptionsBool("-test_expand_points_empty", "For -test_expand_points, rank 0 will have empty input array", "ex18.c", options->testExpandPointsEmpty, &options->testExpandPointsEmpty, NULL));
287c4762a1bSJed Brown   }
2885f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsBool("-test_heavy", "Run the heavy PointSF test", "ex18.c", options->testHeavy, &options->testHeavy, NULL));
2895f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsBool("-custom_view", "Custom DMPlex view", "ex18.c", options->customView, &options->customView, NULL));
2905f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex18.c", options->dim, &options->dim, &flg1,1,3));
291c4762a1bSJed Brown   dim = 3;
2925f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsIntArray("-faces", "Number of faces per dimension", "ex18.c", options->faces, &dim, &flg2));
293c4762a1bSJed Brown   if (flg2) {
2942c71b3e2SJacob Faibussowitsch     PetscCheckFalse(flg1 && dim != options->dim,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);
295c4762a1bSJed Brown     options->dim = dim;
296c4762a1bSJed Brown   }
2975f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsString("-filename", "The mesh file", "ex18.c", options->filename, options->filename, sizeof(options->filename), NULL));
2985f80ce2aSJacob Faibussowitsch   CHKERRQ(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));
2995f80ce2aSJacob Faibussowitsch   CHKERRQ(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));
3002c71b3e2SJacob Faibussowitsch   PetscCheckFalse(flg2 != options->testOrientIF,comm, PETSC_ERR_ARG_OUTOFRANGE, "neither or both -rotate_interface_0 -rotate_interface_1 must be set");
301c4762a1bSJed Brown   if (options->testOrientIF) {
302c4762a1bSJed Brown     PetscInt i;
303c4762a1bSJed Brown     for (i=0; i<2; i++) {
304c4762a1bSJed Brown       if (options->ornt[i] >= 10) options->ornt[i] = -(options->ornt[i]-10);  /* 11 12 13 become -1 -2 -3 */
305c4762a1bSJed Brown     }
306c4762a1bSJed Brown     options->filename[0]  = 0;
307c4762a1bSJed Brown     options->useGenerator = PETSC_FALSE;
308c4762a1bSJed Brown     options->dim          = 3;
309c4762a1bSJed Brown     options->cellSimplex  = PETSC_TRUE;
310c4762a1bSJed Brown     options->interpolate  = CREATE;
311c4762a1bSJed Brown     options->distribute   = PETSC_FALSE;
312c4762a1bSJed Brown   }
3131e1ea65dSPierre Jolivet   ierr = PetscOptionsEnd();CHKERRQ(ierr);
314c4762a1bSJed Brown   PetscFunctionReturn(0);
315c4762a1bSJed Brown }
316c4762a1bSJed Brown 
317c4762a1bSJed Brown static PetscErrorCode CreateMesh_1D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
318c4762a1bSJed Brown {
319c4762a1bSJed Brown   PetscInt       testNum = user->testNum;
320c4762a1bSJed Brown   PetscMPIInt    rank,size;
321c4762a1bSJed Brown   PetscInt       numCorners=2,i;
322c4762a1bSJed Brown   PetscInt       numCells,numVertices,network;
323a4a685f2SJacob Faibussowitsch   PetscInt       *cells;
324c4762a1bSJed Brown   PetscReal      *coords;
325c4762a1bSJed Brown 
326c4762a1bSJed Brown   PetscFunctionBegin;
3275f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
3285f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(comm, &size));
3292c71b3e2SJacob Faibussowitsch   PetscCheckFalse(size > 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for <=2 processes",testNum);
330c4762a1bSJed Brown 
331c4762a1bSJed Brown   numCells = 3;
3325f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetInt(NULL, NULL, "-ncells", &numCells, NULL));
3332c71b3e2SJacob Faibussowitsch   PetscCheckFalse(numCells < 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test ncells must >=3",numCells);
334c4762a1bSJed Brown 
335c4762a1bSJed Brown   if (size == 1) {
336c4762a1bSJed Brown     numVertices = numCells + 1;
3375f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscMalloc2(2*numCells,&cells,2*numVertices,&coords));
338c4762a1bSJed Brown     for (i=0; i<numCells; i++) {
339c4762a1bSJed Brown       cells[2*i] = i; cells[2*i+1] = i + 1;
34071bbd31fSVaclav Hapla       coords[2*i] = i; coords[2*i+1] = i + 1;
341c4762a1bSJed Brown     }
342c4762a1bSJed Brown 
3435f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexCreateFromCellListPetsc(comm, user->dim, numCells, numVertices, numCorners, PETSC_FALSE, cells, user->dim, coords, dm));
3445f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscFree2(cells,coords));
345c4762a1bSJed Brown     PetscFunctionReturn(0);
346c4762a1bSJed Brown   }
347c4762a1bSJed Brown 
348c4762a1bSJed Brown   network = 0;
3495f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetInt(NULL, NULL, "-network_case", &network, NULL));
350c4762a1bSJed Brown   if (network == 0) {
351c4762a1bSJed Brown     switch (rank) {
352c4762a1bSJed Brown     case 0:
353c4762a1bSJed Brown     {
354c4762a1bSJed Brown       numCells    = 2;
355c4762a1bSJed Brown       numVertices = numCells;
3565f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscMalloc2(2*numCells,&cells,2*numCells,&coords));
357c4762a1bSJed Brown       cells[0] = 0; cells[1] = 1;
358c4762a1bSJed Brown       cells[2] = 1; cells[3] = 2;
359c4762a1bSJed Brown       coords[0] = 0.; coords[1] = 1.;
360c4762a1bSJed Brown       coords[2] = 1.; coords[3] = 2.;
361c4762a1bSJed Brown     }
362c4762a1bSJed Brown     break;
363c4762a1bSJed Brown     case 1:
364c4762a1bSJed Brown     {
365c4762a1bSJed Brown       numCells    -= 2;
366c4762a1bSJed Brown       numVertices = numCells + 1;
3675f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscMalloc2(2*numCells,&cells,2*numCells,&coords));
368c4762a1bSJed Brown       for (i=0; i<numCells; i++) {
369c4762a1bSJed Brown         cells[2*i] = 2+i; cells[2*i+1] = 2 + i + 1;
370c4762a1bSJed Brown         coords[2*i] = 2+i; coords[2*i+1] = 2 + i + 1;
371c4762a1bSJed Brown       }
372c4762a1bSJed Brown     }
373c4762a1bSJed Brown     break;
37498921bdaSJacob Faibussowitsch     default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
375c4762a1bSJed Brown     }
376c4762a1bSJed Brown   } else { /* network_case = 1 */
377c4762a1bSJed Brown     /* ----------------------- */
378c4762a1bSJed Brown     switch (rank) {
379c4762a1bSJed Brown     case 0:
380c4762a1bSJed Brown     {
381c4762a1bSJed Brown       numCells    = 2;
382c4762a1bSJed Brown       numVertices = 3;
3835f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscMalloc2(2*numCells,&cells,2*numCells,&coords));
384c4762a1bSJed Brown       cells[0] = 0; cells[1] = 3;
385c4762a1bSJed Brown       cells[2] = 3; cells[3] = 1;
386c4762a1bSJed Brown     }
387c4762a1bSJed Brown     break;
388c4762a1bSJed Brown     case 1:
389c4762a1bSJed Brown     {
390c4762a1bSJed Brown       numCells    = 1;
391c4762a1bSJed Brown       numVertices = 1;
3925f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscMalloc2(2*numCells,&cells,2*numCells,&coords));
393c4762a1bSJed Brown       cells[0] = 3; cells[1] = 2;
394c4762a1bSJed Brown     }
395c4762a1bSJed Brown     break;
39698921bdaSJacob Faibussowitsch     default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
397c4762a1bSJed Brown     }
398c4762a1bSJed Brown   }
3995f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, PETSC_FALSE, cells, user->dim, coords, NULL, NULL, dm));
4005f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscFree2(cells,coords));
401c4762a1bSJed Brown   PetscFunctionReturn(0);
402c4762a1bSJed Brown }
403c4762a1bSJed Brown 
404c4762a1bSJed Brown static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
405c4762a1bSJed Brown {
406c4762a1bSJed Brown   PetscInt       testNum = user->testNum, p;
407c4762a1bSJed Brown   PetscMPIInt    rank, size;
408c4762a1bSJed Brown 
409c4762a1bSJed Brown   PetscFunctionBegin;
4105f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
4115f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(comm, &size));
412c4762a1bSJed Brown   switch (testNum) {
413c4762a1bSJed Brown   case 0:
4142c71b3e2SJacob Faibussowitsch     PetscCheckFalse(size != 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum);
415c4762a1bSJed Brown     switch (rank) {
416c4762a1bSJed Brown       case 0:
417c4762a1bSJed Brown       {
418c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
419a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {0, 1, 2};
420c4762a1bSJed Brown         PetscReal      coords[4] = {-0.5, 0.5, 0.0, 0.0};
421c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
422c4762a1bSJed Brown 
4235f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
4245f80ce2aSJacob Faibussowitsch         for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
425c4762a1bSJed Brown       }
426c4762a1bSJed Brown       break;
427c4762a1bSJed Brown       case 1:
428c4762a1bSJed Brown       {
429c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
430a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {1, 3, 2};
431c4762a1bSJed Brown         PetscReal      coords[4] = {0.0, 1.0, 0.5, 0.5};
432c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
433c4762a1bSJed Brown 
4345f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
4355f80ce2aSJacob Faibussowitsch         for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
436c4762a1bSJed Brown       }
437c4762a1bSJed Brown       break;
43898921bdaSJacob Faibussowitsch       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
439c4762a1bSJed Brown     }
440c4762a1bSJed Brown     break;
441c4762a1bSJed Brown   case 1:
4422c71b3e2SJacob Faibussowitsch     PetscCheckFalse(size != 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum);
443c4762a1bSJed Brown     switch (rank) {
444c4762a1bSJed Brown       case 0:
445c4762a1bSJed Brown       {
446c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
447a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {0, 1, 2};
448c4762a1bSJed Brown         PetscReal      coords[4] = {0.0, 1.0, 0.0, 0.0};
449c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
450c4762a1bSJed Brown 
4515f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
4525f80ce2aSJacob Faibussowitsch         for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
453c4762a1bSJed Brown       }
454c4762a1bSJed Brown       break;
455c4762a1bSJed Brown       case 1:
456c4762a1bSJed Brown       {
457c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
458a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {0, 2, 3};
459c4762a1bSJed Brown         PetscReal      coords[4] = {0.5, 0.5, 1.0, 1.0};
460c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
461c4762a1bSJed Brown 
4625f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
4635f80ce2aSJacob Faibussowitsch         for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
464c4762a1bSJed Brown       }
465c4762a1bSJed Brown       break;
466c4762a1bSJed Brown       case 2:
467c4762a1bSJed Brown       {
468c4762a1bSJed Brown         const PetscInt numCells  = 2, numVertices = 1, numCorners = 3;
469a4a685f2SJacob Faibussowitsch         const PetscInt cells[6]  = {2, 4, 3, 2, 1, 4};
470c4762a1bSJed Brown         PetscReal      coords[2] = {1.0, 0.0};
471c4762a1bSJed Brown         PetscInt       markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
472c4762a1bSJed Brown 
4735f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
4745f80ce2aSJacob Faibussowitsch         for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
475c4762a1bSJed Brown       }
476c4762a1bSJed Brown       break;
47798921bdaSJacob Faibussowitsch       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
478c4762a1bSJed Brown     }
479c4762a1bSJed Brown     break;
480c4762a1bSJed Brown   case 2:
4812c71b3e2SJacob Faibussowitsch     PetscCheckFalse(size != 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum);
482c4762a1bSJed Brown     switch (rank) {
483c4762a1bSJed Brown       case 0:
484c4762a1bSJed Brown       {
485c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
486a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {1, 2, 0};
487c4762a1bSJed Brown         PetscReal      coords[4] = {0.5, 0.5, 0.0, 1.0};
488c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
489c4762a1bSJed Brown 
4905f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
4915f80ce2aSJacob Faibussowitsch         for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
492c4762a1bSJed Brown       }
493c4762a1bSJed Brown       break;
494c4762a1bSJed Brown       case 1:
495c4762a1bSJed Brown       {
496c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 3;
497a4a685f2SJacob Faibussowitsch         const PetscInt cells[3]  = {1, 0, 3};
498c4762a1bSJed Brown         PetscReal      coords[4] = {0.0, 0.0, 1.0, 1.0};
499c4762a1bSJed Brown         PetscInt       markerPoints[6] = {1, 1, 2, 1, 3, 1};
500c4762a1bSJed Brown 
5015f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
5025f80ce2aSJacob Faibussowitsch         for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
503c4762a1bSJed Brown       }
504c4762a1bSJed Brown       break;
505c4762a1bSJed Brown       case 2:
506c4762a1bSJed Brown       {
507c4762a1bSJed Brown         const PetscInt numCells  = 2, numVertices = 1, numCorners = 3;
508a4a685f2SJacob Faibussowitsch         const PetscInt cells[6]  = {0, 4, 3, 0, 2, 4};
509c4762a1bSJed Brown         PetscReal      coords[2] = {1.0, 0.0};
510c4762a1bSJed Brown         PetscInt       markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
511c4762a1bSJed Brown 
5125f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
5135f80ce2aSJacob Faibussowitsch         for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
514c4762a1bSJed Brown       }
515c4762a1bSJed Brown       break;
51698921bdaSJacob Faibussowitsch       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
517c4762a1bSJed Brown     }
518c4762a1bSJed Brown     break;
51998921bdaSJacob Faibussowitsch   default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum);
520c4762a1bSJed Brown   }
521c4762a1bSJed Brown   PetscFunctionReturn(0);
522c4762a1bSJed Brown }
523c4762a1bSJed Brown 
524c4762a1bSJed Brown static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
525c4762a1bSJed Brown {
526c4762a1bSJed Brown   PetscInt       testNum = user->testNum, p;
527c4762a1bSJed Brown   PetscMPIInt    rank, size;
528c4762a1bSJed Brown 
529c4762a1bSJed Brown   PetscFunctionBegin;
5305f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
5315f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(comm, &size));
532c4762a1bSJed Brown   switch (testNum) {
533c4762a1bSJed Brown   case 0:
5342c71b3e2SJacob Faibussowitsch     PetscCheckFalse(size != 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum);
535c4762a1bSJed Brown     switch (rank) {
536c4762a1bSJed Brown       case 0:
537c4762a1bSJed Brown       {
538c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 2, numCorners = 4;
539a4a685f2SJacob Faibussowitsch         const PetscInt cells[4]  = {0, 2, 1, 3};
540c4762a1bSJed Brown         PetscReal      coords[6] = {0.0, 0.0, -0.5,  0.0, -0.5, 0.0};
541c4762a1bSJed Brown         PetscInt       markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1};
542c4762a1bSJed Brown 
5435f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
5445f80ce2aSJacob Faibussowitsch         for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
545c4762a1bSJed Brown       }
546c4762a1bSJed Brown       break;
547c4762a1bSJed Brown       case 1:
548c4762a1bSJed Brown       {
549c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 3, numCorners = 4;
550a4a685f2SJacob Faibussowitsch         const PetscInt cells[4]  = {1, 2, 4, 3};
551c4762a1bSJed Brown         PetscReal      coords[9] = {1.0, 0.0, 0.0,  0.0, 0.5, 0.0,  0.0, 0.0, 0.5};
552c4762a1bSJed Brown         PetscInt       markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1};
553c4762a1bSJed Brown 
5545f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
5555f80ce2aSJacob Faibussowitsch         for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
556c4762a1bSJed Brown       }
557c4762a1bSJed Brown       break;
55898921bdaSJacob Faibussowitsch       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
559c4762a1bSJed Brown     }
560c4762a1bSJed Brown     break;
56198921bdaSJacob Faibussowitsch   default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum);
562c4762a1bSJed Brown   }
563c4762a1bSJed Brown   if (user->testOrientIF) {
564c4762a1bSJed Brown     PetscInt ifp[] = {8, 6};
565c4762a1bSJed Brown 
5665f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscObjectSetName((PetscObject) *dm, "Mesh before orientation"));
5675f80ce2aSJacob Faibussowitsch     CHKERRQ(DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view"));
568c4762a1bSJed Brown     /* rotate interface face ifp[rank] by given orientation ornt[rank] */
5695f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexOrientPoint(*dm, ifp[rank], user->ornt[rank]));
5705f80ce2aSJacob Faibussowitsch     CHKERRQ(DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view"));
5715f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexCheckFaces(*dm, 0));
5725f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexOrientInterface_Internal(*dm));
5735f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscPrintf(comm, "Orientation test PASSED\n"));
574c4762a1bSJed Brown   }
575c4762a1bSJed Brown   PetscFunctionReturn(0);
576c4762a1bSJed Brown }
577c4762a1bSJed Brown 
578c4762a1bSJed Brown static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
579c4762a1bSJed Brown {
580c4762a1bSJed Brown   PetscInt       testNum = user->testNum, p;
581c4762a1bSJed Brown   PetscMPIInt    rank, size;
582c4762a1bSJed Brown 
583c4762a1bSJed Brown   PetscFunctionBegin;
5845f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
5855f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(comm, &size));
586c4762a1bSJed Brown   switch (testNum) {
587c4762a1bSJed Brown   case 0:
5882c71b3e2SJacob Faibussowitsch     PetscCheckFalse(size != 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum);
589c4762a1bSJed Brown     switch (rank) {
590c4762a1bSJed Brown       case 0:
591c4762a1bSJed Brown       {
592c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 3, numCorners = 4;
593a4a685f2SJacob Faibussowitsch         const PetscInt cells[4]  = {0, 1, 2, 3};
594c4762a1bSJed Brown         PetscReal      coords[6] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0};
595c4762a1bSJed Brown         PetscInt       markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1};
596c4762a1bSJed Brown 
5975f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
5985f80ce2aSJacob Faibussowitsch         for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
599c4762a1bSJed Brown       }
600c4762a1bSJed Brown       break;
601c4762a1bSJed Brown       case 1:
602c4762a1bSJed Brown       {
603c4762a1bSJed Brown         const PetscInt numCells  = 1, numVertices = 3, numCorners = 4;
604a4a685f2SJacob Faibussowitsch         const PetscInt cells[4]  = {1, 4, 5, 2};
605c4762a1bSJed Brown         PetscReal      coords[6] = {-0.5, 1.0, 0.5, 0.0, 0.5, 1.0};
606c4762a1bSJed Brown         PetscInt       markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1};
607c4762a1bSJed Brown 
6085f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
6095f80ce2aSJacob Faibussowitsch         for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
610c4762a1bSJed Brown       }
611c4762a1bSJed Brown       break;
61298921bdaSJacob Faibussowitsch       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
613c4762a1bSJed Brown     }
614c4762a1bSJed Brown     break;
61598921bdaSJacob Faibussowitsch   default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum);
616c4762a1bSJed Brown   }
617c4762a1bSJed Brown   PetscFunctionReturn(0);
618c4762a1bSJed Brown }
619c4762a1bSJed Brown 
620c4762a1bSJed Brown static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm)
621c4762a1bSJed Brown {
622c4762a1bSJed Brown   PetscInt       testNum = user->testNum, p;
623c4762a1bSJed Brown   PetscMPIInt    rank, size;
624c4762a1bSJed Brown 
625c4762a1bSJed Brown   PetscFunctionBegin;
6265f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
6275f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(comm, &size));
628c4762a1bSJed Brown   switch (testNum) {
629c4762a1bSJed Brown   case 0:
6302c71b3e2SJacob Faibussowitsch     PetscCheckFalse(size != 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum);
631c4762a1bSJed Brown     switch (rank) {
632c4762a1bSJed Brown     case 0:
633c4762a1bSJed Brown     {
634c4762a1bSJed Brown       const PetscInt numCells  = 1, numVertices = 6, numCorners = 8;
635a4a685f2SJacob Faibussowitsch       const PetscInt cells[8]  = {0, 3, 2, 1, 4, 5, 6, 7};
636c4762a1bSJed 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};
637c4762a1bSJed Brown       PetscInt       markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1};
638c4762a1bSJed Brown 
6395f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
6405f80ce2aSJacob Faibussowitsch       for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
641c4762a1bSJed Brown     }
642c4762a1bSJed Brown     break;
643c4762a1bSJed Brown     case 1:
644c4762a1bSJed Brown     {
645c4762a1bSJed Brown       const PetscInt numCells  = 1, numVertices = 6, numCorners = 8;
646a4a685f2SJacob Faibussowitsch       const PetscInt cells[8]  = {1, 2, 9, 8, 5, 10, 11, 6};
647c4762a1bSJed 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};
648c4762a1bSJed Brown       PetscInt       markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1};
649c4762a1bSJed Brown 
6505f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm));
6515f80ce2aSJacob Faibussowitsch       for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]));
652c4762a1bSJed Brown     }
653c4762a1bSJed Brown     break;
65498921bdaSJacob Faibussowitsch     default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank);
655c4762a1bSJed Brown     }
656c4762a1bSJed Brown   break;
65798921bdaSJacob Faibussowitsch   default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum);
658c4762a1bSJed Brown   }
659c4762a1bSJed Brown   PetscFunctionReturn(0);
660c4762a1bSJed Brown }
661c4762a1bSJed Brown 
662c4762a1bSJed Brown static PetscErrorCode CustomView(DM dm, PetscViewer v)
663c4762a1bSJed Brown {
664c4762a1bSJed Brown   DMPlexInterpolatedFlag interpolated;
665c4762a1bSJed Brown   PetscBool              distributed;
666c4762a1bSJed Brown 
667c4762a1bSJed Brown   PetscFunctionBegin;
6685f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexIsDistributed(dm, &distributed));
6695f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexIsInterpolatedCollective(dm, &interpolated));
6705f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPrintf(v, "DMPlexIsDistributed: %s\n", PetscBools[distributed]));
6715f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPrintf(v, "DMPlexIsInterpolatedCollective: %s\n", DMPlexInterpolatedFlags[interpolated]));
672c4762a1bSJed Brown   PetscFunctionReturn(0);
673c4762a1bSJed Brown }
674c4762a1bSJed Brown 
675c4762a1bSJed Brown static PetscErrorCode CreateMeshFromFile(MPI_Comm comm, AppCtx *user, DM *dm, DM *serialDM)
676c4762a1bSJed Brown {
677c4762a1bSJed Brown   const char    *filename       = user->filename;
678c4762a1bSJed Brown   PetscBool      testHeavy      = user->testHeavy;
679c4762a1bSJed Brown   PetscBool      interpCreate   = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE;
680c4762a1bSJed Brown   PetscBool      distributed    = PETSC_FALSE;
681c4762a1bSJed Brown 
682c4762a1bSJed Brown   PetscFunctionBegin;
683c4762a1bSJed Brown   *serialDM = NULL;
6845f80ce2aSJacob Faibussowitsch   if (testHeavy && interpCreate) CHKERRQ(DMPlexSetOrientInterface_Private(NULL, PETSC_FALSE));
6855f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscLogStagePush(stage[0]));
6865f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexCreateFromFile(comm, filename, "ex18_plex", interpCreate, dm)); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */
6875f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscLogStagePop());
6885f80ce2aSJacob Faibussowitsch   if (testHeavy && interpCreate) CHKERRQ(DMPlexSetOrientInterface_Private(NULL, PETSC_TRUE));
6895f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexIsDistributed(*dm, &distributed));
6905f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(comm, "DMPlexCreateFromFile produced %s mesh.\n", distributed ? "distributed" : "serial"));
691c4762a1bSJed Brown   if (testHeavy && distributed) {
6925f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscOptionsSetValue(NULL, "-dm_plex_hdf5_force_sequential", NULL));
6935f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexCreateFromFile(comm, filename, "ex18_plex", interpCreate, serialDM));
6945f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexIsDistributed(*serialDM, &distributed));
69528b400f6SJacob Faibussowitsch     PetscCheck(!distributed,comm, PETSC_ERR_PLIB, "unable to create a serial DM from file");
696c4762a1bSJed Brown   }
6975f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetDimension(*dm, &user->dim));
698c4762a1bSJed Brown   PetscFunctionReturn(0);
699c4762a1bSJed Brown }
700c4762a1bSJed Brown 
701c4762a1bSJed Brown static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
702c4762a1bSJed Brown {
703c4762a1bSJed Brown   PetscPartitioner part;
704c4762a1bSJed Brown   PortableBoundary boundary     = NULL;
705c4762a1bSJed Brown   DM             serialDM       = NULL;
706c4762a1bSJed Brown   PetscBool      cellSimplex    = user->cellSimplex;
707c4762a1bSJed Brown   PetscBool      useGenerator   = user->useGenerator;
708c4762a1bSJed Brown   PetscBool      interpCreate   = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE;
709c4762a1bSJed Brown   PetscBool      interpSerial   = user->interpolate == AFTER_CREATE ? PETSC_TRUE : PETSC_FALSE;
710c4762a1bSJed Brown   PetscBool      interpParallel = user->interpolate == AFTER_DISTRIBUTE ? PETSC_TRUE : PETSC_FALSE;
711c4762a1bSJed Brown   PetscBool      testHeavy      = user->testHeavy;
712c4762a1bSJed Brown   PetscMPIInt    rank;
713c4762a1bSJed Brown 
714c4762a1bSJed Brown   PetscFunctionBegin;
7155f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
716c4762a1bSJed Brown   if (user->filename[0]) {
7175f80ce2aSJacob Faibussowitsch     CHKERRQ(CreateMeshFromFile(comm, user, dm, &serialDM));
718c4762a1bSJed Brown   } else if (useGenerator) {
7195f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscLogStagePush(stage[0]));
7205f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexCreateBoxMesh(comm, user->dim, cellSimplex, user->faces, NULL, NULL, NULL, interpCreate, dm));
7215f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscLogStagePop());
722c4762a1bSJed Brown   } else {
7235f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscLogStagePush(stage[0]));
724c4762a1bSJed Brown     switch (user->dim) {
725c4762a1bSJed Brown     case 1:
7265f80ce2aSJacob Faibussowitsch       CHKERRQ(CreateMesh_1D(comm, interpCreate, user, dm));
727c4762a1bSJed Brown       break;
728c4762a1bSJed Brown     case 2:
729c4762a1bSJed Brown       if (cellSimplex) {
7305f80ce2aSJacob Faibussowitsch         CHKERRQ(CreateSimplex_2D(comm, interpCreate, user, dm));
731c4762a1bSJed Brown       } else {
7325f80ce2aSJacob Faibussowitsch         CHKERRQ(CreateQuad_2D(comm, interpCreate, user, dm));
733c4762a1bSJed Brown       }
734c4762a1bSJed Brown       break;
735c4762a1bSJed Brown     case 3:
736c4762a1bSJed Brown       if (cellSimplex) {
7375f80ce2aSJacob Faibussowitsch         CHKERRQ(CreateSimplex_3D(comm, interpCreate, user, dm));
738c4762a1bSJed Brown       } else {
7395f80ce2aSJacob Faibussowitsch         CHKERRQ(CreateHex_3D(comm, interpCreate, user, dm));
740c4762a1bSJed Brown       }
741c4762a1bSJed Brown       break;
742c4762a1bSJed Brown     default:
74398921bdaSJacob Faibussowitsch       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", user->dim);
744c4762a1bSJed Brown     }
7455f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscLogStagePop());
746c4762a1bSJed Brown   }
7472c71b3e2SJacob Faibussowitsch   PetscCheckFalse(user->ncoords % user->dim,comm, PETSC_ERR_ARG_OUTOFRANGE, "length of coordinates array %D must be divisable by spatial dimension %D", user->ncoords, user->dim);
7485f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectSetName((PetscObject) *dm, "Original Mesh"));
7495f80ce2aSJacob Faibussowitsch   CHKERRQ(DMViewFromOptions(*dm, NULL, "-orig_dm_view"));
750c4762a1bSJed Brown 
751c4762a1bSJed Brown   if (interpSerial) {
752c4762a1bSJed Brown     DM idm;
753c4762a1bSJed Brown 
7545f80ce2aSJacob Faibussowitsch     if (testHeavy) CHKERRQ(DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE));
7555f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscLogStagePush(stage[2]));
7565f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexInterpolate(*dm, &idm)); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */
7575f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscLogStagePop());
7585f80ce2aSJacob Faibussowitsch     if (testHeavy) CHKERRQ(DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE));
7595f80ce2aSJacob Faibussowitsch     CHKERRQ(DMDestroy(dm));
760c4762a1bSJed Brown     *dm = idm;
7615f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscObjectSetName((PetscObject) *dm, "Interpolated Mesh"));
7625f80ce2aSJacob Faibussowitsch     CHKERRQ(DMViewFromOptions(*dm, NULL, "-intp_dm_view"));
763c4762a1bSJed Brown   }
764c4762a1bSJed Brown 
765c4762a1bSJed Brown   /* Set partitioner options */
7665f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexGetPartitioner(*dm, &part));
767c4762a1bSJed Brown   if (part) {
7685f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE));
7695f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscPartitionerSetFromOptions(part));
770c4762a1bSJed Brown   }
771c4762a1bSJed Brown 
7725f80ce2aSJacob Faibussowitsch   if (user->customView) CHKERRQ(CustomView(*dm, PETSC_VIEWER_STDOUT_(comm)));
773c4762a1bSJed Brown   if (testHeavy) {
774c4762a1bSJed Brown     PetscBool distributed;
775c4762a1bSJed Brown 
7765f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexIsDistributed(*dm, &distributed));
777c4762a1bSJed Brown     if (!serialDM && !distributed) {
778c4762a1bSJed Brown       serialDM = *dm;
7795f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscObjectReference((PetscObject)*dm));
780c4762a1bSJed Brown     }
781c4762a1bSJed Brown     if (serialDM) {
7825f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexGetExpandedBoundary_Private(serialDM, &boundary));
783c4762a1bSJed Brown     }
784c4762a1bSJed Brown     if (boundary) {
785c4762a1bSJed Brown       /* check DM which has been created in parallel and already interpolated */
7865f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexCheckPointSFHeavy(*dm, boundary));
787c4762a1bSJed Brown     }
788c4762a1bSJed Brown     /* Orient interface because it could be deliberately skipped above. It is idempotent. */
7895f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexOrientInterface_Internal(*dm));
790c4762a1bSJed Brown   }
791c4762a1bSJed Brown   if (user->distribute) {
792c4762a1bSJed Brown     DM               pdm = NULL;
793c4762a1bSJed Brown 
794c4762a1bSJed Brown     /* Redistribute mesh over processes using that partitioner */
7955f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscLogStagePush(stage[1]));
7965f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexDistribute(*dm, 0, NULL, &pdm));
7975f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscLogStagePop());
798c4762a1bSJed Brown     if (pdm) {
7995f80ce2aSJacob Faibussowitsch       CHKERRQ(DMDestroy(dm));
800c4762a1bSJed Brown       *dm  = pdm;
8015f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscObjectSetName((PetscObject) *dm, "Redistributed Mesh"));
8025f80ce2aSJacob Faibussowitsch       CHKERRQ(DMViewFromOptions(*dm, NULL, "-dist_dm_view"));
803c4762a1bSJed Brown     }
804c4762a1bSJed Brown 
805c4762a1bSJed Brown     if (interpParallel) {
806c4762a1bSJed Brown       DM idm;
807c4762a1bSJed Brown 
8085f80ce2aSJacob Faibussowitsch       if (testHeavy) CHKERRQ(DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE));
8095f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscLogStagePush(stage[2]));
8105f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexInterpolate(*dm, &idm)); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */
8115f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscLogStagePop());
8125f80ce2aSJacob Faibussowitsch       if (testHeavy) CHKERRQ(DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE));
8135f80ce2aSJacob Faibussowitsch       CHKERRQ(DMDestroy(dm));
814c4762a1bSJed Brown       *dm = idm;
8155f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscObjectSetName((PetscObject) *dm, "Interpolated Redistributed Mesh"));
8165f80ce2aSJacob Faibussowitsch       CHKERRQ(DMViewFromOptions(*dm, NULL, "-intp_dm_view"));
817c4762a1bSJed Brown     }
818c4762a1bSJed Brown   }
819c4762a1bSJed Brown   if (testHeavy) {
820c4762a1bSJed Brown     if (boundary) {
8215f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexCheckPointSFHeavy(*dm, boundary));
822c4762a1bSJed Brown     }
823c4762a1bSJed Brown     /* Orient interface because it could be deliberately skipped above. It is idempotent. */
8245f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexOrientInterface_Internal(*dm));
825c4762a1bSJed Brown   }
826c4762a1bSJed Brown 
8275f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectSetName((PetscObject) *dm, "Parallel Mesh"));
8285f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
8295f80ce2aSJacob Faibussowitsch   CHKERRQ(DMSetFromOptions(*dm));
8305f80ce2aSJacob Faibussowitsch   CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view"));
831c4762a1bSJed Brown 
8325f80ce2aSJacob Faibussowitsch   if (user->customView) CHKERRQ(CustomView(*dm, PETSC_VIEWER_STDOUT_(comm)));
8335f80ce2aSJacob Faibussowitsch   CHKERRQ(DMDestroy(&serialDM));
8345f80ce2aSJacob Faibussowitsch   CHKERRQ(PortableBoundaryDestroy(&boundary));
835c4762a1bSJed Brown   PetscFunctionReturn(0);
836c4762a1bSJed Brown }
837c4762a1bSJed Brown 
838d3e1f4ccSVaclav Hapla #define ps2d(number) ((double) PetscRealPart(number))
8399fbee547SJacob Faibussowitsch static inline PetscErrorCode coord2str(char buf[], size_t len, PetscInt dim, const PetscScalar coords[], PetscReal tol)
840c4762a1bSJed Brown {
841c4762a1bSJed Brown   PetscFunctionBegin;
8422c71b3e2SJacob Faibussowitsch   PetscCheckFalse(dim > 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "dim must be less than or equal 3");
843c4762a1bSJed Brown   if (tol >= 1e-3) {
844c4762a1bSJed Brown     switch (dim) {
8455f80ce2aSJacob Faibussowitsch       case 1: CHKERRQ(PetscSNPrintf(buf,len,"(%12.3f)",ps2d(coords[0])));
8465f80ce2aSJacob Faibussowitsch       case 2: CHKERRQ(PetscSNPrintf(buf,len,"(%12.3f, %12.3f)",ps2d(coords[0]),ps2d(coords[1])));
8475f80ce2aSJacob Faibussowitsch       case 3: CHKERRQ(PetscSNPrintf(buf,len,"(%12.3f, %12.3f, %12.3f)",ps2d(coords[0]),ps2d(coords[1]),ps2d(coords[2])));
848c4762a1bSJed Brown     }
849c4762a1bSJed Brown   } else {
850c4762a1bSJed Brown     switch (dim) {
8515f80ce2aSJacob Faibussowitsch       case 1: CHKERRQ(PetscSNPrintf(buf,len,"(%12.6f)",ps2d(coords[0])));
8525f80ce2aSJacob Faibussowitsch       case 2: CHKERRQ(PetscSNPrintf(buf,len,"(%12.6f, %12.6f)",ps2d(coords[0]),ps2d(coords[1])));
8535f80ce2aSJacob Faibussowitsch       case 3: CHKERRQ(PetscSNPrintf(buf,len,"(%12.6f, %12.6f, %12.6f)",ps2d(coords[0]),ps2d(coords[1]),ps2d(coords[2])));
854c4762a1bSJed Brown     }
855c4762a1bSJed Brown   }
856c4762a1bSJed Brown   PetscFunctionReturn(0);
857c4762a1bSJed Brown }
858c4762a1bSJed Brown 
859d3e1f4ccSVaclav Hapla static PetscErrorCode ViewVerticesFromCoords(DM dm, Vec coordsVec, PetscReal tol, PetscViewer viewer)
860c4762a1bSJed Brown {
861d3e1f4ccSVaclav Hapla   PetscInt       dim, i, npoints;
862d3e1f4ccSVaclav Hapla   IS             pointsIS;
863d3e1f4ccSVaclav Hapla   const PetscInt *points;
864d3e1f4ccSVaclav Hapla   const PetscScalar *coords;
865c4762a1bSJed Brown   char           coordstr[128];
866c4762a1bSJed Brown   MPI_Comm       comm;
867c4762a1bSJed Brown   PetscMPIInt    rank;
868c4762a1bSJed Brown 
869c4762a1bSJed Brown   PetscFunctionBegin;
8705f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm));
8715f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
8725f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetDimension(dm, &dim));
8735f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPushSynchronized(viewer));
8745f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexFindVertices(dm, coordsVec, tol, &pointsIS));
8755f80ce2aSJacob Faibussowitsch   CHKERRQ(ISGetIndices(pointsIS, &points));
8765f80ce2aSJacob Faibussowitsch   CHKERRQ(ISGetLocalSize(pointsIS, &npoints));
8775f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(coordsVec, &coords));
878c4762a1bSJed Brown   for (i=0; i < npoints; i++) {
8795f80ce2aSJacob Faibussowitsch     CHKERRQ(coord2str(coordstr, sizeof(coordstr), dim, &coords[i*dim], tol));
8805f80ce2aSJacob Faibussowitsch     if (rank == 0 && i) CHKERRQ(PetscViewerASCIISynchronizedPrintf(viewer, "-----\n"));
8815f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %s --> points[%D] = %D\n", rank, coordstr, i, points[i]));
8825f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscViewerFlush(viewer));
883c4762a1bSJed Brown   }
8845f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPopSynchronized(viewer));
8855f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(coordsVec, &coords));
8865f80ce2aSJacob Faibussowitsch   CHKERRQ(ISRestoreIndices(pointsIS, &points));
8875f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&pointsIS));
888c4762a1bSJed Brown   PetscFunctionReturn(0);
889c4762a1bSJed Brown }
890c4762a1bSJed Brown 
891c4762a1bSJed Brown static PetscErrorCode TestExpandPoints(DM dm, AppCtx *user)
892c4762a1bSJed Brown {
893c4762a1bSJed Brown   IS                is;
894c4762a1bSJed Brown   PetscSection      *sects;
895c4762a1bSJed Brown   IS                *iss;
896c4762a1bSJed Brown   PetscInt          d,depth;
897c4762a1bSJed Brown   PetscMPIInt       rank;
898c4762a1bSJed Brown   PetscViewer       viewer=PETSC_VIEWER_STDOUT_WORLD, sviewer;
899c4762a1bSJed Brown 
900c4762a1bSJed Brown   PetscFunctionBegin;
9015f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank));
902dd400576SPatrick Sanan   if (user->testExpandPointsEmpty && rank == 0) {
9035f80ce2aSJacob Faibussowitsch     CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &is));
904c4762a1bSJed Brown   } else {
9055f80ce2aSJacob Faibussowitsch     CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, user->nPointsToExpand, user->pointsToExpand, PETSC_USE_POINTER, &is));
906c4762a1bSJed Brown   }
9075f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexGetConeRecursive(dm, is, &depth, &iss, &sects));
9085f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer));
9095f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPrintf(sviewer, "[%d] ==========================\n",rank));
910c4762a1bSJed Brown   for (d=depth-1; d>=0; d--) {
911c4762a1bSJed Brown     IS          checkIS;
912c4762a1bSJed Brown     PetscBool   flg;
913c4762a1bSJed Brown 
9145f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscViewerASCIIPrintf(sviewer, "depth %D ---------------\n",d));
9155f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscSectionView(sects[d], sviewer));
9165f80ce2aSJacob Faibussowitsch     CHKERRQ(ISView(iss[d], sviewer));
917c4762a1bSJed Brown     /* check reverse operation */
918c4762a1bSJed Brown     if (d < depth-1) {
9195f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexExpandedConesToFaces_Private(dm, iss[d], sects[d], &checkIS));
9205f80ce2aSJacob Faibussowitsch       CHKERRQ(ISEqualUnsorted(checkIS, iss[d+1], &flg));
92128b400f6SJacob Faibussowitsch       PetscCheck(flg,PetscObjectComm((PetscObject) checkIS), PETSC_ERR_PLIB, "DMPlexExpandedConesToFaces_Private produced wrong IS");
9225f80ce2aSJacob Faibussowitsch       CHKERRQ(ISDestroy(&checkIS));
923c4762a1bSJed Brown     }
924c4762a1bSJed Brown   }
9255f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer));
9265f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerFlush(viewer));
9275f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexRestoreConeRecursive(dm, is, &depth, &iss, &sects));
9285f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&is));
929c4762a1bSJed Brown   PetscFunctionReturn(0);
930c4762a1bSJed Brown }
931c4762a1bSJed Brown 
932c4762a1bSJed Brown static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM dm, IS is, PetscSection section, IS *newis)
933c4762a1bSJed Brown {
934c4762a1bSJed Brown   PetscInt          n,n1,ncone,numCoveredPoints,o,p,q,start,end;
935c4762a1bSJed Brown   const PetscInt    *coveredPoints;
936c4762a1bSJed Brown   const PetscInt    *arr, *cone;
937c4762a1bSJed Brown   PetscInt          *newarr;
938c4762a1bSJed Brown 
939c4762a1bSJed Brown   PetscFunctionBegin;
9405f80ce2aSJacob Faibussowitsch   CHKERRQ(ISGetLocalSize(is, &n));
9415f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionGetStorageSize(section, &n1));
9425f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionGetChart(section, &start, &end));
9432c71b3e2SJacob Faibussowitsch   PetscCheckFalse(n != n1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "IS size = %D != %D = section storage size", n, n1);
9445f80ce2aSJacob Faibussowitsch   CHKERRQ(ISGetIndices(is, &arr));
9455f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscMalloc1(end-start, &newarr));
946c4762a1bSJed Brown   for (q=start; q<end; q++) {
9475f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscSectionGetDof(section, q, &ncone));
9485f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscSectionGetOffset(section, q, &o));
949c4762a1bSJed Brown     cone = &arr[o];
950c4762a1bSJed Brown     if (ncone == 1) {
951c4762a1bSJed Brown       numCoveredPoints = 1;
952c4762a1bSJed Brown       p = cone[0];
953c4762a1bSJed Brown     } else {
954c4762a1bSJed Brown       PetscInt i;
955c4762a1bSJed Brown       p = PETSC_MAX_INT;
956c4762a1bSJed Brown       for (i=0; i<ncone; i++) if (cone[i] < 0) {p = -1; break;}
957c4762a1bSJed Brown       if (p >= 0) {
9585f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexGetJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints));
9592c71b3e2SJacob Faibussowitsch         PetscCheckFalse(numCoveredPoints > 1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "more than one covered points for section point %D",q);
960c4762a1bSJed Brown         if (numCoveredPoints) p = coveredPoints[0];
961c4762a1bSJed Brown         else                  p = -2;
9625f80ce2aSJacob Faibussowitsch         CHKERRQ(DMPlexRestoreJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints));
963c4762a1bSJed Brown       }
964c4762a1bSJed Brown     }
965c4762a1bSJed Brown     newarr[q-start] = p;
966c4762a1bSJed Brown   }
9675f80ce2aSJacob Faibussowitsch   CHKERRQ(ISRestoreIndices(is, &arr));
9685f80ce2aSJacob Faibussowitsch   CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, end-start, newarr, PETSC_OWN_POINTER, newis));
969c4762a1bSJed Brown   PetscFunctionReturn(0);
970c4762a1bSJed Brown }
971c4762a1bSJed Brown 
972c4762a1bSJed Brown static PetscErrorCode DMPlexExpandedVerticesToFaces_Private(DM dm, IS boundary_expanded_is, PetscInt depth, PetscSection sections[], IS *boundary_is)
973c4762a1bSJed Brown {
974c4762a1bSJed Brown   PetscInt          d;
975c4762a1bSJed Brown   IS                is,newis;
976c4762a1bSJed Brown 
977c4762a1bSJed Brown   PetscFunctionBegin;
978c4762a1bSJed Brown   is = boundary_expanded_is;
9795f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectReference((PetscObject)is));
980c4762a1bSJed Brown   for (d = 0; d < depth-1; ++d) {
9815f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexExpandedConesToFaces_Private(dm, is, sections[d], &newis));
9825f80ce2aSJacob Faibussowitsch     CHKERRQ(ISDestroy(&is));
983c4762a1bSJed Brown     is = newis;
984c4762a1bSJed Brown   }
985c4762a1bSJed Brown   *boundary_is = is;
986c4762a1bSJed Brown   PetscFunctionReturn(0);
987c4762a1bSJed Brown }
988c4762a1bSJed Brown 
989c4762a1bSJed Brown #define CHKERRQI(incall,ierr) if (ierr) {incall = PETSC_FALSE; }
990c4762a1bSJed Brown 
991c4762a1bSJed Brown static PetscErrorCode DMLabelViewFromOptionsOnComm_Private(DMLabel label, const char optionname[], MPI_Comm comm)
992c4762a1bSJed Brown {
993c4762a1bSJed Brown   PetscViewer       viewer;
994c4762a1bSJed Brown   PetscBool         flg;
995c4762a1bSJed Brown   static PetscBool  incall = PETSC_FALSE;
996c4762a1bSJed Brown   PetscViewerFormat format;
997c4762a1bSJed Brown 
998c4762a1bSJed Brown   PetscFunctionBegin;
999c4762a1bSJed Brown   if (incall) PetscFunctionReturn(0);
1000c4762a1bSJed Brown   incall = PETSC_TRUE;
10015f80ce2aSJacob Faibussowitsch   CHKERRQI(incall,PetscOptionsGetViewer(comm,((PetscObject)label)->options,((PetscObject)label)->prefix,optionname,&viewer,&format,&flg));
1002c4762a1bSJed Brown   if (flg) {
10035f80ce2aSJacob Faibussowitsch     CHKERRQI(incall,PetscViewerPushFormat(viewer,format));
10045f80ce2aSJacob Faibussowitsch     CHKERRQI(incall,DMLabelView(label, viewer));
10055f80ce2aSJacob Faibussowitsch     CHKERRQI(incall,PetscViewerPopFormat(viewer));
10065f80ce2aSJacob Faibussowitsch     CHKERRQI(incall,PetscViewerDestroy(&viewer));
1007c4762a1bSJed Brown   }
1008c4762a1bSJed Brown   incall = PETSC_FALSE;
1009c4762a1bSJed Brown   PetscFunctionReturn(0);
1010c4762a1bSJed Brown }
1011c4762a1bSJed Brown 
1012c4762a1bSJed Brown /* TODO: this is hotfixing DMLabelGetStratumIS() - it should be fixed systematically instead */
10139fbee547SJacob Faibussowitsch static inline PetscErrorCode DMLabelGetStratumISOnComm_Private(DMLabel label, PetscInt value, MPI_Comm comm, IS *is)
1014c4762a1bSJed Brown {
1015c4762a1bSJed Brown   IS                tmpis;
1016c4762a1bSJed Brown 
1017c4762a1bSJed Brown   PetscFunctionBegin;
10185f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelGetStratumIS(label, value, &tmpis));
10195f80ce2aSJacob Faibussowitsch   if (!tmpis) CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &tmpis));
10205f80ce2aSJacob Faibussowitsch   CHKERRQ(ISOnComm(tmpis, comm, PETSC_COPY_VALUES, is));
10215f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&tmpis));
1022c4762a1bSJed Brown   PetscFunctionReturn(0);
1023c4762a1bSJed Brown }
1024c4762a1bSJed Brown 
1025c4762a1bSJed Brown /* currently only for simple PetscSection without fields or constraints */
1026c4762a1bSJed Brown static PetscErrorCode PetscSectionReplicate_Private(MPI_Comm comm, PetscMPIInt rootrank, PetscSection sec0, PetscSection *secout)
1027c4762a1bSJed Brown {
1028c4762a1bSJed Brown   PetscSection      sec;
1029c4762a1bSJed Brown   PetscInt          chart[2], p;
1030c4762a1bSJed Brown   PetscInt          *dofarr;
1031c4762a1bSJed Brown   PetscMPIInt       rank;
1032c4762a1bSJed Brown 
1033c4762a1bSJed Brown   PetscFunctionBegin;
10345f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
1035c4762a1bSJed Brown   if (rank == rootrank) {
10365f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscSectionGetChart(sec0, &chart[0], &chart[1]));
1037c4762a1bSJed Brown   }
10385f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Bcast(chart, 2, MPIU_INT, rootrank, comm));
10395f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscMalloc1(chart[1]-chart[0], &dofarr));
1040c4762a1bSJed Brown   if (rank == rootrank) {
1041c4762a1bSJed Brown     for (p = chart[0]; p < chart[1]; p++) {
10425f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscSectionGetDof(sec0, p, &dofarr[p-chart[0]]));
1043c4762a1bSJed Brown     }
1044c4762a1bSJed Brown   }
10455f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Bcast(dofarr, chart[1]-chart[0], MPIU_INT, rootrank, comm));
10465f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionCreate(comm, &sec));
10475f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionSetChart(sec, chart[0], chart[1]));
1048c4762a1bSJed Brown   for (p = chart[0]; p < chart[1]; p++) {
10495f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscSectionSetDof(sec, p, dofarr[p-chart[0]]));
1050c4762a1bSJed Brown   }
10515f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionSetUp(sec));
10525f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscFree(dofarr));
1053c4762a1bSJed Brown   *secout = sec;
1054c4762a1bSJed Brown   PetscFunctionReturn(0);
1055c4762a1bSJed Brown }
1056c4762a1bSJed Brown 
1057c4762a1bSJed Brown static PetscErrorCode DMPlexExpandedVerticesCoordinatesToFaces_Private(DM ipdm, PortableBoundary bnd, IS *face_is)
1058c4762a1bSJed Brown {
1059c4762a1bSJed Brown   IS                  faces_expanded_is;
1060c4762a1bSJed Brown 
1061c4762a1bSJed Brown   PetscFunctionBegin;
10625f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexFindVertices(ipdm, bnd->coordinates, 0.0, &faces_expanded_is));
10635f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexExpandedVerticesToFaces_Private(ipdm, faces_expanded_is, bnd->depth, bnd->sections, face_is));
10645f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&faces_expanded_is));
1065c4762a1bSJed Brown   PetscFunctionReturn(0);
1066c4762a1bSJed Brown }
1067c4762a1bSJed Brown 
1068c4762a1bSJed Brown /* hack disabling DMPlexOrientInterface() call in DMPlexInterpolate() via -dm_plex_interpolate_orient_interfaces option */
1069c4762a1bSJed Brown static PetscErrorCode DMPlexSetOrientInterface_Private(DM dm, PetscBool enable)
1070c4762a1bSJed Brown {
1071c4762a1bSJed Brown   PetscOptions      options = NULL;
1072c4762a1bSJed Brown   const char        *prefix = NULL;
1073c4762a1bSJed Brown   const char        opt[] = "-dm_plex_interpolate_orient_interfaces";
1074c4762a1bSJed Brown   char              prefix_opt[512];
1075c4762a1bSJed Brown   PetscBool         flg, set;
1076c4762a1bSJed Brown   static PetscBool  wasSetTrue = PETSC_FALSE;
1077c4762a1bSJed Brown 
1078c4762a1bSJed Brown   PetscFunctionBegin;
1079c4762a1bSJed Brown   if (dm) {
10805f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1081c4762a1bSJed Brown     options = ((PetscObject)dm)->options;
1082c4762a1bSJed Brown   }
10835f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscStrcpy(prefix_opt, "-"));
10845f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscStrlcat(prefix_opt, prefix, sizeof(prefix_opt)));
10855f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscStrlcat(prefix_opt, &opt[1], sizeof(prefix_opt)));
10865f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetBool(options, prefix, opt, &flg, &set));
1087c4762a1bSJed Brown   if (!enable) {
1088c4762a1bSJed Brown     if (set && flg) wasSetTrue = PETSC_TRUE;
10895f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscOptionsSetValue(options, prefix_opt, "0"));
1090c4762a1bSJed Brown   } else if (set && !flg) {
1091c4762a1bSJed Brown     if (wasSetTrue) {
10925f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscOptionsSetValue(options, prefix_opt, "1"));
1093c4762a1bSJed Brown     } else {
1094c4762a1bSJed Brown       /* default is PETSC_TRUE */
10955f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscOptionsClearValue(options, prefix_opt));
1096c4762a1bSJed Brown     }
1097c4762a1bSJed Brown     wasSetTrue = PETSC_FALSE;
1098c4762a1bSJed Brown   }
109976bd3646SJed Brown   if (PetscDefined(USE_DEBUG)) {
11005f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscOptionsGetBool(options, prefix, opt, &flg, &set));
11012c71b3e2SJacob Faibussowitsch     PetscCheckFalse(set && flg != enable,PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "PetscOptionsSetValue did not have the desired effect");
1102c4762a1bSJed Brown   }
1103c4762a1bSJed Brown   PetscFunctionReturn(0);
1104c4762a1bSJed Brown }
1105c4762a1bSJed Brown 
1106c4762a1bSJed Brown /* get coordinate description of the whole-domain boundary */
1107c4762a1bSJed Brown static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM dm, PortableBoundary *boundary)
1108c4762a1bSJed Brown {
1109c4762a1bSJed Brown   PortableBoundary       bnd0, bnd;
1110c4762a1bSJed Brown   MPI_Comm               comm;
1111c4762a1bSJed Brown   DM                     idm;
1112c4762a1bSJed Brown   DMLabel                label;
1113c4762a1bSJed Brown   PetscInt               d;
1114c4762a1bSJed Brown   const char             boundaryName[] = "DMPlexDistributeInterpolateMarkInterface_boundary";
1115c4762a1bSJed Brown   IS                     boundary_is;
1116c4762a1bSJed Brown   IS                     *boundary_expanded_iss;
1117c4762a1bSJed Brown   PetscMPIInt            rootrank = 0;
1118c4762a1bSJed Brown   PetscMPIInt            rank, size;
1119c4762a1bSJed Brown   PetscInt               value = 1;
1120c4762a1bSJed Brown   DMPlexInterpolatedFlag intp;
1121c4762a1bSJed Brown   PetscBool              flg;
1122c4762a1bSJed Brown 
1123c4762a1bSJed Brown   PetscFunctionBegin;
11245f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscNew(&bnd));
11255f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm));
11265f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm, &rank));
11275f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(comm, &size));
11285f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexIsDistributed(dm, &flg));
112928b400f6SJacob Faibussowitsch   PetscCheck(!flg,comm, PETSC_ERR_ARG_WRONG, "serial DM (all points on one rank) needed");
1130c4762a1bSJed Brown 
1131c4762a1bSJed Brown   /* interpolate serial DM if not yet interpolated */
11325f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexIsInterpolatedCollective(dm, &intp));
1133c4762a1bSJed Brown   if (intp == DMPLEX_INTERPOLATED_FULL) {
1134c4762a1bSJed Brown     idm = dm;
11355f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscObjectReference((PetscObject)dm));
1136c4762a1bSJed Brown   } else {
11375f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexInterpolate(dm, &idm));
11385f80ce2aSJacob Faibussowitsch     CHKERRQ(DMViewFromOptions(idm, NULL, "-idm_view"));
1139c4762a1bSJed Brown   }
1140c4762a1bSJed Brown 
1141c4762a1bSJed Brown   /* mark whole-domain boundary of the serial DM */
11425f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, boundaryName, &label));
11435f80ce2aSJacob Faibussowitsch   CHKERRQ(DMAddLabel(idm, label));
11445f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexMarkBoundaryFaces(idm, value, label));
11455f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelViewFromOptionsOnComm_Private(label, "-idm_boundary_view", comm));
11465f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelGetStratumIS(label, value, &boundary_is));
1147c4762a1bSJed Brown 
1148c4762a1bSJed Brown   /* translate to coordinates */
11495f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscNew(&bnd0));
11505f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetCoordinatesLocalSetUp(idm));
1151c4762a1bSJed Brown   if (rank == rootrank) {
11525f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexGetConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections));
11535f80ce2aSJacob Faibussowitsch     CHKERRQ(DMGetCoordinatesLocalTuple(dm, boundary_expanded_iss[0], NULL, &bnd0->coordinates));
1154c4762a1bSJed Brown     /* self-check */
1155c4762a1bSJed Brown     {
1156c4762a1bSJed Brown       IS is0;
11575f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexExpandedVerticesCoordinatesToFaces_Private(idm, bnd0, &is0));
11585f80ce2aSJacob Faibussowitsch       CHKERRQ(ISEqual(is0, boundary_is, &flg));
115928b400f6SJacob Faibussowitsch       PetscCheck(flg,PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMPlexExpandedVerticesCoordinatesToFaces_Private produced a wrong IS");
11605f80ce2aSJacob Faibussowitsch       CHKERRQ(ISDestroy(&is0));
1161c4762a1bSJed Brown     }
1162c4762a1bSJed Brown   } else {
11635f80ce2aSJacob Faibussowitsch     CHKERRQ(VecCreateSeq(PETSC_COMM_SELF, 0, &bnd0->coordinates));
1164c4762a1bSJed Brown   }
1165c4762a1bSJed Brown 
1166c4762a1bSJed Brown   {
1167c4762a1bSJed Brown     Vec         tmp;
1168c4762a1bSJed Brown     VecScatter  sc;
1169c4762a1bSJed Brown     IS          xis;
1170c4762a1bSJed Brown     PetscInt    n;
1171c4762a1bSJed Brown 
1172c4762a1bSJed Brown     /* just convert seq vectors to mpi vector */
11735f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetLocalSize(bnd0->coordinates, &n));
11745f80ce2aSJacob Faibussowitsch     CHKERRMPI(MPI_Bcast(&n, 1, MPIU_INT, rootrank, comm));
1175c4762a1bSJed Brown     if (rank == rootrank) {
11765f80ce2aSJacob Faibussowitsch       CHKERRQ(VecCreateMPI(comm, n, n, &tmp));
1177c4762a1bSJed Brown     } else {
11785f80ce2aSJacob Faibussowitsch       CHKERRQ(VecCreateMPI(comm, 0, n, &tmp));
1179c4762a1bSJed Brown     }
11805f80ce2aSJacob Faibussowitsch     CHKERRQ(VecCopy(bnd0->coordinates, tmp));
11815f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&bnd0->coordinates));
1182c4762a1bSJed Brown     bnd0->coordinates = tmp;
1183c4762a1bSJed Brown 
1184c4762a1bSJed Brown     /* replicate coordinates from root rank to all ranks */
11855f80ce2aSJacob Faibussowitsch     CHKERRQ(VecCreateMPI(comm, n, n*size, &bnd->coordinates));
11865f80ce2aSJacob Faibussowitsch     CHKERRQ(ISCreateStride(comm, n, 0, 1, &xis));
11875f80ce2aSJacob Faibussowitsch     CHKERRQ(VecScatterCreate(bnd0->coordinates, xis, bnd->coordinates, NULL, &sc));
11885f80ce2aSJacob Faibussowitsch     CHKERRQ(VecScatterBegin(sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD));
11895f80ce2aSJacob Faibussowitsch     CHKERRQ(VecScatterEnd(  sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD));
11905f80ce2aSJacob Faibussowitsch     CHKERRQ(VecScatterDestroy(&sc));
11915f80ce2aSJacob Faibussowitsch     CHKERRQ(ISDestroy(&xis));
1192c4762a1bSJed Brown   }
1193c4762a1bSJed Brown   bnd->depth = bnd0->depth;
11945f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Bcast(&bnd->depth, 1, MPIU_INT, rootrank, comm));
11955f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscMalloc1(bnd->depth, &bnd->sections));
1196c4762a1bSJed Brown   for (d=0; d<bnd->depth; d++) {
11975f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscSectionReplicate_Private(comm, rootrank, (rank == rootrank) ? bnd0->sections[d] : NULL, &bnd->sections[d]));
1198c4762a1bSJed Brown   }
1199c4762a1bSJed Brown 
1200c4762a1bSJed Brown   if (rank == rootrank) {
12015f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexRestoreConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections));
1202c4762a1bSJed Brown   }
12035f80ce2aSJacob Faibussowitsch   CHKERRQ(PortableBoundaryDestroy(&bnd0));
12045f80ce2aSJacob Faibussowitsch   CHKERRQ(DMRemoveLabelBySelf(idm, &label, PETSC_TRUE));
12055f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelDestroy(&label));
12065f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&boundary_is));
12075f80ce2aSJacob Faibussowitsch   CHKERRQ(DMDestroy(&idm));
1208c4762a1bSJed Brown   *boundary = bnd;
1209c4762a1bSJed Brown   PetscFunctionReturn(0);
1210c4762a1bSJed Brown }
1211c4762a1bSJed Brown 
1212c4762a1bSJed Brown /* get faces of inter-partition interface */
1213c4762a1bSJed Brown static PetscErrorCode DMPlexGetInterfaceFaces_Private(DM ipdm, IS boundary_faces_is, IS *interface_faces_is)
1214c4762a1bSJed Brown {
1215c4762a1bSJed Brown   MPI_Comm               comm;
1216c4762a1bSJed Brown   DMLabel                label;
1217c4762a1bSJed Brown   IS                     part_boundary_faces_is;
1218c4762a1bSJed Brown   const char             partBoundaryName[] = "DMPlexDistributeInterpolateMarkInterface_partBoundary";
1219c4762a1bSJed Brown   PetscInt               value = 1;
1220c4762a1bSJed Brown   DMPlexInterpolatedFlag intp;
1221c4762a1bSJed Brown 
1222c4762a1bSJed Brown   PetscFunctionBegin;
12235f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)ipdm, &comm));
12245f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexIsInterpolatedCollective(ipdm, &intp));
12252c71b3e2SJacob Faibussowitsch   PetscCheckFalse(intp != DMPLEX_INTERPOLATED_FULL,comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex");
1226c4762a1bSJed Brown 
1227c4762a1bSJed Brown   /* get ipdm partition boundary (partBoundary) */
12285f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, partBoundaryName, &label));
12295f80ce2aSJacob Faibussowitsch   CHKERRQ(DMAddLabel(ipdm, label));
12305f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexMarkBoundaryFaces(ipdm, value, label));
12315f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelViewFromOptionsOnComm_Private(label, "-ipdm_part_boundary_view", comm));
12325f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelGetStratumISOnComm_Private(label, value, comm, &part_boundary_faces_is));
12335f80ce2aSJacob Faibussowitsch   CHKERRQ(DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE));
12345f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelDestroy(&label));
1235c4762a1bSJed Brown 
1236c4762a1bSJed Brown   /* remove ipdm whole-domain boundary (boundary_faces_is) from ipdm partition boundary (part_boundary_faces_is), resulting just in inter-partition interface */
12375f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDifference(part_boundary_faces_is,boundary_faces_is,interface_faces_is));
12385f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&part_boundary_faces_is));
1239c4762a1bSJed Brown   PetscFunctionReturn(0);
1240c4762a1bSJed Brown }
1241c4762a1bSJed Brown 
1242c4762a1bSJed Brown /* compute inter-partition interface including edges and vertices */
1243c4762a1bSJed Brown static PetscErrorCode DMPlexComputeCompleteInterface_Private(DM ipdm, IS interface_faces_is, IS *interface_is)
1244c4762a1bSJed Brown {
1245c4762a1bSJed Brown   DMLabel                label;
1246c4762a1bSJed Brown   PetscInt               value = 1;
1247c4762a1bSJed Brown   const char             interfaceName[] = "DMPlexDistributeInterpolateMarkInterface_interface";
1248c4762a1bSJed Brown   DMPlexInterpolatedFlag intp;
1249c4762a1bSJed Brown   MPI_Comm               comm;
1250c4762a1bSJed Brown 
1251c4762a1bSJed Brown   PetscFunctionBegin;
12525f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)ipdm, &comm));
12535f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexIsInterpolatedCollective(ipdm, &intp));
12542c71b3e2SJacob Faibussowitsch   PetscCheckFalse(intp != DMPLEX_INTERPOLATED_FULL,comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex");
1255c4762a1bSJed Brown 
12565f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, interfaceName, &label));
12575f80ce2aSJacob Faibussowitsch   CHKERRQ(DMAddLabel(ipdm, label));
12585f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelSetStratumIS(label, value, interface_faces_is));
12595f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelViewFromOptionsOnComm_Private(label, "-interface_faces_view", comm));
12605f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexLabelComplete(ipdm, label));
12615f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelViewFromOptionsOnComm_Private(label, "-interface_view", comm));
12625f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelGetStratumISOnComm_Private(label, value, comm, interface_is));
12635f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectSetName((PetscObject)*interface_is, "interface_is"));
12645f80ce2aSJacob Faibussowitsch   CHKERRQ(ISViewFromOptions(*interface_is, NULL, "-interface_is_view"));
12655f80ce2aSJacob Faibussowitsch   CHKERRQ(DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE));
12665f80ce2aSJacob Faibussowitsch   CHKERRQ(DMLabelDestroy(&label));
1267c4762a1bSJed Brown   PetscFunctionReturn(0);
1268c4762a1bSJed Brown }
1269c4762a1bSJed Brown 
1270c4762a1bSJed Brown static PetscErrorCode PointSFGetOutwardInterfacePoints(PetscSF sf, IS *is)
1271c4762a1bSJed Brown {
1272c4762a1bSJed Brown   PetscInt        n;
1273c4762a1bSJed Brown   const PetscInt  *arr;
1274c4762a1bSJed Brown 
1275c4762a1bSJed Brown   PetscFunctionBegin;
12765f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSFGetGraph(sf, NULL, &n, &arr, NULL));
12775f80ce2aSJacob Faibussowitsch   CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_USE_POINTER, is));
1278c4762a1bSJed Brown   PetscFunctionReturn(0);
1279c4762a1bSJed Brown }
1280c4762a1bSJed Brown 
1281c4762a1bSJed Brown static PetscErrorCode PointSFGetInwardInterfacePoints(PetscSF sf, IS *is)
1282c4762a1bSJed Brown {
1283c4762a1bSJed Brown   PetscInt        n;
1284c4762a1bSJed Brown   const PetscInt  *rootdegree;
1285c4762a1bSJed Brown   PetscInt        *arr;
1286c4762a1bSJed Brown 
1287c4762a1bSJed Brown   PetscFunctionBegin;
12885f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSFSetUp(sf));
12895f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSFComputeDegreeBegin(sf, &rootdegree));
12905f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSFComputeDegreeEnd(sf, &rootdegree));
12915f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSFComputeMultiRootOriginalNumbering(sf, rootdegree, &n, &arr));
12925f80ce2aSJacob Faibussowitsch   CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_OWN_POINTER, is));
1293c4762a1bSJed Brown   PetscFunctionReturn(0);
1294c4762a1bSJed Brown }
1295c4762a1bSJed Brown 
1296c4762a1bSJed Brown static PetscErrorCode PointSFGetInterfacePoints_Private(PetscSF pointSF, IS *is)
1297c4762a1bSJed Brown {
1298c4762a1bSJed Brown   IS pointSF_out_is, pointSF_in_is;
1299c4762a1bSJed Brown 
1300c4762a1bSJed Brown   PetscFunctionBegin;
13015f80ce2aSJacob Faibussowitsch   CHKERRQ(PointSFGetOutwardInterfacePoints(pointSF, &pointSF_out_is));
13025f80ce2aSJacob Faibussowitsch   CHKERRQ(PointSFGetInwardInterfacePoints(pointSF, &pointSF_in_is));
13035f80ce2aSJacob Faibussowitsch   CHKERRQ(ISExpand(pointSF_out_is, pointSF_in_is, is));
13045f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&pointSF_out_is));
13055f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&pointSF_in_is));
1306c4762a1bSJed Brown   PetscFunctionReturn(0);
1307c4762a1bSJed Brown }
1308c4762a1bSJed Brown 
13095f80ce2aSJacob Faibussowitsch #define CHKERRMY(ierr) PetscCheck(!ierr,PETSC_COMM_SELF, PETSC_ERR_PLIB, "PointSF is wrong. Unable to show details!")
1310c4762a1bSJed Brown 
1311c4762a1bSJed Brown static PetscErrorCode ViewPointsWithType_Internal(DM dm, IS pointsIS, PetscViewer v)
1312c4762a1bSJed Brown {
1313c4762a1bSJed Brown   DMLabel         label;
1314c4762a1bSJed Brown   PetscSection    coordsSection;
1315c4762a1bSJed Brown   Vec             coordsVec;
1316c4762a1bSJed Brown   PetscScalar     *coordsScalar;
1317c4762a1bSJed Brown   PetscInt        coneSize, depth, dim, i, p, npoints;
1318c4762a1bSJed Brown   const PetscInt  *points;
1319c4762a1bSJed Brown 
1320c4762a1bSJed Brown   PetscFunctionBegin;
13215f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetDimension(dm, &dim));
13225f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetCoordinateSection(dm, &coordsSection));
13235f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetCoordinatesLocal(dm, &coordsVec));
13245f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArray(coordsVec, &coordsScalar));
13255f80ce2aSJacob Faibussowitsch   CHKERRQ(ISGetLocalSize(pointsIS, &npoints));
13265f80ce2aSJacob Faibussowitsch   CHKERRQ(ISGetIndices(pointsIS, &points));
13275f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexGetDepthLabel(dm, &label));
13285f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPushTab(v));
1329c4762a1bSJed Brown   for (i=0; i<npoints; i++) {
1330c4762a1bSJed Brown     p = points[i];
13315f80ce2aSJacob Faibussowitsch     CHKERRQ(DMLabelGetValue(label, p, &depth));
1332c4762a1bSJed Brown     if (!depth) {
1333d3e1f4ccSVaclav Hapla       PetscInt        n, o;
1334c4762a1bSJed Brown       char            coordstr[128];
1335c4762a1bSJed Brown 
13365f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscSectionGetDof(coordsSection, p, &n));
13375f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscSectionGetOffset(coordsSection, p, &o));
13385f80ce2aSJacob Faibussowitsch       CHKERRQ(coord2str(coordstr, sizeof(coordstr), n, &coordsScalar[o], 1.0));
13395f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscViewerASCIISynchronizedPrintf(v, "vertex %D w/ coordinates %s\n", p, coordstr));
1340c4762a1bSJed Brown     } else {
1341c4762a1bSJed Brown       char            entityType[16];
1342c4762a1bSJed Brown 
1343c4762a1bSJed Brown       switch (depth) {
13445f80ce2aSJacob Faibussowitsch         case 1: CHKERRQ(PetscStrcpy(entityType, "edge")); break;
13455f80ce2aSJacob Faibussowitsch         case 2: CHKERRQ(PetscStrcpy(entityType, "face")); break;
13465f80ce2aSJacob Faibussowitsch         case 3: CHKERRQ(PetscStrcpy(entityType, "cell")); break;
13472a27bf02SStefano Zampini         default: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for depth <= 3");
1348c4762a1bSJed Brown       }
1349c4762a1bSJed Brown       if (depth == dim && dim < 3) {
13505f80ce2aSJacob Faibussowitsch         CHKERRQ(PetscStrlcat(entityType, " (cell)", sizeof(entityType)));
1351c4762a1bSJed Brown       }
13525f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscViewerASCIISynchronizedPrintf(v, "%s %D\n", entityType, p));
1353c4762a1bSJed Brown     }
13545f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexGetConeSize(dm, p, &coneSize));
1355c4762a1bSJed Brown     if (coneSize) {
1356c4762a1bSJed Brown       const PetscInt *cone;
1357c4762a1bSJed Brown       IS             coneIS;
1358c4762a1bSJed Brown 
13595f80ce2aSJacob Faibussowitsch       CHKERRQ(DMPlexGetCone(dm, p, &cone));
13605f80ce2aSJacob Faibussowitsch       CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, coneSize, cone, PETSC_USE_POINTER, &coneIS));
13615f80ce2aSJacob Faibussowitsch       CHKERRQ(ViewPointsWithType_Internal(dm, coneIS, v));
13625f80ce2aSJacob Faibussowitsch       CHKERRQ(ISDestroy(&coneIS));
1363c4762a1bSJed Brown     }
1364c4762a1bSJed Brown   }
13655f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPopTab(v));
13665f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArray(coordsVec, &coordsScalar));
13675f80ce2aSJacob Faibussowitsch   CHKERRQ(ISRestoreIndices(pointsIS, &points));
1368c4762a1bSJed Brown   PetscFunctionReturn(0);
1369c4762a1bSJed Brown }
1370c4762a1bSJed Brown 
1371c4762a1bSJed Brown static PetscErrorCode ViewPointsWithType(DM dm, IS points, PetscViewer v)
1372c4762a1bSJed Brown {
1373c4762a1bSJed Brown   PetscBool       flg;
1374c4762a1bSJed Brown   PetscInt        npoints;
1375c4762a1bSJed Brown   PetscMPIInt     rank;
1376c4762a1bSJed Brown 
1377c4762a1bSJed Brown   PetscFunctionBegin;
13785f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &flg));
137928b400f6SJacob Faibussowitsch   PetscCheck(flg,PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for ASCII viewer");
13805f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));
13815f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPushSynchronized(v));
13825f80ce2aSJacob Faibussowitsch   CHKERRQ(ISGetLocalSize(points, &npoints));
1383c4762a1bSJed Brown   if (npoints) {
13845f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", rank));
13855f80ce2aSJacob Faibussowitsch     CHKERRQ(ViewPointsWithType_Internal(dm, points, v));
1386c4762a1bSJed Brown   }
13875f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerFlush(v));
13885f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerASCIIPopSynchronized(v));
1389c4762a1bSJed Brown   PetscFunctionReturn(0);
1390c4762a1bSJed Brown }
1391c4762a1bSJed Brown 
1392c4762a1bSJed Brown static PetscErrorCode DMPlexComparePointSFWithInterface_Private(DM ipdm, IS interface_is)
1393c4762a1bSJed Brown {
1394c4762a1bSJed Brown   PetscSF         pointsf;
1395c4762a1bSJed Brown   IS              pointsf_is;
1396c4762a1bSJed Brown   PetscBool       flg;
1397c4762a1bSJed Brown   MPI_Comm        comm;
1398c4762a1bSJed Brown   PetscMPIInt     size;
1399c4762a1bSJed Brown 
1400c4762a1bSJed Brown   PetscFunctionBegin;
14015f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)ipdm, &comm));
14025f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(comm, &size));
14035f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetPointSF(ipdm, &pointsf));
1404c4762a1bSJed Brown   if (pointsf) {
1405c4762a1bSJed Brown     PetscInt nroots;
14065f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscSFGetGraph(pointsf, &nroots, NULL, NULL, NULL));
1407c4762a1bSJed Brown     if (nroots < 0) pointsf = NULL; /* uninitialized SF */
1408c4762a1bSJed Brown   }
1409c4762a1bSJed Brown   if (!pointsf) {
1410c4762a1bSJed Brown     PetscInt N=0;
14115f80ce2aSJacob Faibussowitsch     if (interface_is) CHKERRQ(ISGetSize(interface_is, &N));
141228b400f6SJacob Faibussowitsch     PetscCheck(!N,comm, PETSC_ERR_PLIB, "interface_is should be NULL or empty for PointSF being NULL");
1413c4762a1bSJed Brown     PetscFunctionReturn(0);
1414c4762a1bSJed Brown   }
1415c4762a1bSJed Brown 
1416c4762a1bSJed Brown   /* get PointSF points as IS pointsf_is */
14175f80ce2aSJacob Faibussowitsch   CHKERRQ(PointSFGetInterfacePoints_Private(pointsf, &pointsf_is));
1418c4762a1bSJed Brown 
1419c4762a1bSJed Brown   /* compare pointsf_is with interface_is */
14205f80ce2aSJacob Faibussowitsch   CHKERRQ(ISEqual(interface_is, pointsf_is, &flg));
14215f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Allreduce(MPI_IN_PLACE,&flg,1,MPIU_BOOL,MPI_LAND,comm));
1422c4762a1bSJed Brown   if (!flg) {
1423c4762a1bSJed Brown     IS pointsf_extra_is, pointsf_missing_is;
1424c4762a1bSJed Brown     PetscViewer errv = PETSC_VIEWER_STDERR_(comm);
14255f80ce2aSJacob Faibussowitsch     CHKERRMY(ISDifference(interface_is, pointsf_is, &pointsf_missing_is));
14265f80ce2aSJacob Faibussowitsch     CHKERRMY(ISDifference(pointsf_is, interface_is, &pointsf_extra_is));
14275f80ce2aSJacob Faibussowitsch     CHKERRMY(PetscViewerASCIIPrintf(errv, "Points missing in PointSF:\n"));
14285f80ce2aSJacob Faibussowitsch     CHKERRMY(ViewPointsWithType(ipdm, pointsf_missing_is, errv));
14295f80ce2aSJacob Faibussowitsch     CHKERRMY(PetscViewerASCIIPrintf(errv, "Extra points in PointSF:\n"));
14305f80ce2aSJacob Faibussowitsch     CHKERRMY(ViewPointsWithType(ipdm, pointsf_extra_is, errv));
14315f80ce2aSJacob Faibussowitsch     CHKERRMY(ISDestroy(&pointsf_extra_is));
14325f80ce2aSJacob Faibussowitsch     CHKERRMY(ISDestroy(&pointsf_missing_is));
1433c4762a1bSJed Brown     SETERRQ(comm, PETSC_ERR_PLIB, "PointSF is wrong! See details above.");
1434c4762a1bSJed Brown   }
14355f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&pointsf_is));
1436c4762a1bSJed Brown   PetscFunctionReturn(0);
1437c4762a1bSJed Brown }
1438c4762a1bSJed Brown 
1439c4762a1bSJed Brown /* remove faces & edges from label, leave just vertices */
1440c4762a1bSJed Brown static PetscErrorCode DMPlexISFilterVertices_Private(DM dm, IS points)
1441c4762a1bSJed Brown {
1442c4762a1bSJed Brown   PetscInt        vStart, vEnd;
1443c4762a1bSJed Brown   MPI_Comm        comm;
1444c4762a1bSJed Brown 
1445c4762a1bSJed Brown   PetscFunctionBegin;
14465f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm));
14475f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
14485f80ce2aSJacob Faibussowitsch   CHKERRQ(ISGeneralFilter(points, vStart, vEnd));
1449c4762a1bSJed Brown   PetscFunctionReturn(0);
1450c4762a1bSJed Brown }
1451c4762a1bSJed Brown 
1452c4762a1bSJed Brown /*
14532e58f0efSBarry Smith   DMPlexCheckPointSFHeavy - Thoroughly test that the PointSF after parallel DMPlexInterpolate() includes exactly all interface points.
1454c4762a1bSJed Brown 
1455c4762a1bSJed Brown   Collective
1456c4762a1bSJed Brown 
1457c4762a1bSJed Brown   Input Parameters:
1458c4762a1bSJed Brown . dm - The DMPlex object
1459c4762a1bSJed Brown 
1460c4762a1bSJed Brown   Notes:
1461c4762a1bSJed Brown   The input DMPlex must be serial (one partition has all points, the other partitions have no points).
1462c4762a1bSJed 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).
1463c4762a1bSJed Brown   This is mainly intended for debugging/testing purposes.
1464c4762a1bSJed Brown 
1465c4762a1bSJed Brown   Algorithm:
1466c4762a1bSJed Brown   1. boundary faces of the serial version of the whole mesh are found using DMPlexMarkBoundaryFaces()
1467c4762a1bSJed 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
1468c4762a1bSJed Brown   3. the mesh is distributed or loaded in parallel
1469c4762a1bSJed Brown   4. boundary faces of the distributed mesh are reconstructed from PortableBoundary using DMPlexFindVertices()
1470c4762a1bSJed Brown   5. partition boundary faces of the parallel mesh are found using DMPlexMarkBoundaryFaces()
1471c4762a1bSJed Brown   6. partition interfaces are computed as set difference of partition boundary faces minus the reconstructed boundary
1472c4762a1bSJed 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
1473c4762a1bSJed Brown 
1474c4762a1bSJed Brown   Level: developer
1475c4762a1bSJed Brown 
1476c4762a1bSJed Brown .seealso: DMGetPointSF(), DMPlexCheckSymmetry(), DMPlexCheckSkeleton(), DMPlexCheckFaces()
1477c4762a1bSJed Brown */
1478c4762a1bSJed Brown static PetscErrorCode DMPlexCheckPointSFHeavy(DM dm, PortableBoundary bnd)
1479c4762a1bSJed Brown {
1480c4762a1bSJed Brown   DM                     ipdm=NULL;
1481c4762a1bSJed Brown   IS                     boundary_faces_is, interface_faces_is, interface_is;
1482c4762a1bSJed Brown   DMPlexInterpolatedFlag intp;
1483c4762a1bSJed Brown   MPI_Comm               comm;
1484c4762a1bSJed Brown 
1485c4762a1bSJed Brown   PetscFunctionBegin;
14865f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm));
1487c4762a1bSJed Brown 
14885f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexIsInterpolatedCollective(dm, &intp));
1489c4762a1bSJed Brown   if (intp == DMPLEX_INTERPOLATED_FULL) {
1490c4762a1bSJed Brown     ipdm = dm;
1491c4762a1bSJed Brown   } else {
1492c4762a1bSJed Brown     /* create temporary interpolated DM if input DM is not interpolated */
14935f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexSetOrientInterface_Private(dm, PETSC_FALSE));
14945f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexInterpolate(dm, &ipdm)); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexComparePointSFWithInterface_Private() below */
14955f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexSetOrientInterface_Private(dm, PETSC_TRUE));
1496c4762a1bSJed Brown   }
14975f80ce2aSJacob Faibussowitsch   CHKERRQ(DMViewFromOptions(ipdm, NULL, "-ipdm_view"));
1498c4762a1bSJed Brown 
1499c4762a1bSJed Brown   /* recover ipdm whole-domain boundary faces from the expanded vertices coordinates */
15005f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexExpandedVerticesCoordinatesToFaces_Private(ipdm, bnd, &boundary_faces_is));
1501c4762a1bSJed Brown   /* get inter-partition interface faces (interface_faces_is)*/
15025f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexGetInterfaceFaces_Private(ipdm, boundary_faces_is, &interface_faces_is));
1503c4762a1bSJed Brown   /* compute inter-partition interface including edges and vertices (interface_is) */
15045f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexComputeCompleteInterface_Private(ipdm, interface_faces_is, &interface_is));
1505c4762a1bSJed Brown   /* destroy immediate ISs */
15065f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&boundary_faces_is));
15075f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&interface_faces_is));
1508c4762a1bSJed Brown 
1509c4762a1bSJed Brown   /* for uninterpolated case, keep just vertices in interface */
1510c4762a1bSJed Brown   if (!intp) {
15115f80ce2aSJacob Faibussowitsch     CHKERRQ(DMPlexISFilterVertices_Private(ipdm, interface_is));
15125f80ce2aSJacob Faibussowitsch     CHKERRQ(DMDestroy(&ipdm));
1513c4762a1bSJed Brown   }
1514c4762a1bSJed Brown 
1515c4762a1bSJed Brown   /* compare PointSF with the boundary reconstructed from coordinates */
15165f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexComparePointSFWithInterface_Private(dm, interface_is));
15175f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(comm, "DMPlexCheckPointSFHeavy PASSED\n"));
15185f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&interface_is));
1519c4762a1bSJed Brown   PetscFunctionReturn(0);
1520c4762a1bSJed Brown }
1521c4762a1bSJed Brown 
1522c4762a1bSJed Brown int main(int argc, char **argv)
1523c4762a1bSJed Brown {
1524c4762a1bSJed Brown   DM             dm;
1525c4762a1bSJed Brown   AppCtx         user;
1526c4762a1bSJed Brown 
1527*b122ec5aSJacob Faibussowitsch   CHKERRQ(PetscInitialize(&argc, &argv, NULL, help));
15285f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscLogStageRegister("create",&stage[0]));
15295f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscLogStageRegister("distribute",&stage[1]));
15305f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscLogStageRegister("interpolate",&stage[2]));
15315f80ce2aSJacob Faibussowitsch   CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user));
15325f80ce2aSJacob Faibussowitsch   CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1533c4762a1bSJed Brown   if (user.nPointsToExpand) {
15345f80ce2aSJacob Faibussowitsch     CHKERRQ(TestExpandPoints(dm, &user));
1535c4762a1bSJed Brown   }
1536c4762a1bSJed Brown   if (user.ncoords) {
1537d3e1f4ccSVaclav Hapla     Vec coords;
1538d3e1f4ccSVaclav Hapla 
15395f80ce2aSJacob Faibussowitsch     CHKERRQ(VecCreateSeqWithArray(PETSC_COMM_SELF, user.ncoords, user.ncoords, user.coords, &coords));
15405f80ce2aSJacob Faibussowitsch     CHKERRQ(ViewVerticesFromCoords(dm, coords, user.coordsTol, PETSC_VIEWER_STDOUT_WORLD));
15415f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&coords));
1542c4762a1bSJed Brown   }
15435f80ce2aSJacob Faibussowitsch   CHKERRQ(DMDestroy(&dm));
1544*b122ec5aSJacob Faibussowitsch   CHKERRQ(PetscFinalize());
1545*b122ec5aSJacob Faibussowitsch   return 0;
1546c4762a1bSJed Brown }
1547c4762a1bSJed Brown 
1548c4762a1bSJed Brown /*TEST
1549c4762a1bSJed Brown 
1550c4762a1bSJed Brown   testset:
1551c4762a1bSJed Brown     nsize: 2
1552c4762a1bSJed Brown     args: -dm_view ascii::ascii_info_detail
1553c4762a1bSJed Brown     args: -dm_plex_check_all
1554c4762a1bSJed Brown     test:
1555c4762a1bSJed Brown       suffix: 1_tri_dist0
1556c4762a1bSJed Brown       args: -distribute 0 -interpolate {{none create}separate output}
1557c4762a1bSJed Brown     test:
1558c4762a1bSJed Brown       suffix: 1_tri_dist1
1559c4762a1bSJed Brown       args: -distribute 1 -interpolate {{none create after_distribute}separate output}
1560c4762a1bSJed Brown     test:
1561c4762a1bSJed Brown       suffix: 1_quad_dist0
1562c4762a1bSJed Brown       args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output}
1563c4762a1bSJed Brown     test:
1564c4762a1bSJed Brown       suffix: 1_quad_dist1
1565c4762a1bSJed Brown       args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output}
1566c4762a1bSJed Brown     test:
1567c4762a1bSJed Brown       suffix: 1_1d_dist1
1568c4762a1bSJed Brown       args: -dim 1 -distribute 1
1569c4762a1bSJed Brown 
1570c4762a1bSJed Brown   testset:
1571c4762a1bSJed Brown     nsize: 3
1572c4762a1bSJed Brown     args: -testnum 1 -interpolate create
1573c4762a1bSJed Brown     args: -dm_plex_check_all
1574c4762a1bSJed Brown     test:
1575c4762a1bSJed Brown       suffix: 2
1576c4762a1bSJed Brown       args: -dm_view ascii::ascii_info_detail
1577c4762a1bSJed Brown     test:
1578c4762a1bSJed Brown       suffix: 2a
1579c4762a1bSJed Brown       args: -dm_plex_check_cones_conform_on_interfaces_verbose
1580c4762a1bSJed Brown     test:
1581c4762a1bSJed Brown       suffix: 2b
1582c4762a1bSJed Brown       args: -test_expand_points 0,1,2,5,6
1583c4762a1bSJed Brown     test:
1584c4762a1bSJed Brown       suffix: 2c
1585c4762a1bSJed Brown       args: -test_expand_points 0,1,2,5,6 -test_expand_points_empty
1586c4762a1bSJed Brown 
1587c4762a1bSJed Brown   testset:
1588c4762a1bSJed Brown     # the same as 1% for 3D
1589c4762a1bSJed Brown     nsize: 2
1590c4762a1bSJed Brown     args: -dim 3 -dm_view ascii::ascii_info_detail
1591c4762a1bSJed Brown     args: -dm_plex_check_all
1592c4762a1bSJed Brown     test:
1593c4762a1bSJed Brown       suffix: 4_tet_dist0
1594c4762a1bSJed Brown       args: -distribute 0 -interpolate {{none create}separate output}
1595c4762a1bSJed Brown     test:
1596c4762a1bSJed Brown       suffix: 4_tet_dist1
1597c4762a1bSJed Brown       args: -distribute 1 -interpolate {{none create after_distribute}separate output}
1598c4762a1bSJed Brown     test:
1599c4762a1bSJed Brown       suffix: 4_hex_dist0
1600c4762a1bSJed Brown       args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output}
1601c4762a1bSJed Brown     test:
1602c4762a1bSJed Brown       suffix: 4_hex_dist1
1603c4762a1bSJed Brown       args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output}
1604c4762a1bSJed Brown 
1605c4762a1bSJed Brown   test:
1606c4762a1bSJed Brown     # the same as 4_tet_dist0 but test different initial orientations
1607c4762a1bSJed Brown     suffix: 4_tet_test_orient
1608c4762a1bSJed Brown     nsize: 2
1609c4762a1bSJed Brown     args: -dim 3 -distribute 0
1610c4762a1bSJed Brown     args: -dm_plex_check_all
1611c4762a1bSJed Brown     args: -rotate_interface_0 {{0 1 2 11 12 13}}
1612c4762a1bSJed Brown     args: -rotate_interface_1 {{0 1 2 11 12 13}}
1613c4762a1bSJed Brown 
1614c4762a1bSJed Brown   testset:
1615c4762a1bSJed Brown     requires: exodusii
1616c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/TwoQuads.exo
1617c4762a1bSJed Brown     args: -dm_view ascii::ascii_info_detail
1618c4762a1bSJed Brown     args: -dm_plex_check_all
1619c4762a1bSJed Brown     args: -custom_view
1620c4762a1bSJed Brown     test:
1621c4762a1bSJed Brown       suffix: 5_seq
1622c4762a1bSJed Brown       nsize: 1
1623c4762a1bSJed Brown       args: -distribute 0 -interpolate {{none create}separate output}
1624c4762a1bSJed Brown     test:
162530602db0SMatthew G. Knepley       # Detail viewing in a non-distributed mesh is broken because the DMLabelView() is collective, but the label is not shared
1626c4762a1bSJed Brown       suffix: 5_dist0
1627c4762a1bSJed Brown       nsize: 2
162830602db0SMatthew G. Knepley       args: -distribute 0 -interpolate {{none create}separate output} -dm_view
1629c4762a1bSJed Brown     test:
1630c4762a1bSJed Brown       suffix: 5_dist1
1631c4762a1bSJed Brown       nsize: 2
1632c4762a1bSJed Brown       args: -distribute 1 -interpolate {{none create after_distribute}separate output}
1633c4762a1bSJed Brown 
1634c4762a1bSJed Brown   testset:
1635c4762a1bSJed Brown     nsize: {{1 2 4}}
1636c4762a1bSJed Brown     args: -use_generator
1637c4762a1bSJed Brown     args: -dm_plex_check_all
1638c4762a1bSJed Brown     args: -distribute -interpolate none
1639c4762a1bSJed Brown     test:
1640c4762a1bSJed Brown       suffix: 6_tri
1641c4762a1bSJed Brown       requires: triangle
1642c0517cd5SMatthew G. Knepley       args: -faces {{2,2  1,3  7,4}} -cell_simplex 1 -dm_generator triangle
1643c4762a1bSJed Brown     test:
1644c4762a1bSJed Brown       suffix: 6_quad
1645c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 0
1646c4762a1bSJed Brown     test:
1647c4762a1bSJed Brown       suffix: 6_tet
1648c4762a1bSJed Brown       requires: ctetgen
1649c0517cd5SMatthew G. Knepley       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 1 -dm_generator ctetgen
1650c4762a1bSJed Brown     test:
1651c4762a1bSJed Brown       suffix: 6_hex
1652c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 0
1653c4762a1bSJed Brown   testset:
1654c4762a1bSJed Brown     nsize: {{1 2 4}}
1655c4762a1bSJed Brown     args: -use_generator
1656c4762a1bSJed Brown     args: -dm_plex_check_all
1657c4762a1bSJed Brown     args: -distribute -interpolate create
1658c4762a1bSJed Brown     test:
1659c4762a1bSJed Brown       suffix: 6_int_tri
1660c4762a1bSJed Brown       requires: triangle
1661c0517cd5SMatthew G. Knepley       args: -faces {{2,2  1,3  7,4}} -cell_simplex 1 -dm_generator triangle
1662c4762a1bSJed Brown     test:
1663c4762a1bSJed Brown       suffix: 6_int_quad
1664c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 0
1665c4762a1bSJed Brown     test:
1666c4762a1bSJed Brown       suffix: 6_int_tet
1667c4762a1bSJed Brown       requires: ctetgen
1668c0517cd5SMatthew G. Knepley       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 1 -dm_generator ctetgen
1669c4762a1bSJed Brown     test:
1670c4762a1bSJed Brown       suffix: 6_int_hex
1671c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 0
1672c4762a1bSJed Brown   testset:
1673c4762a1bSJed Brown     nsize: {{2 4}}
1674c4762a1bSJed Brown     args: -use_generator
1675c4762a1bSJed Brown     args: -dm_plex_check_all
1676c4762a1bSJed Brown     args: -distribute -interpolate after_distribute
1677c4762a1bSJed Brown     test:
1678c4762a1bSJed Brown       suffix: 6_parint_tri
1679c4762a1bSJed Brown       requires: triangle
1680c0517cd5SMatthew G. Knepley       args: -faces {{2,2  1,3  7,4}} -cell_simplex 1 -dm_generator triangle
1681c4762a1bSJed Brown     test:
1682c4762a1bSJed Brown       suffix: 6_parint_quad
1683c4762a1bSJed Brown       args: -faces {{2,2  1,3  7,4}} -cell_simplex 0
1684c4762a1bSJed Brown     test:
1685c4762a1bSJed Brown       suffix: 6_parint_tet
1686c4762a1bSJed Brown       requires: ctetgen
1687c0517cd5SMatthew G. Knepley       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 1 -dm_generator ctetgen
1688c4762a1bSJed Brown     test:
1689c4762a1bSJed Brown       suffix: 6_parint_hex
1690c4762a1bSJed Brown       args: -faces {{2,2,2  1,3,5  3,4,7}} -cell_simplex 0
1691c4762a1bSJed Brown 
1692c4762a1bSJed Brown   testset: # 7 EXODUS
1693c4762a1bSJed Brown     requires: exodusii
1694c4762a1bSJed Brown     args: -dm_plex_check_all
1695c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo
1696c4762a1bSJed Brown     args: -distribute
1697c4762a1bSJed Brown     test: # seq load, simple partitioner
1698c4762a1bSJed Brown       suffix: 7_exo
1699c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1700c4762a1bSJed Brown       args: -interpolate none
1701c4762a1bSJed Brown     test: # seq load, seq interpolation, simple partitioner
1702c4762a1bSJed Brown       suffix: 7_exo_int_simple
1703c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1704c4762a1bSJed Brown       args: -interpolate create
1705c4762a1bSJed Brown     test: # seq load, seq interpolation, metis partitioner
1706c4762a1bSJed Brown       suffix: 7_exo_int_metis
1707c4762a1bSJed Brown       requires: parmetis
1708c4762a1bSJed Brown       nsize: {{2 4 5}}
1709c4762a1bSJed Brown       args: -interpolate create
1710c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1711c4762a1bSJed Brown     test: # seq load, simple partitioner, par interpolation
1712c4762a1bSJed Brown       suffix: 7_exo_simple_int
1713c4762a1bSJed Brown       nsize: {{2 4 5}}
1714c4762a1bSJed Brown       args: -interpolate after_distribute
1715c4762a1bSJed Brown     test: # seq load, metis partitioner, par interpolation
1716c4762a1bSJed Brown       suffix: 7_exo_metis_int
1717c4762a1bSJed Brown       requires: parmetis
1718c4762a1bSJed Brown       nsize: {{2 4 5}}
1719c4762a1bSJed Brown       args: -interpolate after_distribute
1720c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1721c4762a1bSJed Brown 
1722c4762a1bSJed Brown   testset: # 7 HDF5 SEQUANTIAL LOAD
1723c4762a1bSJed Brown     requires: hdf5 !complex
1724c4762a1bSJed Brown     args: -dm_plex_check_all
1725c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf
1726c4762a1bSJed Brown     args: -dm_plex_hdf5_force_sequential
1727c4762a1bSJed Brown     args: -distribute
1728c4762a1bSJed Brown     test: # seq load, simple partitioner
1729c4762a1bSJed Brown       suffix: 7_seq_hdf5_simple
1730c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1731c4762a1bSJed Brown       args: -interpolate none
1732c4762a1bSJed Brown     test: # seq load, seq interpolation, simple partitioner
1733c4762a1bSJed Brown       suffix: 7_seq_hdf5_int_simple
1734c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1735c4762a1bSJed Brown       args: -interpolate after_create
1736c4762a1bSJed Brown     test: # seq load, seq interpolation, metis partitioner
1737c4762a1bSJed Brown       nsize: {{2 4 5}}
1738c4762a1bSJed Brown       suffix: 7_seq_hdf5_int_metis
1739c4762a1bSJed Brown       requires: parmetis
1740c4762a1bSJed Brown       args: -interpolate after_create
1741c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1742c4762a1bSJed Brown     test: # seq load, simple partitioner, par interpolation
1743c4762a1bSJed Brown       suffix: 7_seq_hdf5_simple_int
1744c4762a1bSJed Brown       nsize: {{2 4 5}}
1745c4762a1bSJed Brown       args: -interpolate after_distribute
1746c4762a1bSJed Brown     test: # seq load, metis partitioner, par interpolation
1747c4762a1bSJed Brown       nsize: {{2 4 5}}
1748c4762a1bSJed Brown       suffix: 7_seq_hdf5_metis_int
1749c4762a1bSJed Brown       requires: parmetis
1750c4762a1bSJed Brown       args: -interpolate after_distribute
1751c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1752c4762a1bSJed Brown 
1753c4762a1bSJed Brown   testset: # 7 HDF5 PARALLEL LOAD
1754c4762a1bSJed Brown     requires: hdf5 !complex
1755c4762a1bSJed Brown     nsize: {{2 4 5}}
1756c4762a1bSJed Brown     args: -dm_plex_check_all
1757c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf
1758c4762a1bSJed Brown     test: # par load
1759c4762a1bSJed Brown       suffix: 7_par_hdf5
1760c4762a1bSJed Brown       args: -interpolate none
1761c4762a1bSJed Brown     test: # par load, par interpolation
1762c4762a1bSJed Brown       suffix: 7_par_hdf5_int
1763c4762a1bSJed Brown       args: -interpolate after_create
1764c4762a1bSJed Brown     test: # par load, parmetis repartitioner
1765c4762a1bSJed Brown       TODO: Parallel partitioning of uninterpolated meshes not supported
1766c4762a1bSJed Brown       suffix: 7_par_hdf5_parmetis
1767c4762a1bSJed Brown       requires: parmetis
1768c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1769c4762a1bSJed Brown       args: -interpolate none
1770c4762a1bSJed Brown     test: # par load, par interpolation, parmetis repartitioner
1771c4762a1bSJed Brown       suffix: 7_par_hdf5_int_parmetis
1772c4762a1bSJed Brown       requires: parmetis
1773c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1774c4762a1bSJed Brown       args: -interpolate after_create
1775c4762a1bSJed Brown     test: # par load, parmetis partitioner, par interpolation
1776c4762a1bSJed Brown       TODO: Parallel partitioning of uninterpolated meshes not supported
1777c4762a1bSJed Brown       suffix: 7_par_hdf5_parmetis_int
1778c4762a1bSJed Brown       requires: parmetis
1779c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1780c4762a1bSJed Brown       args: -interpolate after_distribute
1781c4762a1bSJed Brown 
1782c4762a1bSJed Brown     test:
1783c4762a1bSJed Brown       suffix: 7_hdf5_hierarch
1784c4762a1bSJed Brown       requires: hdf5 ptscotch !complex
1785c4762a1bSJed Brown       nsize: {{2 3 4}separate output}
1786c4762a1bSJed Brown       args: -distribute
1787c4762a1bSJed Brown       args: -interpolate after_create
1788c4762a1bSJed Brown       args: -petscpartitioner_type matpartitioning -petscpartitioner_view ::ascii_info
1789c4762a1bSJed Brown       args: -mat_partitioning_type hierarch -mat_partitioning_hierarchical_nfineparts 2
1790c4762a1bSJed Brown       args: -mat_partitioning_hierarchical_coarseparttype ptscotch -mat_partitioning_hierarchical_fineparttype ptscotch
1791c4762a1bSJed Brown 
1792c4762a1bSJed Brown   test:
1793c4762a1bSJed Brown     suffix: 8
1794c4762a1bSJed Brown     requires: hdf5 !complex
1795c4762a1bSJed Brown     nsize: 4
1796c4762a1bSJed Brown     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf
1797c4762a1bSJed Brown     args: -distribute 0 -interpolate after_create
1798c4762a1bSJed 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
1799c4762a1bSJed Brown     args: -dm_plex_check_all
1800c4762a1bSJed Brown     args: -custom_view
1801c4762a1bSJed Brown 
1802c4762a1bSJed Brown   testset: # 9 HDF5 SEQUANTIAL LOAD
1803c4762a1bSJed Brown     requires: hdf5 !complex datafilespath
1804c4762a1bSJed Brown     args: -dm_plex_check_all
1805c4762a1bSJed 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
1806c4762a1bSJed Brown     args: -dm_plex_hdf5_force_sequential
1807c4762a1bSJed Brown     args: -distribute
1808c4762a1bSJed Brown     test: # seq load, simple partitioner
1809c4762a1bSJed Brown       suffix: 9_seq_hdf5_simple
1810c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1811c4762a1bSJed Brown       args: -interpolate none
1812c4762a1bSJed Brown     test: # seq load, seq interpolation, simple partitioner
1813c4762a1bSJed Brown       suffix: 9_seq_hdf5_int_simple
1814c4762a1bSJed Brown       nsize: {{1 2 4 5}}
1815c4762a1bSJed Brown       args: -interpolate after_create
1816c4762a1bSJed Brown     test: # seq load, seq interpolation, metis partitioner
1817c4762a1bSJed Brown       nsize: {{2 4 5}}
1818c4762a1bSJed Brown       suffix: 9_seq_hdf5_int_metis
1819c4762a1bSJed Brown       requires: parmetis
1820c4762a1bSJed Brown       args: -interpolate after_create
1821c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1822c4762a1bSJed Brown     test: # seq load, simple partitioner, par interpolation
1823c4762a1bSJed Brown       suffix: 9_seq_hdf5_simple_int
1824c4762a1bSJed Brown       nsize: {{2 4 5}}
1825c4762a1bSJed Brown       args: -interpolate after_distribute
1826c4762a1bSJed Brown     test: # seq load, simple partitioner, par interpolation
1827c4762a1bSJed Brown       # This is like 9_seq_hdf5_simple_int but testing error output of DMPlexCheckPointSFHeavy().
1828c4762a1bSJed Brown       # Once 9_seq_hdf5_simple_int gets fixed, this one gets broken.
1829c4762a1bSJed Brown       # We can then provide an intentionally broken mesh instead.
1830c4762a1bSJed Brown       TODO: This test is broken because PointSF is fixed.
1831c4762a1bSJed Brown       suffix: 9_seq_hdf5_simple_int_err
1832c4762a1bSJed Brown       nsize: 4
1833c4762a1bSJed Brown       args: -interpolate after_distribute
1834c4762a1bSJed Brown       filter: sed -e "/PETSC ERROR/,$$d"
1835c4762a1bSJed Brown     test: # seq load, metis partitioner, par interpolation
1836c4762a1bSJed Brown       nsize: {{2 4 5}}
1837c4762a1bSJed Brown       suffix: 9_seq_hdf5_metis_int
1838c4762a1bSJed Brown       requires: parmetis
1839c4762a1bSJed Brown       args: -interpolate after_distribute
1840c4762a1bSJed Brown       args: -petscpartitioner_type parmetis
1841c4762a1bSJed Brown 
1842c4762a1bSJed Brown   testset: # 9 HDF5 PARALLEL LOAD
1843c4762a1bSJed Brown     requires: hdf5 !complex datafilespath
1844c4762a1bSJed Brown     nsize: {{2 4 5}}
1845c4762a1bSJed Brown     args: -dm_plex_check_all
1846c4762a1bSJed 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
1847c4762a1bSJed Brown     test: # par load
1848c4762a1bSJed Brown       suffix: 9_par_hdf5
1849c4762a1bSJed Brown       args: -interpolate none
1850c4762a1bSJed Brown     test: # par load, par interpolation
1851c4762a1bSJed Brown       suffix: 9_par_hdf5_int
1852c4762a1bSJed Brown       args: -interpolate after_create
1853c4762a1bSJed Brown     test: # par load, parmetis repartitioner
1854c4762a1bSJed Brown       TODO: Parallel partitioning of uninterpolated meshes not supported
1855c4762a1bSJed Brown       suffix: 9_par_hdf5_parmetis
1856c4762a1bSJed Brown       requires: parmetis
1857c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1858c4762a1bSJed Brown       args: -interpolate none
1859c4762a1bSJed Brown     test: # par load, par interpolation, parmetis repartitioner
1860c4762a1bSJed Brown       suffix: 9_par_hdf5_int_parmetis
1861c4762a1bSJed Brown       requires: parmetis
1862c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1863c4762a1bSJed Brown       args: -interpolate after_create
1864c4762a1bSJed Brown     test: # par load, parmetis partitioner, par interpolation
1865c4762a1bSJed Brown       TODO: Parallel partitioning of uninterpolated meshes not supported
1866c4762a1bSJed Brown       suffix: 9_par_hdf5_parmetis_int
1867c4762a1bSJed Brown       requires: parmetis
1868c4762a1bSJed Brown       args: -distribute -petscpartitioner_type parmetis
1869c4762a1bSJed Brown       args: -interpolate after_distribute
1870c4762a1bSJed Brown 
1871c4762a1bSJed Brown   testset: # 10 HDF5 PARALLEL LOAD
1872c4762a1bSJed Brown     requires: hdf5 !complex datafilespath
1873c4762a1bSJed Brown     nsize: {{2 4 7}}
1874c4762a1bSJed Brown     args: -dm_plex_check_all
1875c4762a1bSJed 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
1876c4762a1bSJed Brown     test: # par load, par interpolation
1877c4762a1bSJed Brown       suffix: 10_par_hdf5_int
1878c4762a1bSJed Brown       args: -interpolate after_create
1879c4762a1bSJed Brown TEST*/
1880