Actual source code: plexcreate.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  1: #define PETSCDM_DLL
  2:  #include <petsc/private/dmpleximpl.h>
  3:  #include <petsc/private/hashseti.h>
  4:  #include <petscsf.h>

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

  9:   Collective on MPI_Comm

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

 19:   Output Parameter:
 20: . dm - The DM object

 22:   Level: beginner

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

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

 65:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
 66:         for (p = 0; p < 4; ++p) {DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);}
 67:       } else {
 68:         PetscInt    numPoints[2]        = {6, 2};
 69:         PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
 70:         PetscInt    cones[8]            = {2, 3, 4, 5,  3, 6, 7, 4};
 71:         PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
 72:         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};

 74:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
 75:       }
 76:       break;
 77:     case 3:
 78:       if (simplex) {
 79:         PetscInt    numPoints[2]        = {5, 2};
 80:         PetscInt    coneSize[7]         = {4, 4, 0, 0, 0, 0, 0};
 81:         PetscInt    cones[8]            = {4, 3, 5, 2,  5, 3, 4, 6};
 82:         PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
 83:         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};
 84:         PetscInt    markerPoints[10]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};

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

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

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

120:     DMPlexInterpolate(*newdm, &idm);
121:     DMDestroy(newdm);
122:     *newdm = idm;
123:   }
124:   {
125:     DM refinedMesh     = NULL;
126:     DM distributedMesh = NULL;

128:     /* Distribute mesh over processes */
129:     DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);
130:     if (distributedMesh) {
131:       DMDestroy(newdm);
132:       *newdm = distributedMesh;
133:     }
134:     if (refinementUniform) {
135:       DMPlexSetRefinementUniform(*newdm, refinementUniform);
136:       DMRefine(*newdm, comm, &refinedMesh);
137:       if (refinedMesh) {
138:         DMDestroy(newdm);
139:         *newdm = refinedMesh;
140:       }
141:     }
142:   }
143:   return(0);
144: }

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

149:   Collective on MPI_Comm

151:   Input Parameters:
152: + comm  - The communicator for the DM object
153: . lower - The lower left corner coordinates
154: . upper - The upper right corner coordinates
155: - edges - The number of cells in each direction

157:   Output Parameter:
158: . dm  - The DM object

160:   Note: Here is the numbering returned for 2 cells in each direction:
161: $ 18--5-17--4--16
162: $  |     |     |
163: $  6    10     3
164: $  |     |     |
165: $ 19-11-20--9--15
166: $  |     |     |
167: $  7     8     2
168: $  |     |     |
169: $ 12--0-13--1--14

171:   Level: beginner

173: .keywords: DM, create
174: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate()
175: @*/
176: PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
177: {
178:   const PetscInt numVertices    = (edges[0]+1)*(edges[1]+1);
179:   const PetscInt numEdges       = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
180:   PetscInt       markerTop      = 1;
181:   PetscInt       markerBottom   = 1;
182:   PetscInt       markerRight    = 1;
183:   PetscInt       markerLeft     = 1;
184:   PetscBool      markerSeparate = PETSC_FALSE;
185:   Vec            coordinates;
186:   PetscSection   coordSection;
187:   PetscScalar    *coords;
188:   PetscInt       coordSize;
189:   PetscMPIInt    rank;
190:   PetscInt       v, vx, vy;

194:   PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
195:   if (markerSeparate) {
196:     markerTop    = 3;
197:     markerBottom = 1;
198:     markerRight  = 2;
199:     markerLeft   = 4;
200:   }
201:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
202:   if (!rank) {
203:     PetscInt e, ex, ey;

205:     DMPlexSetChart(dm, 0, numEdges+numVertices);
206:     for (e = 0; e < numEdges; ++e) {
207:       DMPlexSetConeSize(dm, e, 2);
208:     }
209:     DMSetUp(dm); /* Allocate space for cones */
210:     for (vx = 0; vx <= edges[0]; vx++) {
211:       for (ey = 0; ey < edges[1]; ey++) {
212:         PetscInt edge   = vx*edges[1] + ey + edges[0]*(edges[1]+1);
213:         PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
214:         PetscInt cone[2];

216:         cone[0] = vertex; cone[1] = vertex+edges[0]+1;
217:         DMPlexSetCone(dm, edge, cone);
218:         if (vx == edges[0]) {
219:           DMSetLabelValue(dm, "marker", edge,    markerRight);
220:           DMSetLabelValue(dm, "marker", cone[0], markerRight);
221:           if (ey == edges[1]-1) {
222:             DMSetLabelValue(dm, "marker", cone[1], markerRight);
223:             DMSetLabelValue(dm, "Face Sets", cone[1], markerRight);
224:           }
225:         } else if (vx == 0) {
226:           DMSetLabelValue(dm, "marker", edge,    markerLeft);
227:           DMSetLabelValue(dm, "marker", cone[0], markerLeft);
228:           if (ey == edges[1]-1) {
229:             DMSetLabelValue(dm, "marker", cone[1], markerLeft);
230:             DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft);
231:           }
232:         }
233:       }
234:     }
235:     for (vy = 0; vy <= edges[1]; vy++) {
236:       for (ex = 0; ex < edges[0]; ex++) {
237:         PetscInt edge   = vy*edges[0]     + ex;
238:         PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
239:         PetscInt cone[2];

241:         cone[0] = vertex; cone[1] = vertex+1;
242:         DMPlexSetCone(dm, edge, cone);
243:         if (vy == edges[1]) {
244:           DMSetLabelValue(dm, "marker", edge,    markerTop);
245:           DMSetLabelValue(dm, "marker", cone[0], markerTop);
246:           if (ex == edges[0]-1) {
247:             DMSetLabelValue(dm, "marker", cone[1], markerTop);
248:             DMSetLabelValue(dm, "Face Sets", cone[1], markerTop);
249:           }
250:         } else if (vy == 0) {
251:           DMSetLabelValue(dm, "marker", edge,    markerBottom);
252:           DMSetLabelValue(dm, "marker", cone[0], markerBottom);
253:           if (ex == edges[0]-1) {
254:             DMSetLabelValue(dm, "marker", cone[1], markerBottom);
255:             DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom);
256:           }
257:         }
258:       }
259:     }
260:   }
261:   DMPlexSymmetrize(dm);
262:   DMPlexStratify(dm);
263:   /* Build coordinates */
264:   DMSetCoordinateDim(dm, 2);
265:   DMGetCoordinateSection(dm, &coordSection);
266:   PetscSectionSetNumFields(coordSection, 1);
267:   PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);
268:   PetscSectionSetFieldComponents(coordSection, 0, 2);
269:   for (v = numEdges; v < numEdges+numVertices; ++v) {
270:     PetscSectionSetDof(coordSection, v, 2);
271:     PetscSectionSetFieldDof(coordSection, v, 0, 2);
272:   }
273:   PetscSectionSetUp(coordSection);
274:   PetscSectionGetStorageSize(coordSection, &coordSize);
275:   VecCreate(PETSC_COMM_SELF, &coordinates);
276:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
277:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
278:   VecSetBlockSize(coordinates, 2);
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: }

293: /*@
294:   DMPlexCreateCubeBoundary - Creates a 2D mesh that is the boundary of a cubic lattice.

296:   Collective on MPI_Comm

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

304:   Output Parameter:
305: . dm  - The DM object

307:   Level: beginner

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

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

333:     DMPlexSetChart(dm, 0, numFaces+numVertices);
334:     for (f = 0; f < numFaces; ++f) {
335:       DMPlexSetConeSize(dm, f, 4);
336:     }
337:     DMSetUp(dm); /* Allocate space for cones */

339:     /* Side 0 (Top) */
340:     for (vy = 0; vy < faces[1]; vy++) {
341:       for (vx = 0; vx < faces[0]; vx++) {
342:         voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx;
343:         cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0];
344:         DMPlexSetCone(dm, iface, cone);
345:         DMSetLabelValue(dm, "marker", iface, 1);
346:         DMSetLabelValue(dm, "marker", voffset+0, 1);
347:         DMSetLabelValue(dm, "marker", voffset+1, 1);
348:         DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);
349:         DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);
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:         DMSetLabelValue(dm, "marker", iface, 1);
361:         DMSetLabelValue(dm, "marker", voffset+0, 1);
362:         DMSetLabelValue(dm, "marker", voffset+1, 1);
363:         DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);
364:         DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);
365:         iface++;
366:       }
367:     }

369:     /* Side 2 (Front) */
370:     for (vz = 0; vz < faces[2]; vz++) {
371:       for (vx = 0; vx < faces[0]; vx++) {
372:         voffset = numFaces + vz*vertices[0]*vertices[1] + vx;
373:         cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1];
374:         DMPlexSetCone(dm, iface, cone);
375:         DMSetLabelValue(dm, "marker", iface, 1);
376:         DMSetLabelValue(dm, "marker", voffset+0, 1);
377:         DMSetLabelValue(dm, "marker", voffset+1, 1);
378:         DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);
379:         DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);
380:         iface++;
381:       }
382:     }

384:     /* Side 3 (Back) */
385:     for (vz = 0; vz < faces[2]; vz++) {
386:       for (vx = 0; vx < faces[0]; vx++) {
387:         voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx;
388:         cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1;
389:         cone[2] = voffset+1; cone[3] = voffset;
390:         DMPlexSetCone(dm, iface, cone);
391:         DMSetLabelValue(dm, "marker", iface, 1);
392:         DMSetLabelValue(dm, "marker", voffset+0, 1);
393:         DMSetLabelValue(dm, "marker", voffset+1, 1);
394:         DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);
395:         DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);
396:         iface++;
397:       }
398:     }

400:     /* Side 4 (Left) */
401:     for (vz = 0; vz < faces[2]; vz++) {
402:       for (vy = 0; vy < faces[1]; vy++) {
403:         voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0];
404:         cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1];
405:         cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0];
406:         DMPlexSetCone(dm, iface, cone);
407:         DMSetLabelValue(dm, "marker", iface, 1);
408:         DMSetLabelValue(dm, "marker", voffset+0, 1);
409:         DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);
410:         DMSetLabelValue(dm, "marker", voffset+vertices[1]+0, 1);
411:         DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);
412:         iface++;
413:       }
414:     }

416:     /* Side 5 (Right) */
417:     for (vz = 0; vz < faces[2]; vz++) {
418:       for (vy = 0; vy < faces[1]; vy++) {
419:         voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + faces[0];
420:         cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset;
421:         cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0];
422:         DMPlexSetCone(dm, iface, cone);
423:         DMSetLabelValue(dm, "marker", iface, 1);
424:         DMSetLabelValue(dm, "marker", voffset+0, 1);
425:         DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);
426:         DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);
427:         DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);
428:         iface++;
429:       }
430:     }
431:   }
432:   DMPlexSymmetrize(dm);
433:   DMPlexStratify(dm);
434:   /* Build coordinates */
435:   DMSetCoordinateDim(dm, 3);
436:   DMGetCoordinateSection(dm, &coordSection);
437:   PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);
438:   for (v = numFaces; v < numFaces+numVertices; ++v) {
439:     PetscSectionSetDof(coordSection, v, 3);
440:   }
441:   PetscSectionSetUp(coordSection);
442:   PetscSectionGetStorageSize(coordSection, &coordSize);
443:   VecCreate(PETSC_COMM_SELF, &coordinates);
444:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
445:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
446:   VecSetBlockSize(coordinates, 3);
447:   VecSetType(coordinates,VECSTANDARD);
448:   VecGetArray(coordinates, &coords);
449:   for (vz = 0; vz <= faces[2]; ++vz) {
450:     for (vy = 0; vy <= faces[1]; ++vy) {
451:       for (vx = 0; vx <= faces[0]; ++vx) {
452:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
453:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
454:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
455:       }
456:     }
457:   }
458:   VecRestoreArray(coordinates, &coords);
459:   DMSetCoordinatesLocal(dm, coordinates);
460:   VecDestroy(&coordinates);
461:   return(0);
462: }

464: static PetscErrorCode DMPlexCreateLineMesh_Internal(MPI_Comm comm,PetscInt segments,PetscReal lower,PetscReal upper,DMBoundaryType bd,DM *dm)
465: {
466:   PetscInt       i,fStart,fEnd,numCells = 0,numVerts = 0;
467:   PetscInt       numPoints[2],*coneSize,*cones,*coneOrientations;
468:   PetscScalar    *vertexCoords;
469:   PetscReal      L,maxCell;
470:   PetscBool      markerSeparate = PETSC_FALSE;
471:   PetscInt       markerLeft  = 1, faceMarkerLeft  = 1;
472:   PetscInt       markerRight = 1, faceMarkerRight = 2;
473:   PetscBool      wrap = (bd == DM_BOUNDARY_PERIODIC || bd == DM_BOUNDARY_TWIST) ? PETSC_TRUE : PETSC_FALSE;
474:   PetscMPIInt    rank;


480:   DMCreate(comm,dm);
481:   DMSetType(*dm,DMPLEX);
482:   DMSetDimension(*dm,1);
483:   DMCreateLabel(*dm,"marker");
484:   DMCreateLabel(*dm,"Face Sets");

486:   MPI_Comm_rank(comm,&rank);
487:   if (!rank) numCells = segments;
488:   if (!rank) numVerts = segments + (wrap ? 0 : 1);

490:   numPoints[0] = numVerts ; numPoints[1] = numCells;
491:   PetscMalloc4(numCells+numVerts,&coneSize,numCells*2,&cones,numCells+numVerts,&coneOrientations,numVerts,&vertexCoords);
492:   PetscMemzero(coneOrientations,(numCells+numVerts)*sizeof(PetscInt));
493:   for (i = 0; i < numCells; ++i) { coneSize[i] = 2; }
494:   for (i = 0; i < numVerts; ++i) { coneSize[numCells+i] = 0; }
495:   for (i = 0; i < numCells; ++i) { cones[2*i] = numCells + i%numVerts; cones[2*i+1] = numCells + (i+1)%numVerts; }
496:   for (i = 0; i < numVerts; ++i) { vertexCoords[i] = lower + (upper-lower)*((PetscReal)i/(PetscReal)numCells); }
497:   DMPlexCreateFromDAG(*dm,1,numPoints,coneSize,cones,coneOrientations,vertexCoords);
498:   PetscFree4(coneSize,cones,coneOrientations,vertexCoords);

500:   PetscOptionsGetBool(((PetscObject)*dm)->options,((PetscObject)*dm)->prefix,"-dm_plex_separate_marker",&markerSeparate,NULL);
501:   if (markerSeparate) { markerLeft = faceMarkerLeft; markerRight = faceMarkerRight;}
502:   if (!wrap && !rank) {
503:     DMPlexGetHeightStratum(*dm,1,&fStart,&fEnd);
504:     DMSetLabelValue(*dm,"marker",fStart,markerLeft);
505:     DMSetLabelValue(*dm,"marker",fEnd-1,markerRight);
506:     DMSetLabelValue(*dm,"Face Sets",fStart,faceMarkerLeft);
507:     DMSetLabelValue(*dm,"Face Sets",fEnd-1,faceMarkerRight);
508:   }
509:   if (wrap) {
510:     L       = upper - lower;
511:     maxCell = (PetscReal)1.1*(L/(PetscReal)PetscMax(1,segments));
512:     DMSetPeriodicity(*dm,PETSC_TRUE,&maxCell,&L,&bd);
513:   }
514:   return(0);
515: }

