Actual source code: plexcreate.c

petsc-dev 2014-08-26
Report Typos and Errors
  1: #define PETSCDM_DLL
  2: #include <petsc-private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
  3: #include <petscdmda.h>
  4: #include <petscsf.h>

  8: /*@
  9:   DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement.

 11:   Collective on MPI_Comm

 13:   Input Parameters:
 14: + comm - The communicator for the DM object
 15: . dim - The spatial dimension
 16: . simplex - Flag for simplicial cells, otherwise they are tensor product cells
 17: . interpolate - Flag to create intermediate mesh pieces (edges, faces)
 18: . refinementUniform - Flag for uniform parallel refinement
 19: - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell

 21:   Output Parameter:
 22: . dm  - The DM object

 24:   Level: beginner

 26: .keywords: DM, create
 27: .seealso: DMSetType(), DMCreate()
 28: @*/
 29: PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm)
 30: {
 31:   DM             dm;
 32:   PetscInt       p;
 33:   PetscMPIInt    rank;

 37:   DMCreate(comm, &dm);
 38:   DMSetType(dm, DMPLEX);
 39:   DMSetDimension(dm, dim);
 40:   MPI_Comm_rank(comm, &rank);
 41:   switch (dim) {
 42:   case 2:
 43:     if (simplex) {PetscObjectSetName((PetscObject) dm, "triangular");}
 44:     else         {PetscObjectSetName((PetscObject) dm, "quadrilateral");}
 45:     break;
 46:   case 3:
 47:     if (simplex) {PetscObjectSetName((PetscObject) dm, "tetrahedral");}
 48:     else         {PetscObjectSetName((PetscObject) dm, "hexahedral");}
 49:     break;
 50:   default:
 51:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
 52:   }
 53:   if (rank) {
 54:     PetscInt numPoints[2] = {0, 0};
 55:     DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);
 56:   } else {
 57:     switch (dim) {
 58:     case 2:
 59:       if (simplex) {
 60:         PetscInt    numPoints[2]        = {4, 2};
 61:         PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
 62:         PetscInt    cones[6]            = {2, 3, 4,  5, 4, 3};
 63:         PetscInt    coneOrientations[6] = {0, 0, 0,  0, 0, 0};
 64:         PetscScalar vertexCoords[8]     = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
 65:         PetscInt    markerPoints[8]     = {2, 1, 3, 1, 4, 1, 5, 1};

 67:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
 68:         for (p = 0; p < 4; ++p) {DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);}
 69:       } else {
 70:         PetscInt    numPoints[2]        = {6, 2};
 71:         PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
 72:         PetscInt    cones[8]            = {2, 3, 4, 5,  3, 6, 7, 4};
 73:         PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
 74:         PetscScalar vertexCoords[12]    = {-1.0, -0.5,  0.0, -0.5,  0.0, 0.5,  -1.0, 0.5,  1.0, -0.5,  1.0, 0.5};

 76:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
 77:       }
 78:       break;
 79:     case 3:
 80:       if (simplex) {
 81:         PetscInt    numPoints[2]        = {5, 2};
 82:         PetscInt    coneSize[7]         = {4, 4, 0, 0, 0, 0, 0};
 83:         PetscInt    cones[8]            = {4, 3, 5, 2,  5, 3, 4, 6};
 84:         PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
 85:         PetscScalar vertexCoords[15]    = {-1.0, 0.0, 0.0,  0.0, -1.0, 0.0,  0.0, 0.0, 1.0,  0.0, 1.0, 0.0,  1.0, 0.0, 0.0};
 86:         PetscInt    markerPoints[10]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};

 88:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
 89:         for (p = 0; p < 5; ++p) {DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);}
 90:       } else {
 91:         PetscInt    numPoints[2]         = {12, 2};
 92:         PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 93:         PetscInt    cones[16]            = {2, 3, 4, 5, 6, 7, 8, 9,  5, 4, 10, 11, 7, 12, 13, 8};
 94:         PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0,  0, 0,  0,  0, 0,  0,  0, 0};
 95:         PetscScalar vertexCoords[36]     = {-1.0, -0.5, -0.5,  -1.0,  0.5, -0.5,  0.0,  0.5, -0.5,   0.0, -0.5, -0.5,
 96:                                             -1.0, -0.5,  0.5,   0.0, -0.5,  0.5,  0.0,  0.5,  0.5,  -1.0,  0.5,  0.5,
 97:                                              1.0,  0.5, -0.5,   1.0, -0.5, -0.5,  1.0, -0.5,  0.5,   1.0,  0.5,  0.5};

 99:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
100:       }
101:       break;
102:     default:
103:       SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
104:     }
105:   }
106:   *newdm = dm;
107:   if (refinementLimit > 0.0) {
108:     DM rdm;
109:     const char *name;

111:     DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);
112:     DMPlexSetRefinementLimit(*newdm, refinementLimit);
113:     DMRefine(*newdm, comm, &rdm);
114:     PetscObjectGetName((PetscObject) *newdm, &name);
115:     PetscObjectSetName((PetscObject)    rdm,  name);
116:     DMDestroy(newdm);
117:     *newdm = rdm;
118:   }
119:   if (interpolate) {
120:     DM idm = NULL;
121:     const char *name;

123:     DMPlexInterpolate(*newdm, &idm);
124:     PetscObjectGetName((PetscObject) *newdm, &name);
125:     PetscObjectSetName((PetscObject)    idm,  name);
126:     DMPlexCopyCoordinates(*newdm, idm);
127:     DMPlexCopyLabels(*newdm, idm);
128:     DMDestroy(newdm);
129:     *newdm = idm;
130:   }
131:   {
132:     DM refinedMesh     = NULL;
133:     DM distributedMesh = NULL;

135:     /* Distribute mesh over processes */
136:     DMPlexDistribute(*newdm, NULL, 0, NULL, &distributedMesh);
137:     if (distributedMesh) {
138:       DMDestroy(newdm);
139:       *newdm = distributedMesh;
140:     }
141:     if (refinementUniform) {
142:       DMPlexSetRefinementUniform(*newdm, refinementUniform);
143:       DMRefine(*newdm, comm, &refinedMesh);
144:       if (refinedMesh) {
145:         DMDestroy(newdm);
146:         *newdm = refinedMesh;
147:       }
148:     }
149:   }
150:   return(0);
151: }

