Actual source code: plexcreate.c

petsc-3.5.1 2014-07-24
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:   DMPlexSetDimension(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:   DMPlexSetDimension(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:   DMPlexSetDimension(*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:   PetscSection   coordSection;
768:   Vec            coords;

772:   DMGetPointSF(dmNew, &sf);
773:   DMSetPointSF(dm, sf);
774:   DMGetCoordinateSection(dmNew, &coordSection);
775:   DMGetCoordinatesLocal(dmNew, &coords);
776:   DMSetCoordinateSection(dm, coordSection);
777:   DMSetCoordinatesLocal(dm, coords);
778:   DMDestroy_Plex(dm);
779:   dm->data = dmNew->data;
780:   ((DM_Plex *) dmNew->data)->refct++;
781:   return(0);
782: }

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

798:   DMGetCoordinateDM(dmA, &coordDMA);
799:   DMGetCoordinateDM(dmB, &coordDMB);
800:   PetscObjectReference((PetscObject) coordDMA);
801:   DMSetCoordinateDM(dmA, coordDMB);
802:   DMSetCoordinateDM(dmB, coordDMA);
803:   PetscObjectDereference((PetscObject) coordDMA);

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

819: PetscErrorCode  DMSetFromOptions_NonRefinement_Plex(DM dm)
820: {
821:   DM_Plex       *mesh = (DM_Plex*) dm->data;
822:   DMBoundary     b;

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

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

839:       DMPlexGetLabel(dm, b->labelname, &label);
840:       for (i = 0; i < len; ++i) {
841:         PetscBool has;

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

862: PetscErrorCode  DMSetFromOptions_Plex(DM dm)
863: {
864:   PetscInt       refine = 0, r;
865:   PetscBool      isHierarchy;

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

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

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

912: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
913: {

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

926: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
927: {

931:   DMCreateLocalVector_Section_Private(dm,vec);
932:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);
933:   VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);
934:   return(0);
935: }

939: PetscErrorCode DMInitialize_Plex(DM dm)
940: {
942:   dm->ops->view                            = DMView_Plex;
943:   dm->ops->load                            = DMLoad_Plex;
944:   dm->ops->setfromoptions                  = DMSetFromOptions_Plex;
945:   dm->ops->clone                           = DMClone_Plex;
946:   dm->ops->setup                           = DMSetUp_Plex;
947:   dm->ops->createdefaultsection            = DMCreateDefaultSection_Plex;
948:   dm->ops->createglobalvector              = DMCreateGlobalVector_Plex;
949:   dm->ops->createlocalvector               = DMCreateLocalVector_Plex;
950:   dm->ops->getlocaltoglobalmapping         = NULL;
951:   dm->ops->createfieldis                   = NULL;
952:   dm->ops->createcoordinatedm              = DMCreateCoordinateDM_Plex;
953:   dm->ops->getcoloring                     = NULL;
954:   dm->ops->creatematrix                    = DMCreateMatrix_Plex;
955:   dm->ops->createinterpolation             = DMCreateInterpolation_Plex;
956:   dm->ops->getaggregates                   = NULL;
957:   dm->ops->getinjection                    = DMCreateInjection_Plex;
958:   dm->ops->refine                          = DMRefine_Plex;
959:   dm->ops->coarsen                         = DMCoarsen_Plex;
960:   dm->ops->refinehierarchy                 = DMRefineHierarchy_Plex;
961:   dm->ops->coarsenhierarchy                = NULL;
962:   dm->ops->globaltolocalbegin              = NULL;
963:   dm->ops->globaltolocalend                = NULL;
964:   dm->ops->localtoglobalbegin              = NULL;
965:   dm->ops->localtoglobalend                = NULL;
966:   dm->ops->destroy                         = DMDestroy_Plex;
967:   dm->ops->createsubdm                     = DMCreateSubDM_Plex;
968:   dm->ops->locatepoints                    = DMLocatePoints_Plex;
969:   return(0);
970: }

974: PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
975: {
976:   DM_Plex        *mesh = (DM_Plex *) dm->data;

980:   mesh->refct++;
981:   (*newdm)->data = mesh;
982:   PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
983:   DMInitialize_Plex(*newdm);
984:   return(0);
985: }

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

993:   Level: intermediate

995: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
996: M*/

1000: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
1001: {
1002:   DM_Plex        *mesh;
1003:   PetscInt       unit, d;

1008:   PetscNewLog(dm,&mesh);
1009:   dm->data = mesh;

1011:   mesh->refct             = 1;
1012:   mesh->dim               = 0;
1013:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
1014:   mesh->maxConeSize       = 0;
1015:   mesh->cones             = NULL;
1016:   mesh->coneOrientations  = NULL;
1017:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
1018:   mesh->maxSupportSize    = 0;
1019:   mesh->supports          = NULL;
1020:   mesh->refinementUniform = PETSC_TRUE;
1021:   mesh->refinementLimit   = -1.0;

1023:   mesh->facesTmp = NULL;

1025:   mesh->subpointMap = NULL;

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

1029:   mesh->labels              = NULL;
1030:   mesh->depthLabel          = NULL;
1031:   mesh->globalVertexNumbers = NULL;
1032:   mesh->globalCellNumbers   = NULL;
1033:   for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
1034:   mesh->vtkCellHeight       = 0;
1035:   mesh->useCone             = PETSC_FALSE;
1036:   mesh->useClosure          = PETSC_TRUE;

1038:   mesh->printSetValues = PETSC_FALSE;
1039:   mesh->printFEM       = 0;
1040:   mesh->printTol       = 1.0e-10;

1042:   DMInitialize_Plex(dm);
1043:   return(0);
1044: }

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

1051:   Collective on MPI_Comm

1053:   Input Parameter:
1054: . comm - The communicator for the DMPlex object

1056:   Output Parameter:
1057: . mesh  - The DMPlex object

1059:   Level: beginner

1061: .keywords: DMPlex, create
1062: @*/
1063: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
1064: {

1069:   DMCreate(comm, mesh);
1070:   DMSetType(*mesh, DMPLEX);
1071:   return(0);
1072: }

1076: /*
1077:   This takes as input the common mesh generator output, a list of the vertices for each cell
1078: */
1079: PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[])
1080: {
1081:   PetscInt      *cone, c, p;

1085:   DMPlexSetChart(dm, 0, numCells+numVertices);
1086:   for (c = 0; c < numCells; ++c) {
1087:     DMPlexSetConeSize(dm, c, numCorners);
1088:   }
1089:   DMSetUp(dm);
1090:   DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);
1091:   for (c = 0; c < numCells; ++c) {
1092:     for (p = 0; p < numCorners; ++p) {
1093:       cone[p] = cells[c*numCorners+p]+numCells;
1094:     }
1095:     DMPlexSetCone(dm, c, cone);
1096:   }
1097:   DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);
1098:   DMPlexSymmetrize(dm);
1099:   DMPlexStratify(dm);
1100:   return(0);
1101: }

1105: /*
1106:   This takes as input the coordinates for each vertex
1107: */
1108: PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
1109: {
1110:   PetscSection   coordSection;
1111:   Vec            coordinates;
1112:   PetscScalar   *coords;
1113:   PetscInt       coordSize, v, d;

1117:   DMGetCoordinateSection(dm, &coordSection);
1118:   PetscSectionSetNumFields(coordSection, 1);
1119:   PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
1120:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1121:   for (v = numCells; v < numCells+numVertices; ++v) {
1122:     PetscSectionSetDof(coordSection, v, spaceDim);
1123:     PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
1124:   }
1125:   PetscSectionSetUp(coordSection);
1126:   PetscSectionGetStorageSize(coordSection, &coordSize);
1127:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1128:   VecSetBlockSize(coordinates, spaceDim);
1129:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
1130:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1131:   VecSetType(coordinates,VECSTANDARD);
1132:   VecGetArray(coordinates, &coords);
1133:   for (v = 0; v < numVertices; ++v) {
1134:     for (d = 0; d < spaceDim; ++d) {
1135:       coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
1136:     }
1137:   }
1138:   VecRestoreArray(coordinates, &coords);
1139:   DMSetCoordinatesLocal(dm, coordinates);
1140:   VecDestroy(&coordinates);
1141:   return(0);
1142: }

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

1149:   Input Parameters:
1150: + comm - The communicator
1151: . dim - The topological dimension of the mesh
1152: . numCells - The number of cells
1153: . numVertices - The number of vertices
1154: . numCorners - The number of vertices for each cell
1155: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
1156: . cells - An array of numCells*numCorners numbers, the vertices for each cell
1157: . spaceDim - The spatial dimension used for coordinates
1158: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex

1160:   Output Parameter:
1161: . dm - The DM

1163:   Note: Two triangles sharing a face
1164: $
1165: $        2
1166: $      / | \
1167: $     /  |  \
1168: $    /   |   \
1169: $   0  0 | 1  3
1170: $    \   |   /
1171: $     \  |  /
1172: $      \ | /
1173: $        1
1174: would have input
1175: $  numCells = 2, numVertices = 4
1176: $  cells = [0 1 2  1 3 2]
1177: $
1178: which would result in the DMPlex
1179: $
1180: $        4
1181: $      / | \
1182: $     /  |  \
1183: $    /   |   \
1184: $   2  0 | 1  5
1185: $    \   |   /
1186: $     \  |  /
1187: $      \ | /
1188: $        3

1190:   Level: beginner

1192: .seealso: DMPlexCreateFromDAG(), DMPlexCreate()
1193: @*/
1194: 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)
1195: {

1199:   DMCreate(comm, dm);
1200:   DMSetType(*dm, DMPLEX);
1201:   DMPlexSetDimension(*dm, dim);
1202:   DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);
1203:   if (interpolate) {
1204:     DM idm = NULL;

1206:     DMPlexInterpolate(*dm, &idm);
1207:     DMDestroy(dm);
1208:     *dm  = idm;
1209:   }
1210:   DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);
1211:   return(0);
1212: }

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