517: static PetscErrorCode DMPlexCreateBoxMesh_Simplex_Internal(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm)
518: {
519:   DM             boundary;
520:   PetscInt       i;

525:   for (i = 0; i < dim; ++i) if (periodicity[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity is not supported for simplex meshes");
526:   DMCreate(comm, &boundary);
528:   DMSetType(boundary, DMPLEX);
529:   DMSetDimension(boundary, dim-1);
530:   DMSetCoordinateDim(boundary, dim);
531:   switch (dim) {
532:   case 2: DMPlexCreateSquareBoundary(boundary, lower, upper, faces);break;
533:   case 3: DMPlexCreateCubeBoundary(boundary, lower, upper, faces);break;
534:   default: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
535:   }
536:   DMPlexGenerate(boundary, NULL, interpolate, dm);
537:   DMDestroy(&boundary);
538:   return(0);
539: }

541: static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ)
542: {
543:   DMLabel        cutLabel = NULL;
544:   PetscInt       markerTop      = 1, faceMarkerTop      = 1;
545:   PetscInt       markerBottom   = 1, faceMarkerBottom   = 1;
546:   PetscInt       markerFront    = 1, faceMarkerFront    = 1;
547:   PetscInt       markerBack     = 1, faceMarkerBack     = 1;
548:   PetscInt       markerRight    = 1, faceMarkerRight    = 1;
549:   PetscInt       markerLeft     = 1, faceMarkerLeft     = 1;
550:   PetscInt       dim;
551:   PetscBool      markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE;
552:   PetscMPIInt    rank;

556:   DMGetDimension(dm,&dim);
557:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
558:   DMCreateLabel(dm,"marker");
559:   DMCreateLabel(dm,"Face Sets");
560:   PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);
561:   if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ||
562:       bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ||
563:       bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) {

565:     if (cutMarker) {DMCreateLabel(dm, "periodic_cut"); DMGetLabel(dm, "periodic_cut", &cutLabel);}
566:   }
567:   switch (dim) {
568:   case 2:
569:     faceMarkerTop    = 3;
570:     faceMarkerBottom = 1;
571:     faceMarkerRight  = 2;
572:     faceMarkerLeft   = 4;
573:     break;
574:   case 3:
575:     faceMarkerBottom = 1;
576:     faceMarkerTop    = 2;
577:     faceMarkerFront  = 3;
578:     faceMarkerBack   = 4;
579:     faceMarkerRight  = 5;
580:     faceMarkerLeft   = 6;
581:     break;
582:   default:
583:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim);
584:     break;
585:   }
586:   PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
587:   if (markerSeparate) {
588:     markerBottom = faceMarkerBottom;
589:     markerTop    = faceMarkerTop;
590:     markerFront  = faceMarkerFront;
591:     markerBack   = faceMarkerBack;
592:     markerRight  = faceMarkerRight;
593:     markerLeft   = faceMarkerLeft;
594:   }
595:   {
596:     const PetscInt numXEdges    = !rank ? edges[0] : 0;
597:     const PetscInt numYEdges    = !rank ? edges[1] : 0;
598:     const PetscInt numZEdges    = !rank ? edges[2] : 0;
599:     const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0;
600:     const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0;
601:     const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0;
602:     const PetscInt numCells     = numXEdges*numYEdges*numZEdges;
603:     const PetscInt numXFaces    = numYEdges*numZEdges;
604:     const PetscInt numYFaces    = numXEdges*numZEdges;
605:     const PetscInt numZFaces    = numXEdges*numYEdges;
606:     const PetscInt numTotXFaces = numXVertices*numXFaces;
607:     const PetscInt numTotYFaces = numYVertices*numYFaces;
608:     const PetscInt numTotZFaces = numZVertices*numZFaces;
609:     const PetscInt numFaces     = numTotXFaces + numTotYFaces + numTotZFaces;
610:     const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices;
611:     const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices;
612:     const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices;
613:     const PetscInt numVertices  = numXVertices*numYVertices*numZVertices;
614:     const PetscInt numEdges     = numTotXEdges + numTotYEdges + numTotZEdges;
615:     const PetscInt firstVertex  = (dim == 2) ? numFaces : numCells;
616:     const PetscInt firstXFace   = (dim == 2) ? 0 : numCells + numVertices;
617:     const PetscInt firstYFace   = firstXFace + numTotXFaces;
618:     const PetscInt firstZFace   = firstYFace + numTotYFaces;
619:     const PetscInt firstXEdge   = numCells + numFaces + numVertices;
620:     const PetscInt firstYEdge   = firstXEdge + numTotXEdges;
621:     const PetscInt firstZEdge   = firstYEdge + numTotYEdges;
622:     Vec            coordinates;
623:     PetscSection   coordSection;
624:     PetscScalar   *coords;
625:     PetscInt       coordSize;
626:     PetscInt       v, vx, vy, vz;
627:     PetscInt       c, f, fx, fy, fz, e, ex, ey, ez;

629:     DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);
630:     for (c = 0; c < numCells; c++) {
631:       DMPlexSetConeSize(dm, c, 6);
632:     }
633:     for (f = firstXFace; f < firstXFace+numFaces; ++f) {
634:       DMPlexSetConeSize(dm, f, 4);
635:     }
636:     for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
637:       DMPlexSetConeSize(dm, e, 2);
638:     }
639:     DMSetUp(dm); /* Allocate space for cones */
640:     /* Build cells */
641:     for (fz = 0; fz < numZEdges; ++fz) {
642:       for (fy = 0; fy < numYEdges; ++fy) {
643:         for (fx = 0; fx < numXEdges; ++fx) {
644:           PetscInt cell    = (fz*numYEdges + fy)*numXEdges + fx;
645:           PetscInt faceB   = firstZFace + (fy*numXEdges+fx)*numZVertices +   fz;
646:           PetscInt faceT   = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices);
647:           PetscInt faceF   = firstYFace + (fz*numXEdges+fx)*numYVertices +   fy;
648:           PetscInt faceK   = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices);
649:           PetscInt faceL   = firstXFace + (fz*numYEdges+fy)*numXVertices +   fx;
650:           PetscInt faceR   = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices);
651:                             /* B,  T,  F,  K,  R,  L */
652:           PetscInt ornt[6] = {-4,  0,  0, -1,  0, -4}; /* ??? */
653:           PetscInt cone[6];

655:           /* no boundary twisting in 3D */
656:           cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL;
657:           DMPlexSetCone(dm, cell, cone);
658:           DMPlexSetConeOrientation(dm, cell, ornt);
659:           if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, cell, 2);}
660:           if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, cell, 2);}
661:           if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, cell, 2);}
662:         }
663:       }
664:     }
665:     /* Build x faces */
666:     for (fz = 0; fz < numZEdges; ++fz) {
667:       for (fy = 0; fy < numYEdges; ++fy) {
668:         for (fx = 0; fx < numXVertices; ++fx) {
669:           PetscInt face    = firstXFace + (fz*numYEdges+fy)*numXVertices + fx;
670:           PetscInt edgeL   = firstZEdge + (  fy*                 numXVertices+fx)*numZEdges + fz;
671:           PetscInt edgeR   = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz;
672:           PetscInt edgeB   = firstYEdge + (  fz*                 numXVertices+fx)*numYEdges + fy;
673:           PetscInt edgeT   = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy;
674:           PetscInt ornt[4] = {0, 0, -2, -2};
675:           PetscInt cone[4];

677:           if (dim == 3) {
678:             /* markers */
679:             if (bdX != DM_BOUNDARY_PERIODIC) {
680:               if (fx == numXVertices-1) {
681:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);
682:                 DMSetLabelValue(dm, "marker", face, markerRight);
683:               }
684:               else if (fx == 0) {
685:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);
686:                 DMSetLabelValue(dm, "marker", face, markerLeft);
687:               }
688:             }
689:           }
690:           cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
691:           DMPlexSetCone(dm, face, cone);
692:           DMPlexSetConeOrientation(dm, face, ornt);
693:         }
694:       }
695:     }
696:     /* Build y faces */
697:     for (fz = 0; fz < numZEdges; ++fz) {
698:       for (fx = 0; fx < numXEdges; ++fx) {
699:         for (fy = 0; fy < numYVertices; ++fy) {
700:           PetscInt face    = firstYFace + (fz*numXEdges+fx)*numYVertices + fy;
701:           PetscInt edgeL   = firstZEdge + (fy*numXVertices+  fx                 )*numZEdges + fz;
702:           PetscInt edgeR   = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz;
703:           PetscInt edgeB   = firstXEdge + (  fz                 *numYVertices+fy)*numXEdges + fx;
704:           PetscInt edgeT   = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx;
705:           PetscInt ornt[4] = {0, 0, -2, -2};
706:           PetscInt cone[4];

708:           if (dim == 3) {
709:             /* markers */
710:             if (bdY != DM_BOUNDARY_PERIODIC) {
711:               if (fy == numYVertices-1) {
712:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);
713:                 DMSetLabelValue(dm, "marker", face, markerBack);
714:               }
715:               else if (fy == 0) {
716:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);
717:                 DMSetLabelValue(dm, "marker", face, markerFront);
718:               }
719:             }
720:           }
721:           cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
722:           DMPlexSetCone(dm, face, cone);
723:           DMPlexSetConeOrientation(dm, face, ornt);
724:         }
725:       }
726:     }
727:     /* Build z faces */
728:     for (fy = 0; fy < numYEdges; ++fy) {
729:       for (fx = 0; fx < numXEdges; ++fx) {
730:         for (fz = 0; fz < numZVertices; fz++) {
731:           PetscInt face    = firstZFace + (fy*numXEdges+fx)*numZVertices + fz;
732:           PetscInt edgeL   = firstYEdge + (fz*numXVertices+  fx                 )*numYEdges + fy;
733:           PetscInt edgeR   = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy;
734:           PetscInt edgeB   = firstXEdge + (fz*numYVertices+  fy                 )*numXEdges + fx;
735:           PetscInt edgeT   = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx;
736:           PetscInt ornt[4] = {0, 0, -2, -2};
737:           PetscInt cone[4];

739:           if (dim == 2) {
740:             if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;}
741:             if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] =  0;}
742:             if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, face, 2);}
743:             if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, face, 2);}
744:           } else {
745:             /* markers */
746:             if (bdZ != DM_BOUNDARY_PERIODIC) {
747:               if (fz == numZVertices-1) {
748:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);
749:                 DMSetLabelValue(dm, "marker", face, markerTop);
750:               }
751:               else if (fz == 0) {
752:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);
753:                 DMSetLabelValue(dm, "marker", face, markerBottom);
754:               }
755:             }
756:           }
757:           cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
758:           DMPlexSetCone(dm, face, cone);
759:           DMPlexSetConeOrientation(dm, face, ornt);
760:         }
761:       }
762:     }
763:     /* Build Z edges*/
764:     for (vy = 0; vy < numYVertices; vy++) {
765:       for (vx = 0; vx < numXVertices; vx++) {
766:         for (ez = 0; ez < numZEdges; ez++) {
767:           const PetscInt edge    = firstZEdge  + (vy*numXVertices+vx)*numZEdges + ez;
768:           const PetscInt vertexB = firstVertex + (  ez                 *numYVertices+vy)*numXVertices + vx;
769:           const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx;
770:           PetscInt       cone[2];

772:           if (dim == 3) {
773:             if (bdX != DM_BOUNDARY_PERIODIC) {
774:               if (vx == numXVertices-1) {
775:                 DMSetLabelValue(dm, "marker", edge, markerRight);
776:               }
777:               else if (vx == 0) {
778:                 DMSetLabelValue(dm, "marker", edge, markerLeft);
779:               }
780:             }
781:             if (bdY != DM_BOUNDARY_PERIODIC) {
782:               if (vy == numYVertices-1) {
783:                 DMSetLabelValue(dm, "marker", edge, markerBack);
784:               }
785:               else if (vy == 0) {
786:                 DMSetLabelValue(dm, "marker", edge, markerFront);
787:               }
788:             }
789:           }
790:           cone[0] = vertexB; cone[1] = vertexT;
791:           DMPlexSetCone(dm, edge, cone);
792:         }
793:       }
794:     }
795:     /* Build Y edges*/
796:     for (vz = 0; vz < numZVertices; vz++) {
797:       for (vx = 0; vx < numXVertices; vx++) {
798:         for (ey = 0; ey < numYEdges; ey++) {
799:           const PetscInt nextv   = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx;
800:           const PetscInt edge    = firstYEdge  + (vz*numXVertices+vx)*numYEdges + ey;
801:           const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx;
802:           const PetscInt vertexK = firstVertex + nextv;
803:           PetscInt       cone[2];

805:           cone[0] = vertexF; cone[1] = vertexK;
806:           DMPlexSetCone(dm, edge, cone);
807:           if (dim == 2) {
808:             if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) {
809:               if (vx == numXVertices-1) {
810:                 DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);
811:                 DMSetLabelValue(dm, "marker", edge,    markerRight);
812:                 DMSetLabelValue(dm, "marker", cone[0], markerRight);
813:                 if (ey == numYEdges-1) {
814:                   DMSetLabelValue(dm, "marker", cone[1], markerRight);
815:                 }
816:               } else if (vx == 0) {
817:                 DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);
818:                 DMSetLabelValue(dm, "marker", edge,    markerLeft);
819:                 DMSetLabelValue(dm, "marker", cone[0], markerLeft);
820:                 if (ey == numYEdges-1) {
821:                   DMSetLabelValue(dm, "marker", cone[1], markerLeft);
822:                 }
823:               }
824:             } else {
825:               if (vx == 0 && cutLabel) {
826:                 DMLabelSetValue(cutLabel, edge,    1);
827:                 DMLabelSetValue(cutLabel, cone[0], 1);
828:                 if (ey == numYEdges-1) {
829:                   DMLabelSetValue(cutLabel, cone[1], 1);
830:                 }
831:               }
832:             }
833:           } else {
834:             if (bdX != DM_BOUNDARY_PERIODIC) {
835:               if (vx == numXVertices-1) {
836:                 DMSetLabelValue(dm, "marker", edge, markerRight);
837:               } else if (vx == 0) {
838:                 DMSetLabelValue(dm, "marker", edge, markerLeft);
839:               }
840:             }
841:             if (bdZ != DM_BOUNDARY_PERIODIC) {
842:               if (vz == numZVertices-1) {
843:                 DMSetLabelValue(dm, "marker", edge, markerTop);
844:               } else if (vz == 0) {
845:                 DMSetLabelValue(dm, "marker", edge, markerBottom);
846:               }
847:             }
848:           }
849:         }
850:       }
851:     }
852:     /* Build X edges*/
853:     for (vz = 0; vz < numZVertices; vz++) {
854:       for (vy = 0; vy < numYVertices; vy++) {
855:         for (ex = 0; ex < numXEdges; ex++) {
856:           const PetscInt nextv   = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices;
857:           const PetscInt edge    = firstXEdge  + (vz*numYVertices+vy)*numXEdges + ex;
858:           const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex;
859:           const PetscInt vertexR = firstVertex + nextv;
860:           PetscInt       cone[2];

862:           cone[0] = vertexL; cone[1] = vertexR;
863:           DMPlexSetCone(dm, edge, cone);
864:           if (dim == 2) {
865:             if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) {
866:               if (vy == numYVertices-1) {
867:                 DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);
868:                 DMSetLabelValue(dm, "marker", edge,    markerTop);
869:                 DMSetLabelValue(dm, "marker", cone[0], markerTop);
870:                 if (ex == numXEdges-1) {
871:                   DMSetLabelValue(dm, "marker", cone[1], markerTop);
872:                 }
873:               } else if (vy == 0) {
874:                 DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);
875:                 DMSetLabelValue(dm, "marker", edge,    markerBottom);
876:                 DMSetLabelValue(dm, "marker", cone[0], markerBottom);
877:                 if (ex == numXEdges-1) {
878:                   DMSetLabelValue(dm, "marker", cone[1], markerBottom);
879:                 }
880:               }
881:             } else {
882:               if (vy == 0 && cutLabel) {
883:                 DMLabelSetValue(cutLabel, edge,    1);
884:                 DMLabelSetValue(cutLabel, cone[0], 1);
885:                 if (ex == numXEdges-1) {
886:                   DMLabelSetValue(cutLabel, cone[1], 1);
887:                 }
888:               }
889:             }
890:           } else {
891:             if (bdY != DM_BOUNDARY_PERIODIC) {
892:               if (vy == numYVertices-1) {
893:                 DMSetLabelValue(dm, "marker", edge, markerBack);
894:               }
895:               else if (vy == 0) {
896:                 DMSetLabelValue(dm, "marker", edge, markerFront);
897:               }
898:             }
899:             if (bdZ != DM_BOUNDARY_PERIODIC) {
900:               if (vz == numZVertices-1) {
901:                 DMSetLabelValue(dm, "marker", edge, markerTop);
902:               }
903:               else if (vz == 0) {
904:                 DMSetLabelValue(dm, "marker", edge, markerBottom);
905:               }
906:             }
907:           }
908:         }
909:       }
910:     }
911:     DMPlexSymmetrize(dm);
912:     DMPlexStratify(dm);
913:     /* Build coordinates */
914:     DMGetCoordinateSection(dm, &coordSection);
915:     PetscSectionSetNumFields(coordSection, 1);
916:     PetscSectionSetFieldComponents(coordSection, 0, dim);
917:     PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);
918:     for (v = firstVertex; v < firstVertex+numVertices; ++v) {
919:       PetscSectionSetDof(coordSection, v, dim);
920:       PetscSectionSetFieldDof(coordSection, v, 0, dim);
921:     }
922:     PetscSectionSetUp(coordSection);
923:     PetscSectionGetStorageSize(coordSection, &coordSize);
924:     VecCreate(PETSC_COMM_SELF, &coordinates);
925:     PetscObjectSetName((PetscObject) coordinates, "coordinates");
926:     VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
927:     VecSetBlockSize(coordinates, dim);
928:     VecSetType(coordinates,VECSTANDARD);
929:     VecGetArray(coordinates, &coords);
930:     for (vz = 0; vz < numZVertices; ++vz) {
931:       for (vy = 0; vy < numYVertices; ++vy) {
932:         for (vx = 0; vx < numXVertices; ++vx) {
933:           coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
934:           coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
935:           if (dim == 3) {
936:             coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz;
937:           }
938:         }
939:       }
940:     }
941:     VecRestoreArray(coordinates, &coords);
942:     DMSetCoordinatesLocal(dm, coordinates);
943:     VecDestroy(&coordinates);
944:   }
945:   return(0);
946: }