155: /*@
156:   DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice.

158:   Collective on MPI_Comm

160:   Input Parameters:
161: + comm  - The communicator for the DM object
162: . lower - The lower left corner coordinates
163: . upper - The upper right corner coordinates
164: - edges - The number of cells in each direction

166:   Output Parameter:
167: . dm  - The DM object

169:   Note: Here is the numbering returned for 2 cells in each direction:
170: $ 18--5-17--4--16
171: $  |     |     |
172: $  6    10     3
173: $  |     |     |
174: $ 19-11-20--9--15
175: $  |     |     |
176: $  7     8     2
177: $  |     |     |
178: $ 12--0-13--1--14

180:   Level: beginner

182: .keywords: DM, create
183: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate()
184: @*/
185: PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
186: {
187:   PetscInt       numVertices    = (edges[0]+1)*(edges[1]+1);
188:   PetscInt       numEdges       = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
189:   PetscInt       markerTop      = 1;
190:   PetscInt       markerBottom   = 1;
191:   PetscInt       markerRight    = 1;
192:   PetscInt       markerLeft     = 1;
193:   PetscBool      markerSeparate = PETSC_FALSE;
194:   Vec            coordinates;
195:   PetscSection   coordSection;
196:   PetscScalar    *coords;
197:   PetscInt       coordSize;
198:   PetscMPIInt    rank;
199:   PetscInt       v, vx, vy;

203:   PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
204:   if (markerSeparate) {
205:     markerTop    = 1;
206:     markerBottom = 0;
207:     markerRight  = 0;
208:     markerLeft   = 0;
209:   }
210:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
211:   if (!rank) {
212:     PetscInt e, ex, ey;

214:     DMPlexSetChart(dm, 0, numEdges+numVertices);
215:     for (e = 0; e < numEdges; ++e) {
216:       DMPlexSetConeSize(dm, e, 2);
217:     }
218:     DMSetUp(dm); /* Allocate space for cones */
219:     for (vx = 0; vx <= edges[0]; vx++) {
220:       for (ey = 0; ey < edges[1]; ey++) {
221:         PetscInt edge   = vx*edges[1] + ey + edges[0]*(edges[1]+1);
222:         PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
223:         PetscInt cone[2];

225:         cone[0] = vertex; cone[1] = vertex+edges[0]+1;
226:         DMPlexSetCone(dm, edge, cone);
227:         if (vx == edges[0]) {
228:           DMPlexSetLabelValue(dm, "marker", edge,    markerRight);
229:           DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
230:           if (ey == edges[1]-1) {
231:             DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
232:           }
233:         } else if (vx == 0) {
234:           DMPlexSetLabelValue(dm, "marker", edge,    markerLeft);
235:           DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
236:           if (ey == edges[1]-1) {
237:             DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
238:           }
239:         }
240:       }
241:     }
242:     for (vy = 0; vy <= edges[1]; vy++) {
243:       for (ex = 0; ex < edges[0]; ex++) {
244:         PetscInt edge   = vy*edges[0]     + ex;
245:         PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
246:         PetscInt cone[2];

248:         cone[0] = vertex; cone[1] = vertex+1;
249:         DMPlexSetCone(dm, edge, cone);
250:         if (vy == edges[1]) {
251:           DMPlexSetLabelValue(dm, "marker", edge,    markerTop);
252:           DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
253:           if (ex == edges[0]-1) {
254:             DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
255:           }
256:         } else if (vy == 0) {
257:           DMPlexSetLabelValue(dm, "marker", edge,    markerBottom);
258:           DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
259:           if (ex == edges[0]-1) {
260:             DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
261:           }
262:         }
263:       }
264:     }
265:   }
266:   DMPlexSymmetrize(dm);
267:   DMPlexStratify(dm);
268:   /* Build coordinates */
269:   DMGetCoordinateSection(dm, &coordSection);
270:   PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);
271:   for (v = numEdges; v < numEdges+numVertices; ++v) {
272:     PetscSectionSetDof(coordSection, v, 2);
273:   }
274:   PetscSectionSetUp(coordSection);
275:   PetscSectionGetStorageSize(coordSection, &coordSize);
276:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
277:   VecSetBlockSize(coordinates, 2);
278:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
279:   VecSetType(coordinates,VECSTANDARD);
280:   VecGetArray(coordinates, &coords);
281:   for (vy = 0; vy <= edges[1]; ++vy) {
282:     for (vx = 0; vx <= edges[0]; ++vx) {
283:       coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
284:       coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
285:     }
286:   }
287:   VecRestoreArray(coordinates, &coords);
288:   DMSetCoordinatesLocal(dm, coordinates);
289:   VecDestroy(&coordinates);
290:   return(0);
291: }

