xref: /petsc/src/dm/impls/plex/plexcreate.c (revision 6113b4542af7371eb509ce9949a149bedab1a29f) 
1 #define PETSCDM_DLL
2 #include <petsc-private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
3 #include <petscdmda.h>
4 #include <petscsf.h>
5 
6 #undef __FUNCT__
7 #define __FUNCT__ "DMSetFromOptions_Plex"
8 PetscErrorCode  DMSetFromOptions_Plex(DM dm)
9 {
10   DM_Plex        *mesh = (DM_Plex*) dm->data;
11   PetscErrorCode ierr;
12 
13   PetscFunctionBegin;
14   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
15   ierr = PetscOptionsHead("DMPlex Options");CHKERRQ(ierr);
16   /* Handle DMPlex refinement */
17   /* Handle associated vectors */
18   /* Handle viewing */
19   ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr);
20   ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr);
21   ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", PETSC_FALSE, &mesh->printTol, NULL);CHKERRQ(ierr);
22   ierr = PetscOptionsTail();CHKERRQ(ierr);
23   PetscFunctionReturn(0);
24 }
25 
26 #undef __FUNCT__
27 #define __FUNCT__ "DMPlexCreateSquareBoundary"
28 /*
29  Simple square boundary:
30 
31  18--5-17--4--16
32   |     |     |
33   6    10     3
34   |     |     |
35  19-11-20--9--15
36   |     |     |
37   7     8     2
38   |     |     |
39  12--0-13--1--14
40 */
41 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
42 {
43   PetscInt       numVertices    = (edges[0]+1)*(edges[1]+1);
44   PetscInt       numEdges       = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
45   PetscInt       markerTop      = 1;
46   PetscInt       markerBottom   = 1;
47   PetscInt       markerRight    = 1;
48   PetscInt       markerLeft     = 1;
49   PetscBool      markerSeparate = PETSC_FALSE;
50   Vec            coordinates;
51   PetscSection   coordSection;
52   PetscScalar    *coords;
53   PetscInt       coordSize;
54   PetscMPIInt    rank;
55   PetscInt       v, vx, vy;
56   PetscErrorCode ierr;
57 
58   PetscFunctionBegin;
59   ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr);
60   if (markerSeparate) {
61     markerTop    = 1;
62     markerBottom = 0;
63     markerRight  = 0;
64     markerLeft   = 0;
65   }
66   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
67   if (!rank) {
68     PetscInt e, ex, ey;
69 
70     ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr);
71     for (e = 0; e < numEdges; ++e) {
72       ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr);
73     }
74     ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */
75     for (vx = 0; vx <= edges[0]; vx++) {
76       for (ey = 0; ey < edges[1]; ey++) {
77         PetscInt edge   = vx*edges[1] + ey + edges[0]*(edges[1]+1);
78         PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
79         PetscInt cone[2];
80 
81         cone[0] = vertex; cone[1] = vertex+edges[0]+1;
82         ierr    = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
83         if (vx == edges[0]) {
84           ierr = DMPlexSetLabelValue(dm, "marker", edge,    markerRight);CHKERRQ(ierr);
85           ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr);
86           if (ey == edges[1]-1) {
87             ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr);
88           }
89         } else if (vx == 0) {
90           ierr = DMPlexSetLabelValue(dm, "marker", edge,    markerLeft);CHKERRQ(ierr);
91           ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr);
92           if (ey == edges[1]-1) {
93             ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr);
94           }
95         }
96       }
97     }
98     for (vy = 0; vy <= edges[1]; vy++) {
99       for (ex = 0; ex < edges[0]; ex++) {
100         PetscInt edge   = vy*edges[0]     + ex;
101         PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
102         PetscInt cone[2];
103 
104         cone[0] = vertex; cone[1] = vertex+1;
105         ierr    = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
106         if (vy == edges[1]) {
107           ierr = DMPlexSetLabelValue(dm, "marker", edge,    markerTop);CHKERRQ(ierr);
108           ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr);
109           if (ex == edges[0]-1) {
110             ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr);
111           }
112         } else if (vy == 0) {
113           ierr = DMPlexSetLabelValue(dm, "marker", edge,    markerBottom);CHKERRQ(ierr);
114           ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr);
115           if (ex == edges[0]-1) {
116             ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr);
117           }
118         }
119       }
120     }
121   }
122   ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
123   ierr = DMPlexStratify(dm);CHKERRQ(ierr);
124   /* Build coordinates */
125   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
126   ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr);
127   for (v = numEdges; v < numEdges+numVertices; ++v) {
128     ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr);
129   }
130   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
131   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
132   ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
133   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
134   ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
135   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
136   for (vy = 0; vy <= edges[1]; ++vy) {
137     for (vx = 0; vx <= edges[0]; ++vx) {
138       coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
139       coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
140     }
141   }
142   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
143   ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
144   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
145   PetscFunctionReturn(0);
146 }
147 
148 #undef __FUNCT__
149 #define __FUNCT__ "DMPlexCreateCubeBoundary"
150 /*
151  Simple cubic boundary:
152 
153      2-------3
154     /|      /|
155    6-------7 |
156    | |     | |
157    | 0-----|-1
158    |/      |/
159    4-------5
160 */
161 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
162 {
163   PetscInt       numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
164   PetscInt       numFaces    = 6;
165   Vec            coordinates;
166   PetscSection   coordSection;
167   PetscScalar    *coords;
168   PetscInt       coordSize;
169   PetscMPIInt    rank;
170   PetscInt       v, vx, vy, vz;
171   PetscErrorCode ierr;
172 
173   PetscFunctionBegin;
174   if ((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Must have at least 1 face per side");
175   if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Currently can't handle more than 1 face per side");
176   ierr = PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);CHKERRQ(ierr);
177   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
178   if (!rank) {
179     PetscInt f;
180 
181     ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr);
182     for (f = 0; f < numFaces; ++f) {
183       ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr);
184     }
185     ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */
186     for (v = 0; v < numFaces+numVertices; ++v) {
187       ierr = DMPlexSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr);
188     }
189     { /* Side 0 (Front) */
190       PetscInt cone[4];
191       cone[0] = numFaces+4; cone[1] = numFaces+5; cone[2] = numFaces+7; cone[3] = numFaces+6;
192       ierr    = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr);
193     }
194     { /* Side 1 (Back) */
195       PetscInt cone[4];
196       cone[0] = numFaces+1; cone[1] = numFaces+0; cone[2] = numFaces+2; cone[3] = numFaces+3;
197       ierr    = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr);
198     }
199     { /* Side 2 (Bottom) */
200       PetscInt cone[4];
201       cone[0] = numFaces+0; cone[1] = numFaces+1; cone[2] = numFaces+5; cone[3] = numFaces+4;
202       ierr    = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr);
203     }
204     { /* Side 3 (Top) */
205       PetscInt cone[4];
206       cone[0] = numFaces+6; cone[1] = numFaces+7; cone[2] = numFaces+3; cone[3] = numFaces+2;
207       ierr    = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr);
208     }
209     { /* Side 4 (Left) */
210       PetscInt cone[4];
211       cone[0] = numFaces+0; cone[1] = numFaces+4; cone[2] = numFaces+6; cone[3] = numFaces+2;
212       ierr    = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr);
213     }
214     { /* Side 5 (Right) */
215       PetscInt cone[4];
216       cone[0] = numFaces+5; cone[1] = numFaces+1; cone[2] = numFaces+3; cone[3] = numFaces+7;
217       ierr    = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr);
218     }
219   }
220   ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
221   ierr = DMPlexStratify(dm);CHKERRQ(ierr);
222   /* Build coordinates */
223   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
224   ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr);
225   for (v = numFaces; v < numFaces+numVertices; ++v) {
226     ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr);
227   }
228   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
229   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
230   ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
231   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
232   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
233   ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
234   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
235   for (vz = 0; vz <= faces[2]; ++vz) {
236     for (vy = 0; vy <= faces[1]; ++vy) {
237       for (vx = 0; vx <= faces[0]; ++vx) {
238         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
239         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
240         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
241       }
242     }
243   }
244   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
245   ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
246   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
247   PetscFunctionReturn(0);
248 }
249 
250 #undef __FUNCT__
251 #define __FUNCT__ "DMPlexCreateSquareMesh"
252 /*
253  Simple square mesh:
254 
255  22--8-23--9--24
256   |     |     |
257  13  2 14  3  15
258   |     |     |
259  19--6-20--7--21
260   |     |     |
261  10  0 11  1 12
262   |     |     |
263  16--4-17--5--18
264 */
265 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
266 {
267   PetscInt       markerTop      = 1;
268   PetscInt       markerBottom   = 1;
269   PetscInt       markerRight    = 1;
270   PetscInt       markerLeft     = 1;
271   PetscBool      markerSeparate = PETSC_FALSE;
272   PetscMPIInt    rank;
273   PetscErrorCode ierr;
274 
275   PetscFunctionBegin;
276   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
277   ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr);
278   if (markerSeparate) {
279     markerTop    = 3;
280     markerBottom = 1;
281     markerRight  = 2;
282     markerLeft   = 4;
283   }
284   {
285     const PetscInt numXEdges    = !rank ? edges[0]   : 0;
286     const PetscInt numYEdges    = !rank ? edges[1]   : 0;
287     const PetscInt numXVertices = !rank ? edges[0]+1 : 0;
288     const PetscInt numYVertices = !rank ? edges[1]+1 : 0;
289     const PetscInt numTotXEdges = numXEdges*numYVertices;
290     const PetscInt numTotYEdges = numYEdges*numXVertices;
291     const PetscInt numVertices  = numXVertices*numYVertices;
292     const PetscInt numEdges     = numTotXEdges + numTotYEdges;
293     const PetscInt numFaces     = numXEdges*numYEdges;
294     const PetscInt firstVertex  = numFaces;
295     const PetscInt firstXEdge   = numFaces + numVertices;
296     const PetscInt firstYEdge   = numFaces + numVertices + numTotXEdges;
297     Vec            coordinates;
298     PetscSection   coordSection;
299     PetscScalar    *coords;
300     PetscInt       coordSize;
301     PetscInt       v, vx, vy;
302     PetscInt       f, fx, fy, e, ex, ey;
303 
304     ierr = DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);CHKERRQ(ierr);
305     for (f = 0; f < numFaces; ++f) {
306       ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr);
307     }
308     for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
309       ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr);
310     }
311     ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */
312     /* Build faces */
313     for (fy = 0; fy < numYEdges; fy++) {
314       for (fx = 0; fx < numXEdges; fx++) {
315         const PetscInt face    = fy*numXEdges + fx;
316         const PetscInt edgeL   = firstYEdge + fx*numYEdges + fy;
317         const PetscInt edgeB   = firstXEdge + fy*numXEdges + fx;
318         const PetscInt ornt[4] = {0,     0,               -2,              -2};
319         PetscInt       cone[4];
320 
321         cone[0] = edgeB; cone[1] = edgeL+numYEdges; cone[2] = edgeB+numXEdges; cone[3] = edgeL;
322         ierr    = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr);
323         ierr    = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr);
324       }
325     }
326     /* Build Y edges*/
327     for (vx = 0; vx < numXVertices; vx++) {
328       for (ey = 0; ey < numYEdges; ey++) {
329         const PetscInt edge   = firstYEdge  + vx*numYEdges + ey;
330         const PetscInt vertex = firstVertex + ey*numXVertices + vx;
331         PetscInt       cone[2];
332 
333         cone[0] = vertex; cone[1] = vertex+numXVertices;
334         ierr    = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
335         if (vx == numXVertices-1) {
336           ierr = DMPlexSetLabelValue(dm, "marker", edge,    markerRight);CHKERRQ(ierr);
337           ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr);
338           if (ey == numYEdges-1) {
339             ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr);
340           }
341         } else if (vx == 0) {
342           ierr = DMPlexSetLabelValue(dm, "marker", edge,    markerLeft);CHKERRQ(ierr);
343           ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr);
344           if (ey == numYEdges-1) {
345             ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr);
346           }
347         }
348       }
349     }
350     /* Build X edges*/
351     for (vy = 0; vy < numYVertices; vy++) {
352       for (ex = 0; ex < numXEdges; ex++) {
353         const PetscInt edge   = firstXEdge  + vy*numXEdges + ex;
354         const PetscInt vertex = firstVertex + vy*numXVertices + ex;
355         PetscInt       cone[2];
356 
357         cone[0] = vertex; cone[1] = vertex+1;
358         ierr    = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
359         if (vy == numYVertices-1) {
360           ierr = DMPlexSetLabelValue(dm, "marker", edge,    markerTop);CHKERRQ(ierr);
361           ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr);
362           if (ex == numXEdges-1) {
363             ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr);
364           }
365         } else if (vy == 0) {
366           ierr = DMPlexSetLabelValue(dm, "marker", edge,    markerBottom);CHKERRQ(ierr);
367           ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr);
368           if (ex == numXEdges-1) {
369             ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr);
370           }
371         }
372       }
373     }
374     ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
375     ierr = DMPlexStratify(dm);CHKERRQ(ierr);
376     /* Build coordinates */
377     ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
378     ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr);
379     for (v = firstVertex; v < firstVertex+numVertices; ++v) {
380       ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr);
381     }
382     ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
383     ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
384     ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
385     ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
386     ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
387     ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
388     for (vy = 0; vy < numYVertices; ++vy) {
389       for (vx = 0; vx < numXVertices; ++vx) {
390         coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
391         coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
392       }
393     }
394     ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
395     ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
396     ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
397   }
398   PetscFunctionReturn(0);
399 }
400 
401 #undef __FUNCT__
402 #define __FUNCT__ "DMPlexCreateBoxMesh"
403 /*@
404   DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box).
405 
406   Collective on MPI_Comm
407 
408   Input Parameters:
409 + comm - The communicator for the DM object
410 . dim - The spatial dimension
411 - interpolate - Flag to create intermediate mesh pieces (edges, faces)
412 
413   Output Parameter:
414 . dm  - The DM object
415 
416   Level: beginner
417 
418 .keywords: DM, create
419 .seealso: DMSetType(), DMCreate()
420 @*/
421 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm)
422 {
423   DM             boundary;
424   PetscErrorCode ierr;
425 
426   PetscFunctionBegin;
427   PetscValidPointer(dm, 4);
428   ierr = DMCreate(comm, &boundary);CHKERRQ(ierr);
429   PetscValidLogicalCollectiveInt(boundary,dim,2);
430   ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr);
431   ierr = DMPlexSetDimension(boundary, dim-1);CHKERRQ(ierr);
432   switch (dim) {
433   case 2:
434   {
435     PetscReal lower[2] = {0.0, 0.0};
436     PetscReal upper[2] = {1.0, 1.0};
437     PetscInt  edges[2] = {2, 2};
438 
439     ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr);
440     break;
441   }
442   case 3:
443   {
444     PetscReal lower[3] = {0.0, 0.0, 0.0};
445     PetscReal upper[3] = {1.0, 1.0, 1.0};
446     PetscInt  faces[3] = {1, 1, 1};
447 
448     ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr);
449     break;
450   }
451   default:
452     SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
453   }
454   ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr);
455   ierr = DMDestroy(&boundary);CHKERRQ(ierr);
456   PetscFunctionReturn(0);
457 }
458 
459 #undef __FUNCT__
460 #define __FUNCT__ "DMPlexCreateHexBoxMesh"
461 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm)
462 {
463   PetscErrorCode ierr;
464 
465   PetscFunctionBegin;
466   PetscValidPointer(dm, 4);
467   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
468   PetscValidLogicalCollectiveInt(*dm,dim,2);
469   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
470   ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr);
471   switch (dim) {
472   case 2:
473   {
474     PetscReal lower[2] = {0.0, 0.0};
475     PetscReal upper[2] = {1.0, 1.0};
476 
477     ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells);CHKERRQ(ierr);
478     break;
479   }
480 #if 0
481   case 3:
482   {
483     PetscReal lower[3] = {0.0, 0.0, 0.0};
484     PetscReal upper[3] = {1.0, 1.0, 1.0};
485 
486     ierr = DMPlexCreateCubeMesh(boundary, lower, upper, cells);CHKERRQ(ierr);
487     break;
488   }
489 #endif
490   default:
491     SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
492   }
493   PetscFunctionReturn(0);
494 }
495 
496 /* External function declarations here */
497 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
498 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J);
499 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
500 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined);
501 extern PetscErrorCode DMSetUp_Plex(DM dm);
502 extern PetscErrorCode DMDestroy_Plex(DM dm);
503 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
504 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
505 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS);
506 
507 #undef __FUNCT__
508 #define __FUNCT__ "DMCreateGlobalVector_Plex"
509 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
510 {
511   PetscErrorCode ierr;
512 
513   PetscFunctionBegin;
514   ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr);
515   /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */
516   ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);CHKERRQ(ierr);
517   PetscFunctionReturn(0);
518 }
519 
520 #undef __FUNCT__
521 #define __FUNCT__ "DMCreateLocalVector_Plex"
522 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
523 {
524   PetscErrorCode ierr;
525 
526   PetscFunctionBegin;
527   ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr);
528   ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);CHKERRQ(ierr);
529   PetscFunctionReturn(0);
530 }
531 
532 
533 #undef __FUNCT__
534 #define __FUNCT__ "DMInitialize_Plex"
535 PetscErrorCode DMInitialize_Plex(DM dm)
536 {
537   PetscErrorCode ierr;
538 
539   PetscFunctionBegin;
540   ierr = PetscStrallocpy(VECSTANDARD, (char**)&dm->vectype);CHKERRQ(ierr);
541 
542   dm->ops->view                            = DMView_Plex;
543   dm->ops->setfromoptions                  = DMSetFromOptions_Plex;
544   dm->ops->setup                           = DMSetUp_Plex;
545   dm->ops->createglobalvector              = DMCreateGlobalVector_Plex;
546   dm->ops->createlocalvector               = DMCreateLocalVector_Plex;
547   dm->ops->getlocaltoglobalmapping         = NULL;
548   dm->ops->getlocaltoglobalmappingblock    = NULL;
549   dm->ops->createfieldis                   = NULL;
550   dm->ops->createcoordinatedm              = DMCreateCoordinateDM_Plex;
551   dm->ops->getcoloring                     = 0;
552   dm->ops->creatematrix                    = DMCreateMatrix_Plex;
553   dm->ops->createinterpolation             = 0;
554   dm->ops->getaggregates                   = 0;
555   dm->ops->getinjection                    = 0;
556   dm->ops->refine                          = DMRefine_Plex;
557   dm->ops->coarsen                         = 0;
558   dm->ops->refinehierarchy                 = 0;
559   dm->ops->coarsenhierarchy                = 0;
560   dm->ops->globaltolocalbegin              = NULL;
561   dm->ops->globaltolocalend                = NULL;
562   dm->ops->localtoglobalbegin              = NULL;
563   dm->ops->localtoglobalend                = NULL;
564   dm->ops->destroy                         = DMDestroy_Plex;
565   dm->ops->createsubdm                     = DMCreateSubDM_Plex;
566   dm->ops->locatepoints                    = DMLocatePoints_Plex;
567   PetscFunctionReturn(0);
568 }
569 
570 /*MC
571   DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram.
572                     In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is
573                     specified by a PetscSection object. Ownership in the global representation is determined by
574                     ownership of the underlying DMPlex points. This is specified by another PetscSection object.
575 
576   Level: intermediate
577 
578 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
579 M*/
580 
581 #undef __FUNCT__
582 #define __FUNCT__ "DMCreate_Plex"
583 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
584 {
585   DM_Plex        *mesh;
586   PetscInt       unit, d;
587   PetscErrorCode ierr;
588 
589   PetscFunctionBegin;
590   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
591   ierr     = PetscNewLog(dm, DM_Plex, &mesh);CHKERRQ(ierr);
592   dm->data = mesh;
593 
594   mesh->refct             = 1;
595   mesh->dim               = 0;
596   ierr                    = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr);
597   mesh->maxConeSize       = 0;
598   mesh->cones             = NULL;
599   mesh->coneOrientations  = NULL;
600   ierr                    = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr);
601   mesh->maxSupportSize    = 0;
602   mesh->supports          = NULL;
603   mesh->refinementUniform = PETSC_TRUE;
604   mesh->refinementLimit   = -1.0;
605 
606   mesh->facesTmp = NULL;
607 
608   mesh->subpointMap = NULL;
609 
610   for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0;
611 
612   mesh->labels              = NULL;
613   mesh->globalVertexNumbers = NULL;
614   mesh->globalCellNumbers   = NULL;
615   for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
616   mesh->vtkCellHeight     = 0;
617   mesh->preallocCenterDim = -1;
618 
619   mesh->printSetValues = PETSC_FALSE;
620   mesh->printFEM       = 0;
621   mesh->printTol       = 1.0e-10;
622 
623   ierr = DMInitialize_Plex(dm);CHKERRQ(ierr);
624   PetscFunctionReturn(0);
625 }
626 
627 #undef __FUNCT__
628 #define __FUNCT__ "DMPlexCreate"
629 /*@
630   DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.
631 
632   Collective on MPI_Comm
633 
634   Input Parameter:
635 . comm - The communicator for the DMPlex object
636 
637   Output Parameter:
638 . mesh  - The DMPlex object
639 
640   Level: beginner
641 
642 .keywords: DMPlex, create
643 @*/
644 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
645 {
646   PetscErrorCode ierr;
647 
648   PetscFunctionBegin;
649   PetscValidPointer(mesh,2);
650   ierr = DMCreate(comm, mesh);CHKERRQ(ierr);
651   ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr);
652   PetscFunctionReturn(0);
653 }
654 
655 #undef __FUNCT__
656 #define __FUNCT__ "DMPlexClone"
657 /*@
658   DMPlexClone - Creates a DMPlex object with the same mesh as the original.
659 
660   Collective on MPI_Comm
661 
662   Input Parameter:
663 . dm - The original DMPlex object
664 
665   Output Parameter:
666 . newdm  - The new DMPlex object
667 
668   Level: beginner
669 
670 .keywords: DMPlex, create
671 @*/
672 PetscErrorCode DMPlexClone(DM dm, DM *newdm)
673 {
674   DM_Plex        *mesh;
675   Vec             coords;
676   void           *ctx;
677   PetscErrorCode ierr;
678 
679   PetscFunctionBegin;
680   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
681   PetscValidPointer(newdm,2);
682   ierr         = DMCreate(PetscObjectComm((PetscObject)dm), newdm);CHKERRQ(ierr);
683   ierr         = PetscSFDestroy(&(*newdm)->sf);CHKERRQ(ierr);
684   ierr         = PetscObjectReference((PetscObject) dm->sf);CHKERRQ(ierr);
685   (*newdm)->sf = dm->sf;
686   mesh         = (DM_Plex*) dm->data;
687   mesh->refct++;
688   (*newdm)->data = mesh;
689   ierr           = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr);
690   ierr           = DMInitialize_Plex(*newdm);CHKERRQ(ierr);
691   ierr           = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr);
692   ierr           = DMSetApplicationContext(*newdm, ctx);CHKERRQ(ierr);
693   ierr           = DMGetCoordinatesLocal(dm, &coords);CHKERRQ(ierr);
694   if (coords) {
695     ierr         = DMSetCoordinatesLocal(*newdm, coords);CHKERRQ(ierr);
696   } else {
697     ierr         = DMGetCoordinates(dm, &coords);CHKERRQ(ierr);
698     if (coords) {ierr = DMSetCoordinates(*newdm, coords);CHKERRQ(ierr);}
699   }
700   PetscFunctionReturn(0);
701 }
702 
703 #undef __FUNCT__
704 #define __FUNCT__ "DMPlexBuildFromCellList_Private"
705 /*
706   This takes as input the common mesh generator output, a list of the vertices for each cell
707 */
708 PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[])
709 {
710   PetscInt      *cone, c, p;
711   PetscErrorCode ierr;
712 
713   PetscFunctionBegin;
714   ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr);
715   for (c = 0; c < numCells; ++c) {
716     ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);
717   }
718   ierr = DMSetUp(dm);CHKERRQ(ierr);
719   ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr);
720   for (c = 0; c < numCells; ++c) {
721     for (p = 0; p < numCorners; ++p) {
722       cone[p] = cells[c*numCorners+p]+numCells;
723     }
724     ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr);
725   }
726   ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr);
727   ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
728   ierr = DMPlexStratify(dm);CHKERRQ(ierr);
729   PetscFunctionReturn(0);
730 }
731 
732 #undef __FUNCT__
733 #define __FUNCT__ "DMPlexBuildCoordinates_Private"
734 /*
735   This takes as input the coordinates for each vertex
736 */
737 PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
738 {
739   PetscSection   coordSection;
740   Vec            coordinates;
741   PetscScalar   *coords;
742   PetscInt       coordSize, v, d;
743   PetscErrorCode ierr;
744 
745   PetscFunctionBegin;
746   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
747   ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
748   ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr);
749   ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr);
750   for (v = numCells; v < numCells+numVertices; ++v) {
751     ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr);
752     ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr);
753   }
754   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
755   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
756   ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
757   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
758   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
759   ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
760   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
761   for (v = 0; v < numVertices; ++v) {
762     for (d = 0; d < spaceDim; ++d) {
763       coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
764     }
765   }
766   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
767   ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
768   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
769   PetscFunctionReturn(0);
770 }
771 
772 #undef __FUNCT__
773 #define __FUNCT__ "DMPlexCreateFromCellList"
774 /*@C
775   DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM
776 
777   Input Parameters:
778 + comm - The communicator
779 . dim - The topological dimension of the mesh
780 . numCells - The number of cells
781 . numVertices - The number of vertices
782 . numCorners - The number of vertices for each cell
783 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
784 . cells - An array of numCells*numCorners numbers, the vertices for each cell
785 . spaceDim - The spatial dimension used for coordinates
786 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
787 
788   Output Parameter:
789 . dm - The DM
790 
791   Note: Two triangles sharing a face
792 $
793 $        2
794 $      / | \
795 $     /  |  \
796 $    /   |   \
797 $   0  0 | 1  3
798 $    \   |   /
799 $     \  |  /
800 $      \ | /
801 $        1
802 would have input
803 $  numCells = 2, numVertices = 4
804 $  cells = [0 1 2  1 3 2]
805 $
806 which would result in the DMPlex
807 $
808 $        4
809 $      / | \
810 $     /  |  \
811 $    /   |   \
812 $   2  0 | 1  5
813 $    \   |   /
814 $     \  |  /
815 $      \ | /
816 $        3
817 
818   Level: beginner
819 
820 .seealso: DMPlexCreate()
821 @*/
822 PetscErrorCode DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm)
823 {
824   PetscErrorCode ierr;
825 
826   PetscFunctionBegin;
827   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
828   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
829   ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr);
830   ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr);
831   if (interpolate) {
832     DM idm;
833 
834     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
835     ierr = DMDestroy(dm);CHKERRQ(ierr);
836     *dm  = idm;
837   }
838   ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr);
839   PetscFunctionReturn(0);
840 }
841 
842 #undef __FUNCT__
843 #define __FUNCT__ "DMPlexCreateFromDAG"
844 /*
845   This takes as input the raw Hasse Diagram data
846 */
847 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
848 {
849   Vec            coordinates;
850   PetscSection   coordSection;
851   PetscScalar    *coords;
852   PetscInt       coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off;
853   PetscErrorCode ierr;
854 
855   PetscFunctionBegin;
856   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
857   for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
858   ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr);
859   for (p = pStart; p < pEnd; ++p) {
860     ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr);
861   }
862   ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */
863   for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
864     ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr);
865     ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr);
866   }
867   ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
868   ierr = DMPlexStratify(dm);CHKERRQ(ierr);
869   /* Build coordinates */
870   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
871   ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
872   ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
873   ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr);
874   for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
875     ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
876     ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
877   }
878   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
879   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
880   ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
881   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
882   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
883   ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
884   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
885   for (v = 0; v < numPoints[0]; ++v) {
886     PetscInt off;
887 
888     ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr);
889     for (d = 0; d < dim; ++d) {
890       coords[off+d] = vertexCoords[v*dim+d];
891     }
892   }
893   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
894   ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
895   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
896   PetscFunctionReturn(0);
897 }
898