948: static PetscErrorCode DMPlexCreateBoxMesh_Tensor_Internal(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm)
949: {
950:   PetscInt       i;

955:   DMCreate(comm, dm);
957:   DMSetType(*dm, DMPLEX);
958:   DMSetDimension(*dm, dim);
959:   switch (dim) {
960:   case 2: {DMPlexCreateCubeMesh_Internal(*dm, lower, upper, faces, periodicity[0], periodicity[1], DM_BOUNDARY_NONE);break;}
961:   case 3: {DMPlexCreateCubeMesh_Internal(*dm, lower, upper, faces, periodicity[0], periodicity[1], periodicity[2]);break;}
962:   default: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
963:   }
964:   if (periodicity[0] == DM_BOUNDARY_PERIODIC || periodicity[0] == DM_BOUNDARY_TWIST ||
965:       periodicity[1] == DM_BOUNDARY_PERIODIC || periodicity[1] == DM_BOUNDARY_TWIST ||
966:       (dim > 2 && (periodicity[2] == DM_BOUNDARY_PERIODIC || periodicity[2] == DM_BOUNDARY_TWIST))) {
967:     PetscReal L[3];
968:     PetscReal maxCell[3];

970:     for (i = 0; i < dim; i++) {
971:       L[i]       = upper[i] - lower[i];
972:       maxCell[i] = 1.1 * (L[i] / PetscMax(1,faces[i]));
973:     }
974:     DMSetPeriodicity(*dm,PETSC_TRUE,maxCell,L,periodicity);
975:   }
976:   if (!interpolate) {
977:     DM udm;

979:     DMPlexUninterpolate(*dm, &udm);
980:     DMPlexCopyCoordinates(*dm, udm);
981:     DMDestroy(dm);
982:     *dm  = udm;
983:   }
984:   return(0);
985: }

987: /*@C
988:   DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices or tensor cells (hexahedra).

990:   Collective on MPI_Comm

992:   Input Parameters:
993: + comm        - The communicator for the DM object
994: . dim         - The spatial dimension
995: . simplex     - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells
996: . faces       - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D
997: . lower       - The lower left corner, or NULL for (0, 0, 0)
998: . upper       - The upper right corner, or NULL for (1, 1, 1)
999: . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE
1000: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

1002:   Output Parameter:
1003: . dm  - The DM object

1005:   Options Database Keys:
1006: + -dm_plex_box_lower <x,y,z> - Specify lower-left-bottom coordinates for the box
1007: - -dm_plex_box_upper <x,y,z> - Specify upper-right-top coordinates for the box

1009:   Note: Here is the numbering returned for 2 faces in each direction for tensor cells:
1010: $ 10---17---11---18----12
1011: $  |         |         |
1012: $  |         |         |
1013: $ 20    2   22    3    24
1014: $  |         |         |
1015: $  |         |         |
1016: $  7---15----8---16----9
1017: $  |         |         |
1018: $  |         |         |
1019: $ 19    0   21    1   23
1020: $  |         |         |
1021: $  |         |         |
1022: $  4---13----5---14----6

1024: and for simplicial cells

1026: $ 14----8---15----9----16
1027: $  |\     5  |\      7 |
1028: $  | \       | \       |
1029: $ 13   2    14    3    15
1030: $  | 4   \   | 6   \   |
1031: $  |       \ |       \ |
1032: $ 11----6---12----7----13
1033: $  |\        |\        |
1034: $  | \    1  | \     3 |
1035: $ 10   0    11    1    12
1036: $  | 0   \   | 2   \   |
1037: $  |       \ |       \ |
1038: $  8----4----9----5----10

1040:   Level: beginner

1042: .keywords: DM, create
1043: .seealso: DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate()
1044: @*/
1045: PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm)
1046: {
1047:   PetscInt       fac[3] = {0, 0, 0};
1048:   PetscReal      low[3] = {0, 0, 0};
1049:   PetscReal      upp[3] = {1, 1, 1};
1050:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
1051:   PetscInt       i, n;
1052:   PetscBool      flg;

1056:   for (i = 0; i < dim; ++i) fac[i] = faces ? faces[i] : (dim == 1 ? 1 : 4-dim);
1057:   if (lower) for (i = 0; i < dim; ++i) low[i] = lower[i];
1058:   if (upper) for (i = 0; i < dim; ++i) upp[i] = upper[i];
1059:   if (periodicity) for (i = 0; i < dim; ++i) bdt[i] = periodicity[i];
1060:   /* Allow bounds to be specified from the command line */
1061:   n    = 3;
1062:   PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_lower", low, &n, &flg);
1063:   if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim);
1064:   n    = 3;
1065:   PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_upper", upp, &n, &flg);
1066:   if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim);

1068:   if (dim == 1)      {DMPlexCreateLineMesh_Internal(comm, fac[0], low[0], upp[0], bdt[0], dm);}
1069:   else if (simplex)  {DMPlexCreateBoxMesh_Simplex_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);}
1070:   else               {DMPlexCreateBoxMesh_Tensor_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);}
1071:   return(0);
1072: }

1074: /*@
1075:   DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells.

1077:   Collective on MPI_Comm

1079:   Input Parameters:
1080: + comm        - The communicator for the DM object
1081: . faces       - Number of faces per dimension, or NULL for (1, 1, 1)
1082: . lower       - The lower left corner, or NULL for (0, 0, 0)
1083: . upper       - The upper right corner, or NULL for (1, 1, 1)
1084: . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE
1085: . ordExt      - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first
1086: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

1088:   Output Parameter:
1089: . dm  - The DM object

1091:   Level: beginner

1093: .keywords: DM, create
1094: .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1095: @*/
1096: PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool ordExt, PetscBool interpolate, DM *dm)
1097: {
1098:   DM             bdm, botdm;
1099:   PetscInt       i;
1100:   PetscInt       fac[3] = {0, 0, 0};
1101:   PetscReal      low[3] = {0, 0, 0};
1102:   PetscReal      upp[3] = {1, 1, 1};
1103:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};

1107:   for (i = 0; i < 3; ++i) fac[i] = faces ? (faces[i] > 0 ? faces[i] : 1) : 1;
1108:   if (lower) for (i = 0; i < 3; ++i) low[i] = lower[i];
1109:   if (upper) for (i = 0; i < 3; ++i) upp[i] = upper[i];
1110:   if (periodicity) for (i = 0; i < 3; ++i) bdt[i] = periodicity[i];
1111:   for (i = 0; i < 3; ++i) if (bdt[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity not yet supported");

1113:   DMCreate(comm, &bdm);
1114:   DMSetType(bdm, DMPLEX);
1115:   DMSetDimension(bdm, 1);
1116:   DMSetCoordinateDim(bdm, 2);
1117:   DMPlexCreateSquareBoundary(bdm, low, upp, fac);
1118:   DMPlexGenerate(bdm, NULL, PETSC_FALSE, &botdm);
1119:   DMDestroy(&bdm);
1120:   DMPlexExtrude(botdm, fac[2], upp[2] - low[2], ordExt, interpolate, dm);
1121:   if (low[2] != 0.0) {
1122:     Vec         v;
1123:     PetscScalar *x;
1124:     PetscInt    cDim, n;

1126:     DMGetCoordinatesLocal(*dm, &v);
1127:     VecGetBlockSize(v, &cDim);
1128:     VecGetLocalSize(v, &n);
1129:     VecGetArray(v, &x);
1130:     x   += cDim;
1131:     for (i=0; i<n; i+=cDim) x[i] += low[2];
1132:     VecRestoreArray(v,&x);
1133:     DMSetCoordinatesLocal(*dm, v);
1134:   }
1135:   DMDestroy(&botdm);
1136:   return(0);
1137: }

1139: /*@
1140:   DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells.

1142:   Collective on idm

1144:   Input Parameters:
1145: + idm - The mesh to be extruted
1146: . layers - The number of layers
1147: . height - The height of the extruded layer
1148: . ordExt - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first
1149: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

1151:   Output Parameter:
1152: . dm  - The DM object

1154:   Notes: The object created is an hybrid mesh, the vertex ordering in the cone of the cell is that of the prismatic cells

1156:   Level: advanced

1158: .keywords: DM, create
1159: .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMPlexSetHybridBounds(), DMSetType(), DMCreate()
1160: @*/
1161: PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool ordExt, PetscBool interpolate, DM* dm)
1162: {
1163:   PetscScalar       *coordsB;
1164:   const PetscScalar *coordsA;
1165:   PetscReal         *normals = NULL;
1166:   Vec               coordinatesA, coordinatesB;
1167:   PetscSection      coordSectionA, coordSectionB;
1168:   PetscInt          dim, cDim, cDimB, c, l, v, coordSize, *newCone;
1169:   PetscInt          cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices;
1170:   PetscErrorCode    ierr;

1177:   DMGetDimension(idm, &dim);
1178:   if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim);

1180:   DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);
1181:   DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);
1182:   numCells = (cEnd - cStart)*layers;
1183:   numVertices = (vEnd - vStart)*(layers+1);
1184:   DMCreate(PetscObjectComm((PetscObject)idm), dm);
1185:   DMSetType(*dm, DMPLEX);
1186:   DMSetDimension(*dm, dim+1);
1187:   DMPlexSetChart(*dm, 0, numCells+numVertices);
1188:   for (c = cStart, cellV = 0; c < cEnd; ++c) {
1189:     PetscInt *closure = NULL;
1190:     PetscInt closureSize, numCorners = 0;

1192:     DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);
1193:     for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++;
1194:     DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);
1195:     for (l = 0; l < layers; ++l) {
1196:       DMPlexSetConeSize(*dm, ordExt ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, 2*numCorners);
1197:     }
1198:     cellV = PetscMax(numCorners,cellV);
1199:   }
1200:   DMPlexSetHybridBounds(*dm, 0, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1201:   DMSetUp(*dm);

1203:   DMGetCoordinateDim(idm, &cDim);
1204:   if (dim != cDim) {
1205:     PetscCalloc1(cDim*(vEnd - vStart), &normals);
1206:   }
1207:   PetscMalloc1(3*cellV,&newCone);
1208:   for (c = cStart; c < cEnd; ++c) {
1209:     PetscInt *closure = NULL;
1210:     PetscInt closureSize, numCorners = 0, l;
1211:     PetscReal normal[3] = {0, 0, 0};

1213:     if (normals) {
1214:       DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);
1215:     }
1216:     DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);
1217:     for (v = 0; v < closureSize*2; v += 2) {
1218:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
1219:         PetscInt d;

1221:         newCone[numCorners++] = closure[v] - vStart;
1222:         if (normals) { for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d]; }
1223:       }
1224:     }
1225:     switch (numCorners) {
1226:     case 4: /* do nothing */
1227:     case 2: /* do nothing */
1228:       break;
1229:     case 3: /* from counter-clockwise to wedge ordering */
1230:       l = newCone[1];
1231:       newCone[1] = newCone[2];
1232:       newCone[2] = l;
1233:       break;
1234:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported number of corners: %D", numCorners);
1235:     }
1236:     DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);
1237:     for (l = 0; l < layers; ++l) {
1238:       PetscInt i;

1240:       for (i = 0; i < numCorners; ++i) {
1241:         newCone[  numCorners + i] = ordExt ? (layers+1)*newCone[i] + l     + numCells :     l*(vEnd - vStart) + newCone[i] + numCells;
1242:         newCone[2*numCorners + i] = ordExt ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells;
1243:       }
1244:       DMPlexSetCone(*dm, ordExt ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);
1245:     }
1246:   }
1247:   DMPlexSymmetrize(*dm);
1248:   DMPlexStratify(*dm);
1249:   PetscFree(newCone);

1251:   cDimB = cDim == dim ? cDim+1 : cDim;
1252:   DMGetCoordinateSection(*dm, &coordSectionB);
1253:   PetscSectionSetNumFields(coordSectionB, 1);
1254:   PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);
1255:   PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);
1256:   for (v = numCells; v < numCells+numVertices; ++v) {
1257:     PetscSectionSetDof(coordSectionB, v, cDimB);
1258:     PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);
1259:   }
1260:   PetscSectionSetUp(coordSectionB);
1261:   PetscSectionGetStorageSize(coordSectionB, &coordSize);
1262:   VecCreate(PETSC_COMM_SELF, &coordinatesB);
1263:   PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1264:   VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);
1265:   VecSetBlockSize(coordinatesB, cDimB);
1266:   VecSetType(coordinatesB,VECSTANDARD);

1268:   DMGetCoordinateSection(idm, &coordSectionA);
1269:   DMGetCoordinatesLocal(idm, &coordinatesA);
1270:   VecGetArray(coordinatesB, &coordsB);
1271:   VecGetArrayRead(coordinatesA, &coordsA);
1272:   for (v = vStart; v < vEnd; ++v) {
1273:     const PetscScalar *cptr;
1274:     PetscReal         ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.};
1275:     PetscReal         *normal, norm, h = height/layers;
1276:     PetscInt          offA, d, cDimA = cDim;

1278:     normal = normals ? normals + cDimB*(v - vStart) : (cDim > 1 ? ones3 : ones2);
1279:     if (normals) {
1280:       for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d];
1281:       for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm);
1282:     }

1284:     PetscSectionGetOffset(coordSectionA, v, &offA);
1285:     cptr = coordsA + offA;
1286:     for (l = 0; l < layers+1; ++l) {
1287:       PetscInt offB, d, newV;

1289:       newV = ordExt ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells;
1290:       PetscSectionGetOffset(coordSectionB, newV, &offB);
1291:       for (d = 0; d < cDimA; ++d) { coordsB[offB+d]  = cptr[d]; }
1292:       for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*h : 0.0; }
1293:       cptr    = coordsB + offB;
1294:       cDimA   = cDimB;
1295:     }
1296:   }
1297:   VecRestoreArrayRead(coordinatesA, &coordsA);
1298:   VecRestoreArray(coordinatesB, &coordsB);
1299:   DMSetCoordinatesLocal(*dm, coordinatesB);
1300:   VecDestroy(&coordinatesB);
1301:   PetscFree(normals);
1302:   if (interpolate) {
1303:     DM idm;

1305:     DMPlexInterpolate(*dm, &idm);
1306:     DMPlexCopyCoordinates(*dm, idm);
1307:     DMDestroy(dm);
1308:     *dm  = idm;
1309:   }
1310:   return(0);
1311: }

1313: /*@C
1314:   DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database.

1316:   Logically Collective on DM

1318:   Input Parameters:
1319: + dm - the DM context
1320: - prefix - the prefix to prepend to all option names

1322:   Notes:
1323:   A hyphen (-) must NOT be given at the beginning of the prefix name.
1324:   The first character of all runtime options is AUTOMATICALLY the hyphen.

1326:   Level: advanced

1328: .seealso: SNESSetFromOptions()
1329: @*/
1330: PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[])
1331: {
1332:   DM_Plex       *mesh = (DM_Plex *) dm->data;

1337:   PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);
1338:   PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);
1339:   return(0);
1340: }