295: /*@
296:   DMPlexCreateCubeBoundary - Creates a 2D mesh the is the boundary of a cubic lattice.

298:   Collective on MPI_Comm

300:   Input Parameters:
301: + comm  - The communicator for the DM object
302: . lower - The lower left front corner coordinates
303: . upper - The upper right back corner coordinates
304: - edges - The number of cells in each direction

306:   Output Parameter:
307: . dm  - The DM object

309:   Level: beginner

311: .keywords: DM, create
312: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate()
313: @*/
314: PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
315: {
316:   PetscInt       vertices[3], numVertices;
317:   PetscInt       numFaces    = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2];
318:   Vec            coordinates;
319:   PetscSection   coordSection;
320:   PetscScalar    *coords;
321:   PetscInt       coordSize;
322:   PetscMPIInt    rank;
323:   PetscInt       v, vx, vy, vz;
324:   PetscInt       voffset, iface=0, cone[4];

328:   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");
329:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
330:   vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1;
331:   numVertices = vertices[0]*vertices[1]*vertices[2];
332:   if (!rank) {
333:     PetscInt f;

335:     DMPlexSetChart(dm, 0, numFaces+numVertices);
336:     for (f = 0; f < numFaces; ++f) {
337:       DMPlexSetConeSize(dm, f, 4);
338:     }
339:     DMSetUp(dm); /* Allocate space for cones */
340:     for (v = 0; v < numFaces+numVertices; ++v) {
341:       DMPlexSetLabelValue(dm, "marker", v, 1);
342:     }

344:     /* Side 0 (Top) */
345:     for (vy = 0; vy < faces[1]; vy++) {
346:       for (vx = 0; vx < faces[0]; vx++) {
347:         voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx;
348:         cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0];
349:         DMPlexSetCone(dm, iface, cone);
350:         iface++;
351:       }
352:     }

354:     /* Side 1 (Bottom) */
355:     for (vy = 0; vy < faces[1]; vy++) {
356:       for (vx = 0; vx < faces[0]; vx++) {
357:         voffset = numFaces + vy*(faces[0]+1) + vx;
358:         cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1;
359:         DMPlexSetCone(dm, iface, cone);
360:         iface++;
361:       }
362:     }

364:     /* Side 2 (Front) */
365:     for (vz = 0; vz < faces[2]; vz++) {
366:       for (vx = 0; vx < faces[0]; vx++) {
367:         voffset = numFaces + vz*vertices[0]*vertices[1] + vx;
368:         cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1];
369:         DMPlexSetCone(dm, iface, cone);
370:         iface++;
371:       }
372:     }

374:     /* Side 3 (Back) */
375:     for (vz = 0; vz < faces[2]; vz++) {
376:       for (vx = 0; vx < faces[0]; vx++) {
377:         voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx;
378:         cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1;
379:         cone[2] = voffset+1; cone[3] = voffset;
380:         DMPlexSetCone(dm, iface, cone);
381:         iface++;
382:       }
383:     }

385:     /* Side 4 (Left) */
386:     for (vz = 0; vz < faces[2]; vz++) {
387:       for (vy = 0; vy < faces[1]; vy++) {
388:         voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0];
389:         cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1];
390:         cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0];
391:         DMPlexSetCone(dm, iface, cone);
392:         iface++;
393:       }
394:     }

396:     /* Side 5 (Right) */
397:     for (vz = 0; vz < faces[2]; vz++) {
398:       for (vy = 0; vy < faces[1]; vy++) {
399:         voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + vx;
400:         cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset;
401:         cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0];
402:         DMPlexSetCone(dm, iface, cone);
403:         iface++;
404:       }
405:     }
406:   }
407:   DMPlexSymmetrize(dm);
408:   DMPlexStratify(dm);
409:   /* Build coordinates */
410:   DMGetCoordinateSection(dm, &coordSection);
411:   PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);
412:   for (v = numFaces; v < numFaces+numVertices; ++v) {
413:     PetscSectionSetDof(coordSection, v, 3);
414:   }
415:   PetscSectionSetUp(coordSection);
416:   PetscSectionGetStorageSize(coordSection, &coordSize);
417:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
418:   VecSetBlockSize(coordinates, 3);
419:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
420:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
421:   VecSetType(coordinates,VECSTANDARD);
422:   VecGetArray(coordinates, &coords);
423:   for (vz = 0; vz <= faces[2]; ++vz) {
424:     for (vy = 0; vy <= faces[1]; ++vy) {
425:       for (vx = 0; vx <= faces[0]; ++vx) {
426:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
427:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
428:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
429:       }
430:     }
431:   }
432:   VecRestoreArray(coordinates, &coords);
433:   DMSetCoordinatesLocal(dm, coordinates);
434:   VecDestroy(&coordinates);
435:   return(0);
436: }

