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