1342: /*@
1343:   DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra.

1345:   Collective on MPI_Comm

1347:   Input Parameters:
1348: + comm      - The communicator for the DM object
1349: . numRefine - The number of regular refinements to the basic 5 cell structure
1350: - periodicZ - The boundary type for the Z direction

1352:   Output Parameter:
1353: . dm  - The DM object

1355:   Note: Here is the output numbering looking from the bottom of the cylinder:
1356: $       17-----14
1357: $        |     |
1358: $        |  2  |
1359: $        |     |
1360: $ 17-----8-----7-----14
1361: $  |     |     |     |
1362: $  |  3  |  0  |  1  |
1363: $  |     |     |     |
1364: $ 19-----5-----6-----13
1365: $        |     |
1366: $        |  4  |
1367: $        |     |
1368: $       19-----13
1369: $
1370: $ and up through the top
1371: $
1372: $       18-----16
1373: $        |     |
1374: $        |  2  |
1375: $        |     |
1376: $ 18----10----11-----16
1377: $  |     |     |     |
1378: $  |  3  |  0  |  1  |
1379: $  |     |     |     |
1380: $ 20-----9----12-----15
1381: $        |     |
1382: $        |  4  |
1383: $        |     |
1384: $       20-----15

1386:   Level: beginner

1388: .keywords: DM, create
1389: .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1390: @*/
1391: PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm)
1392: {
1393:   const PetscInt dim = 3;
1394:   PetscInt       numCells, numVertices, r;
1395:   PetscMPIInt    rank;

1400:   MPI_Comm_rank(comm, &rank);
1401:   if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine);
1402:   DMCreate(comm, dm);
1403:   DMSetType(*dm, DMPLEX);
1404:   DMSetDimension(*dm, dim);
1405:   /* Create topology */
1406:   {
1407:     PetscInt cone[8], c;

1409:     numCells    = !rank ?  5 : 0;
1410:     numVertices = !rank ? 16 : 0;
1411:     if (periodicZ == DM_BOUNDARY_PERIODIC) {
1412:       numCells   *= 3;
1413:       numVertices = !rank ? 24 : 0;
1414:     }
1415:     DMPlexSetChart(*dm, 0, numCells+numVertices);
1416:     for (c = 0; c < numCells; c++) {DMPlexSetConeSize(*dm, c, 8);}
1417:     DMSetUp(*dm);
1418:     if (!rank) {
1419:       if (periodicZ == DM_BOUNDARY_PERIODIC) {
1420:         cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16;
1421:         cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34;
1422:         DMPlexSetCone(*dm, 0, cone);
1423:         cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23;
1424:         cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */
1425:         DMPlexSetCone(*dm, 1, cone);
1426:         cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17;
1427:         cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38;
1428:         DMPlexSetCone(*dm, 2, cone);
1429:         cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15;
1430:         cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38;
1431:         DMPlexSetCone(*dm, 3, cone);
1432:         cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23;
1433:         cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31;
1434:         DMPlexSetCone(*dm, 4, cone);

1436:         cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32;
1437:         cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20;
1438:         DMPlexSetCone(*dm, 5, cone);
1439:         cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36;
1440:         cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21;
1441:         DMPlexSetCone(*dm, 6, cone);
1442:         cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33;
1443:         cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28;
1444:         DMPlexSetCone(*dm, 7, cone);
1445:         cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31;
1446:         cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28;
1447:         DMPlexSetCone(*dm, 8, cone);
1448:         cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36;
1449:         cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19;
1450:         DMPlexSetCone(*dm, 9, cone);

1452:         cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22;
1453:         cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18;
1454:         DMPlexSetCone(*dm, 10, cone);
1455:         cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25;
1456:         cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17;
1457:         DMPlexSetCone(*dm, 11, cone);
1458:         cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21;
1459:         cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27;
1460:         DMPlexSetCone(*dm, 12, cone);
1461:         cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19;
1462:         cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27;
1463:         DMPlexSetCone(*dm, 13, cone);
1464:         cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25;
1465:         cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15;
1466:         DMPlexSetCone(*dm, 14, cone);
1467:       } else {
1468:         cone[0] =  5; cone[1] =  8; cone[2] =  7; cone[3] =  6;
1469:         cone[4] =  9; cone[5] = 12; cone[6] = 11; cone[7] = 10;
1470:         DMPlexSetCone(*dm, 0, cone);
1471:         cone[0] =  6; cone[1] =  7; cone[2] = 14; cone[3] = 13;
1472:         cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11;
1473:         DMPlexSetCone(*dm, 1, cone);
1474:         cone[0] =  8; cone[1] = 17; cone[2] = 14; cone[3] =  7;
1475:         cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18;
1476:         DMPlexSetCone(*dm, 2, cone);
1477:         cone[0] = 19; cone[1] = 17; cone[2] =  8; cone[3] =  5;
1478:         cone[4] = 20; cone[5] =  9; cone[6] = 10; cone[7] = 18;
1479:         DMPlexSetCone(*dm, 3, cone);
1480:         cone[0] = 19; cone[1] =  5; cone[2] =  6; cone[3] = 13;
1481:         cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] =  9;
1482:         DMPlexSetCone(*dm, 4, cone);
1483:       }
1484:     }
1485:     DMPlexSymmetrize(*dm);
1486:     DMPlexStratify(*dm);
1487:   }
1488:   /* Interpolate */
1489:   {
1490:     DM idm;

1492:     DMPlexInterpolate(*dm, &idm);
1493:     DMDestroy(dm);
1494:     *dm  = idm;
1495:   }
1496:   /* Create cube geometry */
1497:   {
1498:     Vec             coordinates;
1499:     PetscSection    coordSection;
1500:     PetscScalar    *coords;
1501:     PetscInt        coordSize, v;
1502:     const PetscReal dis = 1.0/PetscSqrtReal(2.0);
1503:     const PetscReal ds2 = dis/2.0;

1505:     /* Build coordinates */
1506:     DMGetCoordinateSection(*dm, &coordSection);
1507:     PetscSectionSetNumFields(coordSection, 1);
1508:     PetscSectionSetFieldComponents(coordSection, 0, dim);
1509:     PetscSectionSetChart(coordSection, numCells, numCells+numVertices);
1510:     for (v = numCells; v < numCells+numVertices; ++v) {
1511:       PetscSectionSetDof(coordSection, v, dim);
1512:       PetscSectionSetFieldDof(coordSection, v, 0, dim);
1513:     }
1514:     PetscSectionSetUp(coordSection);
1515:     PetscSectionGetStorageSize(coordSection, &coordSize);
1516:     VecCreate(PETSC_COMM_SELF, &coordinates);
1517:     PetscObjectSetName((PetscObject) coordinates, "coordinates");
1518:     VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1519:     VecSetBlockSize(coordinates, dim);
1520:     VecSetType(coordinates,VECSTANDARD);
1521:     VecGetArray(coordinates, &coords);
1522:     if (!rank) {
1523:       coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0;
1524:       coords[1*dim+0] =  ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0;
1525:       coords[2*dim+0] =  ds2; coords[2*dim+1] =  ds2; coords[2*dim+2] = 0.0;
1526:       coords[3*dim+0] = -ds2; coords[3*dim+1] =  ds2; coords[3*dim+2] = 0.0;
1527:       coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0;
1528:       coords[5*dim+0] = -ds2; coords[5*dim+1] =  ds2; coords[5*dim+2] = 1.0;
1529:       coords[6*dim+0] =  ds2; coords[6*dim+1] =  ds2; coords[6*dim+2] = 1.0;
1530:       coords[7*dim+0] =  ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0;
1531:       coords[ 8*dim+0] =  dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0;
1532:       coords[ 9*dim+0] =  dis; coords[ 9*dim+1] =  dis; coords[ 9*dim+2] = 0.0;
1533:       coords[10*dim+0] =  dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0;
1534:       coords[11*dim+0] =  dis; coords[11*dim+1] =  dis; coords[11*dim+2] = 1.0;
1535:       coords[12*dim+0] = -dis; coords[12*dim+1] =  dis; coords[12*dim+2] = 0.0;
1536:       coords[13*dim+0] = -dis; coords[13*dim+1] =  dis; coords[13*dim+2] = 1.0;
1537:       coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0;
1538:       coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0;
1539:       if (periodicZ == DM_BOUNDARY_PERIODIC) {
1540:         /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5;
1541:         /* 16 32 22 */ coords[17*dim+0] =  ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5;
1542:         /* 17 33 21 */ coords[18*dim+0] =  ds2; coords[18*dim+1] =  ds2; coords[18*dim+2] = 0.5;
1543:         /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] =  ds2; coords[19*dim+2] = 0.5;
1544:         /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5;
1545:         /* 23 36 25 */ coords[21*dim+0] =  dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5;
1546:         /* 24 37 26 */ coords[22*dim+0] =  dis; coords[22*dim+1] =  dis; coords[22*dim+2] = 0.5;
1547:         /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] =  dis; coords[23*dim+2] = 0.5;
1548:       }
1549:     }
1550:     VecRestoreArray(coordinates, &coords);
1551:     DMSetCoordinatesLocal(*dm, coordinates);
1552:     VecDestroy(&coordinates);
1553:   }
1554:   /* Create periodicity */
1555:   if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) {
1556:     PetscReal      L[3];
1557:     PetscReal      maxCell[3];
1558:     DMBoundaryType bdType[3];
1559:     PetscReal      lower[3] = {0.0, 0.0, 0.0};
1560:     PetscReal      upper[3] = {1.0, 1.0, 1.5};
1561:     PetscInt       i, numZCells = 3;

1563:     bdType[0] = DM_BOUNDARY_NONE;
1564:     bdType[1] = DM_BOUNDARY_NONE;
1565:     bdType[2] = periodicZ;
1566:     for (i = 0; i < dim; i++) {
1567:       L[i]       = upper[i] - lower[i];
1568:       maxCell[i] = 1.1 * (L[i] / numZCells);
1569:     }
1570:     DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);
1571:   }
1572:   /* Refine topology */
1573:   for (r = 0; r < numRefine; ++r) {
1574:     DM rdm = NULL;

1576:     DMRefine(*dm, comm, &rdm);
1577:     DMDestroy(dm);
1578:     *dm  = rdm;
1579:   }
1580:   /* Remap geometry to cylinder
1581:        Interior square: Linear interpolation is correct
1582:        The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing
1583:        such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance

1585:          phi     = arctan(y/x)
1586:          d_close = sqrt(1/8 + 1/4 sin^2(phi))
1587:          d_far   = sqrt(1/2 + sin^2(phi))

1589:        so we remap them using

1591:          x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close)
1592:          y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close)

1594:        If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y.
1595:   */
1596:   {
1597:     Vec           coordinates;
1598:     PetscSection  coordSection;
1599:     PetscScalar  *coords;
1600:     PetscInt      vStart, vEnd, v;
1601:     const PetscReal dis = 1.0/PetscSqrtReal(2.0);
1602:     const PetscReal ds2 = 0.5*dis;

1604:     DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
1605:     DMGetCoordinateSection(*dm, &coordSection);
1606:     DMGetCoordinatesLocal(*dm, &coordinates);
1607:     VecGetArray(coordinates, &coords);
1608:     for (v = vStart; v < vEnd; ++v) {
1609:       PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc;
1610:       PetscInt  off;

1612:       PetscSectionGetOffset(coordSection, v, &off);
1613:       if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue;
1614:       x    = PetscRealPart(coords[off]);
1615:       y    = PetscRealPart(coords[off+1]);
1616:       phi  = PetscAtan2Real(y, x);
1617:       sinp = PetscSinReal(phi);
1618:       cosp = PetscCosReal(phi);
1619:       if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) {
1620:         dc = PetscAbsReal(ds2/sinp);
1621:         df = PetscAbsReal(dis/sinp);
1622:         xc = ds2*x/PetscAbsReal(y);
1623:         yc = ds2*PetscSignReal(y);
1624:       } else {
1625:         dc = PetscAbsReal(ds2/cosp);
1626:         df = PetscAbsReal(dis/cosp);
1627:         xc = ds2*PetscSignReal(x);
1628:         yc = ds2*y/PetscAbsReal(x);
1629:       }
1630:       coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc);
1631:       coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc);
1632:     }
1633:     VecRestoreArray(coordinates, &coords);
1634:     if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) {
1635:       DMLocalizeCoordinates(*dm);
1636:     }
1637:   }
1638:   return(0);
1639: }

1641: /*@
1642:   DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges.

1644:   Collective on MPI_Comm

1646:   Input Parameters:
1647: + comm - The communicator for the DM object
1648: . n    - The number of wedges around the origin
1649: - interpolate - Create edges and faces

1651:   Output Parameter:
1652: . dm  - The DM object

1654:   Level: beginner

1656: .keywords: DM, create
1657: .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1658: @*/
1659: PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm)
1660: {
1661:   const PetscInt dim = 3;
1662:   PetscInt       numCells, numVertices;
1663:   PetscMPIInt    rank;

1668:   MPI_Comm_rank(comm, &rank);
1669:   if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n);
1670:   DMCreate(comm, dm);
1671:   DMSetType(*dm, DMPLEX);
1672:   DMSetDimension(*dm, dim);
1673:   /* Create topology */
1674:   {
1675:     PetscInt cone[6], c;

1677:     numCells    = !rank ?        n : 0;
1678:     numVertices = !rank ?  2*(n+1) : 0;
1679:     DMPlexSetChart(*dm, 0, numCells+numVertices);
1680:     DMPlexSetHybridBounds(*dm, 0, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1681:     for (c = 0; c < numCells; c++) {DMPlexSetConeSize(*dm, c, 6);}
1682:     DMSetUp(*dm);
1683:     for (c = 0; c < numCells; c++) {
1684:       cone[0] =  c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n;
1685:       cone[3] =  c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n;
1686:       DMPlexSetCone(*dm, c, cone);
1687:     }
1688:     DMPlexSymmetrize(*dm);
1689:     DMPlexStratify(*dm);
1690:   }
1691:   /* Interpolate */
1692:   if (interpolate) {
1693:     DM idm;

1695:     DMPlexInterpolate(*dm, &idm);
1696:     DMDestroy(dm);
1697:     *dm  = idm;
1698:   }
1699:   /* Create cylinder geometry */
1700:   {
1701:     Vec          coordinates;
1702:     PetscSection coordSection;
1703:     PetscScalar *coords;
1704:     PetscInt     coordSize, v, c;

1706:     /* Build coordinates */
1707:     DMGetCoordinateSection(*dm, &coordSection);
1708:     PetscSectionSetNumFields(coordSection, 1);
1709:     PetscSectionSetFieldComponents(coordSection, 0, dim);
1710:     PetscSectionSetChart(coordSection, numCells, numCells+numVertices);
1711:     for (v = numCells; v < numCells+numVertices; ++v) {
1712:       PetscSectionSetDof(coordSection, v, dim);
1713:       PetscSectionSetFieldDof(coordSection, v, 0, dim);
1714:     }
1715:     PetscSectionSetUp(coordSection);
1716:     PetscSectionGetStorageSize(coordSection, &coordSize);
1717:     VecCreate(PETSC_COMM_SELF, &coordinates);
1718:     PetscObjectSetName((PetscObject) coordinates, "coordinates");
1719:     VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1720:     VecSetBlockSize(coordinates, dim);
1721:     VecSetType(coordinates,VECSTANDARD);
1722:     VecGetArray(coordinates, &coords);
1723:     for (c = 0; c < numCells; c++) {
1724:       coords[(c+0*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+2] = 1.0;
1725:       coords[(c+1*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+2] = 0.0;
1726:     }
1727:     if (!rank) {
1728:       coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0;
1729:       coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0;
1730:     }
1731:     VecRestoreArray(coordinates, &coords);
1732:     DMSetCoordinatesLocal(*dm, coordinates);
1733:     VecDestroy(&coordinates);
1734:   }
1735:   return(0);
1736: }

1738: PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[])
1739: {
1740:   PetscReal prod = 0.0;
1741:   PetscInt  i;
1742:   for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]);
1743:   return PetscSqrtReal(prod);
1744: }
1745: PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[])
1746: {
1747:   PetscReal prod = 0.0;
1748:   PetscInt  i;
1749:   for (i = 0; i < dim; ++i) prod += x[i]*y[i];
1750:   return prod;
1751: }