440: /*@
441:   DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice.

443:   Collective on MPI_Comm

445:   Input Parameters:
446: + comm  - The communicator for the DM object
447: . lower - The lower left corner coordinates
448: . upper - The upper right corner coordinates
449: . edges - The number of cells in each direction
450: . bdX   - The boundary type for the X direction
451: - bdY   - The boundary type for the Y direction

453:   Output Parameter:
454: . dm  - The DM object

456:   Note: Here is the numbering returned for 2 cells in each direction:
457: $ 22--8-23--9--24
458: $  |     |     |
459: $ 13  2 14  3  15
460: $  |     |     |
461: $ 19--6-20--7--21
462: $  |     |     |
463: $ 10  0 11  1 12
464: $  |     |     |
465: $ 16--4-17--5--18

467:   Level: beginner

469: .keywords: DM, create
470: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate()
471: @*/
472: PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY)
473: {
474:   PetscInt       markerTop      = 1;
475:   PetscInt       markerBottom   = 1;
476:   PetscInt       markerRight    = 1;
477:   PetscInt       markerLeft     = 1;
478:   PetscBool      markerSeparate = PETSC_FALSE;
479:   PetscMPIInt    rank;

483:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
484:   PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
485:   if (markerSeparate) {
486:     markerTop    = 3;
487:     markerBottom = 1;
488:     markerRight  = 2;
489:     markerLeft   = 4;
490:   }
491:   {
492:     const PetscInt numXEdges    = !rank ? edges[0]   : 0;
493:     const PetscInt numYEdges    = !rank ? edges[1]   : 0;
494:     const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0;
495:     const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0;
496:     const PetscInt numTotXEdges = numXEdges*numYVertices;
497:     const PetscInt numTotYEdges = numYEdges*numXVertices;
498:     const PetscInt numVertices  = numXVertices*numYVertices;
499:     const PetscInt numEdges     = numTotXEdges + numTotYEdges;
500:     const PetscInt numFaces     = numXEdges*numYEdges;
501:     const PetscInt firstVertex  = numFaces;
502:     const PetscInt firstXEdge   = numFaces + numVertices;
503:     const PetscInt firstYEdge   = numFaces + numVertices + numTotXEdges;
504:     Vec            coordinates;
505:     PetscSection   coordSection;
506:     PetscScalar   *coords;
507:     PetscInt       coordSize;
508:     PetscInt       v, vx, vy;
509:     PetscInt       f, fx, fy, e, ex, ey;

511:     DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);
512:     for (f = 0; f < numFaces; ++f) {
513:       DMPlexSetConeSize(dm, f, 4);
514:     }
515:     for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
516:       DMPlexSetConeSize(dm, e, 2);
517:     }
518:     DMSetUp(dm); /* Allocate space for cones */
519:     /* Build faces */
520:     for (fy = 0; fy < numYEdges; ++fy) {
521:       for (fx = 0; fx < numXEdges; ++fx) {
522:         PetscInt face    = fy*numXEdges + fx;
523:         PetscInt edgeL   = firstYEdge +   fx                 *numYEdges + fy;
524:         PetscInt edgeR   = firstYEdge + ((fx+1)%numXVertices)*numYEdges + fy;
525:         PetscInt edgeB   = firstXEdge +   fy                 *numXEdges + fx;
526:         PetscInt edgeT   = firstXEdge + ((fy+1)%numYVertices)*numXEdges + fx;
527:         PetscInt ornt[4] = {0, 0, -2, -2};
528:         PetscInt cone[4];

530:         if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;}
531:         if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] =  0;}
532:         cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
533:         DMPlexSetCone(dm, face, cone);
534:         DMPlexSetConeOrientation(dm, face, ornt);
535:       }
536:     }
537:     /* Build Y edges*/
538:     for (vx = 0; vx < numXVertices; vx++) {
539:       for (ey = 0; ey < numYEdges; ey++) {
540:         const PetscInt nextv   = (bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : ((ey+1)%numYVertices)*numXVertices + vx;
541:         const PetscInt edge    = firstYEdge  + vx*numYEdges + ey;
542:         const PetscInt vertexB = firstVertex + ey*numXVertices + vx;
543:         const PetscInt vertexT = firstVertex + nextv;
544:         PetscInt       cone[2];

546:         cone[0] = vertexB; cone[1] = vertexT;
547:         DMPlexSetCone(dm, edge, cone);
548:         if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) {
549:           if (vx == numXVertices-1) {
550:             DMPlexSetLabelValue(dm, "marker", edge,    markerRight);
551:             DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
552:             if (ey == numYEdges-1) {
553:               DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
554:             }
555:           } else if (vx == 0) {
556:             DMPlexSetLabelValue(dm, "marker", edge,    markerLeft);
557:             DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
558:             if (ey == numYEdges-1) {
559:               DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
560:             }
561:           }
562:         }
563:       }
564:     }
565:     /* Build X edges*/
566:     for (vy = 0; vy < numYVertices; vy++) {
567:       for (ex = 0; ex < numXEdges; ex++) {
568:         const PetscInt nextv   = (bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : vy*numXVertices + (ex+1)%numXVertices;
569:         const PetscInt edge    = firstXEdge  + vy*numXEdges + ex;
570:         const PetscInt vertexL = firstVertex + vy*numXVertices + ex;
571:         const PetscInt vertexR = firstVertex + nextv;
572:         PetscInt       cone[2];

574:         cone[0] = vertexL; cone[1] = vertexR;
575:         DMPlexSetCone(dm, edge, cone);
576:         if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) {
577:           if (vy == numYVertices-1) {
578:             DMPlexSetLabelValue(dm, "marker", edge,    markerTop);
579:             DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
580:             if (ex == numXEdges-1) {
581:               DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
582:             }
583:           } else if (vy == 0) {
584:             DMPlexSetLabelValue(dm, "marker", edge,    markerBottom);
585:             DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
586:             if (ex == numXEdges-1) {
587:               DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
588:             }
589:           }
590:         }
591:       }
592:     }
593:     DMPlexSymmetrize(dm);
594:     DMPlexStratify(dm);
595:     /* Build coordinates */
596:     DMGetCoordinateSection(dm, &coordSection);
597:     PetscSectionSetNumFields(coordSection, 1);
598:     PetscSectionSetFieldComponents(coordSection, 0, 2);
599:     PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);
600:     for (v = firstVertex; v < firstVertex+numVertices; ++v) {
601:       PetscSectionSetDof(coordSection, v, 2);
602:       PetscSectionSetFieldDof(coordSection, v, 0, 2);
603:     }
604:     PetscSectionSetUp(coordSection);
605:     PetscSectionGetStorageSize(coordSection, &coordSize);
606:     VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
607:     VecSetBlockSize(coordinates, 2);
608:     VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
609:     VecSetType(coordinates,VECSTANDARD);
610:     VecGetArray(coordinates, &coords);
611:     for (vy = 0; vy < numYVertices; ++vy) {
612:       for (vx = 0; vx < numXVertices; ++vx) {
613:         coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
614:         coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
615:       }
616:     }
617:     VecRestoreArray(coordinates, &coords);
618:     DMSetCoordinatesLocal(dm, coordinates);
619:     VecDestroy(&coordinates);
620:   }
621:   return(0);
622: }

