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