1753: /*@
1754:   DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d.

1756:   Collective on MPI_Comm

1758:   Input Parameters:
1759: . comm  - The communicator for the DM object
1760: . dim   - The dimension
1761: - simplex - Use simplices, or tensor product cells

1763:   Output Parameter:
1764: . dm  - The DM object

1766:   Level: beginner

1768: .keywords: DM, create
1769: .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1770: @*/
1771: PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *dm)
1772: {
1773:   const PetscInt  embedDim = dim+1;
1774:   PetscSection    coordSection;
1775:   Vec             coordinates;
1776:   PetscScalar    *coords;
1777:   PetscReal      *coordsIn;
1778:   PetscInt        numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e;
1779:   PetscMPIInt     rank;
1780:   PetscErrorCode  ierr;

1784:   DMCreate(comm, dm);
1785:   DMSetType(*dm, DMPLEX);
1786:   DMSetDimension(*dm, dim);
1787:   DMSetCoordinateDim(*dm, dim+1);
1788:   MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);
1789:   switch (dim) {
1790:   case 2:
1791:     if (simplex) {
1792:       DM              idm;
1793:       const PetscReal edgeLen     = 2.0/(1.0 + PETSC_PHI);
1794:       const PetscReal vertex[3]   = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)};
1795:       const PetscInt  degree      = 5;
1796:       PetscInt        s[3]        = {1, 1, 1};
1797:       PetscInt        cone[3];
1798:       PetscInt       *graph, p, i, j, k;

1800:       numCells    = !rank ? 20 : 0;
1801:       numVerts    = !rank ? 12 : 0;
1802:       firstVertex = numCells;
1803:       /* Use icosahedron, which for a unit sphere has coordinates which are all cyclic permutations of

1805:            (0, \pm 1/\phi+1, \pm \phi/\phi+1)

1807:          where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge
1808:          length is then given by 2/\phi = 2 * 2.73606 = 5.47214.
1809:        */
1810:       /* Construct vertices */
1811:       PetscCalloc1(numVerts * embedDim, &coordsIn);
1812:       for (p = 0, i = 0; p < embedDim; ++p) {
1813:         for (s[1] = -1; s[1] < 2; s[1] += 2) {
1814:           for (s[2] = -1; s[2] < 2; s[2] += 2) {
1815:             for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim];
1816:             ++i;
1817:           }
1818:         }
1819:       }
1820:       /* Construct graph */
1821:       PetscCalloc1(numVerts * numVerts, &graph);
1822:       for (i = 0; i < numVerts; ++i) {
1823:         for (j = 0, k = 0; j < numVerts; ++j) {
1824:           if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;}
1825:         }
1826:         if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree);
1827:       }
1828:       /* Build Topology */
1829:       DMPlexSetChart(*dm, 0, numCells+numVerts);
1830:       for (c = 0; c < numCells; c++) {
1831:         DMPlexSetConeSize(*dm, c, embedDim);
1832:       }
1833:       DMSetUp(*dm); /* Allocate space for cones */
1834:       /* Cells */
1835:       for (i = 0, c = 0; i < numVerts; ++i) {
1836:         for (j = 0; j < i; ++j) {
1837:           for (k = 0; k < j; ++k) {
1838:             if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) {
1839:               cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k;
1840:               /* Check orientation */
1841:               {
1842:                 const PetscInt epsilon[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
1843:                 PetscReal normal[3];
1844:                 PetscInt  e, f;

1846:                 for (d = 0; d < embedDim; ++d) {
1847:                   normal[d] = 0.0;
1848:                   for (e = 0; e < embedDim; ++e) {
1849:                     for (f = 0; f < embedDim; ++f) {
1850:                       normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]);
1851:                     }
1852:                   }
1853:                 }
1854:                 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;}
1855:               }
1856:               DMPlexSetCone(*dm, c++, cone);
1857:             }
1858:           }
1859:         }
1860:       }
1861:       DMPlexSymmetrize(*dm);
1862:       DMPlexStratify(*dm);
1863:       PetscFree(graph);
1864:       /* Interpolate mesh */
1865:       DMPlexInterpolate(*dm, &idm);
1866:       DMDestroy(dm);
1867:       *dm  = idm;
1868:     } else {
1869:       /*
1870:         12-21--13
1871:          |     |
1872:         25  4  24
1873:          |     |
1874:   12-25--9-16--8-24--13
1875:    |     |     |     |
1876:   23  5 17  0 15  3  22
1877:    |     |     |     |
1878:   10-20--6-14--7-19--11
1879:          |     |
1880:         20  1  19
1881:          |     |
1882:         10-18--11
1883:          |     |
1884:         23  2  22
1885:          |     |
1886:         12-21--13
1887:        */
1888:       const PetscReal dist = 1.0/PetscSqrtReal(3.0);
1889:       PetscInt        cone[4], ornt[4];

1891:       numCells    = !rank ?  6 : 0;
1892:       numEdges    = !rank ? 12 : 0;
1893:       numVerts    = !rank ?  8 : 0;
1894:       firstVertex = numCells;
1895:       firstEdge   = numCells + numVerts;
1896:       /* Build Topology */
1897:       DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);
1898:       for (c = 0; c < numCells; c++) {
1899:         DMPlexSetConeSize(*dm, c, 4);
1900:       }
1901:       for (e = firstEdge; e < firstEdge+numEdges; ++e) {
1902:         DMPlexSetConeSize(*dm, e, 2);
1903:       }
1904:       DMSetUp(*dm); /* Allocate space for cones */
1905:       /* Cell 0 */
1906:       cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17;
1907:       DMPlexSetCone(*dm, 0, cone);
1908:       ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0;
1909:       DMPlexSetConeOrientation(*dm, 0, ornt);
1910:       /* Cell 1 */
1911:       cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20;
1912:       DMPlexSetCone(*dm, 1, cone);
1913:       ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0;
1914:       DMPlexSetConeOrientation(*dm, 1, ornt);
1915:       /* Cell 2 */
1916:       cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23;
1917:       DMPlexSetCone(*dm, 2, cone);
1918:       ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0;
1919:       DMPlexSetConeOrientation(*dm, 2, ornt);
1920:       /* Cell 3 */
1921:       cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15;
1922:       DMPlexSetCone(*dm, 3, cone);
1923:       ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2;
1924:       DMPlexSetConeOrientation(*dm, 3, ornt);
1925:       /* Cell 4 */
1926:       cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25;
1927:       DMPlexSetCone(*dm, 4, cone);
1928:       ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0;
1929:       DMPlexSetConeOrientation(*dm, 4, ornt);
1930:       /* Cell 5 */
1931:       cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23;
1932:       DMPlexSetCone(*dm, 5, cone);
1933:       ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2;
1934:       DMPlexSetConeOrientation(*dm, 5, ornt);
1935:       /* Edges */
1936:       cone[0] =  6; cone[1] =  7;
1937:       DMPlexSetCone(*dm, 14, cone);
1938:       cone[0] =  7; cone[1] =  8;
1939:       DMPlexSetCone(*dm, 15, cone);
1940:       cone[0] =  8; cone[1] =  9;
1941:       DMPlexSetCone(*dm, 16, cone);
1942:       cone[0] =  9; cone[1] =  6;
1943:       DMPlexSetCone(*dm, 17, cone);
1944:       cone[0] = 10; cone[1] = 11;
1945:       DMPlexSetCone(*dm, 18, cone);
1946:       cone[0] = 11; cone[1] =  7;
1947:       DMPlexSetCone(*dm, 19, cone);
1948:       cone[0] =  6; cone[1] = 10;
1949:       DMPlexSetCone(*dm, 20, cone);
1950:       cone[0] = 12; cone[1] = 13;
1951:       DMPlexSetCone(*dm, 21, cone);
1952:       cone[0] = 13; cone[1] = 11;
1953:       DMPlexSetCone(*dm, 22, cone);
1954:       cone[0] = 10; cone[1] = 12;
1955:       DMPlexSetCone(*dm, 23, cone);
1956:       cone[0] = 13; cone[1] =  8;
1957:       DMPlexSetCone(*dm, 24, cone);
1958:       cone[0] = 12; cone[1] =  9;
1959:       DMPlexSetCone(*dm, 25, cone);
1960:       DMPlexSymmetrize(*dm);
1961:       DMPlexStratify(*dm);
1962:       /* Build coordinates */
1963:       PetscCalloc1(numVerts * embedDim, &coordsIn);
1964:       coordsIn[0*embedDim+0] = -dist; coordsIn[0*embedDim+1] =  dist; coordsIn[0*embedDim+2] = -dist;
1965:       coordsIn[1*embedDim+0] =  dist; coordsIn[1*embedDim+1] =  dist; coordsIn[1*embedDim+2] = -dist;
1966:       coordsIn[2*embedDim+0] =  dist; coordsIn[2*embedDim+1] = -dist; coordsIn[2*embedDim+2] = -dist;
1967:       coordsIn[3*embedDim+0] = -dist; coordsIn[3*embedDim+1] = -dist; coordsIn[3*embedDim+2] = -dist;
1968:       coordsIn[4*embedDim+0] = -dist; coordsIn[4*embedDim+1] =  dist; coordsIn[4*embedDim+2] =  dist;
1969:       coordsIn[5*embedDim+0] =  dist; coordsIn[5*embedDim+1] =  dist; coordsIn[5*embedDim+2] =  dist;
1970:       coordsIn[6*embedDim+0] = -dist; coordsIn[6*embedDim+1] = -dist; coordsIn[6*embedDim+2] =  dist;
1971:       coordsIn[7*embedDim+0] =  dist; coordsIn[7*embedDim+1] = -dist; coordsIn[7*embedDim+2] =  dist;
1972:     }
1973:     break;
1974:   case 3:
1975:     if (simplex) {
1976:       DM              idm;
1977:       const PetscReal edgeLen         = 1.0/PETSC_PHI;
1978:       const PetscReal vertexA[4]      = {0.5, 0.5, 0.5, 0.5};
1979:       const PetscReal vertexB[4]      = {1.0, 0.0, 0.0, 0.0};
1980:       const PetscReal vertexC[4]      = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0};
1981:       const PetscInt  degree          = 12;
1982:       PetscInt        s[4]            = {1, 1, 1};
1983:       PetscInt        evenPerm[12][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 0, 3, 2}, {1, 2, 0, 3}, {1, 3, 2, 0},
1984:                                          {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}};
1985:       PetscInt        cone[4];
1986:       PetscInt       *graph, p, i, j, k, l;

1988:       numCells    = !rank ? 600 : 0;
1989:       numVerts    = !rank ? 120 : 0;
1990:       firstVertex = numCells;
1991:       /* Use the 600-cell, which for a unit sphere has coordinates which are

1993:            1/2 (\pm 1, \pm 1,    \pm 1, \pm 1)                          16
1994:                (\pm 1,    0,       0,      0)  all cyclic permutations   8
1995:            1/2 (\pm 1, \pm phi, \pm 1/phi, 0)  all even permutations    96

1997:          where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge
1998:          length is then given by 1/\phi = 2.73606.

2000:          http://buzzard.pugetsound.edu/sage-practice/ch03s03.html
2001:          http://mathworld.wolfram.com/600-Cell.html
2002:        */
2003:       /* Construct vertices */
2004:       PetscCalloc1(numVerts * embedDim, &coordsIn);
2005:       i    = 0;
2006:       for (s[0] = -1; s[0] < 2; s[0] += 2) {
2007:         for (s[1] = -1; s[1] < 2; s[1] += 2) {
2008:           for (s[2] = -1; s[2] < 2; s[2] += 2) {
2009:             for (s[3] = -1; s[3] < 2; s[3] += 2) {
2010:               for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d];
2011:               ++i;
2012:             }
2013:           }
2014:         }
2015:       }
2016:       for (p = 0; p < embedDim; ++p) {
2017:         s[1] = s[2] = s[3] = 1;
2018:         for (s[0] = -1; s[0] < 2; s[0] += 2) {
2019:           for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim];
2020:           ++i;
2021:         }
2022:       }
2023:       for (p = 0; p < 12; ++p) {
2024:         s[3] = 1;
2025:         for (s[0] = -1; s[0] < 2; s[0] += 2) {
2026:           for (s[1] = -1; s[1] < 2; s[1] += 2) {
2027:             for (s[2] = -1; s[2] < 2; s[2] += 2) {
2028:               for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]];
2029:               ++i;
2030:             }
2031:           }
2032:         }
2033:       }
2034:       if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts);
2035:       /* Construct graph */
2036:       PetscCalloc1(numVerts * numVerts, &graph);
2037:       for (i = 0; i < numVerts; ++i) {
2038:         for (j = 0, k = 0; j < numVerts; ++j) {
2039:           if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;}
2040:         }
2041:         if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree);
2042:       }
2043:       /* Build Topology */
2044:       DMPlexSetChart(*dm, 0, numCells+numVerts);
2045:       for (c = 0; c < numCells; c++) {
2046:         DMPlexSetConeSize(*dm, c, embedDim);
2047:       }
2048:       DMSetUp(*dm); /* Allocate space for cones */
2049:       /* Cells */
2050:       for (i = 0, c = 0; i < numVerts; ++i) {
2051:         for (j = 0; j < i; ++j) {
2052:           for (k = 0; k < j; ++k) {
2053:             for (l = 0; l < k; ++l) {
2054:               if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] &&
2055:                   graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) {
2056:                 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l;
2057:                 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */
2058:                 {
2059:                   const PetscInt epsilon[4][4][4][4] = {{{{0,  0,  0,  0}, { 0, 0,  0,  0}, { 0,  0, 0,  0}, { 0,  0,  0, 0}},
2060:                                                          {{0,  0,  0,  0}, { 0, 0,  0,  0}, { 0,  0, 0,  1}, { 0,  0, -1, 0}},
2061:                                                          {{0,  0,  0,  0}, { 0, 0,  0, -1}, { 0,  0, 0,  0}, { 0,  1,  0, 0}},
2062:                                                          {{0,  0,  0,  0}, { 0, 0,  1,  0}, { 0, -1, 0,  0}, { 0,  0,  0, 0}}},

2064:                                                         {{{0,  0,  0,  0}, { 0, 0,  0,  0}, { 0,  0, 0, -1}, { 0,  0,  1, 0}},
2065:                                                          {{0,  0,  0,  0}, { 0, 0,  0,  0}, { 0,  0, 0,  0}, { 0,  0,  0, 0}},
2066:                                                          {{0,  0,  0,  1}, { 0, 0,  0,  0}, { 0,  0, 0,  0}, {-1,  0,  0, 0}},
2067:                                                          {{0,  0, -1,  0}, { 0, 0,  0,  0}, { 1,  0, 0,  0}, { 0,  0,  0, 0}}},

2069:                                                         {{{0,  0,  0,  0}, { 0, 0,  0,  1}, { 0,  0, 0,  0}, { 0, -1,  0, 0}},
2070:                                                          {{0,  0,  0, -1}, { 0, 0,  0,  0}, { 0,  0, 0,  0}, { 1,  0,  0, 0}},
2071:                                                          {{0,  0,  0,  0}, { 0, 0,  0,  0}, { 0,  0, 0,  0}, { 0,  0,  0, 0}},
2072:                                                          {{0,  1,  0,  0}, {-1, 0,  0,  0}, { 0,  0, 0,  0}, { 0,  0,  0, 0}}},

2074:                                                         {{{0,  0,  0,  0}, { 0, 0, -1,  0}, { 0,  1, 0,  0}, { 0,  0,  0, 0}},
2075:                                                          {{0,  0,  1,  0}, { 0, 0,  0,  0}, {-1,  0, 0,  0}, { 0,  0,  0, 0}},
2076:                                                          {{0, -1,  0,  0}, { 1, 0,  0,  0}, { 0,  0, 0,  0}, { 0,  0,  0, 0}},
2077:                                                          {{0,  0,  0,  0}, { 0, 0,  0,  0}, { 0,  0, 0,  0}, { 0,  0,  0, 0}}}};
2078:                   PetscReal normal[4];
2079:                   PetscInt  e, f, g;

2081:                   for (d = 0; d < embedDim; ++d) {
2082:                     normal[d] = 0.0;
2083:                     for (e = 0; e < embedDim; ++e) {
2084:                       for (f = 0; f < embedDim; ++f) {
2085:                         for (g = 0; g < embedDim; ++g) {
2086:                           normal[d] += epsilon[d][e][f][g]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f])*(coordsIn[l*embedDim+f] - coordsIn[i*embedDim+f]);
2087:                         }
2088:                       }
2089:                     }
2090:                   }
2091:                   if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;}
2092:                 }
2093:                 DMPlexSetCone(*dm, c++, cone);
2094:               }
2095:             }
2096:           }
2097:         }
2098:       }
2099:       DMPlexSymmetrize(*dm);
2100:       DMPlexStratify(*dm);
2101:       PetscFree(graph);
2102:       /* Interpolate mesh */
2103:       DMPlexInterpolate(*dm, &idm);
2104:       DMDestroy(dm);
2105:       *dm  = idm;
2106:       break;
2107:     }
2108:   default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim);
2109:   }
2110:   /* Create coordinates */
2111:   DMGetCoordinateSection(*dm, &coordSection);
2112:   PetscSectionSetNumFields(coordSection, 1);
2113:   PetscSectionSetFieldComponents(coordSection, 0, embedDim);
2114:   PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);
2115:   for (v = firstVertex; v < firstVertex+numVerts; ++v) {
2116:     PetscSectionSetDof(coordSection, v, embedDim);
2117:     PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
2118:   }
2119:   PetscSectionSetUp(coordSection);
2120:   PetscSectionGetStorageSize(coordSection, &coordSize);
2121:   VecCreate(PETSC_COMM_SELF, &coordinates);
2122:   VecSetBlockSize(coordinates, embedDim);
2123:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
2124:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
2125:   VecSetType(coordinates,VECSTANDARD);
2126:   VecGetArray(coordinates, &coords);
2127:   for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];}
2128:   VecRestoreArray(coordinates, &coords);
2129:   DMSetCoordinatesLocal(*dm, coordinates);
2130:   VecDestroy(&coordinates);
2131:   PetscFree(coordsIn);
2132:   return(0);
2133: }