626: /*@
627:   DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices.

629:   Collective on MPI_Comm

631:   Input Parameters:
632: + comm - The communicator for the DM object
633: . dim - The spatial dimension
634: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

636:   Output Parameter:
637: . dm  - The DM object

639:   Level: beginner

641: .keywords: DM, create
642: .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate()
643: @*/
644: PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm)
645: {
646:   DM             boundary;

651:   DMCreate(comm, &boundary);
653:   DMSetType(boundary, DMPLEX);
654:   DMSetDimension(boundary, dim-1);
655:   switch (dim) {
656:   case 2:
657:   {
658:     PetscReal lower[2] = {0.0, 0.0};
659:     PetscReal upper[2] = {1.0, 1.0};
660:     PetscInt  edges[2] = {2, 2};

662:     DMPlexCreateSquareBoundary(boundary, lower, upper, edges);
663:     break;
664:   }
665:   case 3:
666:   {
667:     PetscReal lower[3] = {0.0, 0.0, 0.0};
668:     PetscReal upper[3] = {1.0, 1.0, 1.0};
669:     PetscInt  faces[3] = {1, 1, 1};

671:     DMPlexCreateCubeBoundary(boundary, lower, upper, faces);
672:     break;
673:   }
674:   default:
675:     SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
676:   }
677:   DMPlexGenerate(boundary, NULL, interpolate, dm);
678:   DMDestroy(&boundary);
679:   return(0);
680: }

684: /*@
685:   DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra.

687:   Collective on MPI_Comm

689:   Input Parameters:
690: + comm  - The communicator for the DM object
691: . dim   - The spatial dimension
692: . periodicX - The boundary type for the X direction
693: . periodicY - The boundary type for the Y direction
694: . periodicZ - The boundary type for the Z direction
695: - cells - The number of cells in each direction

697:   Output Parameter:
698: . dm  - The DM object

700:   Level: beginner

702: .keywords: DM, create
703: .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
704: @*/
705: PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm)
706: {

711:   DMCreate(comm, dm);
713:   DMSetType(*dm, DMPLEX);
714:   DMSetDimension(*dm, dim);
715:   switch (dim) {
716:   case 2:
717:   {
718:     PetscReal lower[2] = {0.0, 0.0};
719:     PetscReal upper[2] = {1.0, 1.0};

721:     DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);
722:     break;
723:   }
724: #if 0
725:   case 3:
726:   {
727:     PetscReal lower[3] = {0.0, 0.0, 0.0};
728:     PetscReal upper[3] = {1.0, 1.0, 1.0};

730:     DMPlexCreateCubeMesh(boundary, lower, upper, cells, periodicX, periodicY, periodicZ);
731:     break;
732:   }
733: #endif
734:   default:
735:     SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
736:   }
737:   return(0);
738: }

740: /* External function declarations here */
741: extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
742: extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, VecScatter *ctx);
743: extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm);
744: extern PetscErrorCode DMCreateMatrix_Plex(DM dm,  Mat *J);
745: extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
746: extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined);
747: extern PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened);
748: extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]);
749: extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
750: extern PetscErrorCode DMSetUp_Plex(DM dm);
751: extern PetscErrorCode DMDestroy_Plex(DM dm);
752: extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
753: extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
754: extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
755: extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS);

759: /* Replace dm with the contents of dmNew
760:    - Share the DM_Plex structure
761:    - Share the coordinates
762:    - Share the SF
763: */
764: static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew)
765: {
766:   PetscSF          sf;
767:   DM               coordDM;
768:   Vec              coords;
769:   const PetscReal *maxCell, *L;
770:   PetscErrorCode   ierr;

773:   DMGetPointSF(dmNew, &sf);
774:   DMSetPointSF(dm, sf);
775:   DMGetCoordinateDM(dmNew, &coordDM);
776:   DMGetCoordinatesLocal(dmNew, &coords);
777:   DMSetCoordinateDM(dm, coordDM);
778:   DMSetCoordinatesLocal(dm, coords);
779:   DMGetPeriodicity(dm, &maxCell, &L);
780:   if (L) {DMSetPeriodicity(dmNew, maxCell, L);}
781:   DMDestroy_Plex(dm);
782:   dm->data = dmNew->data;
783:   ((DM_Plex *) dmNew->data)->refct++;
784:   return(0);
785: }