1219:   Input Parameters:
1220: + dm - The empty DM object, usually from DMCreate() and DMPlexSetDimension()
1221: . depth - The depth of the DAG
1222: . numPoints - The number of points at each depth
1223: . coneSize - The cone size of each point
1224: . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
1225: . coneOrientations - The orientation of each cone point
1226: - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex

1228:   Output Parameter:
1229: . dm - The DM

1231:   Note: Two triangles sharing a face would have input
1232: $  depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
1233: $  cones = [2 3 4  3 5 4], coneOrientations = [0 0 0  0 0 0]
1234: $ vertexCoords = [-1.0 0.0  0.0 -1.0  0.0 1.0  1.0 0.0]
1235: $
1236: which would result in the DMPlex
1237: $
1238: $        4
1239: $      / | \
1240: $     /  |  \
1241: $    /   |   \
1242: $   2  0 | 1  5
1243: $    \   |   /
1244: $     \  |  /
1245: $      \ | /
1246: $        3
1247: $
1248: $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList()

1250:   Level: advanced

1252: .seealso: DMPlexCreateFromCellList(), DMPlexCreate()
1253: @*/
1254: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
1255: {
1256:   Vec            coordinates;
1257:   PetscSection   coordSection;
1258:   PetscScalar    *coords;
1259:   PetscInt       coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off;

1263:   DMPlexGetDimension(dm, &dim);
1264:   for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
1265:   DMPlexSetChart(dm, pStart, pEnd);
1266:   for (p = pStart; p < pEnd; ++p) {
1267:     DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
1268:   }
1269:   DMSetUp(dm); /* Allocate space for cones */
1270:   for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
1271:     DMPlexSetCone(dm, p, &cones[off]);
1272:     DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
1273:   }
1274:   DMPlexSymmetrize(dm);
1275:   DMPlexStratify(dm);
1276:   /* Build coordinates */
1277:   DMGetCoordinateSection(dm, &coordSection);
1278:   PetscSectionSetNumFields(coordSection, 1);
1279:   PetscSectionSetFieldComponents(coordSection, 0, dim);
1280:   PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
1281:   for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
1282:     PetscSectionSetDof(coordSection, v, dim);
1283:     PetscSectionSetFieldDof(coordSection, v, 0, dim);
1284:   }
1285:   PetscSectionSetUp(coordSection);
1286:   PetscSectionGetStorageSize(coordSection, &coordSize);
1287:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1288:   VecSetBlockSize(coordinates, dim);
1289:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
1290:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1291:   VecSetType(coordinates,VECSTANDARD);
1292:   VecGetArray(coordinates, &coords);
1293:   for (v = 0; v < numPoints[0]; ++v) {
1294:     PetscInt off;

1296:     PetscSectionGetOffset(coordSection, v+firstVertex, &off);
1297:     for (d = 0; d < dim; ++d) {
1298:       coords[off+d] = vertexCoords[v*dim+d];
1299:     }
1300:   }
1301:   VecRestoreArray(coordinates, &coords);
1302:   DMSetCoordinatesLocal(dm, coordinates);
1303:   VecDestroy(&coordinates);
1304:   return(0);
1305: }