2135: /* External function declarations here */
2136: extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
2137: extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
2138: extern PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
2139: extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm);
2140: extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
2141: extern PetscErrorCode DMCreateMatrix_Plex(DM dm,  Mat *J);
2142: extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
2143: extern PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
2144: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
2145: extern PetscErrorCode DMSetUp_Plex(DM dm);
2146: extern PetscErrorCode DMDestroy_Plex(DM dm);
2147: extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
2148: extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
2149: extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
2150: extern PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);
2151: static PetscErrorCode DMInitialize_Plex(DM dm);

2153: /* Replace dm with the contents of dmNew
2154:    - Share the DM_Plex structure
2155:    - Share the coordinates
2156:    - Share the SF
2157: */
2158: static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew)
2159: {
2160:   PetscSF               sf;
2161:   DM                    coordDM, coarseDM;
2162:   Vec                   coords;
2163:   PetscBool             isper;
2164:   const PetscReal      *maxCell, *L;
2165:   const DMBoundaryType *bd;
2166:   PetscErrorCode        ierr;

2169:   DMGetPointSF(dmNew, &sf);
2170:   DMSetPointSF(dm, sf);
2171:   DMGetCoordinateDM(dmNew, &coordDM);
2172:   DMGetCoordinatesLocal(dmNew, &coords);
2173:   DMSetCoordinateDM(dm, coordDM);
2174:   DMSetCoordinatesLocal(dm, coords);
2175:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
2176:   DMSetPeriodicity(dmNew, isper, maxCell, L, bd);
2177:   DMDestroy_Plex(dm);
2178:   DMInitialize_Plex(dm);
2179:   dm->data = dmNew->data;
2180:   ((DM_Plex *) dmNew->data)->refct++;
2181:   dmNew->labels->refct++;
2182:   if (!--(dm->labels->refct)) {
2183:     DMLabelLink next = dm->labels->next;

2185:     /* destroy the labels */
2186:     while (next) {
2187:       DMLabelLink tmp = next->next;

2189:       DMLabelDestroy(&next->label);
2190:       PetscFree(next);
2191:       next = tmp;
2192:     }
2193:     PetscFree(dm->labels);
2194:   }
2195:   dm->labels = dmNew->labels;
2196:   dm->depthLabel = dmNew->depthLabel;
2197:   DMGetCoarseDM(dmNew,&coarseDM);
2198:   DMSetCoarseDM(dm,coarseDM);
2199:   return(0);
2200: }

2202: /* Swap dm with the contents of dmNew
2203:    - Swap the DM_Plex structure
2204:    - Swap the coordinates
2205:    - Swap the point PetscSF
2206: */
2207: static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB)
2208: {
2209:   DM              coordDMA, coordDMB;
2210:   Vec             coordsA,  coordsB;
2211:   PetscSF         sfA,      sfB;
2212:   void            *tmp;
2213:   DMLabelLinkList listTmp;
2214:   DMLabel         depthTmp;
2215:   PetscInt        tmpI;
2216:   PetscErrorCode  ierr;

2219:   DMGetPointSF(dmA, &sfA);
2220:   DMGetPointSF(dmB, &sfB);
2221:   PetscObjectReference((PetscObject) sfA);
2222:   DMSetPointSF(dmA, sfB);
2223:   DMSetPointSF(dmB, sfA);
2224:   PetscObjectDereference((PetscObject) sfA);

2226:   DMGetCoordinateDM(dmA, &coordDMA);
2227:   DMGetCoordinateDM(dmB, &coordDMB);
2228:   PetscObjectReference((PetscObject) coordDMA);
2229:   DMSetCoordinateDM(dmA, coordDMB);
2230:   DMSetCoordinateDM(dmB, coordDMA);
2231:   PetscObjectDereference((PetscObject) coordDMA);

2233:   DMGetCoordinatesLocal(dmA, &coordsA);
2234:   DMGetCoordinatesLocal(dmB, &coordsB);
2235:   PetscObjectReference((PetscObject) coordsA);
2236:   DMSetCoordinatesLocal(dmA, coordsB);
2237:   DMSetCoordinatesLocal(dmB, coordsA);
2238:   PetscObjectDereference((PetscObject) coordsA);

2240:   tmp       = dmA->data;
2241:   dmA->data = dmB->data;
2242:   dmB->data = tmp;
2243:   listTmp   = dmA->labels;
2244:   dmA->labels = dmB->labels;
2245:   dmB->labels = listTmp;
2246:   depthTmp  = dmA->depthLabel;
2247:   dmA->depthLabel = dmB->depthLabel;
2248:   dmB->depthLabel = depthTmp;
2249:   tmpI         = dmA->levelup;
2250:   dmA->levelup = dmB->levelup;
2251:   dmB->levelup = tmpI;
2252:   return(0);
2253: }

2255: static PetscErrorCode DMPlexIsSimplex_Static(DM dm, PetscBool *isSimplex)
2256: {
2257:   PetscInt       dim, cStart, coneSize;

2261:   DMGetDimension(dm, &dim);
2262:   DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
2263:   DMPlexGetConeSize(dm, cStart, &coneSize);
2264:   *isSimplex = coneSize == dim+1 ? PETSC_TRUE : PETSC_FALSE;
2265:   return(0);
2266: }

2268: PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm)
2269: {
2270:   DM_Plex       *mesh = (DM_Plex*) dm->data;

2274:   /* Handle viewing */
2275:   PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);
2276:   PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL);
2277:   PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);
2278:   PetscOptionsInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL);
2279:   /* Point Location */
2280:   PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);
2281:   /* Partitioning and distribution */
2282:   PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);
2283:   /* Generation and remeshing */
2284:   PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);
2285:   /* Projection behavior */
2286:   PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);
2287:   PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);
2288:   /* Checking structure */
2289:   {
2290:     const char *cellTypes[] = {"simplex", "tensor", "unknown", "DMPlexCellType", "DM_PLEX_CELLTYPE_", NULL};
2291:     PetscEnum   ct;
2292:     PetscBool   flg = PETSC_FALSE, flg2 = PETSC_FALSE;

2294:     PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);
2295:     if (flg && flg2) {DMPlexCheckSymmetry(dm);}
2296:     PetscOptionsEnum("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices", "DMPlexCheckSkeleton", cellTypes, (PetscEnum) 0, &ct, &flg);
2297:     if (flg) {
2298:       PetscBool isSimplex = ct ? PETSC_FALSE : PETSC_TRUE;

2300:       if (ct == (PetscEnum) DM_PLEX_CELLTYPE_UNKNOWN) {DMPlexIsSimplex_Static(dm, &isSimplex);}
2301:       DMPlexCheckSkeleton(dm, isSimplex, 0);
2302:     }
2303:     PetscOptionsEnum("-dm_plex_check_faces", "Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type", "DMPlexCheckFaces", cellTypes, (PetscEnum) 0, &ct, &flg);
2304:     if (flg) {
2305:       PetscBool isSimplex = ct ? PETSC_FALSE : PETSC_TRUE;

2307:       if (ct == (PetscEnum) DM_PLEX_CELLTYPE_UNKNOWN) {DMPlexIsSimplex_Static(dm, &isSimplex);}
2308:       DMPlexCheckFaces(dm, isSimplex, 0);
2309:     }
2310:     PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);
2311:     if (flg && flg2) {DMPlexCheckGeometry(dm);}
2312:   }

2314:   PetscPartitionerSetFromOptions(mesh->partitioner);
2315:   return(0);
2316: }

2318: static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm)
2319: {
2320:   PetscInt       refine = 0, coarsen = 0, r;
2321:   PetscBool      isHierarchy;

2326:   PetscOptionsHead(PetscOptionsObject,"DMPlex Options");
2327:   /* Handle DMPlex refinement */
2328:   PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);
2329:   PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);
2330:   if (refine) {DMPlexSetRefinementUniform(dm, PETSC_TRUE);}
2331:   if (refine && isHierarchy) {
2332:     DM *dms, coarseDM;

2334:     DMGetCoarseDM(dm, &coarseDM);
2335:     PetscObjectReference((PetscObject)coarseDM);
2336:     PetscMalloc1(refine,&dms);
2337:     DMRefineHierarchy(dm, refine, dms);
2338:     /* Total hack since we do not pass in a pointer */
2339:     DMPlexSwap_Static(dm, dms[refine-1]);
2340:     if (refine == 1) {
2341:       DMSetCoarseDM(dm, dms[0]);
2342:       DMPlexSetRegularRefinement(dm, PETSC_TRUE);
2343:     } else {
2344:       DMSetCoarseDM(dm, dms[refine-2]);
2345:       DMPlexSetRegularRefinement(dm, PETSC_TRUE);
2346:       DMSetCoarseDM(dms[0], dms[refine-1]);
2347:       DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);
2348:     }
2349:     DMSetCoarseDM(dms[refine-1], coarseDM);
2350:     PetscObjectDereference((PetscObject)coarseDM);
2351:     /* Free DMs */
2352:     for (r = 0; r < refine; ++r) {
2353:       DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);
2354:       DMDestroy(&dms[r]);
2355:     }
2356:     PetscFree(dms);
2357:   } else {
2358:     for (r = 0; r < refine; ++r) {
2359:       DM refinedMesh;

2361:       DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
2362:       DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);
2363:       /* Total hack since we do not pass in a pointer */
2364:       DMPlexReplace_Static(dm, refinedMesh);
2365:       DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
2366:       DMDestroy(&refinedMesh);
2367:     }
2368:   }
2369:   /* Handle DMPlex coarsening */
2370:   PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);
2371:   PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);
2372:   if (coarsen && isHierarchy) {
2373:     DM *dms;

2375:     PetscMalloc1(coarsen, &dms);
2376:     DMCoarsenHierarchy(dm, coarsen, dms);
2377:     /* Free DMs */
2378:     for (r = 0; r < coarsen; ++r) {
2379:       DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);
2380:       DMDestroy(&dms[r]);
2381:     }
2382:     PetscFree(dms);
2383:   } else {
2384:     for (r = 0; r < coarsen; ++r) {
2385:       DM coarseMesh;

2387:       DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
2388:       DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);
2389:       /* Total hack since we do not pass in a pointer */
2390:       DMPlexReplace_Static(dm, coarseMesh);
2391:       DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
2392:       DMDestroy(&coarseMesh);
2393:     }
2394:   }
2395:   /* Handle */
2396:   DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
2397:   PetscOptionsTail();
2398:   return(0);
2399: }

2401: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
2402: {

2406:   DMCreateGlobalVector_Section_Private(dm,vec);
2407:   /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
2408:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);
2409:   VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);
2410:   VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);
2411:   VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);
2412:   return(0);
2413: }

2415: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
2416: {

2420:   DMCreateLocalVector_Section_Private(dm,vec);
2421:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);
2422:   VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);
2423:   return(0);
2424: }

2426: static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
2427: {
2428:   PetscInt       depth, d;

2432:   DMPlexGetDepth(dm, &depth);
2433:   if (depth == 1) {
2434:     DMGetDimension(dm, &d);
2435:     if (dim == 0)      {DMPlexGetDepthStratum(dm, dim, pStart, pEnd);}
2436:     else if (dim == d) {DMPlexGetDepthStratum(dm, 1, pStart, pEnd);}
2437:     else               {*pStart = 0; *pEnd = 0;}
2438:   } else {
2439:     DMPlexGetDepthStratum(dm, dim, pStart, pEnd);
2440:   }
2441:   return(0);
2442: }

2444: static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[])
2445: {
2446:   PetscSF        sf;

2450:   DMGetPointSF(dm, &sf);
2451:   PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);
2452:   return(0);
2453: }

2455: static PetscErrorCode DMHasCreateInjection_Plex(DM dm, PetscBool *flg)
2456: {
2460:   *flg = PETSC_TRUE;
2461:   return(0);
2462: }

2464: static PetscErrorCode DMInitialize_Plex(DM dm)
2465: {

2469:   dm->ops->view                            = DMView_Plex;
2470:   dm->ops->load                            = DMLoad_Plex;
2471:   dm->ops->setfromoptions                  = DMSetFromOptions_Plex;
2472:   dm->ops->clone                           = DMClone_Plex;
2473:   dm->ops->setup                           = DMSetUp_Plex;
2474:   dm->ops->createdefaultsection            = DMCreateDefaultSection_Plex;
2475:   dm->ops->createdefaultconstraints        = DMCreateDefaultConstraints_Plex;
2476:   dm->ops->createglobalvector              = DMCreateGlobalVector_Plex;
2477:   dm->ops->createlocalvector               = DMCreateLocalVector_Plex;
2478:   dm->ops->getlocaltoglobalmapping         = NULL;
2479:   dm->ops->createfieldis                   = NULL;
2480:   dm->ops->createcoordinatedm              = DMCreateCoordinateDM_Plex;
2481:   dm->ops->createcoordinatefield           = DMCreateCoordinateField_Plex;
2482:   dm->ops->getcoloring                     = NULL;
2483:   dm->ops->creatematrix                    = DMCreateMatrix_Plex;
2484:   dm->ops->createinterpolation             = DMCreateInterpolation_Plex;
2485:   dm->ops->createmassmatrix                = DMCreateMassMatrix_Plex;
2486:   dm->ops->getaggregates                   = NULL;
2487:   dm->ops->getinjection                    = DMCreateInjection_Plex;
2488:   dm->ops->hascreateinjection              = DMHasCreateInjection_Plex;
2489:   dm->ops->refine                          = DMRefine_Plex;
2490:   dm->ops->coarsen                         = DMCoarsen_Plex;
2491:   dm->ops->refinehierarchy                 = DMRefineHierarchy_Plex;
2492:   dm->ops->coarsenhierarchy                = DMCoarsenHierarchy_Plex;
2493:   dm->ops->adaptlabel                      = DMAdaptLabel_Plex;
2494:   dm->ops->adaptmetric                     = DMAdaptMetric_Plex;
2495:   dm->ops->globaltolocalbegin              = NULL;
2496:   dm->ops->globaltolocalend                = NULL;
2497:   dm->ops->localtoglobalbegin              = NULL;
2498:   dm->ops->localtoglobalend                = NULL;
2499:   dm->ops->destroy                         = DMDestroy_Plex;
2500:   dm->ops->createsubdm                     = DMCreateSubDM_Plex;
2501:   dm->ops->createsuperdm                   = DMCreateSuperDM_Plex;
2502:   dm->ops->getdimpoints                    = DMGetDimPoints_Plex;
2503:   dm->ops->locatepoints                    = DMLocatePoints_Plex;
2504:   dm->ops->projectfunctionlocal            = DMProjectFunctionLocal_Plex;
2505:   dm->ops->projectfunctionlabellocal       = DMProjectFunctionLabelLocal_Plex;
2506:   dm->ops->projectfieldlocal               = DMProjectFieldLocal_Plex;
2507:   dm->ops->projectfieldlabellocal          = DMProjectFieldLabelLocal_Plex;
2508:   dm->ops->computel2diff                   = DMComputeL2Diff_Plex;
2509:   dm->ops->computel2gradientdiff           = DMComputeL2GradientDiff_Plex;
2510:   dm->ops->computel2fielddiff              = DMComputeL2FieldDiff_Plex;
2511:   dm->ops->getneighbors                    = DMGetNeighors_Plex;
2512:   PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);
2513:   PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);
2514:   return(0);
2515: }