789: /* Swap dm with the contents of dmNew
790:    - Swap the DM_Plex structure
791:    - Swap the coordinates
792: */
793: static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB)
794: {
795:   DM             coordDMA, coordDMB;
796:   Vec            coordsA,  coordsB;
797:   void          *tmp;

801:   DMGetCoordinateDM(dmA, &coordDMA);
802:   DMGetCoordinateDM(dmB, &coordDMB);
803:   PetscObjectReference((PetscObject) coordDMA);
804:   DMSetCoordinateDM(dmA, coordDMB);
805:   DMSetCoordinateDM(dmB, coordDMA);
806:   PetscObjectDereference((PetscObject) coordDMA);

808:   DMGetCoordinatesLocal(dmA, &coordsA);
809:   DMGetCoordinatesLocal(dmB, &coordsB);
810:   PetscObjectReference((PetscObject) coordsA);
811:   DMSetCoordinatesLocal(dmA, coordsB);
812:   DMSetCoordinatesLocal(dmB, coordsA);
813:   PetscObjectDereference((PetscObject) coordsA);
814:   tmp       = dmA->data;
815:   dmA->data = dmB->data;
816:   dmB->data = tmp;
817:   return(0);
818: }

822: PetscErrorCode  DMSetFromOptions_NonRefinement_Plex(DM dm)
823: {
824:   DM_Plex       *mesh = (DM_Plex*) dm->data;
825:   DMBoundary     b;

829:   /* Handle boundary conditions */
830:   PetscOptionsBegin(PetscObjectComm((PetscObject) dm), NULL, "Boundary condition options", "");
831:   for (b = mesh->boundary; b; b = b->next) {
832:     char      optname[1024];
833:     PetscInt  ids[1024], len = 1024, i;
834:     PetscBool flg;

836:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
837:     PetscMemzero(ids, sizeof(ids));
838:     PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
839:     if (flg) {
840:       DMLabel label;

842:       DMPlexGetLabel(dm, b->labelname, &label);
843:       for (i = 0; i < len; ++i) {
844:         PetscBool has;

846:         DMLabelHasValue(label, ids[i], &has);
847:         if (!has) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Boundary id %D is not present in the label %s", ids[i], b->name);
848:       }
849:       b->numids = len;
850:       PetscFree(b->ids);
851:       PetscMalloc1(len, &b->ids);
852:       PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));
853:     }
854:   }
855:   PetscOptionsEnd();
856:   /* Handle viewing */
857:   PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);
858:   PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);
859:   PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", PETSC_FALSE, &mesh->printTol, NULL);
860:   return(0);
861: }

865: PetscErrorCode  DMSetFromOptions_Plex(DM dm)
866: {
867:   PetscInt       refine = 0, r;
868:   PetscBool      isHierarchy;

873:   PetscOptionsHead("DMPlex Options");
874:   /* Handle DMPlex refinement */
875:   PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);
876:   PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);
877:   if (refine) {DMPlexSetRefinementUniform(dm, PETSC_TRUE);}
878:   if (refine && isHierarchy) {
879:     DM *dms;

881:     PetscMalloc1(refine,&dms);
882:     DMRefineHierarchy(dm, refine, dms);
883:     /* Total hack since we do not pass in a pointer */
884:     DMPlexSwap_Static(dm, dms[refine-1]);
885:     if (refine == 1) {
886:       DMPlexSetCoarseDM(dm, dms[0]);
887:     } else {
888:       DMPlexSetCoarseDM(dm, dms[refine-2]);
889:       DMPlexSetCoarseDM(dms[0], dms[refine-1]);
890:     }
891:     /* Free DMs */
892:     for (r = 0; r < refine; ++r) {
893:       DMSetFromOptions_NonRefinement_Plex(dm);
894:       DMDestroy(&dms[r]);
895:     }
896:     PetscFree(dms);
897:   } else {
898:     for (r = 0; r < refine; ++r) {
899:       DM refinedMesh;

901:       DMSetFromOptions_NonRefinement_Plex(dm);
902:       DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);
903:       /* Total hack since we do not pass in a pointer */
904:       DMPlexReplace_Static(dm, refinedMesh);
905:       DMDestroy(&refinedMesh);
906:     }
907:   }
908:   DMSetFromOptions_NonRefinement_Plex(dm);
909:   PetscOptionsTail();
910:   return(0);
911: }

915: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
916: {

920:   DMCreateGlobalVector_Section_Private(dm,vec);
921:   /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
922:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);
923:   VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);
924:   return(0);
925: }

929: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
930: {

934:   DMCreateLocalVector_Section_Private(dm,vec);
935:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);
936:   VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);
937:   return(0);
938: }

942: static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
943: {
944:   PetscInt       depth, d;

948:   DMPlexGetDepth(dm, &depth);
949:   if (depth == 1) {
950:     DMGetDimension(dm, &d);
951:     if (dim == 0)      {DMPlexGetDepthStratum(dm, dim, pStart, pEnd);}
952:     else if (dim == d) {DMPlexGetDepthStratum(dm, 1, pStart, pEnd);}
953:     else               {*pStart = 0; *pEnd = 0;}
954:   } else {
955:     DMPlexGetDepthStratum(dm, dim, pStart, pEnd);
956:   }
957:   return(0);
958: }

962: PetscErrorCode DMInitialize_Plex(DM dm)
963: {
965:   dm->ops->view                            = DMView_Plex;
966:   dm->ops->load                            = DMLoad_Plex;
967:   dm->ops->setfromoptions                  = DMSetFromOptions_Plex;
968:   dm->ops->clone                           = DMClone_Plex;
969:   dm->ops->setup                           = DMSetUp_Plex;
970:   dm->ops->createdefaultsection            = DMCreateDefaultSection_Plex;
971:   dm->ops->createglobalvector              = DMCreateGlobalVector_Plex;
972:   dm->ops->createlocalvector               = DMCreateLocalVector_Plex;
973:   dm->ops->getlocaltoglobalmapping         = NULL;
974:   dm->ops->createfieldis                   = NULL;
975:   dm->ops->createcoordinatedm              = DMCreateCoordinateDM_Plex;
976:   dm->ops->getcoloring                     = NULL;
977:   dm->ops->creatematrix                    = DMCreateMatrix_Plex;
978:   dm->ops->createinterpolation             = DMCreateInterpolation_Plex;
979:   dm->ops->getaggregates                   = NULL;
980:   dm->ops->getinjection                    = DMCreateInjection_Plex;
981:   dm->ops->refine                          = DMRefine_Plex;
982:   dm->ops->coarsen                         = DMCoarsen_Plex;
983:   dm->ops->refinehierarchy                 = DMRefineHierarchy_Plex;
984:   dm->ops->coarsenhierarchy                = NULL;
985:   dm->ops->globaltolocalbegin              = NULL;
986:   dm->ops->globaltolocalend                = NULL;
987:   dm->ops->localtoglobalbegin              = NULL;
988:   dm->ops->localtoglobalend                = NULL;
989:   dm->ops->destroy                         = DMDestroy_Plex;
990:   dm->ops->createsubdm                     = DMCreateSubDM_Plex;
991:   dm->ops->getdimpoints                    = DMGetDimPoints_Plex;
992:   dm->ops->locatepoints                    = DMLocatePoints_Plex;
993:   return(0);
994: }

998: PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
999: {
1000:   DM_Plex        *mesh = (DM_Plex *) dm->data;

1004:   mesh->refct++;
1005:   (*newdm)->data = mesh;
1006:   PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
1007:   DMInitialize_Plex(*newdm);
1008:   return(0);
1009: }

1011: /*MC
1012:   DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram.
1013:                     In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is
1014:                     specified by a PetscSection object. Ownership in the global representation is determined by
1015:                     ownership of the underlying DMPlex points. This is specified by another PetscSection object.

1017:   Level: intermediate

1019: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
1020: M*/

1024: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
1025: {
1026:   DM_Plex        *mesh;
1027:   PetscInt       unit, d;

1032:   PetscNewLog(dm,&mesh);
1033:   dm->dim  = 0;
1034:   dm->data = mesh;

1036:   mesh->refct             = 1;
1037:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
1038:   mesh->maxConeSize       = 0;
1039:   mesh->cones             = NULL;
1040:   mesh->coneOrientations  = NULL;
1041:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
1042:   mesh->maxSupportSize    = 0;
1043:   mesh->supports          = NULL;
1044:   mesh->refinementUniform = PETSC_TRUE;
1045:   mesh->refinementLimit   = -1.0;

1047:   mesh->facesTmp = NULL;

1049:   mesh->subpointMap = NULL;

1051:   for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0;

1053:   mesh->labels              = NULL;
1054:   mesh->depthLabel          = NULL;
1055:   mesh->globalVertexNumbers = NULL;
1056:   mesh->globalCellNumbers   = NULL;
1057:   for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
1058:   mesh->vtkCellHeight       = 0;
1059:   mesh->useCone             = PETSC_FALSE;
1060:   mesh->useClosure          = PETSC_TRUE;

1062:   mesh->printSetValues = PETSC_FALSE;
1063:   mesh->printFEM       = 0;
1064:   mesh->printTol       = 1.0e-10;

1066:   DMInitialize_Plex(dm);
1067:   return(0);
1068: }

1072: /*@
1073:   DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.

1075:   Collective on MPI_Comm

1077:   Input Parameter:
1078: . comm - The communicator for the DMPlex object

1080:   Output Parameter:
1081: . mesh  - The DMPlex object

1083:   Level: beginner

1085: .keywords: DMPlex, create
1086: @*/
1087: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
1088: {

1093:   DMCreate(comm, mesh);
1094:   DMSetType(*mesh, DMPLEX);
1095:   return(0);
1096: }

1100: /*
1101:   This takes as input the common mesh generator output, a list of the vertices for each cell
1102: */
1103: PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[])
1104: {
1105:   PetscInt      *cone, c, p;

1109:   DMPlexSetChart(dm, 0, numCells+numVertices);
1110:   for (c = 0; c < numCells; ++c) {
1111:     DMPlexSetConeSize(dm, c, numCorners);
1112:   }
1113:   DMSetUp(dm);
1114:   DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);
1115:   for (c = 0; c < numCells; ++c) {
1116:     for (p = 0; p < numCorners; ++p) {
1117:       cone[p] = cells[c*numCorners+p]+numCells;
1118:     }
1119:     DMPlexSetCone(dm, c, cone);
1120:   }
1121:   DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);
1122:   DMPlexSymmetrize(dm);
1123:   DMPlexStratify(dm);
1124:   return(0);
1125: }

1129: /*
1130:   This takes as input the coordinates for each vertex
1131: */
1132: PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
1133: {
1134:   PetscSection   coordSection;
1135:   Vec            coordinates;
1136:   PetscScalar   *coords;
1137:   PetscInt       coordSize, v, d;

1141:   DMGetCoordinateSection(dm, &coordSection);
1142:   PetscSectionSetNumFields(coordSection, 1);
1143:   PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
1144:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1145:   for (v = numCells; v < numCells+numVertices; ++v) {
1146:     PetscSectionSetDof(coordSection, v, spaceDim);
1147:     PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
1148:   }
1149:   PetscSectionSetUp(coordSection);
1150:   PetscSectionGetStorageSize(coordSection, &coordSize);
1151:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1152:   VecSetBlockSize(coordinates, spaceDim);
1153:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
1154:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1155:   VecSetType(coordinates,VECSTANDARD);
1156:   VecGetArray(coordinates, &coords);
1157:   for (v = 0; v < numVertices; ++v) {
1158:     for (d = 0; d < spaceDim; ++d) {
1159:       coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
1160:     }
1161:   }
1162:   VecRestoreArray(coordinates, &coords);
1163:   DMSetCoordinatesLocal(dm, coordinates);
1164:   VecDestroy(&coordinates);
1165:   return(0);
1166: }