2517: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
2518: {
2519:   DM_Plex        *mesh = (DM_Plex *) dm->data;

2523:   mesh->refct++;
2524:   (*newdm)->data = mesh;
2525:   PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
2526:   DMInitialize_Plex(*newdm);
2527:   return(0);
2528: }

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

2536:   Options Database Keys:
2537: + -dm_plex_hash_location             - Use grid hashing for point location
2538: . -dm_plex_partition_balance         - Attempt to evenly divide points on partition boundary between processes
2539: . -dm_plex_remesh_bd                 - Allow changes to the boundary on remeshing
2540: . -dm_plex_max_projection_height     - Maxmimum mesh point height used to project locally
2541: . -dm_plex_regular_refinement        - Use special nested projection algorithm for regular refinement
2542: . -dm_plex_check_symmetry            - Check that the adjacency information in the mesh is symmetric
2543: . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices
2544: . -dm_plex_check_faces <celltype>    - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type
2545: . -dm_plex_check_geometry            - Check that cells have positive volume
2546: . -dm_view :mesh.tex:ascii_latex     - View the mesh in LaTeX/TikZ
2547: . -dm_plex_view_scale <num>          - Scale the TikZ
2548: - -dm_plex_print_fem <num>           - View FEM assembly information, such as element vectors and matrices


2551:   Level: intermediate

2553: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
2554: M*/

2556: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
2557: {
2558:   DM_Plex        *mesh;
2559:   PetscInt       unit, d;

2564:   PetscNewLog(dm,&mesh);
2565:   dm->dim  = 0;
2566:   dm->data = mesh;

2568:   mesh->refct             = 1;
2569:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
2570:   mesh->maxConeSize       = 0;
2571:   mesh->cones             = NULL;
2572:   mesh->coneOrientations  = NULL;
2573:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
2574:   mesh->maxSupportSize    = 0;
2575:   mesh->supports          = NULL;
2576:   mesh->refinementUniform = PETSC_TRUE;
2577:   mesh->refinementLimit   = -1.0;

2579:   mesh->facesTmp = NULL;

2581:   mesh->tetgenOpts   = NULL;
2582:   mesh->triangleOpts = NULL;
2583:   PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);
2584:   mesh->remeshBd     = PETSC_FALSE;

2586:   mesh->subpointMap = NULL;

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

2590:   mesh->regularRefinement   = PETSC_FALSE;
2591:   mesh->depthState          = -1;
2592:   mesh->globalVertexNumbers = NULL;
2593:   mesh->globalCellNumbers   = NULL;
2594:   mesh->anchorSection       = NULL;
2595:   mesh->anchorIS            = NULL;
2596:   mesh->createanchors       = NULL;
2597:   mesh->computeanchormatrix = NULL;
2598:   mesh->parentSection       = NULL;
2599:   mesh->parents             = NULL;
2600:   mesh->childIDs            = NULL;
2601:   mesh->childSection        = NULL;
2602:   mesh->children            = NULL;
2603:   mesh->referenceTree       = NULL;
2604:   mesh->getchildsymmetry    = NULL;
2605:   for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
2606:   mesh->vtkCellHeight       = 0;
2607:   mesh->useAnchors          = PETSC_FALSE;

2609:   mesh->maxProjectionHeight = 0;

2611:   mesh->printSetValues = PETSC_FALSE;
2612:   mesh->printFEM       = 0;
2613:   mesh->printTol       = 1.0e-10;

2615:   DMInitialize_Plex(dm);
2616:   return(0);
2617: }

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

2622:   Collective on MPI_Comm

2624:   Input Parameter:
2625: . comm - The communicator for the DMPlex object

2627:   Output Parameter:
2628: . mesh  - The DMPlex object

2630:   Level: beginner

2632: .keywords: DMPlex, create
2633: @*/
2634: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
2635: {

2640:   DMCreate(comm, mesh);
2641:   DMSetType(*mesh, DMPLEX);
2642:   return(0);
2643: }

2645: /*
2646:   This takes as input the common mesh generator output, a list of the vertices for each cell, but vertex numbers are global and an SF is built for them
2647: */
2648: /* TODO: invertCells and spaceDim arguments could be added also to to DMPlexCreateFromCellListParallel(), DMPlexBuildFromCellList_Internal() and DMPlexCreateFromCellList() */
2649: PetscErrorCode DMPlexBuildFromCellList_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells, PetscSF *sfVert)
2650: {
2651:   PetscSF         sfPoint;
2652:   PetscLayout     vLayout;
2653:   PetscHSetI      vhash;
2654:   PetscSFNode    *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex;
2655:   const PetscInt *vrange;
2656:   PetscInt        numVerticesAdj, off = 0, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g;
2657:   PetscMPIInt     rank, size;
2658:   PetscErrorCode  ierr;

2661:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2662:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
2663:   /* Partition vertices */
2664:   PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);
2665:   PetscLayoutSetLocalSize(vLayout, numVertices);
2666:   PetscLayoutSetBlockSize(vLayout, 1);
2667:   PetscLayoutSetUp(vLayout);
2668:   PetscLayoutGetRanges(vLayout, &vrange);
2669:   /* Count vertices and map them to procs */
2670:   PetscHSetICreate(&vhash);
2671:   for (c = 0; c < numCells; ++c) {
2672:     for (p = 0; p < numCorners; ++p) {
2673:       PetscHSetIAdd(vhash, cells[c*numCorners+p]);
2674:     }
2675:   }
2676:   PetscHSetIGetSize(vhash, &numVerticesAdj);
2677:   PetscMalloc1(numVerticesAdj, &verticesAdj);
2678:   PetscHSetIGetElems(vhash, &off, verticesAdj);
2679:   PetscHSetIDestroy(&vhash);
2680:   if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj);
2681:   PetscSortInt(numVerticesAdj, verticesAdj);
2682:   PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);
2683:   for (v = 0; v < numVerticesAdj; ++v) {
2684:     const PetscInt gv = verticesAdj[v];
2685:     PetscInt       vrank;

2687:     PetscFindInt(gv, size+1, vrange, &vrank);
2688:     vrank = vrank < 0 ? -(vrank+2) : vrank;
2689:     remoteVerticesAdj[v].index = gv - vrange[vrank];
2690:     remoteVerticesAdj[v].rank  = vrank;
2691:   }
2692:   /* Create cones */
2693:   DMPlexSetChart(dm, 0, numCells+numVerticesAdj);
2694:   for (c = 0; c < numCells; ++c) {DMPlexSetConeSize(dm, c, numCorners);}
2695:   DMSetUp(dm);
2696:   DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);
2697:   for (c = 0; c < numCells; ++c) {
2698:     for (p = 0; p < numCorners; ++p) {
2699:       const PetscInt gv = cells[c*numCorners+p];
2700:       PetscInt       lv;

2702:       PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);
2703:       if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv);
2704:       cone[p] = lv+numCells;
2705:     }
2706:     if (invertCells) { DMPlexInvertCell_Internal(spaceDim, numCorners, cone); }
2707:     DMPlexSetCone(dm, c, cone);
2708:   }
2709:   DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);
2710:   /* Create SF for vertices */
2711:   PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);
2712:   PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");
2713:   PetscSFSetFromOptions(*sfVert);
2714:   PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);
2715:   PetscFree(verticesAdj);
2716:   /* Build pointSF */
2717:   PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);
2718:   for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;}
2719:   for (v = 0; v < numVertices;    ++v) {vertexOwner[v].index = -1;         vertexOwner[v].rank = -1;}
2720:   PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);
2721:   PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);
2722:   for (v = 0; v < numVertices;    ++v) if (vertexOwner[v].rank < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global vertex %d on rank %d was unclaimed", v + vrange[rank], rank);
2723:   PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);
2724:   PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);
2725:   for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost;
2726:   PetscMalloc1(numVerticesGhost, &localVertex);
2727:   PetscMalloc1(numVerticesGhost, &remoteVertex);
2728:   for (v = 0, g = 0; v < numVerticesAdj; ++v) {
2729:     if (vertexLocal[v].rank != rank) {
2730:       localVertex[g]        = v+numCells;
2731:       remoteVertex[g].index = vertexLocal[v].index;
2732:       remoteVertex[g].rank  = vertexLocal[v].rank;
2733:       ++g;
2734:     }
2735:   }
2736:   PetscFree2(vertexLocal, vertexOwner);
2737:   if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost);
2738:   DMGetPointSF(dm, &sfPoint);
2739:   PetscObjectSetName((PetscObject) sfPoint, "point SF");
2740:   PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);
2741:   PetscLayoutDestroy(&vLayout);
2742:   /* Fill in the rest of the topology structure */
2743:   DMPlexSymmetrize(dm);
2744:   DMPlexStratify(dm);
2745:   return(0);
2746: }

2748: /*
2749:   This takes as input the coordinates for each owned vertex
2750: */
2751: PetscErrorCode DMPlexBuildCoordinates_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[])
2752: {
2753:   PetscSection   coordSection;
2754:   Vec            coordinates;
2755:   PetscScalar   *coords;
2756:   PetscInt       numVertices, numVerticesAdj, coordSize, v;

2760:   DMSetCoordinateDim(dm, spaceDim);
2761:   PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);
2762:   DMGetCoordinateSection(dm, &coordSection);
2763:   PetscSectionSetNumFields(coordSection, 1);
2764:   PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
2765:   PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);
2766:   for (v = numCells; v < numCells+numVerticesAdj; ++v) {
2767:     PetscSectionSetDof(coordSection, v, spaceDim);
2768:     PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
2769:   }
2770:   PetscSectionSetUp(coordSection);
2771:   PetscSectionGetStorageSize(coordSection, &coordSize);
2772:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
2773:   VecSetBlockSize(coordinates, spaceDim);
2774:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
2775:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
2776:   VecSetType(coordinates,VECSTANDARD);
2777:   VecGetArray(coordinates, &coords);
2778:   {
2779:     MPI_Datatype coordtype;

2781:     /* Need a temp buffer for coords if we have complex/single */
2782:     MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);
2783:     MPI_Type_commit(&coordtype);
2784: #if defined(PETSC_USE_COMPLEX)
2785:     {
2786:     PetscScalar *svertexCoords;
2787:     PetscInt    i;
2788:     PetscMalloc1(numV*spaceDim,&svertexCoords);
2789:     for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i];
2790:     PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);
2791:     PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);
2792:     PetscFree(svertexCoords);
2793:     }
2794: #else
2795:     PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);
2796:     PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);
2797: #endif
2798:     MPI_Type_free(&coordtype);
2799:   }
2800:   VecRestoreArray(coordinates, &coords);
2801:   DMSetCoordinatesLocal(dm, coordinates);
2802:   VecDestroy(&coordinates);
2803:   return(0);
2804: }

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

2809:   Input Parameters:
2810: + comm - The communicator
2811: . dim - The topological dimension of the mesh
2812: . numCells - The number of cells owned by this process
2813: . numVertices - The number of vertices owned by this process
2814: . numCorners - The number of vertices for each cell
2815: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
2816: . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell
2817: . spaceDim - The spatial dimension used for coordinates
2818: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex

2820:   Output Parameter:
2821: + dm - The DM
2822: - vertexSF - Optional, SF describing complete vertex ownership

2824:   Note: Two triangles sharing a face
2825: $
2826: $        2
2827: $      / | \
2828: $     /  |  \
2829: $    /   |   \
2830: $   0  0 | 1  3
2831: $    \   |   /
2832: $     \  |  /
2833: $      \ | /
2834: $        1
2835: would have input
2836: $  numCells = 2, numVertices = 4
2837: $  cells = [0 1 2  1 3 2]
2838: $
2839: which would result in the DMPlex
2840: $
2841: $        4
2842: $      / | \
2843: $     /  |  \
2844: $    /   |   \
2845: $   2  0 | 1  5
2846: $    \   |   /
2847: $     \  |  /
2848: $      \ | /
2849: $        3

2851:   Level: beginner

2853: .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate()
2854: @*/
2855: PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm)
2856: {
2857:   PetscSF        sfVert;

2861:   DMCreate(comm, dm);
2862:   DMSetType(*dm, DMPLEX);
2865:   DMSetDimension(*dm, dim);
2866:   DMPlexBuildFromCellList_Parallel_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE, &sfVert);
2867:   if (interpolate) {
2868:     DM idm;

2870:     DMPlexInterpolate(*dm, &idm);
2871:     DMDestroy(dm);
2872:     *dm  = idm;
2873:   }
2874:   DMPlexBuildCoordinates_Parallel_Internal(*dm, spaceDim, numCells, numVertices, sfVert, vertexCoords);
2875:   if (vertexSF) *vertexSF = sfVert;
2876:   else {PetscSFDestroy(&sfVert);}
2877:   return(0);
2878: }

2880: /*
2881:   This takes as input the common mesh generator output, a list of the vertices for each cell
2882: */
2883: PetscErrorCode DMPlexBuildFromCellList_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells)
2884: {
2885:   PetscInt      *cone, c, p;

2889:   DMPlexSetChart(dm, 0, numCells+numVertices);
2890:   for (c = 0; c < numCells; ++c) {
2891:     DMPlexSetConeSize(dm, c, numCorners);
2892:   }
2893:   DMSetUp(dm);
2894:   DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);
2895:   for (c = 0; c < numCells; ++c) {
2896:     for (p = 0; p < numCorners; ++p) {
2897:       cone[p] = cells[c*numCorners+p]+numCells;
2898:     }
2899:     if (invertCells) { DMPlexInvertCell_Internal(spaceDim, numCorners, cone); }
2900:     DMPlexSetCone(dm, c, cone);
2901:   }
2902:   DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);
2903:   DMPlexSymmetrize(dm);
2904:   DMPlexStratify(dm);
2905:   return(0);
2906: }

2908: /*
2909:   This takes as input the coordinates for each vertex
2910: */
2911: PetscErrorCode DMPlexBuildCoordinates_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
2912: {
2913:   PetscSection   coordSection;
2914:   Vec            coordinates;
2915:   DM             cdm;
2916:   PetscScalar   *coords;
2917:   PetscInt       v, d;

2921:   DMSetCoordinateDim(dm, spaceDim);
2922:   DMGetCoordinateSection(dm, &coordSection);
2923:   PetscSectionSetNumFields(coordSection, 1);
2924:   PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
2925:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
2926:   for (v = numCells; v < numCells+numVertices; ++v) {
2927:     PetscSectionSetDof(coordSection, v, spaceDim);
2928:     PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
2929:   }
2930:   PetscSectionSetUp(coordSection);

2932:   DMGetCoordinateDM(dm, &cdm);
2933:   DMCreateLocalVector(cdm, &coordinates);
2934:   VecSetBlockSize(coordinates, spaceDim);
2935:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
2936:   VecGetArray(coordinates, &coords);
2937:   for (v = 0; v < numVertices; ++v) {
2938:     for (d = 0; d < spaceDim; ++d) {
2939:       coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
2940:     }
2941:   }
2942:   VecRestoreArray(coordinates, &coords);
2943:   DMSetCoordinatesLocal(dm, coordinates);
2944:   VecDestroy(&coordinates);
2945:   return(0);
2946: }

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

2951:   Input Parameters:
2952: + comm - The communicator
2953: . dim - The topological dimension of the mesh
2954: . numCells - The number of cells
2955: . numVertices - The number of vertices
2956: . numCorners - The number of vertices for each cell
2957: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
2958: . cells - An array of numCells*numCorners numbers, the vertices for each cell
2959: . spaceDim - The spatial dimension used for coordinates
2960: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex

2962:   Output Parameter:
2963: . dm - The DM

2965:   Note: Two triangles sharing a face
2966: $
2967: $        2
2968: $      / | \
2969: $     /  |  \
2970: $    /   |   \
2971: $   0  0 | 1  3
2972: $    \   |   /
2973: $     \  |  /
2974: $      \ | /
2975: $        1
2976: would have input
2977: $  numCells = 2, numVertices = 4
2978: $  cells = [0 1 2  1 3 2]
2979: $
2980: which would result in the DMPlex
2981: $
2982: $        4
2983: $      / | \
2984: $     /  |  \
2985: $    /   |   \
2986: $   2  0 | 1  5
2987: $    \   |   /
2988: $     \  |  /
2989: $      \ | /
2990: $        3

2992:   Level: beginner

2994: .seealso: DMPlexCreateFromDAG(), DMPlexCreate()
2995: @*/
2996: 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)
2997: {

3001:   DMCreate(comm, dm);
3002:   DMSetType(*dm, DMPLEX);
3003:   DMSetDimension(*dm, dim);
3004:   DMPlexBuildFromCellList_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE);
3005:   if (interpolate) {
3006:     DM idm;

3008:     DMPlexInterpolate(*dm, &idm);
3009:     DMDestroy(dm);
3010:     *dm  = idm;
3011:   }
3012:   DMPlexBuildCoordinates_Internal(*dm, spaceDim, numCells, numVertices, vertexCoords);
3013:   return(0);
3014: }

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

3019:   Input Parameters:
3020: + dm - The empty DM object, usually from DMCreate() and DMSetDimension()
3021: . depth - The depth of the DAG
3022: . numPoints - Array of size depth + 1 containing the number of points at each depth
3023: . coneSize - The cone size of each point
3024: . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
3025: . coneOrientations - The orientation of each cone point
3026: - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim()

3028:   Output Parameter:
3029: . dm - The DM

3031:   Note: Two triangles sharing a face would have input
3032: $  depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
3033: $  cones = [2 3 4  3 5 4], coneOrientations = [0 0 0  0 0 0]
3034: $ vertexCoords = [-1.0 0.0  0.0 -1.0  0.0 1.0  1.0 0.0]
3035: $
3036: which would result in the DMPlex
3037: $
3038: $        4
3039: $      / | \
3040: $     /  |  \
3041: $    /   |   \
3042: $   2  0 | 1  5
3043: $    \   |   /
3044: $     \  |  /
3045: $      \ | /
3046: $        3
3047: $
3048: $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList()

3050:   Level: advanced

3052: .seealso: DMPlexCreateFromCellList(), DMPlexCreate()
3053: @*/
3054: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
3055: {
3056:   Vec            coordinates;
3057:   PetscSection   coordSection;
3058:   PetscScalar    *coords;
3059:   PetscInt       coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off;

3063:   DMGetDimension(dm, &dim);
3064:   DMGetCoordinateDim(dm, &dimEmbed);
3065:   if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim);
3066:   for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
3067:   DMPlexSetChart(dm, pStart, pEnd);
3068:   for (p = pStart; p < pEnd; ++p) {
3069:     DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
3070:     if (firstVertex < 0 && !coneSize[p - pStart]) {
3071:       firstVertex = p - pStart;
3072:     }
3073:   }
3074:   if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]);
3075:   DMSetUp(dm); /* Allocate space for cones */
3076:   for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
3077:     DMPlexSetCone(dm, p, &cones[off]);
3078:     DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
3079:   }
3080:   DMPlexSymmetrize(dm);
3081:   DMPlexStratify(dm);
3082:   /* Build coordinates */
3083:   DMGetCoordinateSection(dm, &coordSection);
3084:   PetscSectionSetNumFields(coordSection, 1);
3085:   PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
3086:   PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
3087:   for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
3088:     PetscSectionSetDof(coordSection, v, dimEmbed);
3089:     PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
3090:   }
3091:   PetscSectionSetUp(coordSection);
3092:   PetscSectionGetStorageSize(coordSection, &coordSize);
3093:   VecCreate(PETSC_COMM_SELF, &coordinates);
3094:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
3095:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
3096:   VecSetBlockSize(coordinates, dimEmbed);
3097:   VecSetType(coordinates,VECSTANDARD);
3098:   VecGetArray(coordinates, &coords);
3099:   for (v = 0; v < numPoints[0]; ++v) {
3100:     PetscInt off;

3102:     PetscSectionGetOffset(coordSection, v+firstVertex, &off);
3103:     for (d = 0; d < dimEmbed; ++d) {
3104:       coords[off+d] = vertexCoords[v*dimEmbed+d];
3105:     }
3106:   }
3107:   VecRestoreArray(coordinates, &coords);
3108:   DMSetCoordinatesLocal(dm, coordinates);
3109:   VecDestroy(&coordinates);
3110:   return(0);
3111: }

3113: /*@C
3114:   DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file.

3116: + comm        - The MPI communicator
3117: . filename    - Name of the .dat file
3118: - interpolate - Create faces and edges in the mesh

3120:   Output Parameter:
3121: . dm  - The DM object representing the mesh

3123:   Note: The format is the simplest possible:
3124: $ Ne
3125: $ v0 v1 ... vk
3126: $ Nv
3127: $ x y z marker

3129:   Level: beginner

3131: .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
3132: @*/
3133: PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
3134: {
3135:   DMLabel         marker;
3136:   PetscViewer     viewer;
3137:   Vec             coordinates;
3138:   PetscSection    coordSection;
3139:   PetscScalar    *coords;
3140:   char            line[PETSC_MAX_PATH_LEN];
3141:   PetscInt        dim = 3, cdim = 3, coordSize, v, c, d;
3142:   PetscMPIInt     rank;
3143:   int             snum, Nv, Nc;
3144:   PetscErrorCode  ierr;

3147:   MPI_Comm_rank(comm, &rank);
3148:   PetscViewerCreate(comm, &viewer);
3149:   PetscViewerSetType(viewer, PETSCVIEWERASCII);
3150:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
3151:   PetscViewerFileSetName(viewer, filename);
3152:   if (!rank) {
3153:     PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
3154:     snum = sscanf(line, "%d %d", &Nc, &Nv);
3155:     if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
3156:   } else {
3157:     Nc = Nv = 0;
3158:   }
3159:   DMCreate(comm, dm);
3160:   DMSetType(*dm, DMPLEX);
3161:   DMPlexSetChart(*dm, 0, Nc+Nv);
3162:   DMSetDimension(*dm, dim);
3163:   DMSetCoordinateDim(*dm, cdim);
3164:   /* Read topology */
3165:   if (!rank) {
3166:     PetscInt cone[8], corners = 8;
3167:     int      vbuf[8], v;

3169:     for (c = 0; c < Nc; ++c) {DMPlexSetConeSize(*dm, c, corners);}
3170:     DMSetUp(*dm);
3171:     for (c = 0; c < Nc; ++c) {
3172:       PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);
3173:       snum = sscanf(line, "%d %d %d %d %d %d %d %d", &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);
3174:       if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
3175:       for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc;
3176:       /* Hexahedra are inverted */
3177:       {
3178:         PetscInt tmp = cone[1];
3179:         cone[1] = cone[3];
3180:         cone[3] = tmp;
3181:       }
3182:       DMPlexSetCone(*dm, c, cone);
3183:     }
3184:   }
3185:   DMPlexSymmetrize(*dm);
3186:   DMPlexStratify(*dm);
3187:   /* Read coordinates */
3188:   DMGetCoordinateSection(*dm, &coordSection);
3189:   PetscSectionSetNumFields(coordSection, 1);
3190:   PetscSectionSetFieldComponents(coordSection, 0, cdim);
3191:   PetscSectionSetChart(coordSection, Nc, Nc + Nv);
3192:   for (v = Nc; v < Nc+Nv; ++v) {
3193:     PetscSectionSetDof(coordSection, v, cdim);
3194:     PetscSectionSetFieldDof(coordSection, v, 0, cdim);
3195:   }
3196:   PetscSectionSetUp(coordSection);
3197:   PetscSectionGetStorageSize(coordSection, &coordSize);
3198:   VecCreate(PETSC_COMM_SELF, &coordinates);
3199:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
3200:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
3201:   VecSetBlockSize(coordinates, cdim);
3202:   VecSetType(coordinates, VECSTANDARD);
3203:   VecGetArray(coordinates, &coords);
3204:   if (!rank) {
3205:     double x[3];
3206:     int    val;

3208:     DMCreateLabel(*dm, "marker");
3209:     DMGetLabel(*dm, "marker", &marker);
3210:     for (v = 0; v < Nv; ++v) {
3211:       PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);
3212:       snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val);
3213:       if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
3214:       for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d];
3215:       if (val) {DMLabelSetValue(marker, v+Nc, val);}
3216:     }
3217:   }
3218:   VecRestoreArray(coordinates, &coords);
3219:   DMSetCoordinatesLocal(*dm, coordinates);
3220:   VecDestroy(&coordinates);
3221:   PetscViewerDestroy(&viewer);
3222:   if (interpolate) {
3223:     DM      idm;
3224:     DMLabel bdlabel;

3226:     DMPlexInterpolate(*dm, &idm);
3227:     DMDestroy(dm);
3228:     *dm  = idm;

3230:     DMGetLabel(*dm, "marker", &bdlabel);
3231:     DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);
3232:     DMPlexLabelComplete(*dm, bdlabel);
3233:   }
3234:   return(0);
3235: }

3237: /*@C
3238:   DMPlexCreateFromFile - This takes a filename and produces a DM

3240:   Input Parameters:
3241: + comm - The communicator
3242: . filename - A file name
3243: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

3245:   Output Parameter:
3246: . dm - The DM

3248:   Options Database Keys:
3249: . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 

3251:   Level: beginner

3253: .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate()
3254: @*/
3255: PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
3256: {
3257:   const char    *extGmsh    = ".msh";
3258:   const char    *extGmsh2   = ".msh2";
3259:   const char    *extGmsh4   = ".msh4";
3260:   const char    *extCGNS    = ".cgns";
3261:   const char    *extExodus  = ".exo";
3262:   const char    *extGenesis = ".gen";
3263:   const char    *extFluent  = ".cas";
3264:   const char    *extHDF5    = ".h5";
3265:   const char    *extMed     = ".med";
3266:   const char    *extPLY     = ".ply";
3267:   const char    *extCV      = ".dat";
3268:   size_t         len;
3269:   PetscBool      isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV;
3270:   PetscMPIInt    rank;

3276:   MPI_Comm_rank(comm, &rank);
3277:   PetscStrlen(filename, &len);
3278:   if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path");
3279:   PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh,    4, &isGmsh);
3280:   PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2,   5, &isGmsh2);
3281:   PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4,   5, &isGmsh4);
3282:   PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS,    5, &isCGNS);
3283:   PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus,  4, &isExodus);
3284:   PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);
3285:   PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent,  4, &isFluent);
3286:   PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5,    3, &isHDF5);
3287:   PetscStrncmp(&filename[PetscMax(0,len-4)], extMed,     4, &isMed);
3288:   PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY,     4, &isPLY);
3289:   PetscStrncmp(&filename[PetscMax(0,len-4)], extCV,      4, &isCV);
3290:   if (isGmsh || isGmsh2 || isGmsh4) {
3291:     DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);
3292:   } else if (isCGNS) {
3293:     DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);
3294:   } else if (isExodus || isGenesis) {
3295:     DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);
3296:   } else if (isFluent) {
3297:     DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);
3298:   } else if (isHDF5) {
3299:     PetscBool      load_hdf5_xdmf = PETSC_FALSE;
3300:     PetscViewer viewer;

3302:     /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */
3303:     PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf",  5, &load_hdf5_xdmf);
3304:     PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);
3305:     PetscViewerCreate(comm, &viewer);
3306:     PetscViewerSetType(viewer, PETSCVIEWERHDF5);
3307:     PetscViewerFileSetMode(viewer, FILE_MODE_READ);
3308:     PetscViewerFileSetName(viewer, filename);
3309:     DMCreate(comm, dm);
3310:     DMSetType(*dm, DMPLEX);
3311:     if (load_hdf5_xdmf) {PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);}
3312:     DMLoad(*dm, viewer);
3313:     if (load_hdf5_xdmf) {PetscViewerPopFormat(viewer);}
3314:     PetscViewerDestroy(&viewer);

3316:     if (interpolate) {
3317:       DM idm;

3319:       DMPlexInterpolate(*dm, &idm);
3320:       DMDestroy(dm);
3321:       *dm  = idm;
3322:     }
3323:   } else if (isMed) {
3324:     DMPlexCreateMedFromFile(comm, filename, interpolate, dm);
3325:   } else if (isPLY) {
3326:     DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);
3327:   } else if (isCV) {
3328:     DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);
3329:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename);
3330:   return(0);
3331: }

3333: /*@
3334:   DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

3336:   Collective on comm

3338:   Input Parameters:
3339: + comm    - The communicator
3340: . dim     - The spatial dimension
3341: - simplex - Flag for simplex, otherwise use a tensor-product cell

3343:   Output Parameter:
3344: . refdm - The reference cell

3346:   Level: intermediate

3348: .keywords: reference cell
3349: .seealso:
3350: @*/
3351: PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm)
3352: {
3353:   DM             rdm;
3354:   Vec            coords;

3358:   DMCreate(comm, &rdm);
3359:   DMSetType(rdm, DMPLEX);
3360:   DMSetDimension(rdm, dim);
3361:   switch (dim) {
3362:   case 0:
3363:   {
3364:     PetscInt    numPoints[1]        = {1};
3365:     PetscInt    coneSize[1]         = {0};
3366:     PetscInt    cones[1]            = {0};
3367:     PetscInt    coneOrientations[1] = {0};
3368:     PetscScalar vertexCoords[1]     = {0.0};

3370:     DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3371:   }
3372:   break;
3373:   case 1:
3374:   {
3375:     PetscInt    numPoints[2]        = {2, 1};
3376:     PetscInt    coneSize[3]         = {2, 0, 0};
3377:     PetscInt    cones[2]            = {1, 2};
3378:     PetscInt    coneOrientations[2] = {0, 0};
3379:     PetscScalar vertexCoords[2]     = {-1.0,  1.0};

3381:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3382:   }
3383:   break;
3384:   case 2:
3385:     if (simplex) {
3386:       PetscInt    numPoints[2]        = {3, 1};
3387:       PetscInt    coneSize[4]         = {3, 0, 0, 0};
3388:       PetscInt    cones[3]            = {1, 2, 3};
3389:       PetscInt    coneOrientations[3] = {0, 0, 0};
3390:       PetscScalar vertexCoords[6]     = {-1.0, -1.0,  1.0, -1.0,  -1.0, 1.0};

3392:       DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3393:     } else {
3394:       PetscInt    numPoints[2]        = {4, 1};
3395:       PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
3396:       PetscInt    cones[4]            = {1, 2, 3, 4};
3397:       PetscInt    coneOrientations[4] = {0, 0, 0, 0};
3398:       PetscScalar vertexCoords[8]     = {-1.0, -1.0,  1.0, -1.0,  1.0, 1.0,  -1.0, 1.0};

3400:       DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3401:     }
3402:   break;
3403:   case 3:
3404:     if (simplex) {
3405:       PetscInt    numPoints[2]        = {4, 1};
3406:       PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
3407:       PetscInt    cones[4]            = {1, 3, 2, 4};
3408:       PetscInt    coneOrientations[4] = {0, 0, 0, 0};
3409:       PetscScalar vertexCoords[12]    = {-1.0, -1.0, -1.0,  1.0, -1.0, -1.0,  -1.0, 1.0, -1.0,  -1.0, -1.0, 1.0};

3411:       DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3412:     } else {
3413:       PetscInt    numPoints[2]        = {8, 1};
3414:       PetscInt    coneSize[9]         = {8, 0, 0, 0, 0, 0, 0, 0, 0};
3415:       PetscInt    cones[8]            = {1, 4, 3, 2, 5, 6, 7, 8};
3416:       PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
3417:       PetscScalar vertexCoords[24]    = {-1.0, -1.0, -1.0,  1.0, -1.0, -1.0,  1.0, 1.0, -1.0,  -1.0, 1.0, -1.0,
3418:                                          -1.0, -1.0,  1.0,  1.0, -1.0,  1.0,  1.0, 1.0,  1.0,  -1.0, 1.0,  1.0};

3420:       DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3421:     }
3422:   break;
3423:   default:
3424:     SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim);
3425:   }
3426:   DMPlexInterpolate(rdm, refdm);
3427:   if (rdm->coordinateDM) {
3428:     DM           ncdm;
3429:     PetscSection cs;
3430:     PetscInt     pEnd = -1;

3432:     DMGetSection(rdm->coordinateDM, &cs);
3433:     if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
3434:     if (pEnd >= 0) {
3435:       DMClone(rdm->coordinateDM, &ncdm);
3436:       DMCopyDisc(rdm->coordinateDM, ncdm);
3437:       DMSetSection(ncdm, cs);
3438:       DMSetCoordinateDM(*refdm, ncdm);
3439:       DMDestroy(&ncdm);
3440:     }
3441:   }
3442:   DMGetCoordinatesLocal(rdm, &coords);
3443:   if (coords) {
3444:     DMSetCoordinatesLocal(*refdm, coords);
3445:   } else {
3446:     DMGetCoordinates(rdm, &coords);
3447:     if (coords) {DMSetCoordinates(*refdm, coords);}
3448:   }
3449:   DMDestroy(&rdm);
3450:   return(0);
3451: }