1170: /*@C
1171:   DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM

1173:   Input Parameters:
1174: + comm - The communicator
1175: . dim - The topological dimension of the mesh
1176: . numCells - The number of cells
1177: . numVertices - The number of vertices
1178: . numCorners - The number of vertices for each cell
1179: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
1180: . cells - An array of numCells*numCorners numbers, the vertices for each cell
1181: . spaceDim - The spatial dimension used for coordinates
1182: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex

1184:   Output Parameter:
1185: . dm - The DM

1187:   Note: Two triangles sharing a face
1188: $
1189: $        2
1190: $      / | \
1191: $     /  |  \
1192: $    /   |   \
1193: $   0  0 | 1  3
1194: $    \   |   /
1195: $     \  |  /
1196: $      \ | /
1197: $        1
1198: would have input
1199: $  numCells = 2, numVertices = 4
1200: $  cells = [0 1 2  1 3 2]
1201: $
1202: which would result in the DMPlex
1203: $
1204: $        4
1205: $      / | \
1206: $     /  |  \
1207: $    /   |   \
1208: $   2  0 | 1  5
1209: $    \   |   /
1210: $     \  |  /
1211: $      \ | /
1212: $        3

1214:   Level: beginner

1216: .seealso: DMPlexCreateFromDAG(), DMPlexCreate()
1217: @*/
1218: 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)
1219: {

1223:   DMCreate(comm, dm);
1224:   DMSetType(*dm, DMPLEX);
1225:   DMSetDimension(*dm, dim);
1226:   DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);
1227:   if (interpolate) {
1228:     DM idm = NULL;

1230:     DMPlexInterpolate(*dm, &idm);
1231:     DMDestroy(dm);
1232:     *dm  = idm;
1233:   }
1234:   DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);
1235:   return(0);
1236: }

1240: /*@
1241:   DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM

1243:   Input Parameters:
1244: + dm - The empty DM object, usually from DMCreate() and DMSetDimension()
1245: . depth - The depth of the DAG
1246: . numPoints - The number of points at each depth
1247: . coneSize - The cone size of each point
1248: . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
1249: . coneOrientations - The orientation of each cone point
1250: - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex

1252:   Output Parameter:
1253: . dm - The DM

1255:   Note: Two triangles sharing a face would have input
1256: $  depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
1257: $  cones = [2 3 4  3 5 4], coneOrientations = [0 0 0  0 0 0]
1258: $ vertexCoords = [-1.0 0.0  0.0 -1.0  0.0 1.0  1.0 0.0]
1259: $
1260: which would result in the DMPlex
1261: $
1262: $        4
1263: $      / | \
1264: $     /  |  \
1265: $    /   |   \
1266: $   2  0 | 1  5
1267: $    \   |   /
1268: $     \  |  /
1269: $      \ | /
1270: $        3
1271: $
1272: $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList()

1274:   Level: advanced

1276: .seealso: DMPlexCreateFromCellList(), DMPlexCreate()
1277: @*/
1278: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
1279: {
1280:   Vec            coordinates;
1281:   PetscSection   coordSection;
1282:   PetscScalar    *coords;
1283:   PetscInt       coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off;

1287:   DMGetDimension(dm, &dim);
1288:   for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
1289:   DMPlexSetChart(dm, pStart, pEnd);
1290:   for (p = pStart; p < pEnd; ++p) {
1291:     DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
1292:   }
1293:   DMSetUp(dm); /* Allocate space for cones */
1294:   for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
1295:     DMPlexSetCone(dm, p, &cones[off]);
1296:     DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
1297:   }
1298:   DMPlexSymmetrize(dm);
1299:   DMPlexStratify(dm);
1300:   /* Build coordinates */
1301:   DMGetCoordinateSection(dm, &coordSection);
1302:   PetscSectionSetNumFields(coordSection, 1);
1303:   PetscSectionSetFieldComponents(coordSection, 0, dim);
1304:   PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
1305:   for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
1306:     PetscSectionSetDof(coordSection, v, dim);
1307:     PetscSectionSetFieldDof(coordSection, v, 0, dim);
1308:   }
1309:   PetscSectionSetUp(coordSection);
1310:   PetscSectionGetStorageSize(coordSection, &coordSize);
1311:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1312:   VecSetBlockSize(coordinates, dim);
1313:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
1314:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1315:   VecSetType(coordinates,VECSTANDARD);
1316:   VecGetArray(coordinates, &coords);
1317:   for (v = 0; v < numPoints[0]; ++v) {
1318:     PetscInt off;

1320:     PetscSectionGetOffset(coordSection, v+firstVertex, &off);
1321:     for (d = 0; d < dim; ++d) {
1322:       coords[off+d] = vertexCoords[v*dim+d];
1323:     }
1324:   }
1325:   VecRestoreArray(coordinates, &coords);
1326:   DMSetCoordinatesLocal(dm, coordinates);
1327:   VecDestroy(&coordinates);
1328:   return(0);
1329: }