Actual source code: plexcreate.c

  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>
  3: #include <petsc/private/hashseti.h>
  4: #include <petscsf.h>
  5: #include <petscdmplextransform.h>
  6: #include <petsc/private/kernels/blockmatmult.h>
  7: #include <petsc/private/kernels/blockinvert.h>

  9: PetscLogEvent DMPLEX_CreateFromFile, DMPLEX_BuildFromCellList, DMPLEX_BuildCoordinatesFromCellList;

 11: /* External function declarations here */
 12: static PetscErrorCode DMInitialize_Plex(DM dm);

 14: /* This copies internal things in the Plex structure that we generally want when making a new, related Plex */
 15: PetscErrorCode DMPlexCopy_Internal(DM dmin, PetscBool copyPeriodicity, PetscBool copyOverlap, DM dmout)
 16: {
 17:   const PetscReal         *maxCell, *Lstart, *L;
 18:   PetscBool                dist;
 19:   DMPlexReorderDefaultFlag reorder;

 21:   if (copyPeriodicity) {
 22:     DMGetPeriodicity(dmin, &maxCell, &Lstart, &L);
 23:     DMSetPeriodicity(dmout, maxCell, Lstart, L);
 24:   }
 25:   DMPlexDistributeGetDefault(dmin, &dist);
 26:   DMPlexDistributeSetDefault(dmout, dist);
 27:   DMPlexReorderGetDefault(dmin, &reorder);
 28:   DMPlexReorderSetDefault(dmout, reorder);
 29:   ((DM_Plex *)dmout->data)->useHashLocation = ((DM_Plex *)dmin->data)->useHashLocation;
 30:   if (copyOverlap) DMPlexSetOverlap_Plex(dmout, dmin, 0);
 31:   return 0;
 32: }

 34: /* Replace dm with the contents of ndm, and then destroy ndm
 35:    - Share the DM_Plex structure
 36:    - Share the coordinates
 37:    - Share the SF
 38: */
 39: PetscErrorCode DMPlexReplace_Internal(DM dm, DM *ndm)
 40: {
 41:   PetscSF          sf;
 42:   DM               dmNew = *ndm, coordDM, coarseDM;
 43:   Vec              coords;
 44:   const PetscReal *maxCell, *Lstart, *L;
 45:   PetscInt         dim, cdim;

 47:   if (dm == dmNew) {
 48:     DMDestroy(ndm);
 49:     return 0;
 50:   }
 51:   dm->setupcalled = dmNew->setupcalled;
 52:   DMGetDimension(dmNew, &dim);
 53:   DMSetDimension(dm, dim);
 54:   DMGetCoordinateDim(dmNew, &cdim);
 55:   DMSetCoordinateDim(dm, cdim);
 56:   DMGetPointSF(dmNew, &sf);
 57:   DMSetPointSF(dm, sf);
 58:   DMGetCoordinateDM(dmNew, &coordDM);
 59:   DMGetCoordinatesLocal(dmNew, &coords);
 60:   DMSetCoordinateDM(dm, coordDM);
 61:   DMSetCoordinatesLocal(dm, coords);
 62:   DMGetCellCoordinateDM(dmNew, &coordDM);
 63:   DMGetCellCoordinatesLocal(dmNew, &coords);
 64:   DMSetCellCoordinateDM(dm, coordDM);
 65:   DMSetCellCoordinatesLocal(dm, coords);
 66:   /* Do not want to create the coordinate field if it does not already exist, so do not call DMGetCoordinateField() */
 67:   DMFieldDestroy(&dm->coordinates[0].field);
 68:   dm->coordinates[0].field            = dmNew->coordinates[0].field;
 69:   ((DM_Plex *)dmNew->data)->coordFunc = ((DM_Plex *)dm->data)->coordFunc;
 70:   DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L);
 71:   DMSetPeriodicity(dm, maxCell, Lstart, L);
 72:   DMDestroy_Plex(dm);
 73:   DMInitialize_Plex(dm);
 74:   dm->data = dmNew->data;
 75:   ((DM_Plex *)dmNew->data)->refct++;
 76:   DMDestroyLabelLinkList_Internal(dm);
 77:   DMCopyLabels(dmNew, dm, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL);
 78:   DMGetCoarseDM(dmNew, &coarseDM);
 79:   DMSetCoarseDM(dm, coarseDM);
 80:   DMDestroy(ndm);
 81:   return 0;
 82: }

 84: /* Swap dm with the contents of dmNew
 85:    - Swap the DM_Plex structure
 86:    - Swap the coordinates
 87:    - Swap the point PetscSF
 88: */
 89: static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB)
 90: {
 91:   DM          coordDMA, coordDMB;
 92:   Vec         coordsA, coordsB;
 93:   PetscSF     sfA, sfB;
 94:   DMField     fieldTmp;
 95:   void       *tmp;
 96:   DMLabelLink listTmp;
 97:   DMLabel     depthTmp;
 98:   PetscInt    tmpI;

100:   if (dmA == dmB) return 0;
101:   DMGetPointSF(dmA, &sfA);
102:   DMGetPointSF(dmB, &sfB);
103:   PetscObjectReference((PetscObject)sfA);
104:   DMSetPointSF(dmA, sfB);
105:   DMSetPointSF(dmB, sfA);
106:   PetscObjectDereference((PetscObject)sfA);

108:   DMGetCoordinateDM(dmA, &coordDMA);
109:   DMGetCoordinateDM(dmB, &coordDMB);
110:   PetscObjectReference((PetscObject)coordDMA);
111:   DMSetCoordinateDM(dmA, coordDMB);
112:   DMSetCoordinateDM(dmB, coordDMA);
113:   PetscObjectDereference((PetscObject)coordDMA);

115:   DMGetCoordinatesLocal(dmA, &coordsA);
116:   DMGetCoordinatesLocal(dmB, &coordsB);
117:   PetscObjectReference((PetscObject)coordsA);
118:   DMSetCoordinatesLocal(dmA, coordsB);
119:   DMSetCoordinatesLocal(dmB, coordsA);
120:   PetscObjectDereference((PetscObject)coordsA);

122:   DMGetCellCoordinateDM(dmA, &coordDMA);
123:   DMGetCellCoordinateDM(dmB, &coordDMB);
124:   PetscObjectReference((PetscObject)coordDMA);
125:   DMSetCellCoordinateDM(dmA, coordDMB);
126:   DMSetCellCoordinateDM(dmB, coordDMA);
127:   PetscObjectDereference((PetscObject)coordDMA);

129:   DMGetCellCoordinatesLocal(dmA, &coordsA);
130:   DMGetCellCoordinatesLocal(dmB, &coordsB);
131:   PetscObjectReference((PetscObject)coordsA);
132:   DMSetCellCoordinatesLocal(dmA, coordsB);
133:   DMSetCellCoordinatesLocal(dmB, coordsA);
134:   PetscObjectDereference((PetscObject)coordsA);

136:   fieldTmp                  = dmA->coordinates[0].field;
137:   dmA->coordinates[0].field = dmB->coordinates[0].field;
138:   dmB->coordinates[0].field = fieldTmp;
139:   fieldTmp                  = dmA->coordinates[1].field;
140:   dmA->coordinates[1].field = dmB->coordinates[1].field;
141:   dmB->coordinates[1].field = fieldTmp;
142:   tmp                       = dmA->data;
143:   dmA->data                 = dmB->data;
144:   dmB->data                 = tmp;
145:   listTmp                   = dmA->labels;
146:   dmA->labels               = dmB->labels;
147:   dmB->labels               = listTmp;
148:   depthTmp                  = dmA->depthLabel;
149:   dmA->depthLabel           = dmB->depthLabel;
150:   dmB->depthLabel           = depthTmp;
151:   depthTmp                  = dmA->celltypeLabel;
152:   dmA->celltypeLabel        = dmB->celltypeLabel;
153:   dmB->celltypeLabel        = depthTmp;
154:   tmpI                      = dmA->levelup;
155:   dmA->levelup              = dmB->levelup;
156:   dmB->levelup              = tmpI;
157:   return 0;
158: }

160: static PetscErrorCode DMPlexInterpolateInPlace_Internal(DM dm)
161: {
162:   DM idm;

164:   DMPlexInterpolate(dm, &idm);
165:   DMPlexCopyCoordinates(dm, idm);
166:   DMPlexReplace_Internal(dm, &idm);
167:   return 0;
168: }

170: /*@C
171:   DMPlexCreateCoordinateSpace - Creates a finite element space for the coordinates

173:   Collective

175:   Input Parameters:
176: + DM        - The DM
177: . degree    - The degree of the finite element or PETSC_DECIDE
178: - coordFunc - An optional function to map new points from refinement to the surface

180:   Level: advanced

182: .seealso: `PetscFECreateLagrange()`, `DMGetCoordinateDM()`
183: @*/
184: PetscErrorCode DMPlexCreateCoordinateSpace(DM dm, PetscInt degree, PetscPointFunc coordFunc)
185: {
186:   DM_Plex     *mesh = (DM_Plex *)dm->data;
187:   DM           cdm;
188:   PetscDS      cds;
189:   PetscFE      fe;
190:   PetscClassId id;

192:   DMGetCoordinateDM(dm, &cdm);
193:   DMGetDS(cdm, &cds);
194:   PetscDSGetDiscretization(cds, 0, (PetscObject *)&fe);
195:   PetscObjectGetClassId((PetscObject)fe, &id);
196:   if (id != PETSCFE_CLASSID) {
197:     PetscBool simplex;
198:     PetscInt  dim, dE, qorder;

200:     DMGetDimension(dm, &dim);
201:     DMGetCoordinateDim(dm, &dE);
202:     qorder = degree;
203:     PetscObjectOptionsBegin((PetscObject)cdm);
204:     PetscOptionsBoundedInt("-coord_dm_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "DMPlexCreateCoordinateSpace", qorder, &qorder, NULL, 0);
205:     PetscOptionsEnd();
206:     if (degree == PETSC_DECIDE) fe = NULL;
207:     else {
208:       DMPlexIsSimplex(dm, &simplex);
209:       PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, degree, qorder, &fe);
210:     }
211:     DMProjectCoordinates(dm, fe);
212:     PetscFEDestroy(&fe);
213:   }
214:   mesh->coordFunc = coordFunc;
215:   return 0;
216: }

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

221:   Collective

223:   Input Parameters:
224: + comm - The communicator for the DM object
225: . dim - The spatial dimension
226: . simplex - Flag for simplicial cells, otherwise they are tensor product cells
227: . interpolate - Flag to create intermediate mesh pieces (edges, faces)
228: - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell

230:   Output Parameter:
231: . dm - The DM object

233:   Level: beginner

235: .seealso: `DMSetType()`, `DMCreate()`
236: @*/
237: PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscReal refinementLimit, DM *newdm)
238: {
239:   DM          dm;
240:   PetscMPIInt rank;

242:   DMCreate(comm, &dm);
243:   DMSetType(dm, DMPLEX);
244:   DMSetDimension(dm, dim);
245:   MPI_Comm_rank(comm, &rank);
246:   switch (dim) {
247:   case 2:
248:     if (simplex) PetscObjectSetName((PetscObject)dm, "triangular");
249:     else PetscObjectSetName((PetscObject)dm, "quadrilateral");
250:     break;
251:   case 3:
252:     if (simplex) PetscObjectSetName((PetscObject)dm, "tetrahedral");
253:     else PetscObjectSetName((PetscObject)dm, "hexahedral");
254:     break;
255:   default:
256:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %" PetscInt_FMT, dim);
257:   }
258:   if (rank) {
259:     PetscInt numPoints[2] = {0, 0};
260:     DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);
261:   } else {
262:     switch (dim) {
263:     case 2:
264:       if (simplex) {
265:         PetscInt    numPoints[2]        = {4, 2};
266:         PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
267:         PetscInt    cones[6]            = {2, 3, 4, 5, 4, 3};
268:         PetscInt    coneOrientations[6] = {0, 0, 0, 0, 0, 0};
269:         PetscScalar vertexCoords[8]     = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};

271:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
272:       } else {
273:         PetscInt    numPoints[2]        = {6, 2};
274:         PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
275:         PetscInt    cones[8]            = {2, 3, 4, 5, 3, 6, 7, 4};
276:         PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
277:         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};

279:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
280:       }
281:       break;
282:     case 3:
283:       if (simplex) {
284:         PetscInt    numPoints[2]        = {5, 2};
285:         PetscInt    coneSize[7]         = {4, 4, 0, 0, 0, 0, 0};
286:         PetscInt    cones[8]            = {4, 3, 5, 2, 5, 3, 4, 6};
287:         PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
288:         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};

290:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
291:       } else {
292:         PetscInt    numPoints[2]         = {12, 2};
293:         PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
294:         PetscInt    cones[16]            = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8};
295:         PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
296:         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, -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5};

298:         DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
299:       }
300:       break;
301:     default:
302:       SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %" PetscInt_FMT, dim);
303:     }
304:   }
305:   *newdm = dm;
306:   if (refinementLimit > 0.0) {
307:     DM          rdm;
308:     const char *name;

310:     DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);
311:     DMPlexSetRefinementLimit(*newdm, refinementLimit);
312:     DMRefine(*newdm, comm, &rdm);
313:     PetscObjectGetName((PetscObject)*newdm, &name);
314:     PetscObjectSetName((PetscObject)rdm, name);
315:     DMDestroy(newdm);
316:     *newdm = rdm;
317:   }
318:   if (interpolate) {
319:     DM idm;

321:     DMPlexInterpolate(*newdm, &idm);
322:     DMDestroy(newdm);
323:     *newdm = idm;
324:   }
325:   return 0;
326: }

328: static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
329: {
330:   const PetscInt numVertices    = 2;
331:   PetscInt       markerRight    = 1;
332:   PetscInt       markerLeft     = 1;
333:   PetscBool      markerSeparate = PETSC_FALSE;
334:   Vec            coordinates;
335:   PetscSection   coordSection;
336:   PetscScalar   *coords;
337:   PetscInt       coordSize;
338:   PetscMPIInt    rank;
339:   PetscInt       cdim = 1, v;

341:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
342:   if (markerSeparate) {
343:     markerRight = 2;
344:     markerLeft  = 1;
345:   }
346:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
347:   if (rank == 0) {
348:     DMPlexSetChart(dm, 0, numVertices);
349:     DMSetUp(dm); /* Allocate space for cones */
350:     DMSetLabelValue(dm, "marker", 0, markerLeft);
351:     DMSetLabelValue(dm, "marker", 1, markerRight);
352:   }
353:   DMPlexSymmetrize(dm);
354:   DMPlexStratify(dm);
355:   /* Build coordinates */
356:   DMSetCoordinateDim(dm, cdim);
357:   DMGetCoordinateSection(dm, &coordSection);
358:   PetscSectionSetNumFields(coordSection, 1);
359:   PetscSectionSetChart(coordSection, 0, numVertices);
360:   PetscSectionSetFieldComponents(coordSection, 0, cdim);
361:   for (v = 0; v < numVertices; ++v) {
362:     PetscSectionSetDof(coordSection, v, cdim);
363:     PetscSectionSetFieldDof(coordSection, v, 0, cdim);
364:   }
365:   PetscSectionSetUp(coordSection);
366:   PetscSectionGetStorageSize(coordSection, &coordSize);
367:   VecCreate(PETSC_COMM_SELF, &coordinates);
368:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
369:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
370:   VecSetBlockSize(coordinates, cdim);
371:   VecSetType(coordinates, VECSTANDARD);
372:   VecGetArray(coordinates, &coords);
373:   coords[0] = lower[0];
374:   coords[1] = upper[0];
375:   VecRestoreArray(coordinates, &coords);
376:   DMSetCoordinatesLocal(dm, coordinates);
377:   VecDestroy(&coordinates);
378:   return 0;
379: }

381: static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
382: {
383:   const PetscInt numVertices    = (edges[0] + 1) * (edges[1] + 1);
384:   const PetscInt numEdges       = edges[0] * (edges[1] + 1) + (edges[0] + 1) * edges[1];
385:   PetscInt       markerTop      = 1;
386:   PetscInt       markerBottom   = 1;
387:   PetscInt       markerRight    = 1;
388:   PetscInt       markerLeft     = 1;
389:   PetscBool      markerSeparate = PETSC_FALSE;
390:   Vec            coordinates;
391:   PetscSection   coordSection;
392:   PetscScalar   *coords;
393:   PetscInt       coordSize;
394:   PetscMPIInt    rank;
395:   PetscInt       v, vx, vy;

397:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
398:   if (markerSeparate) {
399:     markerTop    = 3;
400:     markerBottom = 1;
401:     markerRight  = 2;
402:     markerLeft   = 4;
403:   }
404:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
405:   if (rank == 0) {
406:     PetscInt e, ex, ey;

408:     DMPlexSetChart(dm, 0, numEdges + numVertices);
409:     for (e = 0; e < numEdges; ++e) DMPlexSetConeSize(dm, e, 2);
410:     DMSetUp(dm); /* Allocate space for cones */
411:     for (vx = 0; vx <= edges[0]; vx++) {
412:       for (ey = 0; ey < edges[1]; ey++) {
413:         PetscInt edge   = vx * edges[1] + ey + edges[0] * (edges[1] + 1);
414:         PetscInt vertex = ey * (edges[0] + 1) + vx + numEdges;
415:         PetscInt cone[2];

417:         cone[0] = vertex;
418:         cone[1] = vertex + edges[0] + 1;
419:         DMPlexSetCone(dm, edge, cone);
420:         if (vx == edges[0]) {
421:           DMSetLabelValue(dm, "marker", edge, markerRight);
422:           DMSetLabelValue(dm, "marker", cone[0], markerRight);
423:           if (ey == edges[1] - 1) {
424:             DMSetLabelValue(dm, "marker", cone[1], markerRight);
425:             DMSetLabelValue(dm, "Face Sets", cone[1], markerRight);
426:           }
427:         } else if (vx == 0) {
428:           DMSetLabelValue(dm, "marker", edge, markerLeft);
429:           DMSetLabelValue(dm, "marker", cone[0], markerLeft);
430:           if (ey == edges[1] - 1) {
431:             DMSetLabelValue(dm, "marker", cone[1], markerLeft);
432:             DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft);
433:           }
434:         }
435:       }
436:     }
437:     for (vy = 0; vy <= edges[1]; vy++) {
438:       for (ex = 0; ex < edges[0]; ex++) {
439:         PetscInt edge   = vy * edges[0] + ex;
440:         PetscInt vertex = vy * (edges[0] + 1) + ex + numEdges;
441:         PetscInt cone[2];

443:         cone[0] = vertex;
444:         cone[1] = vertex + 1;
445:         DMPlexSetCone(dm, edge, cone);
446:         if (vy == edges[1]) {
447:           DMSetLabelValue(dm, "marker", edge, markerTop);
448:           DMSetLabelValue(dm, "marker", cone[0], markerTop);
449:           if (ex == edges[0] - 1) {
450:             DMSetLabelValue(dm, "marker", cone[1], markerTop);
451:             DMSetLabelValue(dm, "Face Sets", cone[1], markerTop);
452:           }
453:         } else if (vy == 0) {
454:           DMSetLabelValue(dm, "marker", edge, markerBottom);
455:           DMSetLabelValue(dm, "marker", cone[0], markerBottom);
456:           if (ex == edges[0] - 1) {
457:             DMSetLabelValue(dm, "marker", cone[1], markerBottom);
458:             DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom);
459:           }
460:         }
461:       }
462:     }
463:   }
464:   DMPlexSymmetrize(dm);
465:   DMPlexStratify(dm);
466:   /* Build coordinates */
467:   DMSetCoordinateDim(dm, 2);
468:   DMGetCoordinateSection(dm, &coordSection);
469:   PetscSectionSetNumFields(coordSection, 1);
470:   PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);
471:   PetscSectionSetFieldComponents(coordSection, 0, 2);
472:   for (v = numEdges; v < numEdges + numVertices; ++v) {
473:     PetscSectionSetDof(coordSection, v, 2);
474:     PetscSectionSetFieldDof(coordSection, v, 0, 2);
475:   }
476:   PetscSectionSetUp(coordSection);
477:   PetscSectionGetStorageSize(coordSection, &coordSize);
478:   VecCreate(PETSC_COMM_SELF, &coordinates);
479:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
480:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
481:   VecSetBlockSize(coordinates, 2);
482:   VecSetType(coordinates, VECSTANDARD);
483:   VecGetArray(coordinates, &coords);
484:   for (vy = 0; vy <= edges[1]; ++vy) {
485:     for (vx = 0; vx <= edges[0]; ++vx) {
486:       coords[(vy * (edges[0] + 1) + vx) * 2 + 0] = lower[0] + ((upper[0] - lower[0]) / edges[0]) * vx;
487:       coords[(vy * (edges[0] + 1) + vx) * 2 + 1] = lower[1] + ((upper[1] - lower[1]) / edges[1]) * vy;
488:     }
489:   }
490:   VecRestoreArray(coordinates, &coords);
491:   DMSetCoordinatesLocal(dm, coordinates);
492:   VecDestroy(&coordinates);
493:   return 0;
494: }

496: static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
497: {
498:   PetscInt     vertices[3], numVertices;
499:   PetscInt     numFaces       = 2 * faces[0] * faces[1] + 2 * faces[1] * faces[2] + 2 * faces[0] * faces[2];
500:   PetscInt     markerTop      = 1;
501:   PetscInt     markerBottom   = 1;
502:   PetscInt     markerFront    = 1;
503:   PetscInt     markerBack     = 1;
504:   PetscInt     markerRight    = 1;
505:   PetscInt     markerLeft     = 1;
506:   PetscBool    markerSeparate = PETSC_FALSE;
507:   Vec          coordinates;
508:   PetscSection coordSection;
509:   PetscScalar *coords;
510:   PetscInt     coordSize;
511:   PetscMPIInt  rank;
512:   PetscInt     v, vx, vy, vz;
513:   PetscInt     voffset, iface = 0, cone[4];

516:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
517:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
518:   if (markerSeparate) {
519:     markerBottom = 1;
520:     markerTop    = 2;
521:     markerFront  = 3;
522:     markerBack   = 4;
523:     markerRight  = 5;
524:     markerLeft   = 6;
525:   }
526:   vertices[0] = faces[0] + 1;
527:   vertices[1] = faces[1] + 1;
528:   vertices[2] = faces[2] + 1;
529:   numVertices = vertices[0] * vertices[1] * vertices[2];
530:   if (rank == 0) {
531:     PetscInt f;

533:     DMPlexSetChart(dm, 0, numFaces + numVertices);
534:     for (f = 0; f < numFaces; ++f) DMPlexSetConeSize(dm, f, 4);
535:     DMSetUp(dm); /* Allocate space for cones */

537:     /* Side 0 (Top) */
538:     for (vy = 0; vy < faces[1]; vy++) {
539:       for (vx = 0; vx < faces[0]; vx++) {
540:         voffset = numFaces + vertices[0] * vertices[1] * (vertices[2] - 1) + vy * vertices[0] + vx;
541:         cone[0] = voffset;
542:         cone[1] = voffset + 1;
543:         cone[2] = voffset + vertices[0] + 1;
544:         cone[3] = voffset + vertices[0];
545:         DMPlexSetCone(dm, iface, cone);
546:         DMSetLabelValue(dm, "marker", iface, markerTop);
547:         DMSetLabelValue(dm, "marker", voffset + 0, markerTop);
548:         DMSetLabelValue(dm, "marker", voffset + 1, markerTop);
549:         DMSetLabelValue(dm, "marker", voffset + vertices[0] + 0, markerTop);
550:         DMSetLabelValue(dm, "marker", voffset + vertices[0] + 1, markerTop);
551:         iface++;
552:       }
553:     }

555:     /* Side 1 (Bottom) */
556:     for (vy = 0; vy < faces[1]; vy++) {
557:       for (vx = 0; vx < faces[0]; vx++) {
558:         voffset = numFaces + vy * (faces[0] + 1) + vx;
559:         cone[0] = voffset + 1;
560:         cone[1] = voffset;
561:         cone[2] = voffset + vertices[0];
562:         cone[3] = voffset + vertices[0] + 1;
563:         DMPlexSetCone(dm, iface, cone);
564:         DMSetLabelValue(dm, "marker", iface, markerBottom);
565:         DMSetLabelValue(dm, "marker", voffset + 0, markerBottom);
566:         DMSetLabelValue(dm, "marker", voffset + 1, markerBottom);
567:         DMSetLabelValue(dm, "marker", voffset + vertices[0] + 0, markerBottom);
568:         DMSetLabelValue(dm, "marker", voffset + vertices[0] + 1, markerBottom);
569:         iface++;
570:       }
571:     }

573:     /* Side 2 (Front) */
574:     for (vz = 0; vz < faces[2]; vz++) {
575:       for (vx = 0; vx < faces[0]; vx++) {
576:         voffset = numFaces + vz * vertices[0] * vertices[1] + vx;
577:         cone[0] = voffset;
578:         cone[1] = voffset + 1;
579:         cone[2] = voffset + vertices[0] * vertices[1] + 1;
580:         cone[3] = voffset + vertices[0] * vertices[1];
581:         DMPlexSetCone(dm, iface, cone);
582:         DMSetLabelValue(dm, "marker", iface, markerFront);
583:         DMSetLabelValue(dm, "marker", voffset + 0, markerFront);
584:         DMSetLabelValue(dm, "marker", voffset + 1, markerFront);
585:         DMSetLabelValue(dm, "marker", voffset + vertices[0] * vertices[1] + 0, markerFront);
586:         DMSetLabelValue(dm, "marker", voffset + vertices[0] * vertices[1] + 1, markerFront);
587:         iface++;
588:       }
589:     }

591:     /* Side 3 (Back) */
592:     for (vz = 0; vz < faces[2]; vz++) {
593:       for (vx = 0; vx < faces[0]; vx++) {
594:         voffset = numFaces + vz * vertices[0] * vertices[1] + vertices[0] * (vertices[1] - 1) + vx;
595:         cone[0] = voffset + vertices[0] * vertices[1];
596:         cone[1] = voffset + vertices[0] * vertices[1] + 1;
597:         cone[2] = voffset + 1;
598:         cone[3] = voffset;
599:         DMPlexSetCone(dm, iface, cone);
600:         DMSetLabelValue(dm, "marker", iface, markerBack);
601:         DMSetLabelValue(dm, "marker", voffset + 0, markerBack);
602:         DMSetLabelValue(dm, "marker", voffset + 1, markerBack);
603:         DMSetLabelValue(dm, "marker", voffset + vertices[0] * vertices[1] + 0, markerBack);
604:         DMSetLabelValue(dm, "marker", voffset + vertices[0] * vertices[1] + 1, markerBack);
605:         iface++;
606:       }
607:     }

609:     /* Side 4 (Left) */
610:     for (vz = 0; vz < faces[2]; vz++) {
611:       for (vy = 0; vy < faces[1]; vy++) {
612:         voffset = numFaces + vz * vertices[0] * vertices[1] + vy * vertices[0];
613:         cone[0] = voffset;
614:         cone[1] = voffset + vertices[0] * vertices[1];
615:         cone[2] = voffset + vertices[0] * vertices[1] + vertices[0];
616:         cone[3] = voffset + vertices[0];
617:         DMPlexSetCone(dm, iface, cone);
618:         DMSetLabelValue(dm, "marker", iface, markerLeft);
619:         DMSetLabelValue(dm, "marker", voffset + 0, markerLeft);
620:         DMSetLabelValue(dm, "marker", voffset + vertices[0] + 0, markerLeft);
621:         DMSetLabelValue(dm, "marker", voffset + vertices[1] + 0, markerLeft);
622:         DMSetLabelValue(dm, "marker", voffset + vertices[0] * vertices[1] + vertices[0], markerLeft);
623:         iface++;
624:       }
625:     }

627:     /* Side 5 (Right) */
628:     for (vz = 0; vz < faces[2]; vz++) {
629:       for (vy = 0; vy < faces[1]; vy++) {
630:         voffset = numFaces + vz * vertices[0] * vertices[1] + vy * vertices[0] + faces[0];
631:         cone[0] = voffset + vertices[0] * vertices[1];
632:         cone[1] = voffset;
633:         cone[2] = voffset + vertices[0];
634:         cone[3] = voffset + vertices[0] * vertices[1] + vertices[0];
635:         DMPlexSetCone(dm, iface, cone);
636:         DMSetLabelValue(dm, "marker", iface, markerRight);
637:         DMSetLabelValue(dm, "marker", voffset + 0, markerRight);
638:         DMSetLabelValue(dm, "marker", voffset + vertices[0] + 0, markerRight);
639:         DMSetLabelValue(dm, "marker", voffset + vertices[0] * vertices[1] + 0, markerRight);
640:         DMSetLabelValue(dm, "marker", voffset + vertices[0] * vertices[1] + vertices[0], markerRight);
641:         iface++;
642:       }
643:     }
644:   }
645:   DMPlexSymmetrize(dm);
646:   DMPlexStratify(dm);
647:   /* Build coordinates */
648:   DMSetCoordinateDim(dm, 3);
649:   DMGetCoordinateSection(dm, &coordSection);
650:   PetscSectionSetNumFields(coordSection, 1);
651:   PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);
652:   PetscSectionSetFieldComponents(coordSection, 0, 3);
653:   for (v = numFaces; v < numFaces + numVertices; ++v) {
654:     PetscSectionSetDof(coordSection, v, 3);
655:     PetscSectionSetFieldDof(coordSection, v, 0, 3);
656:   }
657:   PetscSectionSetUp(coordSection);
658:   PetscSectionGetStorageSize(coordSection, &coordSize);
659:   VecCreate(PETSC_COMM_SELF, &coordinates);
660:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
661:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
662:   VecSetBlockSize(coordinates, 3);
663:   VecSetType(coordinates, VECSTANDARD);
664:   VecGetArray(coordinates, &coords);
665:   for (vz = 0; vz <= faces[2]; ++vz) {
666:     for (vy = 0; vy <= faces[1]; ++vy) {
667:       for (vx = 0; vx <= faces[0]; ++vx) {
668:         coords[((vz * (faces[1] + 1) + vy) * (faces[0] + 1) + vx) * 3 + 0] = lower[0] + ((upper[0] - lower[0]) / faces[0]) * vx;
669:         coords[((vz * (faces[1] + 1) + vy) * (faces[0] + 1) + vx) * 3 + 1] = lower[1] + ((upper[1] - lower[1]) / faces[1]) * vy;
670:         coords[((vz * (faces[1] + 1) + vy) * (faces[0] + 1) + vx) * 3 + 2] = lower[2] + ((upper[2] - lower[2]) / faces[2]) * vz;
671:       }
672:     }
673:   }
674:   VecRestoreArray(coordinates, &coords);
675:   DMSetCoordinatesLocal(dm, coordinates);
676:   VecDestroy(&coordinates);
677:   return 0;
678: }

680: static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate)
681: {
683:   DMSetDimension(dm, dim - 1);
684:   DMSetCoordinateDim(dm, dim);
685:   switch (dim) {
686:   case 1:
687:     DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(dm, lower, upper, faces);
688:     break;
689:   case 2:
690:     DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(dm, lower, upper, faces);
691:     break;
692:   case 3:
693:     DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(dm, lower, upper, faces);
694:     break;
695:   default:
696:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Dimension not supported: %" PetscInt_FMT, dim);
697:   }
698:   if (interpolate) DMPlexInterpolateInPlace_Internal(dm);
699:   return 0;
700: }

702: /*@C
703:   DMPlexCreateBoxSurfaceMesh - Creates a mesh on the surface of the tensor product of unit intervals (box) using tensor cells (hexahedra).

705:   Collective

707:   Input Parameters:
708: + comm        - The communicator for the DM object
709: . dim         - The spatial dimension of the box, so the resulting mesh is has dimension dim-1
710: . faces       - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D
711: . lower       - The lower left corner, or NULL for (0, 0, 0)
712: . upper       - The upper right corner, or NULL for (1, 1, 1)
713: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

715:   Output Parameter:
716: . dm  - The DM object

718:   Level: beginner

720: .seealso: `DMSetFromOptions()`, `DMPlexCreateBoxMesh()`, `DMPlexCreateFromFile()`, `DMSetType()`, `DMCreate()`
721: @*/
722: PetscErrorCode DMPlexCreateBoxSurfaceMesh(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate, DM *dm)
723: {
724:   PetscInt  fac[3] = {1, 1, 1};
725:   PetscReal low[3] = {0, 0, 0};
726:   PetscReal upp[3] = {1, 1, 1};

728:   DMCreate(comm, dm);
729:   DMSetType(*dm, DMPLEX);
730:   DMPlexCreateBoxSurfaceMesh_Internal(*dm, dim, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, interpolate);
731:   return 0;
732: }

734: static PetscErrorCode DMPlexCreateLineMesh_Internal(DM dm, PetscInt segments, PetscReal lower, PetscReal upper, DMBoundaryType bd)
735: {
736:   PetscInt     i, fStart, fEnd, numCells = 0, numVerts = 0;
737:   PetscInt     numPoints[2], *coneSize, *cones, *coneOrientations;
738:   PetscScalar *vertexCoords;
739:   PetscReal    L, maxCell;
740:   PetscBool    markerSeparate = PETSC_FALSE;
741:   PetscInt     markerLeft = 1, faceMarkerLeft = 1;
742:   PetscInt     markerRight = 1, faceMarkerRight = 2;
743:   PetscBool    wrap = (bd == DM_BOUNDARY_PERIODIC || bd == DM_BOUNDARY_TWIST) ? PETSC_TRUE : PETSC_FALSE;
744:   PetscMPIInt  rank;


748:   DMSetDimension(dm, 1);
749:   DMCreateLabel(dm, "marker");
750:   DMCreateLabel(dm, "Face Sets");

752:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
753:   if (rank == 0) numCells = segments;
754:   if (rank == 0) numVerts = segments + (wrap ? 0 : 1);

756:   numPoints[0] = numVerts;
757:   numPoints[1] = numCells;
758:   PetscMalloc4(numCells + numVerts, &coneSize, numCells * 2, &cones, numCells + numVerts, &coneOrientations, numVerts, &vertexCoords);
759:   PetscArrayzero(coneOrientations, numCells + numVerts);
760:   for (i = 0; i < numCells; ++i) coneSize[i] = 2;
761:   for (i = 0; i < numVerts; ++i) coneSize[numCells + i] = 0;
762:   for (i = 0; i < numCells; ++i) {
763:     cones[2 * i]     = numCells + i % numVerts;
764:     cones[2 * i + 1] = numCells + (i + 1) % numVerts;
765:   }
766:   for (i = 0; i < numVerts; ++i) vertexCoords[i] = lower + (upper - lower) * ((PetscReal)i / (PetscReal)numCells);
767:   DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
768:   PetscFree4(coneSize, cones, coneOrientations, vertexCoords);

770:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
771:   if (markerSeparate) {
772:     markerLeft  = faceMarkerLeft;
773:     markerRight = faceMarkerRight;
774:   }
775:   if (!wrap && rank == 0) {
776:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
777:     DMSetLabelValue(dm, "marker", fStart, markerLeft);
778:     DMSetLabelValue(dm, "marker", fEnd - 1, markerRight);
779:     DMSetLabelValue(dm, "Face Sets", fStart, faceMarkerLeft);
780:     DMSetLabelValue(dm, "Face Sets", fEnd - 1, faceMarkerRight);
781:   }
782:   if (wrap) {
783:     L       = upper - lower;
784:     maxCell = (PetscReal)1.1 * (L / (PetscReal)PetscMax(1, segments));
785:     DMSetPeriodicity(dm, &maxCell, &lower, &L);
786:   }
787:   DMPlexSetRefinementUniform(dm, PETSC_TRUE);
788:   return 0;
789: }

791: static PetscErrorCode DMPlexCreateBoxMesh_Simplex_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate)
792: {
793:   DM      boundary, vol;
794:   DMLabel bdlabel;

798:   DMCreate(PetscObjectComm((PetscObject)dm), &boundary);
799:   DMSetType(boundary, DMPLEX);
800:   DMPlexCreateBoxSurfaceMesh_Internal(boundary, dim, faces, lower, upper, PETSC_FALSE);
801:   DMPlexGenerate(boundary, NULL, interpolate, &vol);
802:   DMGetLabel(vol, "marker", &bdlabel);
803:   if (bdlabel) DMPlexLabelComplete(vol, bdlabel);
804:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_FALSE, vol);
805:   DMPlexReplace_Internal(dm, &vol);
806:   DMDestroy(&boundary);
807:   return 0;
808: }

810: static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ)
811: {
812:   DMLabel     cutLabel  = NULL;
813:   PetscInt    markerTop = 1, faceMarkerTop = 1;
814:   PetscInt    markerBottom = 1, faceMarkerBottom = 1;
815:   PetscInt    markerFront = 1, faceMarkerFront = 1;
816:   PetscInt    markerBack = 1, faceMarkerBack = 1;
817:   PetscInt    markerRight = 1, faceMarkerRight = 1;
818:   PetscInt    markerLeft = 1, faceMarkerLeft = 1;
819:   PetscInt    dim;
820:   PetscBool   markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE;
821:   PetscMPIInt rank;

823:   DMGetDimension(dm, &dim);
824:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
825:   DMCreateLabel(dm, "marker");
826:   DMCreateLabel(dm, "Face Sets");
827:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);
828:   if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST || bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST || bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) {
829:     if (cutMarker) {
830:       DMCreateLabel(dm, "periodic_cut");
831:       DMGetLabel(dm, "periodic_cut", &cutLabel);
832:     }
833:   }
834:   switch (dim) {
835:   case 2:
836:     faceMarkerTop    = 3;
837:     faceMarkerBottom = 1;
838:     faceMarkerRight  = 2;
839:     faceMarkerLeft   = 4;
840:     break;
841:   case 3:
842:     faceMarkerBottom = 1;
843:     faceMarkerTop    = 2;
844:     faceMarkerFront  = 3;
845:     faceMarkerBack   = 4;
846:     faceMarkerRight  = 5;
847:     faceMarkerLeft   = 6;
848:     break;
849:   default:
850:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", dim);
851:   }
852:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
853:   if (markerSeparate) {
854:     markerBottom = faceMarkerBottom;
855:     markerTop    = faceMarkerTop;
856:     markerFront  = faceMarkerFront;
857:     markerBack   = faceMarkerBack;
858:     markerRight  = faceMarkerRight;
859:     markerLeft   = faceMarkerLeft;
860:   }
861:   {
862:     const PetscInt numXEdges    = rank == 0 ? edges[0] : 0;
863:     const PetscInt numYEdges    = rank == 0 ? edges[1] : 0;
864:     const PetscInt numZEdges    = rank == 0 ? edges[2] : 0;
865:     const PetscInt numXVertices = rank == 0 ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0] + 1) : 0;
866:     const PetscInt numYVertices = rank == 0 ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1] + 1) : 0;
867:     const PetscInt numZVertices = rank == 0 ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2] + 1) : 0;
868:     const PetscInt numCells     = numXEdges * numYEdges * numZEdges;
869:     const PetscInt numXFaces    = numYEdges * numZEdges;
870:     const PetscInt numYFaces    = numXEdges * numZEdges;
871:     const PetscInt numZFaces    = numXEdges * numYEdges;
872:     const PetscInt numTotXFaces = numXVertices * numXFaces;
873:     const PetscInt numTotYFaces = numYVertices * numYFaces;
874:     const PetscInt numTotZFaces = numZVertices * numZFaces;
875:     const PetscInt numFaces     = numTotXFaces + numTotYFaces + numTotZFaces;
876:     const PetscInt numTotXEdges = numXEdges * numYVertices * numZVertices;
877:     const PetscInt numTotYEdges = numYEdges * numXVertices * numZVertices;
878:     const PetscInt numTotZEdges = numZEdges * numXVertices * numYVertices;
879:     const PetscInt numVertices  = numXVertices * numYVertices * numZVertices;
880:     const PetscInt numEdges     = numTotXEdges + numTotYEdges + numTotZEdges;
881:     const PetscInt firstVertex  = (dim == 2) ? numFaces : numCells;
882:     const PetscInt firstXFace   = (dim == 2) ? 0 : numCells + numVertices;
883:     const PetscInt firstYFace   = firstXFace + numTotXFaces;
884:     const PetscInt firstZFace   = firstYFace + numTotYFaces;
885:     const PetscInt firstXEdge   = numCells + numFaces + numVertices;
886:     const PetscInt firstYEdge   = firstXEdge + numTotXEdges;
887:     const PetscInt firstZEdge   = firstYEdge + numTotYEdges;
888:     Vec            coordinates;
889:     PetscSection   coordSection;
890:     PetscScalar   *coords;
891:     PetscInt       coordSize;
892:     PetscInt       v, vx, vy, vz;
893:     PetscInt       c, f, fx, fy, fz, e, ex, ey, ez;

895:     DMPlexSetChart(dm, 0, numCells + numFaces + numEdges + numVertices);
896:     for (c = 0; c < numCells; c++) DMPlexSetConeSize(dm, c, 6);
897:     for (f = firstXFace; f < firstXFace + numFaces; ++f) DMPlexSetConeSize(dm, f, 4);
898:     for (e = firstXEdge; e < firstXEdge + numEdges; ++e) DMPlexSetConeSize(dm, e, 2);
899:     DMSetUp(dm); /* Allocate space for cones */
900:     /* Build cells */
901:     for (fz = 0; fz < numZEdges; ++fz) {
902:       for (fy = 0; fy < numYEdges; ++fy) {
903:         for (fx = 0; fx < numXEdges; ++fx) {
904:           PetscInt cell  = (fz * numYEdges + fy) * numXEdges + fx;
905:           PetscInt faceB = firstZFace + (fy * numXEdges + fx) * numZVertices + fz;
906:           PetscInt faceT = firstZFace + (fy * numXEdges + fx) * numZVertices + ((fz + 1) % numZVertices);
907:           PetscInt faceF = firstYFace + (fz * numXEdges + fx) * numYVertices + fy;
908:           PetscInt faceK = firstYFace + (fz * numXEdges + fx) * numYVertices + ((fy + 1) % numYVertices);
909:           PetscInt faceL = firstXFace + (fz * numYEdges + fy) * numXVertices + fx;
910:           PetscInt faceR = firstXFace + (fz * numYEdges + fy) * numXVertices + ((fx + 1) % numXVertices);
911:           /* B,  T,  F,  K,  R,  L */
912:           PetscInt ornt[6] = {-2, 0, 0, -3, 0, -2}; /* ??? */
913:           PetscInt cone[6];

915:           /* no boundary twisting in 3D */
916:           cone[0] = faceB;
917:           cone[1] = faceT;
918:           cone[2] = faceF;
919:           cone[3] = faceK;
920:           cone[4] = faceR;
921:           cone[5] = faceL;
922:           DMPlexSetCone(dm, cell, cone);
923:           DMPlexSetConeOrientation(dm, cell, ornt);
924:           if (bdX != DM_BOUNDARY_NONE && fx == numXEdges - 1 && cutLabel) DMLabelSetValue(cutLabel, cell, 2);
925:           if (bdY != DM_BOUNDARY_NONE && fy == numYEdges - 1 && cutLabel) DMLabelSetValue(cutLabel, cell, 2);
926:           if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges - 1 && cutLabel) DMLabelSetValue(cutLabel, cell, 2);
927:         }
928:       }
929:     }
930:     /* Build x faces */
931:     for (fz = 0; fz < numZEdges; ++fz) {
932:       for (fy = 0; fy < numYEdges; ++fy) {
933:         for (fx = 0; fx < numXVertices; ++fx) {
934:           PetscInt face    = firstXFace + (fz * numYEdges + fy) * numXVertices + fx;
935:           PetscInt edgeL   = firstZEdge + (fy * numXVertices + fx) * numZEdges + fz;
936:           PetscInt edgeR   = firstZEdge + (((fy + 1) % numYVertices) * numXVertices + fx) * numZEdges + fz;
937:           PetscInt edgeB   = firstYEdge + (fz * numXVertices + fx) * numYEdges + fy;
938:           PetscInt edgeT   = firstYEdge + (((fz + 1) % numZVertices) * numXVertices + fx) * numYEdges + fy;
939:           PetscInt ornt[4] = {0, 0, -1, -1};
940:           PetscInt cone[4];

942:           if (dim == 3) {
943:             /* markers */
944:             if (bdX != DM_BOUNDARY_PERIODIC) {
945:               if (fx == numXVertices - 1) {
946:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);
947:                 DMSetLabelValue(dm, "marker", face, markerRight);
948:               } else if (fx == 0) {
949:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);
950:                 DMSetLabelValue(dm, "marker", face, markerLeft);
951:               }
952:             }
953:           }
954:           cone[0] = edgeB;
955:           cone[1] = edgeR;
956:           cone[2] = edgeT;
957:           cone[3] = edgeL;
958:           DMPlexSetCone(dm, face, cone);
959:           DMPlexSetConeOrientation(dm, face, ornt);
960:         }
961:       }
962:     }
963:     /* Build y faces */
964:     for (fz = 0; fz < numZEdges; ++fz) {
965:       for (fx = 0; fx < numXEdges; ++fx) {
966:         for (fy = 0; fy < numYVertices; ++fy) {
967:           PetscInt face    = firstYFace + (fz * numXEdges + fx) * numYVertices + fy;
968:           PetscInt edgeL   = firstZEdge + (fy * numXVertices + fx) * numZEdges + fz;
969:           PetscInt edgeR   = firstZEdge + (fy * numXVertices + ((fx + 1) % numXVertices)) * numZEdges + fz;
970:           PetscInt edgeB   = firstXEdge + (fz * numYVertices + fy) * numXEdges + fx;
971:           PetscInt edgeT   = firstXEdge + (((fz + 1) % numZVertices) * numYVertices + fy) * numXEdges + fx;
972:           PetscInt ornt[4] = {0, 0, -1, -1};
973:           PetscInt cone[4];

975:           if (dim == 3) {
976:             /* markers */
977:             if (bdY != DM_BOUNDARY_PERIODIC) {
978:               if (fy == numYVertices - 1) {
979:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);
980:                 DMSetLabelValue(dm, "marker", face, markerBack);
981:               } else if (fy == 0) {
982:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);
983:                 DMSetLabelValue(dm, "marker", face, markerFront);
984:               }
985:             }
986:           }
987:           cone[0] = edgeB;
988:           cone[1] = edgeR;
989:           cone[2] = edgeT;
990:           cone[3] = edgeL;
991:           DMPlexSetCone(dm, face, cone);
992:           DMPlexSetConeOrientation(dm, face, ornt);
993:         }
994:       }
995:     }
996:     /* Build z faces */
997:     for (fy = 0; fy < numYEdges; ++fy) {
998:       for (fx = 0; fx < numXEdges; ++fx) {
999:         for (fz = 0; fz < numZVertices; fz++) {
1000:           PetscInt face    = firstZFace + (fy * numXEdges + fx) * numZVertices + fz;
1001:           PetscInt edgeL   = firstYEdge + (fz * numXVertices + fx) * numYEdges + fy;
1002:           PetscInt edgeR   = firstYEdge + (fz * numXVertices + ((fx + 1) % numXVertices)) * numYEdges + fy;
1003:           PetscInt edgeB   = firstXEdge + (fz * numYVertices + fy) * numXEdges + fx;
1004:           PetscInt edgeT   = firstXEdge + (fz * numYVertices + ((fy + 1) % numYVertices)) * numXEdges + fx;
1005:           PetscInt ornt[4] = {0, 0, -1, -1};
1006:           PetscInt cone[4];

1008:           if (dim == 2) {
1009:             if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges - 1) {
1010:               edgeR += numYEdges - 1 - 2 * fy;
1011:               ornt[1] = -1;
1012:             }
1013:             if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges - 1) {
1014:               edgeT += numXEdges - 1 - 2 * fx;
1015:               ornt[2] = 0;
1016:             }
1017:             if (bdX != DM_BOUNDARY_NONE && fx == numXEdges - 1 && cutLabel) DMLabelSetValue(cutLabel, face, 2);
1018:             if (bdY != DM_BOUNDARY_NONE && fy == numYEdges - 1 && cutLabel) DMLabelSetValue(cutLabel, face, 2);
1019:           } else {
1020:             /* markers */
1021:             if (bdZ != DM_BOUNDARY_PERIODIC) {
1022:               if (fz == numZVertices - 1) {
1023:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);
1024:                 DMSetLabelValue(dm, "marker", face, markerTop);
1025:               } else if (fz == 0) {
1026:                 DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);
1027:                 DMSetLabelValue(dm, "marker", face, markerBottom);
1028:               }
1029:             }
1030:           }
1031:           cone[0] = edgeB;
1032:           cone[1] = edgeR;
1033:           cone[2] = edgeT;
1034:           cone[3] = edgeL;
1035:           DMPlexSetCone(dm, face, cone);
1036:           DMPlexSetConeOrientation(dm, face, ornt);
1037:         }
1038:       }
1039:     }
1040:     /* Build Z edges*/
1041:     for (vy = 0; vy < numYVertices; vy++) {
1042:       for (vx = 0; vx < numXVertices; vx++) {
1043:         for (ez = 0; ez < numZEdges; ez++) {
1044:           const PetscInt edge    = firstZEdge + (vy * numXVertices + vx) * numZEdges + ez;
1045:           const PetscInt vertexB = firstVertex + (ez * numYVertices + vy) * numXVertices + vx;
1046:           const PetscInt vertexT = firstVertex + (((ez + 1) % numZVertices) * numYVertices + vy) * numXVertices + vx;
1047:           PetscInt       cone[2];

1049:           cone[0] = vertexB;
1050:           cone[1] = vertexT;
1051:           DMPlexSetCone(dm, edge, cone);
1052:           if (dim == 3) {
1053:             if (bdX != DM_BOUNDARY_PERIODIC) {
1054:               if (vx == numXVertices - 1) {
1055:                 DMSetLabelValue(dm, "marker", edge, markerRight);
1056:                 DMSetLabelValue(dm, "marker", cone[0], markerRight);
1057:                 if (ez == numZEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerRight);
1058:               } else if (vx == 0) {
1059:                 DMSetLabelValue(dm, "marker", edge, markerLeft);
1060:                 DMSetLabelValue(dm, "marker", cone[0], markerLeft);
1061:                 if (ez == numZEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerLeft);
1062:               }
1063:             }
1064:             if (bdY != DM_BOUNDARY_PERIODIC) {
1065:               if (vy == numYVertices - 1) {
1066:                 DMSetLabelValue(dm, "marker", edge, markerBack);
1067:                 DMSetLabelValue(dm, "marker", cone[0], markerBack);
1068:                 if (ez == numZEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerBack);
1069:               } else if (vy == 0) {
1070:                 DMSetLabelValue(dm, "marker", edge, markerFront);
1071:                 DMSetLabelValue(dm, "marker", cone[0], markerFront);
1072:                 if (ez == numZEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerFront);
1073:               }
1074:             }
1075:           }
1076:         }
1077:       }
1078:     }
1079:     /* Build Y edges*/
1080:     for (vz = 0; vz < numZVertices; vz++) {
1081:       for (vx = 0; vx < numXVertices; vx++) {
1082:         for (ey = 0; ey < numYEdges; ey++) {
1083:           const PetscInt nextv   = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges - 1) ? (numXVertices - vx - 1) : (vz * numYVertices + ((ey + 1) % numYVertices)) * numXVertices + vx;
1084:           const PetscInt edge    = firstYEdge + (vz * numXVertices + vx) * numYEdges + ey;
1085:           const PetscInt vertexF = firstVertex + (vz * numYVertices + ey) * numXVertices + vx;
1086:           const PetscInt vertexK = firstVertex + nextv;
1087:           PetscInt       cone[2];

1089:           cone[0] = vertexF;
1090:           cone[1] = vertexK;
1091:           DMPlexSetCone(dm, edge, cone);
1092:           if (dim == 2) {
1093:             if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) {
1094:               if (vx == numXVertices - 1) {
1095:                 DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);
1096:                 DMSetLabelValue(dm, "marker", edge, markerRight);
1097:                 DMSetLabelValue(dm, "marker", cone[0], markerRight);
1098:                 if (ey == numYEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerRight);
1099:               } else if (vx == 0) {
1100:                 DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);
1101:                 DMSetLabelValue(dm, "marker", edge, markerLeft);
1102:                 DMSetLabelValue(dm, "marker", cone[0], markerLeft);
1103:                 if (ey == numYEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerLeft);
1104:               }
1105:             } else {
1106:               if (vx == 0 && cutLabel) {
1107:                 DMLabelSetValue(cutLabel, edge, 1);
1108:                 DMLabelSetValue(cutLabel, cone[0], 1);
1109:                 if (ey == numYEdges - 1) DMLabelSetValue(cutLabel, cone[1], 1);
1110:               }
1111:             }
1112:           } else {
1113:             if (bdX != DM_BOUNDARY_PERIODIC) {
1114:               if (vx == numXVertices - 1) {
1115:                 DMSetLabelValue(dm, "marker", edge, markerRight);
1116:                 DMSetLabelValue(dm, "marker", cone[0], markerRight);
1117:                 if (ey == numYEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerRight);
1118:               } else if (vx == 0) {
1119:                 DMSetLabelValue(dm, "marker", edge, markerLeft);
1120:                 DMSetLabelValue(dm, "marker", cone[0], markerLeft);
1121:                 if (ey == numYEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerLeft);
1122:               }
1123:             }
1124:             if (bdZ != DM_BOUNDARY_PERIODIC) {
1125:               if (vz == numZVertices - 1) {
1126:                 DMSetLabelValue(dm, "marker", edge, markerTop);
1127:                 DMSetLabelValue(dm, "marker", cone[0], markerTop);
1128:                 if (ey == numYEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerTop);
1129:               } else if (vz == 0) {
1130:                 DMSetLabelValue(dm, "marker", edge, markerBottom);
1131:                 DMSetLabelValue(dm, "marker", cone[0], markerBottom);
1132:                 if (ey == numYEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerBottom);
1133:               }
1134:             }
1135:           }
1136:         }
1137:       }
1138:     }
1139:     /* Build X edges*/
1140:     for (vz = 0; vz < numZVertices; vz++) {
1141:       for (vy = 0; vy < numYVertices; vy++) {
1142:         for (ex = 0; ex < numXEdges; ex++) {
1143:           const PetscInt nextv   = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges - 1) ? (numYVertices - vy - 1) * numXVertices : (vz * numYVertices + vy) * numXVertices + (ex + 1) % numXVertices;
1144:           const PetscInt edge    = firstXEdge + (vz * numYVertices + vy) * numXEdges + ex;
1145:           const PetscInt vertexL = firstVertex + (vz * numYVertices + vy) * numXVertices + ex;
1146:           const PetscInt vertexR = firstVertex + nextv;
1147:           PetscInt       cone[2];

1149:           cone[0] = vertexL;
1150:           cone[1] = vertexR;
1151:           DMPlexSetCone(dm, edge, cone);
1152:           if (dim == 2) {
1153:             if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) {
1154:               if (vy == numYVertices - 1) {
1155:                 DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);
1156:                 DMSetLabelValue(dm, "marker", edge, markerTop);
1157:                 DMSetLabelValue(dm, "marker", cone[0], markerTop);
1158:                 if (ex == numXEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerTop);
1159:               } else if (vy == 0) {
1160:                 DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);
1161:                 DMSetLabelValue(dm, "marker", edge, markerBottom);
1162:                 DMSetLabelValue(dm, "marker", cone[0], markerBottom);
1163:                 if (ex == numXEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerBottom);
1164:               }
1165:             } else {
1166:               if (vy == 0 && cutLabel) {
1167:                 DMLabelSetValue(cutLabel, edge, 1);
1168:                 DMLabelSetValue(cutLabel, cone[0], 1);
1169:                 if (ex == numXEdges - 1) DMLabelSetValue(cutLabel, cone[1], 1);
1170:               }
1171:             }
1172:           } else {
1173:             if (bdY != DM_BOUNDARY_PERIODIC) {
1174:               if (vy == numYVertices - 1) {
1175:                 DMSetLabelValue(dm, "marker", edge, markerBack);
1176:                 DMSetLabelValue(dm, "marker", cone[0], markerBack);
1177:                 if (ex == numXEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerBack);
1178:               } else if (vy == 0) {
1179:                 DMSetLabelValue(dm, "marker", edge, markerFront);
1180:                 DMSetLabelValue(dm, "marker", cone[0], markerFront);
1181:                 if (ex == numXEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerFront);
1182:               }
1183:             }
1184:             if (bdZ != DM_BOUNDARY_PERIODIC) {
1185:               if (vz == numZVertices - 1) {
1186:                 DMSetLabelValue(dm, "marker", edge, markerTop);
1187:                 DMSetLabelValue(dm, "marker", cone[0], markerTop);
1188:                 if (ex == numXEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerTop);
1189:               } else if (vz == 0) {
1190:                 DMSetLabelValue(dm, "marker", edge, markerBottom);
1191:                 DMSetLabelValue(dm, "marker", cone[0], markerBottom);
1192:                 if (ex == numXEdges - 1) DMSetLabelValue(dm, "marker", cone[1], markerBottom);
1193:               }
1194:             }
1195:           }
1196:         }
1197:       }
1198:     }
1199:     DMPlexSymmetrize(dm);
1200:     DMPlexStratify(dm);
1201:     /* Build coordinates */
1202:     DMGetCoordinateSection(dm, &coordSection);
1203:     PetscSectionSetNumFields(coordSection, 1);
1204:     PetscSectionSetFieldComponents(coordSection, 0, dim);
1205:     PetscSectionSetChart(coordSection, firstVertex, firstVertex + numVertices);
1206:     for (v = firstVertex; v < firstVertex + numVertices; ++v) {
1207:       PetscSectionSetDof(coordSection, v, dim);
1208:       PetscSectionSetFieldDof(coordSection, v, 0, dim);
1209:     }
1210:     PetscSectionSetUp(coordSection);
1211:     PetscSectionGetStorageSize(coordSection, &coordSize);
1212:     VecCreate(PETSC_COMM_SELF, &coordinates);
1213:     PetscObjectSetName((PetscObject)coordinates, "coordinates");
1214:     VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1215:     VecSetBlockSize(coordinates, dim);
1216:     VecSetType(coordinates, VECSTANDARD);
1217:     VecGetArray(coordinates, &coords);
1218:     for (vz = 0; vz < numZVertices; ++vz) {
1219:       for (vy = 0; vy < numYVertices; ++vy) {
1220:         for (vx = 0; vx < numXVertices; ++vx) {
1221:           coords[((vz * numYVertices + vy) * numXVertices + vx) * dim + 0] = lower[0] + ((upper[0] - lower[0]) / numXEdges) * vx;
1222:           coords[((vz * numYVertices + vy) * numXVertices + vx) * dim + 1] = lower[1] + ((upper[1] - lower[1]) / numYEdges) * vy;
1223:           if (dim == 3) coords[((vz * numYVertices + vy) * numXVertices + vx) * dim + 2] = lower[2] + ((upper[2] - lower[2]) / numZEdges) * vz;
1224:         }
1225:       }
1226:     }
1227:     VecRestoreArray(coordinates, &coords);
1228:     DMSetCoordinatesLocal(dm, coordinates);
1229:     VecDestroy(&coordinates);
1230:   }
1231:   return 0;
1232: }

1234: static PetscErrorCode DMPlexCreateBoxMesh_Tensor_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[])
1235: {
1236:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
1237:   PetscInt       fac[3] = {0, 0, 0}, d;

1241:   DMSetDimension(dm, dim);
1242:   for (d = 0; d < dim; ++d) {
1243:     fac[d] = faces[d];
1244:     bdt[d] = periodicity[d];
1245:   }
1246:   DMPlexCreateCubeMesh_Internal(dm, lower, upper, fac, bdt[0], bdt[1], bdt[2]);
1247:   if (periodicity[0] == DM_BOUNDARY_PERIODIC || periodicity[0] == DM_BOUNDARY_TWIST || periodicity[1] == DM_BOUNDARY_PERIODIC || periodicity[1] == DM_BOUNDARY_TWIST || (dim > 2 && (periodicity[2] == DM_BOUNDARY_PERIODIC || periodicity[2] == DM_BOUNDARY_TWIST))) {
1248:     PetscReal L[3]       = {-1., -1., 0.};
1249:     PetscReal maxCell[3] = {-1., -1., 0.};

1251:     for (d = 0; d < dim; ++d) {
1252:       if (periodicity[d] != DM_BOUNDARY_NONE) {
1253:         L[d]       = upper[d] - lower[d];
1254:         maxCell[d] = 1.1 * (L[d] / PetscMax(1, faces[d]));
1255:       }
1256:     }
1257:     DMSetPeriodicity(dm, maxCell, lower, L);
1258:   }
1259:   DMPlexSetRefinementUniform(dm, PETSC_TRUE);
1260:   return 0;
1261: }

1263: static PetscErrorCode DMPlexCreateBoxMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate)
1264: {
1265:   if (dim == 1) DMPlexCreateLineMesh_Internal(dm, faces[0], lower[0], upper[0], periodicity[0]);
1266:   else if (simplex) DMPlexCreateBoxMesh_Simplex_Internal(dm, dim, faces, lower, upper, periodicity, interpolate);
1267:   else DMPlexCreateBoxMesh_Tensor_Internal(dm, dim, faces, lower, upper, periodicity);
1268:   if (!interpolate && dim > 1 && !simplex) {
1269:     DM udm;

1271:     DMPlexUninterpolate(dm, &udm);
1272:     DMPlexCopyCoordinates(dm, udm);
1273:     DMPlexReplace_Internal(dm, &udm);
1274:   }
1275:   return 0;
1276: }

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

1281:   Collective

1283:   Input Parameters:
1284: + comm        - The communicator for the DM object
1285: . dim         - The spatial dimension
1286: . simplex     - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells
1287: . faces       - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D
1288: . lower       - The lower left corner, or NULL for (0, 0, 0)
1289: . upper       - The upper right corner, or NULL for (1, 1, 1)
1290: . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE
1291: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

1293:   Output Parameter:
1294: . dm  - The DM object

1296:   Note: If you want to customize this mesh using options, you just need to
1297: $  DMCreate(comm, &dm);
1298: $  DMSetType(dm, DMPLEX);
1299: $  DMSetFromOptions(dm);
1300: and use the options on the DMSetFromOptions() page.

1302:   Here is the numbering returned for 2 faces in each direction for tensor cells:
1303: $ 10---17---11---18----12
1304: $  |         |         |
1305: $  |         |         |
1306: $ 20    2   22    3    24
1307: $  |         |         |
1308: $  |         |         |
1309: $  7---15----8---16----9
1310: $  |         |         |
1311: $  |         |         |
1312: $ 19    0   21    1   23
1313: $  |         |         |
1314: $  |         |         |
1315: $  4---13----5---14----6

1317: and for simplicial cells

1319: $ 14----8---15----9----16
1320: $  |\     5  |\      7 |
1321: $  | \       | \       |
1322: $ 13   2    14    3    15
1323: $  | 4   \   | 6   \   |
1324: $  |       \ |       \ |
1325: $ 11----6---12----7----13
1326: $  |\        |\        |
1327: $  | \    1  | \     3 |
1328: $ 10   0    11    1    12
1329: $  | 0   \   | 2   \   |
1330: $  |       \ |       \ |
1331: $  8----4----9----5----10

1333:   Level: beginner

1335: .seealso: `DMSetFromOptions()`, `DMPlexCreateFromFile()`, `DMPlexCreateHexCylinderMesh()`, `DMSetType()`, `DMCreate()`
1336: @*/
1337: 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)
1338: {
1339:   PetscInt       fac[3] = {1, 1, 1};
1340:   PetscReal      low[3] = {0, 0, 0};
1341:   PetscReal      upp[3] = {1, 1, 1};
1342:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};

1344:   DMCreate(comm, dm);
1345:   DMSetType(*dm, DMPLEX);
1346:   DMPlexCreateBoxMesh_Internal(*dm, dim, simplex, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt, interpolate);
1347:   if (periodicity) DMLocalizeCoordinates(*dm);
1348:   return 0;
1349: }

1351: static PetscErrorCode DMPlexCreateWedgeBoxMesh_Internal(DM dm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[])
1352: {
1353:   DM       bdm, vol;
1354:   PetscInt i;

1357:   DMCreate(PetscObjectComm((PetscObject)dm), &bdm);
1358:   DMSetType(bdm, DMPLEX);
1359:   DMSetDimension(bdm, 2);
1360:   DMPlexCreateBoxMesh_Simplex_Internal(bdm, 2, faces, lower, upper, periodicity, PETSC_TRUE);
1361:   DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], PETSC_TRUE, PETSC_FALSE, NULL, NULL, &vol);
1362:   DMDestroy(&bdm);
1363:   DMPlexReplace_Internal(dm, &vol);
1364:   if (lower[2] != 0.0) {
1365:     Vec          v;
1366:     PetscScalar *x;
1367:     PetscInt     cDim, n;

1369:     DMGetCoordinatesLocal(dm, &v);
1370:     VecGetBlockSize(v, &cDim);
1371:     VecGetLocalSize(v, &n);
1372:     VecGetArray(v, &x);
1373:     x += cDim;
1374:     for (i = 0; i < n; i += cDim) x[i] += lower[2];
1375:     VecRestoreArray(v, &x);
1376:     DMSetCoordinatesLocal(dm, v);
1377:   }
1378:   return 0;
1379: }

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

1384:   Collective

1386:   Input Parameters:
1387: + comm        - The communicator for the DM object
1388: . faces       - Number of faces per dimension, or NULL for (1, 1, 1)
1389: . lower       - The lower left corner, or NULL for (0, 0, 0)
1390: . upper       - The upper right corner, or NULL for (1, 1, 1)
1391: . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE
1392: . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first
1393: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

1395:   Output Parameter:
1396: . dm  - The DM object

1398:   Level: beginner

1400: .seealso: `DMPlexCreateHexCylinderMesh()`, `DMPlexCreateWedgeCylinderMesh()`, `DMExtrude()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()`
1401: @*/
1402: PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm)
1403: {
1404:   PetscInt       fac[3] = {1, 1, 1};
1405:   PetscReal      low[3] = {0, 0, 0};
1406:   PetscReal      upp[3] = {1, 1, 1};
1407:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};

1409:   DMCreate(comm, dm);
1410:   DMSetType(*dm, DMPLEX);
1411:   DMPlexCreateWedgeBoxMesh_Internal(*dm, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt);
1412:   if (!interpolate) {
1413:     DM udm;

1415:     DMPlexUninterpolate(*dm, &udm);
1416:     DMPlexReplace_Internal(*dm, &udm);
1417:   }
1418:   if (periodicity) DMLocalizeCoordinates(*dm);
1419:   return 0;
1420: }

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

1425:   Logically Collective on dm

1427:   Input Parameters:
1428: + dm - the DM context
1429: - prefix - the prefix to prepend to all option names

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

1435:   Level: advanced

1437: .seealso: `SNESSetFromOptions()`
1438: @*/
1439: PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[])
1440: {
1441:   DM_Plex *mesh = (DM_Plex *)dm->data;

1444:   PetscObjectSetOptionsPrefix((PetscObject)dm, prefix);
1445:   PetscObjectSetOptionsPrefix((PetscObject)mesh->partitioner, prefix);
1446:   return 0;
1447: }

1449: /* Remap geometry to cylinder
1450:    TODO: This only works for a single refinement, then it is broken

1452:      Interior square: Linear interpolation is correct
1453:      The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing
1454:      such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance

1456:        phi     = arctan(y/x)
1457:        d_close = sqrt(1/8 + 1/4 sin^2(phi))
1458:        d_far   = sqrt(1/2 + sin^2(phi))

1460:      so we remap them using

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

1465:      If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y.
1466: */
1467: static void snapToCylinder(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1468: {
1469:   const PetscReal dis = 1.0 / PetscSqrtReal(2.0);
1470:   const PetscReal ds2 = 0.5 * dis;

1472:   if ((PetscAbsScalar(u[0]) <= ds2) && (PetscAbsScalar(u[1]) <= ds2)) {
1473:     f0[0] = u[0];
1474:     f0[1] = u[1];
1475:   } else {
1476:     PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc;

1478:     x    = PetscRealPart(u[0]);
1479:     y    = PetscRealPart(u[1]);
1480:     phi  = PetscAtan2Real(y, x);
1481:     sinp = PetscSinReal(phi);
1482:     cosp = PetscCosReal(phi);
1483:     if ((PetscAbsReal(phi) > PETSC_PI / 4.0) && (PetscAbsReal(phi) < 3.0 * PETSC_PI / 4.0)) {
1484:       dc = PetscAbsReal(ds2 / sinp);
1485:       df = PetscAbsReal(dis / sinp);
1486:       xc = ds2 * x / PetscAbsReal(y);
1487:       yc = ds2 * PetscSignReal(y);
1488:     } else {
1489:       dc = PetscAbsReal(ds2 / cosp);
1490:       df = PetscAbsReal(dis / cosp);
1491:       xc = ds2 * PetscSignReal(x);
1492:       yc = ds2 * y / PetscAbsReal(x);
1493:     }
1494:     f0[0] = xc + (u[0] - xc) * (1.0 - dc) / (df - dc);
1495:     f0[1] = yc + (u[1] - yc) * (1.0 - dc) / (df - dc);
1496:   }
1497:   f0[2] = u[2];
1498: }

1500: static PetscErrorCode DMPlexCreateHexCylinderMesh_Internal(DM dm, DMBoundaryType periodicZ)
1501: {
1502:   const PetscInt dim = 3;
1503:   PetscInt       numCells, numVertices;
1504:   PetscMPIInt    rank;

1506:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1507:   DMSetDimension(dm, dim);
1508:   /* Create topology */
1509:   {
1510:     PetscInt cone[8], c;

1512:     numCells    = rank == 0 ? 5 : 0;
1513:     numVertices = rank == 0 ? 16 : 0;
1514:     if (periodicZ == DM_BOUNDARY_PERIODIC) {
1515:       numCells *= 3;
1516:       numVertices = rank == 0 ? 24 : 0;
1517:     }
1518:     DMPlexSetChart(dm, 0, numCells + numVertices);
1519:     for (c = 0; c < numCells; c++) DMPlexSetConeSize(dm, c, 8);
1520:     DMSetUp(dm);
1521:     if (rank == 0) {
1522:       if (periodicZ == DM_BOUNDARY_PERIODIC) {
1523:         cone[0] = 15;
1524:         cone[1] = 18;
1525:         cone[2] = 17;
1526:         cone[3] = 16;
1527:         cone[4] = 31;
1528:         cone[5] = 32;
1529:         cone[6] = 33;
1530:         cone[7] = 34;
1531:         DMPlexSetCone(dm, 0, cone);
1532:         cone[0] = 16;
1533:         cone[1] = 17;
1534:         cone[2] = 24;
1535:         cone[3] = 23;
1536:         cone[4] = 32;
1537:         cone[5] = 36;
1538:         cone[6] = 37;
1539:         cone[7] = 33; /* 22 25 26 21 */
1540:         DMPlexSetCone(dm, 1, cone);
1541:         cone[0] = 18;
1542:         cone[1] = 27;
1543:         cone[2] = 24;
1544:         cone[3] = 17;
1545:         cone[4] = 34;
1546:         cone[5] = 33;
1547:         cone[6] = 37;
1548:         cone[7] = 38;
1549:         DMPlexSetCone(dm, 2, cone);
1550:         cone[0] = 29;
1551:         cone[1] = 27;
1552:         cone[2] = 18;
1553:         cone[3] = 15;
1554:         cone[4] = 35;
1555:         cone[5] = 31;
1556:         cone[6] = 34;
1557:         cone[7] = 38;
1558:         DMPlexSetCone(dm, 3, cone);
1559:         cone[0] = 29;
1560:         cone[1] = 15;
1561:         cone[2] = 16;
1562:         cone[3] = 23;
1563:         cone[4] = 35;
1564:         cone[5] = 36;
1565:         cone[6] = 32;
1566:         cone[7] = 31;
1567:         DMPlexSetCone(dm, 4, cone);

1569:         cone[0] = 31;
1570:         cone[1] = 34;
1571:         cone[2] = 33;
1572:         cone[3] = 32;
1573:         cone[4] = 19;
1574:         cone[5] = 22;
1575:         cone[6] = 21;
1576:         cone[7] = 20;
1577:         DMPlexSetCone(dm, 5, cone);
1578:         cone[0] = 32;
1579:         cone[1] = 33;
1580:         cone[2] = 37;
1581:         cone[3] = 36;
1582:         cone[4] = 22;
1583:         cone[5] = 25;
1584:         cone[6] = 26;
1585:         cone[7] = 21;
1586:         DMPlexSetCone(dm, 6, cone);
1587:         cone[0] = 34;
1588:         cone[1] = 38;
1589:         cone[2] = 37;
1590:         cone[3] = 33;
1591:         cone[4] = 20;
1592:         cone[5] = 21;
1593:         cone[6] = 26;
1594:         cone[7] = 28;
1595:         DMPlexSetCone(dm, 7, cone);
1596:         cone[0] = 35;
1597:         cone[1] = 38;
1598:         cone[2] = 34;
1599:         cone[3] = 31;
1600:         cone[4] = 30;
1601:         cone[5] = 19;
1602:         cone[6] = 20;
1603:         cone[7] = 28;
1604:         DMPlexSetCone(dm, 8, cone);
1605:         cone[0] = 35;
1606:         cone[1] = 31;
1607:         cone[2] = 32;
1608:         cone[3] = 36;
1609:         cone[4] = 30;
1610:         cone[5] = 25;
1611:         cone[6] = 22;
1612:         cone[7] = 19;
1613:         DMPlexSetCone(dm, 9, cone);

1615:         cone[0] = 19;
1616:         cone[1] = 20;
1617:         cone[2] = 21;
1618:         cone[3] = 22;
1619:         cone[4] = 15;
1620:         cone[5] = 16;
1621:         cone[6] = 17;
1622:         cone[7] = 18;
1623:         DMPlexSetCone(dm, 10, cone);
1624:         cone[0] = 22;
1625:         cone[1] = 21;
1626:         cone[2] = 26;
1627:         cone[3] = 25;
1628:         cone[4] = 16;
1629:         cone[5] = 23;
1630:         cone[6] = 24;
1631:         cone[7] = 17;
1632:         DMPlexSetCone(dm, 11, cone);
1633:         cone[0] = 20;
1634:         cone[1] = 28;
1635:         cone[2] = 26;
1636:         cone[3] = 21;
1637:         cone[4] = 18;
1638:         cone[5] = 17;
1639:         cone[6] = 24;
1640:         cone[7] = 27;
1641:         DMPlexSetCone(dm, 12, cone);
1642:         cone[0] = 30;
1643:         cone[1] = 28;
1644:         cone[2] = 20;
1645:         cone[3] = 19;
1646:         cone[4] = 29;
1647:         cone[5] = 15;
1648:         cone[6] = 18;
1649:         cone[7] = 27;
1650:         DMPlexSetCone(dm, 13, cone);
1651:         cone[0] = 30;
1652:         cone[1] = 19;
1653:         cone[2] = 22;
1654:         cone[3] = 25;
1655:         cone[4] = 29;
1656:         cone[5] = 23;
1657:         cone[6] = 16;
1658:         cone[7] = 15;
1659:         DMPlexSetCone(dm, 14, cone);
1660:       } else {
1661:         cone[0] = 5;
1662:         cone[1] = 8;
1663:         cone[2] = 7;
1664:         cone[3] = 6;
1665:         cone[4] = 9;
1666:         cone[5] = 12;
1667:         cone[6] = 11;
1668:         cone[7] = 10;
1669:         DMPlexSetCone(dm, 0, cone);
1670:         cone[0] = 6;
1671:         cone[1] = 7;
1672:         cone[2] = 14;
1673:         cone[3] = 13;
1674:         cone[4] = 12;
1675:         cone[5] = 15;
1676:         cone[6] = 16;
1677:         cone[7] = 11;
1678:         DMPlexSetCone(dm, 1, cone);
1679:         cone[0] = 8;
1680:         cone[1] = 17;
1681:         cone[2] = 14;
1682:         cone[3] = 7;
1683:         cone[4] = 10;
1684:         cone[5] = 11;
1685:         cone[6] = 16;
1686:         cone[7] = 18;
1687:         DMPlexSetCone(dm, 2, cone);
1688:         cone[0] = 19;
1689:         cone[1] = 17;
1690:         cone[2] = 8;
1691:         cone[3] = 5;
1692:         cone[4] = 20;
1693:         cone[5] = 9;
1694:         cone[6] = 10;
1695:         cone[7] = 18;
1696:         DMPlexSetCone(dm, 3, cone);
1697:         cone[0] = 19;
1698:         cone[1] = 5;
1699:         cone[2] = 6;
1700:         cone[3] = 13;
1701:         cone[4] = 20;
1702:         cone[5] = 15;
1703:         cone[6] = 12;
1704:         cone[7] = 9;
1705:         DMPlexSetCone(dm, 4, cone);
1706:       }
1707:     }
1708:     DMPlexSymmetrize(dm);
1709:     DMPlexStratify(dm);
1710:   }
1711:   /* Create cube geometry */
1712:   {
1713:     Vec             coordinates;
1714:     PetscSection    coordSection;
1715:     PetscScalar    *coords;
1716:     PetscInt        coordSize, v;
1717:     const PetscReal dis = 1.0 / PetscSqrtReal(2.0);
1718:     const PetscReal ds2 = dis / 2.0;

1720:     /* Build coordinates */
1721:     DMGetCoordinateSection(dm, &coordSection);
1722:     PetscSectionSetNumFields(coordSection, 1);
1723:     PetscSectionSetFieldComponents(coordSection, 0, dim);
1724:     PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1725:     for (v = numCells; v < numCells + numVertices; ++v) {
1726:       PetscSectionSetDof(coordSection, v, dim);
1727:       PetscSectionSetFieldDof(coordSection, v, 0, dim);
1728:     }
1729:     PetscSectionSetUp(coordSection);
1730:     PetscSectionGetStorageSize(coordSection, &coordSize);
1731:     VecCreate(PETSC_COMM_SELF, &coordinates);
1732:     PetscObjectSetName((PetscObject)coordinates, "coordinates");
1733:     VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1734:     VecSetBlockSize(coordinates, dim);
1735:     VecSetType(coordinates, VECSTANDARD);
1736:     VecGetArray(coordinates, &coords);
1737:     if (rank == 0) {
1738:       coords[0 * dim + 0]  = -ds2;
1739:       coords[0 * dim + 1]  = -ds2;
1740:       coords[0 * dim + 2]  = 0.0;
1741:       coords[1 * dim + 0]  = ds2;
1742:       coords[1 * dim + 1]  = -ds2;
1743:       coords[1 * dim + 2]  = 0.0;
1744:       coords[2 * dim + 0]  = ds2;
1745:       coords[2 * dim + 1]  = ds2;
1746:       coords[2 * dim + 2]  = 0.0;
1747:       coords[3 * dim + 0]  = -ds2;
1748:       coords[3 * dim + 1]  = ds2;
1749:       coords[3 * dim + 2]  = 0.0;
1750:       coords[4 * dim + 0]  = -ds2;
1751:       coords[4 * dim + 1]  = -ds2;
1752:       coords[4 * dim + 2]  = 1.0;
1753:       coords[5 * dim + 0]  = -ds2;
1754:       coords[5 * dim + 1]  = ds2;
1755:       coords[5 * dim + 2]  = 1.0;
1756:       coords[6 * dim + 0]  = ds2;
1757:       coords[6 * dim + 1]  = ds2;
1758:       coords[6 * dim + 2]  = 1.0;
1759:       coords[7 * dim + 0]  = ds2;
1760:       coords[7 * dim + 1]  = -ds2;
1761:       coords[7 * dim + 2]  = 1.0;
1762:       coords[8 * dim + 0]  = dis;
1763:       coords[8 * dim + 1]  = -dis;
1764:       coords[8 * dim + 2]  = 0.0;
1765:       coords[9 * dim + 0]  = dis;
1766:       coords[9 * dim + 1]  = dis;
1767:       coords[9 * dim + 2]  = 0.0;
1768:       coords[10 * dim + 0] = dis;
1769:       coords[10 * dim + 1] = -dis;
1770:       coords[10 * dim + 2] = 1.0;
1771:       coords[11 * dim + 0] = dis;
1772:       coords[11 * dim + 1] = dis;
1773:       coords[11 * dim + 2] = 1.0;
1774:       coords[12 * dim + 0] = -dis;
1775:       coords[12 * dim + 1] = dis;
1776:       coords[12 * dim + 2] = 0.0;
1777:       coords[13 * dim + 0] = -dis;
1778:       coords[13 * dim + 1] = dis;
1779:       coords[13 * dim + 2] = 1.0;
1780:       coords[14 * dim + 0] = -dis;
1781:       coords[14 * dim + 1] = -dis;
1782:       coords[14 * dim + 2] = 0.0;
1783:       coords[15 * dim + 0] = -dis;
1784:       coords[15 * dim + 1] = -dis;
1785:       coords[15 * dim + 2] = 1.0;
1786:       if (periodicZ == DM_BOUNDARY_PERIODIC) {
1787:         /* 15 31 19 */ coords[16 * dim + 0] = -ds2;
1788:         coords[16 * dim + 1]                = -ds2;
1789:         coords[16 * dim + 2]                = 0.5;
1790:         /* 16 32 22 */ coords[17 * dim + 0] = ds2;
1791:         coords[17 * dim + 1]                = -ds2;
1792:         coords[17 * dim + 2]                = 0.5;
1793:         /* 17 33 21 */ coords[18 * dim + 0] = ds2;
1794:         coords[18 * dim + 1]                = ds2;
1795:         coords[18 * dim + 2]                = 0.5;
1796:         /* 18 34 20 */ coords[19 * dim + 0] = -ds2;
1797:         coords[19 * dim + 1]                = ds2;
1798:         coords[19 * dim + 2]                = 0.5;
1799:         /* 29 35 30 */ coords[20 * dim + 0] = -dis;
1800:         coords[20 * dim + 1]                = -dis;
1801:         coords[20 * dim + 2]                = 0.5;
1802:         /* 23 36 25 */ coords[21 * dim + 0] = dis;
1803:         coords[21 * dim + 1]                = -dis;
1804:         coords[21 * dim + 2]                = 0.5;
1805:         /* 24 37 26 */ coords[22 * dim + 0] = dis;
1806:         coords[22 * dim + 1]                = dis;
1807:         coords[22 * dim + 2]                = 0.5;
1808:         /* 27 38 28 */ coords[23 * dim + 0] = -dis;
1809:         coords[23 * dim + 1]                = dis;
1810:         coords[23 * dim + 2]                = 0.5;
1811:       }
1812:     }
1813:     VecRestoreArray(coordinates, &coords);
1814:     DMSetCoordinatesLocal(dm, coordinates);
1815:     VecDestroy(&coordinates);
1816:   }
1817:   /* Create periodicity */
1818:   if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) {
1819:     PetscReal L[3]       = {-1., -1., 0.};
1820:     PetscReal maxCell[3] = {-1., -1., 0.};
1821:     PetscReal lower[3]   = {0.0, 0.0, 0.0};
1822:     PetscReal upper[3]   = {1.0, 1.0, 1.5};
1823:     PetscInt  numZCells  = 3;

1825:     L[2]       = upper[2] - lower[2];
1826:     maxCell[2] = 1.1 * (L[2] / numZCells);
1827:     DMSetPeriodicity(dm, maxCell, lower, L);
1828:   }
1829:   {
1830:     DM          cdm;
1831:     PetscDS     cds;
1832:     PetscScalar c[2] = {1.0, 1.0};

1834:     DMPlexCreateCoordinateSpace(dm, 1, snapToCylinder);
1835:     DMGetCoordinateDM(dm, &cdm);
1836:     DMGetDS(cdm, &cds);
1837:     PetscDSSetConstants(cds, 2, c);
1838:   }
1839:   /* Wait for coordinate creation before doing in-place modification */
1840:   DMPlexInterpolateInPlace_Internal(dm);
1841:   return 0;
1842: }

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

1847:   Collective

1849:   Input Parameters:
1850: + comm      - The communicator for the DM object
1851: - periodicZ - The boundary type for the Z direction

1853:   Output Parameter:
1854: . dm  - The DM object

1856:   Note:
1857:   Here is the output numbering looking from the bottom of the cylinder:
1858: $       17-----14
1859: $        |     |
1860: $        |  2  |
1861: $        |     |
1862: $ 17-----8-----7-----14
1863: $  |     |     |     |
1864: $  |  3  |  0  |  1  |
1865: $  |     |     |     |
1866: $ 19-----5-----6-----13
1867: $        |     |
1868: $        |  4  |
1869: $        |     |
1870: $       19-----13
1871: $
1872: $ and up through the top
1873: $
1874: $       18-----16
1875: $        |     |
1876: $        |  2  |
1877: $        |     |
1878: $ 18----10----11-----16
1879: $  |     |     |     |
1880: $  |  3  |  0  |  1  |
1881: $  |     |     |     |
1882: $ 20-----9----12-----15
1883: $        |     |
1884: $        |  4  |
1885: $        |     |
1886: $       20-----15

1888:   Level: beginner

1890: .seealso: `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()`
1891: @*/
1892: PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, DMBoundaryType periodicZ, DM *dm)
1893: {
1895:   DMCreate(comm, dm);
1896:   DMSetType(*dm, DMPLEX);
1897:   DMPlexCreateHexCylinderMesh_Internal(*dm, periodicZ);
1898:   return 0;
1899: }

1901: static PetscErrorCode DMPlexCreateWedgeCylinderMesh_Internal(DM dm, PetscInt n, PetscBool interpolate)
1902: {
1903:   const PetscInt dim = 3;
1904:   PetscInt       numCells, numVertices, v;
1905:   PetscMPIInt    rank;

1908:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1909:   DMSetDimension(dm, dim);
1910:   /* Must create the celltype label here so that we do not automatically try to compute the types */
1911:   DMCreateLabel(dm, "celltype");
1912:   /* Create topology */
1913:   {
1914:     PetscInt cone[6], c;

1916:     numCells    = rank == 0 ? n : 0;
1917:     numVertices = rank == 0 ? 2 * (n + 1) : 0;
1918:     DMPlexSetChart(dm, 0, numCells + numVertices);
1919:     for (c = 0; c < numCells; c++) DMPlexSetConeSize(dm, c, 6);
1920:     DMSetUp(dm);
1921:     for (c = 0; c < numCells; c++) {
1922:       cone[0] = c + n * 1;
1923:       cone[1] = (c + 1) % n + n * 1;
1924:       cone[2] = 0 + 3 * n;
1925:       cone[3] = c + n * 2;
1926:       cone[4] = (c + 1) % n + n * 2;
1927:       cone[5] = 1 + 3 * n;
1928:       DMPlexSetCone(dm, c, cone);
1929:       DMPlexSetCellType(dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR);
1930:     }
1931:     DMPlexSymmetrize(dm);
1932:     DMPlexStratify(dm);
1933:   }
1934:   for (v = numCells; v < numCells + numVertices; ++v) DMPlexSetCellType(dm, v, DM_POLYTOPE_POINT);
1935:   /* Create cylinder geometry */
1936:   {
1937:     Vec          coordinates;
1938:     PetscSection coordSection;
1939:     PetscScalar *coords;
1940:     PetscInt     coordSize, c;

1942:     /* Build coordinates */
1943:     DMGetCoordinateSection(dm, &coordSection);
1944:     PetscSectionSetNumFields(coordSection, 1);
1945:     PetscSectionSetFieldComponents(coordSection, 0, dim);
1946:     PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1947:     for (v = numCells; v < numCells + numVertices; ++v) {
1948:       PetscSectionSetDof(coordSection, v, dim);
1949:       PetscSectionSetFieldDof(coordSection, v, 0, dim);
1950:     }
1951:     PetscSectionSetUp(coordSection);
1952:     PetscSectionGetStorageSize(coordSection, &coordSize);
1953:     VecCreate(PETSC_COMM_SELF, &coordinates);
1954:     PetscObjectSetName((PetscObject)coordinates, "coordinates");
1955:     VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1956:     VecSetBlockSize(coordinates, dim);
1957:     VecSetType(coordinates, VECSTANDARD);
1958:     VecGetArray(coordinates, &coords);
1959:     for (c = 0; c < numCells; c++) {
1960:       coords[(c + 0 * n) * dim + 0] = PetscCosReal(2.0 * c * PETSC_PI / n);
1961:       coords[(c + 0 * n) * dim + 1] = PetscSinReal(2.0 * c * PETSC_PI / n);
1962:       coords[(c + 0 * n) * dim + 2] = 1.0;
1963:       coords[(c + 1 * n) * dim + 0] = PetscCosReal(2.0 * c * PETSC_PI / n);
1964:       coords[(c + 1 * n) * dim + 1] = PetscSinReal(2.0 * c * PETSC_PI / n);
1965:       coords[(c + 1 * n) * dim + 2] = 0.0;
1966:     }
1967:     if (rank == 0) {
1968:       coords[(2 * n + 0) * dim + 0] = 0.0;
1969:       coords[(2 * n + 0) * dim + 1] = 0.0;
1970:       coords[(2 * n + 0) * dim + 2] = 1.0;
1971:       coords[(2 * n + 1) * dim + 0] = 0.0;
1972:       coords[(2 * n + 1) * dim + 1] = 0.0;
1973:       coords[(2 * n + 1) * dim + 2] = 0.0;
1974:     }
1975:     VecRestoreArray(coordinates, &coords);
1976:     DMSetCoordinatesLocal(dm, coordinates);
1977:     VecDestroy(&coordinates);
1978:   }
1979:   /* Interpolate */
1980:   if (interpolate) DMPlexInterpolateInPlace_Internal(dm);
1981:   return 0;
1982: }

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

1987:   Collective

1989:   Input Parameters:
1990: + comm - The communicator for the DM object
1991: . n    - The number of wedges around the origin
1992: - interpolate - Create edges and faces

1994:   Output Parameter:
1995: . dm  - The DM object

1997:   Level: beginner

1999: .seealso: `DMPlexCreateHexCylinderMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()`
2000: @*/
2001: PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm)
2002: {
2004:   DMCreate(comm, dm);
2005:   DMSetType(*dm, DMPLEX);
2006:   DMPlexCreateWedgeCylinderMesh_Internal(*dm, n, interpolate);
2007:   return 0;
2008: }

2010: static inline PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[])
2011: {
2012:   PetscReal prod = 0.0;
2013:   PetscInt  i;
2014:   for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]);
2015:   return PetscSqrtReal(prod);
2016: }
2017: static inline PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[])
2018: {
2019:   PetscReal prod = 0.0;
2020:   PetscInt  i;
2021:   for (i = 0; i < dim; ++i) prod += x[i] * y[i];
2022:   return prod;
2023: }

2025: /* The first constant is the sphere radius */
2026: static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
2027: {
2028:   PetscReal r     = PetscRealPart(constants[0]);
2029:   PetscReal norm2 = 0.0, fac;
2030:   PetscInt  n     = uOff[1] - uOff[0], d;

2032:   for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d]));
2033:   fac = r / PetscSqrtReal(norm2);
2034:   for (d = 0; d < n; ++d) f0[d] = u[d] * fac;
2035: }

2037: static PetscErrorCode DMPlexCreateSphereMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, PetscReal R)
2038: {
2039:   const PetscInt embedDim = dim + 1;
2040:   PetscSection   coordSection;
2041:   Vec            coordinates;
2042:   PetscScalar   *coords;
2043:   PetscReal     *coordsIn;
2044:   PetscInt       numCells, numEdges, numVerts = 0, firstVertex = 0, v, firstEdge, coordSize, d, c, e;
2045:   PetscMPIInt    rank;

2048:   DMSetDimension(dm, dim);
2049:   DMSetCoordinateDim(dm, dim + 1);
2050:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
2051:   switch (dim) {
2052:   case 2:
2053:     if (simplex) {
2054:       const PetscReal radius    = PetscSqrtReal(1 + PETSC_PHI * PETSC_PHI) / (1.0 + PETSC_PHI);
2055:       const PetscReal edgeLen   = 2.0 / (1.0 + PETSC_PHI) * (R / radius);
2056:       const PetscInt  degree    = 5;
2057:       PetscReal       vertex[3] = {0.0, 1.0 / (1.0 + PETSC_PHI), PETSC_PHI / (1.0 + PETSC_PHI)};
2058:       PetscInt        s[3]      = {1, 1, 1};
2059:       PetscInt        cone[3];
2060:       PetscInt       *graph, p, i, j, k;

2062:       vertex[0] *= R / radius;
2063:       vertex[1] *= R / radius;
2064:       vertex[2] *= R / radius;
2065:       numCells    = rank == 0 ? 20 : 0;
2066:       numVerts    = rank == 0 ? 12 : 0;
2067:       firstVertex = numCells;
2068:       /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of

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

2072:          where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge
2073:          length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393.
2074:       */
2075:       /* Construct vertices */
2076:       PetscCalloc1(numVerts * embedDim, &coordsIn);
2077:       if (rank == 0) {
2078:         for (p = 0, i = 0; p < embedDim; ++p) {
2079:           for (s[1] = -1; s[1] < 2; s[1] += 2) {
2080:             for (s[2] = -1; s[2] < 2; s[2] += 2) {
2081:               for (d = 0; d < embedDim; ++d) coordsIn[i * embedDim + d] = s[(d + p) % embedDim] * vertex[(d + p) % embedDim];
2082:               ++i;
2083:             }
2084:           }
2085:         }
2086:       }
2087:       /* Construct graph */
2088:       PetscCalloc1(numVerts * numVerts, &graph);
2089:       for (i = 0; i < numVerts; ++i) {
2090:         for (j = 0, k = 0; j < numVerts; ++j) {
2091:           if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i * embedDim], &coordsIn[j * embedDim]) - edgeLen) < PETSC_SMALL) {
2092:             graph[i * numVerts + j] = 1;
2093:             ++k;
2094:           }
2095:         }
2097:       }
2098:       /* Build Topology */
2099:       DMPlexSetChart(dm, 0, numCells + numVerts);
2100:       for (c = 0; c < numCells; c++) DMPlexSetConeSize(dm, c, embedDim);
2101:       DMSetUp(dm); /* Allocate space for cones */
2102:       /* Cells */
2103:       for (i = 0, c = 0; i < numVerts; ++i) {
2104:         for (j = 0; j < i; ++j) {
2105:           for (k = 0; k < j; ++k) {
2106:             if (graph[i * numVerts + j] && graph[j * numVerts + k] && graph[k * numVerts + i]) {
2107:               cone[0] = firstVertex + i;
2108:               cone[1] = firstVertex + j;
2109:               cone[2] = firstVertex + k;
2110:               /* Check orientation */
2111:               {
2112:                 const PetscInt epsilon[3][3][3] = {
2113:                   {{0, 0, 0},  {0, 0, 1},  {0, -1, 0}},
2114:                   {{0, 0, -1}, {0, 0, 0},  {1, 0, 0} },
2115:                   {{0, 1, 0},  {-1, 0, 0}, {0, 0, 0} }
2116:                 };
2117:                 PetscReal normal[3];
2118:                 PetscInt  e, f;

2120:                 for (d = 0; d < embedDim; ++d) {
2121:                   normal[d] = 0.0;
2122:                   for (e = 0; e < embedDim; ++e) {
2123:                     for (f = 0; f < embedDim; ++f) normal[d] += epsilon[d][e][f] * (coordsIn[j * embedDim + e] - coordsIn[i * embedDim + e]) * (coordsIn[k * embedDim + f] - coordsIn[i * embedDim + f]);
2124:                   }
2125:                 }
2126:                 if (DotReal(embedDim, normal, &coordsIn[i * embedDim]) < 0) {
2127:                   PetscInt tmp = cone[1];
2128:                   cone[1]      = cone[2];
2129:                   cone[2]      = tmp;
2130:                 }
2131:               }
2132:               DMPlexSetCone(dm, c++, cone);
2133:             }
2134:           }
2135:         }
2136:       }
2137:       DMPlexSymmetrize(dm);
2138:       DMPlexStratify(dm);
2139:       PetscFree(graph);
2140:     } else {
2141:       /*
2142:         12-21--13
2143:          |     |
2144:         25  4  24
2145:          |     |
2146:   12-25--9-16--8-24--13
2147:    |     |     |     |
2148:   23  5 17  0 15  3  22
2149:    |     |     |     |
2150:   10-20--6-14--7-19--11
2151:          |     |
2152:         20  1  19
2153:          |     |
2154:         10-18--11
2155:          |     |
2156:         23  2  22
2157:          |     |
2158:         12-21--13
2159:        */
2160:       PetscInt cone[4], ornt[4];

2162:       numCells    = rank == 0 ? 6 : 0;
2163:       numEdges    = rank == 0 ? 12 : 0;
2164:       numVerts    = rank == 0 ? 8 : 0;
2165:       firstVertex = numCells;
2166:       firstEdge   = numCells + numVerts;
2167:       /* Build Topology */
2168:       DMPlexSetChart(dm, 0, numCells + numEdges + numVerts);
2169:       for (c = 0; c < numCells; c++) DMPlexSetConeSize(dm, c, 4);
2170:       for (e = firstEdge; e < firstEdge + numEdges; ++e) DMPlexSetConeSize(dm, e, 2);
2171:       DMSetUp(dm); /* Allocate space for cones */
2172:       if (rank == 0) {
2173:         /* Cell 0 */
2174:         cone[0] = 14;
2175:         cone[1] = 15;
2176:         cone[2] = 16;
2177:         cone[3] = 17;
2178:         DMPlexSetCone(dm, 0, cone);
2179:         ornt[0] = 0;
2180:         ornt[1] = 0;
2181:         ornt[2] = 0;
2182:         ornt[3] = 0;
2183:         DMPlexSetConeOrientation(dm, 0, ornt);
2184:         /* Cell 1 */
2185:         cone[0] = 18;
2186:         cone[1] = 19;
2187:         cone[2] = 14;
2188:         cone[3] = 20;
2189:         DMPlexSetCone(dm, 1, cone);
2190:         ornt[0] = 0;
2191:         ornt[1] = 0;
2192:         ornt[2] = -1;
2193:         ornt[3] = 0;
2194:         DMPlexSetConeOrientation(dm, 1, ornt);
2195:         /* Cell 2 */
2196:         cone[0] = 21;
2197:         cone[1] = 22;
2198:         cone[2] = 18;
2199:         cone[3] = 23;
2200:         DMPlexSetCone(dm, 2, cone);
2201:         ornt[0] = 0;
2202:         ornt[1] = 0;
2203:         ornt[2] = -1;
2204:         ornt[3] = 0;
2205:         DMPlexSetConeOrientation(dm, 2, ornt);
2206:         /* Cell 3 */
2207:         cone[0] = 19;
2208:         cone[1] = 22;
2209:         cone[2] = 24;
2210:         cone[3] = 15;
2211:         DMPlexSetCone(dm, 3, cone);
2212:         ornt[0] = -1;
2213:         ornt[1] = -1;
2214:         ornt[2] = 0;
2215:         ornt[3] = -1;
2216:         DMPlexSetConeOrientation(dm, 3, ornt);
2217:         /* Cell 4 */
2218:         cone[0] = 16;
2219:         cone[1] = 24;
2220:         cone[2] = 21;
2221:         cone[3] = 25;
2222:         DMPlexSetCone(dm, 4, cone);
2223:         ornt[0] = -1;
2224:         ornt[1] = -1;
2225:         ornt[2] = -1;
2226:         ornt[3] = 0;
2227:         DMPlexSetConeOrientation(dm, 4, ornt);
2228:         /* Cell 5 */
2229:         cone[0] = 20;
2230:         cone[1] = 17;
2231:         cone[2] = 25;
2232:         cone[3] = 23;
2233:         DMPlexSetCone(dm, 5, cone);
2234:         ornt[0] = -1;
2235:         ornt[1] = -1;
2236:         ornt[2] = -1;
2237:         ornt[3] = -1;
2238:         DMPlexSetConeOrientation(dm, 5, ornt);
2239:         /* Edges */
2240:         cone[0] = 6;
2241:         cone[1] = 7;
2242:         DMPlexSetCone(dm, 14, cone);
2243:         cone[0] = 7;
2244:         cone[1] = 8;
2245:         DMPlexSetCone(dm, 15, cone);
2246:         cone[0] = 8;
2247:         cone[1] = 9;
2248:         DMPlexSetCone(dm, 16, cone);
2249:         cone[0] = 9;
2250:         cone[1] = 6;
2251:         DMPlexSetCone(dm, 17, cone);
2252:         cone[0] = 10;
2253:         cone[1] = 11;
2254:         DMPlexSetCone(dm, 18, cone);
2255:         cone[0] = 11;
2256:         cone[1] = 7;
2257:         DMPlexSetCone(dm, 19, cone);
2258:         cone[0] = 6;
2259:         cone[1] = 10;
2260:         DMPlexSetCone(dm, 20, cone);
2261:         cone[0] = 12;
2262:         cone[1] = 13;
2263:         DMPlexSetCone(dm, 21, cone);
2264:         cone[0] = 13;
2265:         cone[1] = 11;
2266:         DMPlexSetCone(dm, 22, cone);
2267:         cone[0] = 10;
2268:         cone[1] = 12;
2269:         DMPlexSetCone(dm, 23, cone);
2270:         cone[0] = 13;
2271:         cone[1] = 8;
2272:         DMPlexSetCone(dm, 24, cone);
2273:         cone[0] = 12;
2274:         cone[1] = 9;
2275:         DMPlexSetCone(dm, 25, cone);
2276:       }
2277:       DMPlexSymmetrize(dm);
2278:       DMPlexStratify(dm);
2279:       /* Build coordinates */
2280:       PetscCalloc1(numVerts * embedDim, &coordsIn);
2281:       if (rank == 0) {
2282:         coordsIn[0 * embedDim + 0] = -R;
2283:         coordsIn[0 * embedDim + 1] = R;
2284:         coordsIn[0 * embedDim + 2] = -R;
2285:         coordsIn[1 * embedDim + 0] = R;
2286:         coordsIn[1 * embedDim + 1] = R;
2287:         coordsIn[1 * embedDim + 2] = -R;
2288:         coordsIn[2 * embedDim + 0] = R;
2289:         coordsIn[2 * embedDim + 1] = -R;
2290:         coordsIn[2 * embedDim + 2] = -R;
2291:         coordsIn[3 * embedDim + 0] = -R;
2292:         coordsIn[3 * embedDim + 1] = -R;
2293:         coordsIn[3 * embedDim + 2] = -R;
2294:         coordsIn[4 * embedDim + 0] = -R;
2295:         coordsIn[4 * embedDim + 1] = R;
2296:         coordsIn[4 * embedDim + 2] = R;
2297:         coordsIn[5 * embedDim + 0] = R;
2298:         coordsIn[5 * embedDim + 1] = R;
2299:         coordsIn[5 * embedDim + 2] = R;
2300:         coordsIn[6 * embedDim + 0] = -R;
2301:         coordsIn[6 * embedDim + 1] = -R;
2302:         coordsIn[6 * embedDim + 2] = R;
2303:         coordsIn[7 * embedDim + 0] = R;
2304:         coordsIn[7 * embedDim + 1] = -R;
2305:         coordsIn[7 * embedDim + 2] = R;
2306:       }
2307:     }
2308:     break;
2309:   case 3:
2310:     if (simplex) {
2311:       const PetscReal edgeLen         = 1.0 / PETSC_PHI;
2312:       PetscReal       vertexA[4]      = {0.5, 0.5, 0.5, 0.5};
2313:       PetscReal       vertexB[4]      = {1.0, 0.0, 0.0, 0.0};
2314:       PetscReal       vertexC[4]      = {0.5, 0.5 * PETSC_PHI, 0.5 / PETSC_PHI, 0.0};
2315:       const PetscInt  degree          = 12;
2316:       PetscInt        s[4]            = {1, 1, 1};
2317:       PetscInt        evenPerm[12][4] = {
2318:         {0, 1, 2, 3},
2319:         {0, 2, 3, 1},
2320:         {0, 3, 1, 2},
2321:         {1, 0, 3, 2},
2322:         {1, 2, 0, 3},
2323:         {1, 3, 2, 0},
2324:         {2, 0, 1, 3},
2325:         {2, 1, 3, 0},
2326:         {2, 3, 0, 1},
2327:         {3, 0, 2, 1},
2328:         {3, 1, 0, 2},
2329:         {3, 2, 1, 0}
2330:       };
2331:       PetscInt  cone[4];
2332:       PetscInt *graph, p, i, j, k, l;

2334:       vertexA[0] *= R;
2335:       vertexA[1] *= R;
2336:       vertexA[2] *= R;
2337:       vertexA[3] *= R;
2338:       vertexB[0] *= R;
2339:       vertexB[1] *= R;
2340:       vertexB[2] *= R;
2341:       vertexB[3] *= R;
2342:       vertexC[0] *= R;
2343:       vertexC[1] *= R;
2344:       vertexC[2] *= R;
2345:       vertexC[3] *= R;
2346:       numCells    = rank == 0 ? 600 : 0;
2347:       numVerts    = rank == 0 ? 120 : 0;
2348:       firstVertex = numCells;
2349:       /* Use the 600-cell, which for a unit sphere has coordinates which are

2351:            1/2 (\pm 1, \pm 1,    \pm 1, \pm 1)                          16
2352:                (\pm 1,    0,       0,      0)  all cyclic permutations   8
2353:            1/2 (\pm 1, \pm phi, \pm 1/phi, 0)  all even permutations    96

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

2358:          http://buzzard.pugetsound.edu/sage-practice/ch03s03.html
2359:          http://mathworld.wolfram.com/600-Cell.html
2360:       */
2361:       /* Construct vertices */
2362:       PetscCalloc1(numVerts * embedDim, &coordsIn);
2363:       i = 0;
2364:       if (rank == 0) {
2365:         for (s[0] = -1; s[0] < 2; s[0] += 2) {
2366:           for (s[1] = -1; s[1] < 2; s[1] += 2) {
2367:             for (s[2] = -1; s[2] < 2; s[2] += 2) {
2368:               for (s[3] = -1; s[3] < 2; s[3] += 2) {
2369:                 for (d = 0; d < embedDim; ++d) coordsIn[i * embedDim + d] = s[d] * vertexA[d];
2370:                 ++i;
2371:               }
2372:             }
2373:           }
2374:         }
2375:         for (p = 0; p < embedDim; ++p) {
2376:           s[1] = s[2] = s[3] = 1;
2377:           for (s[0] = -1; s[0] < 2; s[0] += 2) {
2378:             for (d = 0; d < embedDim; ++d) coordsIn[i * embedDim + d] = s[(d + p) % embedDim] * vertexB[(d + p) % embedDim];
2379:             ++i;
2380:           }
2381:         }
2382:         for (p = 0; p < 12; ++p) {
2383:           s[3] = 1;
2384:           for (s[0] = -1; s[0] < 2; s[0] += 2) {
2385:             for (s[1] = -1; s[1] < 2; s[1] += 2) {
2386:               for (s[2] = -1; s[2] < 2; s[2] += 2) {
2387:                 for (d = 0; d < embedDim; ++d) coordsIn[i * embedDim + d] = s[evenPerm[p][d]] * vertexC[evenPerm[p][d]];
2388:                 ++i;
2389:               }
2390:             }
2391:           }
2392:         }
2393:       }
2395:       /* Construct graph */
2396:       PetscCalloc1(numVerts * numVerts, &graph);
2397:       for (i = 0; i < numVerts; ++i) {
2398:         for (j = 0, k = 0; j < numVerts; ++j) {
2399:           if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i * embedDim], &coordsIn[j * embedDim]) - edgeLen) < PETSC_SMALL) {
2400:             graph[i * numVerts + j] = 1;
2401:             ++k;
2402:           }
2403:         }
2405:       }
2406:       /* Build Topology */
2407:       DMPlexSetChart(dm, 0, numCells + numVerts);
2408:       for (c = 0; c < numCells; c++) DMPlexSetConeSize(dm, c, embedDim);
2409:       DMSetUp(dm); /* Allocate space for cones */
2410:       /* Cells */
2411:       if (rank == 0) {
2412:         for (i = 0, c = 0; i < numVerts; ++i) {
2413:           for (j = 0; j < i; ++j) {
2414:             for (k = 0; k < j; ++k) {
2415:               for (l = 0; l < k; ++l) {
2416:                 if (graph[i * numVerts + j] && graph[j * numVerts + k] && graph[k * numVerts + i] && graph[l * numVerts + i] && graph[l * numVerts + j] && graph[l * numVerts + k]) {
2417:                   cone[0] = firstVertex + i;
2418:                   cone[1] = firstVertex + j;
2419:                   cone[2] = firstVertex + k;
2420:                   cone[3] = firstVertex + l;
2421:                   /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */
2422:                   {
2423:                     const PetscInt epsilon[4][4][4][4] = {
2424:                       {{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}},  {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, -1, 0}}, {{0, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 0, 0}, {0, 1, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 1, 0}, {0, -1, 0, 0}, {0, 0, 0, 0}}},

2426:                       {{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 1, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}},  {{0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}}, {{0, 0, -1, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}}},

2428:                       {{{0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, -1, 0, 0}}, {{0, 0, 0, -1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}},  {{0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}},

2430:                       {{{0, 0, 0, 0}, {0, 0, -1, 0}, {0, 1, 0, 0}, {0, 0, 0, 0}}, {{0, 0, 1, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0}}, {{0, -1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}} }
2431:                     };
2432:                     PetscReal normal[4];
2433:                     PetscInt  e, f, g;

2435:                     for (d = 0; d < embedDim; ++d) {
2436:                       normal[d] = 0.0;
2437:                       for (e = 0; e < embedDim; ++e) {
2438:                         for (f = 0; f < embedDim; ++f) {
2439:                           for (g = 0; g < embedDim; ++g) {
2440:                             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]);
2441:                           }
2442:                         }
2443:                       }
2444:                     }
2445:                     if (DotReal(embedDim, normal, &coordsIn[i * embedDim]) < 0) {
2446:                       PetscInt tmp = cone[1];
2447:                       cone[1]      = cone[2];
2448:                       cone[2]      = tmp;
2449:                     }
2450:                   }
2451:                   DMPlexSetCone(dm, c++, cone);
2452:                 }
2453:               }
2454:             }
2455:           }
2456:         }
2457:       }
2458:       DMPlexSymmetrize(dm);
2459:       DMPlexStratify(dm);
2460:       PetscFree(graph);
2461:     }
2462:     break;
2463:   default:
2464:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension for sphere: %" PetscInt_FMT, dim);
2465:   }
2466:   /* Create coordinates */
2467:   DMGetCoordinateSection(dm, &coordSection);
2468:   PetscSectionSetNumFields(coordSection, 1);
2469:   PetscSectionSetFieldComponents(coordSection, 0, embedDim);
2470:   PetscSectionSetChart(coordSection, firstVertex, firstVertex + numVerts);
2471:   for (v = firstVertex; v < firstVertex + numVerts; ++v) {
2472:     PetscSectionSetDof(coordSection, v, embedDim);
2473:     PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
2474:   }
2475:   PetscSectionSetUp(coordSection);
2476:   PetscSectionGetStorageSize(coordSection, &coordSize);
2477:   VecCreate(PETSC_COMM_SELF, &coordinates);
2478:   VecSetBlockSize(coordinates, embedDim);
2479:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
2480:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
2481:   VecSetType(coordinates, VECSTANDARD);
2482:   VecGetArray(coordinates, &coords);
2483:   for (v = 0; v < numVerts; ++v)
2484:     for (d = 0; d < embedDim; ++d) coords[v * embedDim + d] = coordsIn[v * embedDim + d];
2485:   VecRestoreArray(coordinates, &coords);
2486:   DMSetCoordinatesLocal(dm, coordinates);
2487:   VecDestroy(&coordinates);
2488:   PetscFree(coordsIn);
2489:   {
2490:     DM          cdm;
2491:     PetscDS     cds;
2492:     PetscScalar c = R;

2494:     DMPlexCreateCoordinateSpace(dm, 1, snapToSphere);
2495:     DMGetCoordinateDM(dm, &cdm);
2496:     DMGetDS(cdm, &cds);
2497:     PetscDSSetConstants(cds, 1, &c);
2498:   }
2499:   /* Wait for coordinate creation before doing in-place modification */
2500:   if (simplex) DMPlexInterpolateInPlace_Internal(dm);
2501:   return 0;
2502: }

2504: typedef void (*TPSEvaluateFunc)(const PetscReal[], PetscReal *, PetscReal[], PetscReal (*)[3]);

2506: /*
2507:  The Schwarz P implicit surface is

2509:      f(x) = cos(x0) + cos(x1) + cos(x2) = 0
2510: */
2511: static void TPSEvaluate_SchwarzP(const PetscReal y[3], PetscReal *f, PetscReal grad[], PetscReal (*hess)[3])
2512: {
2513:   PetscReal c[3] = {PetscCosReal(y[0] * PETSC_PI), PetscCosReal(y[1] * PETSC_PI), PetscCosReal(y[2] * PETSC_PI)};
2514:   PetscReal g[3] = {-PetscSinReal(y[0] * PETSC_PI), -PetscSinReal(y[1] * PETSC_PI), -PetscSinReal(y[2] * PETSC_PI)};
2515:   f[0]           = c[0] + c[1] + c[2];
2516:   for (PetscInt i = 0; i < 3; i++) {
2517:     grad[i] = PETSC_PI * g[i];
2518:     for (PetscInt j = 0; j < 3; j++) hess[i][j] = (i == j) ? -PetscSqr(PETSC_PI) * c[i] : 0.;
2519:   }
2520: }

2522: // u[] is a tentative normal on input. Replace with the implicit function gradient in the same direction
2523: static PetscErrorCode TPSExtrudeNormalFunc_SchwarzP(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx)
2524: {
2525:   for (PetscInt i = 0; i < 3; i++) u[i] = -PETSC_PI * PetscSinReal(x[i] * PETSC_PI);
2526:   return 0;
2527: }

2529: /*
2530:  The Gyroid implicit surface is

2532:  f(x,y,z) = sin(pi * x) * cos (pi * (y + 1/2))  + sin(pi * (y + 1/2)) * cos(pi * (z + 1/4)) + sin(pi * (z + 1/4)) * cos(pi * x)

2534: */
2535: static void TPSEvaluate_Gyroid(const PetscReal y[3], PetscReal *f, PetscReal grad[], PetscReal (*hess)[3])
2536: {
2537:   PetscReal s[3] = {PetscSinReal(PETSC_PI * y[0]), PetscSinReal(PETSC_PI * (y[1] + .5)), PetscSinReal(PETSC_PI * (y[2] + .25))};
2538:   PetscReal c[3] = {PetscCosReal(PETSC_PI * y[0]), PetscCosReal(PETSC_PI * (y[1] + .5)), PetscCosReal(PETSC_PI * (y[2] + .25))};
2539:   f[0]           = s[0] * c[1] + s[1] * c[2] + s[2] * c[0];
2540:   grad[0]        = PETSC_PI * (c[0] * c[1] - s[2] * s[0]);
2541:   grad[1]        = PETSC_PI * (c[1] * c[2] - s[0] * s[1]);
2542:   grad[2]        = PETSC_PI * (c[2] * c[0] - s[1] * s[2]);
2543:   hess[0][0]     = -PetscSqr(PETSC_PI) * (s[0] * c[1] + s[2] * c[0]);
2544:   hess[0][1]     = -PetscSqr(PETSC_PI) * (c[0] * s[1]);
2545:   hess[0][2]     = -PetscSqr(PETSC_PI) * (c[2] * s[0]);
2546:   hess[1][0]     = -PetscSqr(PETSC_PI) * (s[1] * c[2] + s[0] * c[1]);
2547:   hess[1][1]     = -PetscSqr(PETSC_PI) * (c[1] * s[2]);
2548:   hess[2][2]     = -PetscSqr(PETSC_PI) * (c[0] * s[1]);
2549:   hess[2][0]     = -PetscSqr(PETSC_PI) * (s[2] * c[0] + s[1] * c[2]);
2550:   hess[2][1]     = -PetscSqr(PETSC_PI) * (c[2] * s[0]);
2551:   hess[2][2]     = -PetscSqr(PETSC_PI) * (c[1] * s[2]);
2552: }

2554: // u[] is a tentative normal on input. Replace with the implicit function gradient in the same direction
2555: static PetscErrorCode TPSExtrudeNormalFunc_Gyroid(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx)
2556: {
2557:   PetscReal s[3] = {PetscSinReal(PETSC_PI * x[0]), PetscSinReal(PETSC_PI * (x[1] + .5)), PetscSinReal(PETSC_PI * (x[2] + .25))};
2558:   PetscReal c[3] = {PetscCosReal(PETSC_PI * x[0]), PetscCosReal(PETSC_PI * (x[1] + .5)), PetscCosReal(PETSC_PI * (x[2] + .25))};
2559:   u[0]           = PETSC_PI * (c[0] * c[1] - s[2] * s[0]);
2560:   u[1]           = PETSC_PI * (c[1] * c[2] - s[0] * s[1]);
2561:   u[2]           = PETSC_PI * (c[2] * c[0] - s[1] * s[2]);
2562:   return 0;
2563: }

2565: /*
2566:    We wish to solve

2568:          min_y || y - x ||^2  subject to f(y) = 0

2570:    Let g(y) = grad(f).  The minimization problem is equivalent to asking to satisfy
2571:    f(y) = 0 and (y-x) is parallel to g(y).  We do this by using Householder QR to obtain a basis for the
2572:    tangent space and ask for both components in the tangent space to be zero.

2574:    Take g to be a column vector and compute the "full QR" factorization Q R = g,
2575:    where Q = I - 2 n n^T is a symmetric orthogonal matrix.
2576:    The first column of Q is parallel to g so the remaining two columns span the null space.
2577:    Let Qn = Q[:,1:] be those remaining columns.  Then Qn Qn^T is an orthogonal projector into the tangent space.
2578:    Since Q is symmetric, this is equivalent to multipyling by Q and taking the last two entries.
2579:    In total, we have a system of 3 equations in 3 unknowns:

2581:      f(y) = 0                       1 equation
2582:      Qn^T (y - x) = 0               2 equations

2584:    Here, we compute the residual and Jacobian of this system.
2585: */
2586: static void TPSNearestPointResJac(TPSEvaluateFunc feval, const PetscScalar x[], const PetscScalar y[], PetscScalar res[], PetscScalar J[])
2587: {
2588:   PetscReal yreal[3] = {PetscRealPart(y[0]), PetscRealPart(y[1]), PetscRealPart(y[2])};
2589:   PetscReal d[3]     = {PetscRealPart(y[0] - x[0]), PetscRealPart(y[1] - x[1]), PetscRealPart(y[2] - x[2])};
2590:   PetscReal f, grad[3], n[3], norm, norm_y[3], nd, nd_y[3], sign;
2591:   PetscReal n_y[3][3] = {
2592:     {0, 0, 0},
2593:     {0, 0, 0},
2594:     {0, 0, 0}
2595:   };

2597:   feval(yreal, &f, grad, n_y);

2599:   for (PetscInt i = 0; i < 3; i++) n[i] = grad[i];
2600:   norm = PetscSqrtReal(PetscSqr(n[0]) + PetscSqr(n[1]) + PetscSqr(n[2]));
2601:   for (PetscInt i = 0; i < 3; i++) norm_y[i] = 1. / norm * n[i] * n_y[i][i];

2603:   // Define the Householder reflector
2604:   sign = n[0] >= 0 ? 1. : -1.;
2605:   n[0] += norm * sign;
2606:   for (PetscInt i = 0; i < 3; i++) n_y[0][i] += norm_y[i] * sign;

2608:   norm      = PetscSqrtReal(PetscSqr(n[0]) + PetscSqr(n[1]) + PetscSqr(n[2]));
2609:   norm_y[0] = 1. / norm * (n[0] * n_y[0][0]);
2610:   norm_y[1] = 1. / norm * (n[0] * n_y[0][1] + n[1] * n_y[1][1]);
2611:   norm_y[2] = 1. / norm * (n[0] * n_y[0][2] + n[2] * n_y[2][2]);

2613:   for (PetscInt i = 0; i < 3; i++) {
2614:     n[i] /= norm;
2615:     for (PetscInt j = 0; j < 3; j++) {
2616:       // note that n[i] is n_old[i]/norm when executing the code below
2617:       n_y[i][j] = n_y[i][j] / norm - n[i] / norm * norm_y[j];
2618:     }
2619:   }

2621:   nd = n[0] * d[0] + n[1] * d[1] + n[2] * d[2];
2622:   for (PetscInt i = 0; i < 3; i++) nd_y[i] = n[i] + n_y[0][i] * d[0] + n_y[1][i] * d[1] + n_y[2][i] * d[2];

2624:   res[0] = f;
2625:   res[1] = d[1] - 2 * n[1] * nd;
2626:   res[2] = d[2] - 2 * n[2] * nd;
2627:   // J[j][i] is J_{ij} (column major)
2628:   for (PetscInt j = 0; j < 3; j++) {
2629:     J[0 + j * 3] = grad[j];
2630:     J[1 + j * 3] = (j == 1) * 1. - 2 * (n_y[1][j] * nd + n[1] * nd_y[j]);
2631:     J[2 + j * 3] = (j == 2) * 1. - 2 * (n_y[2][j] * nd + n[2] * nd_y[j]);
2632:   }
2633: }

2635: /*
2636:    Project x to the nearest point on the implicit surface using Newton's method.
2637: */
2638: static PetscErrorCode TPSNearestPoint(TPSEvaluateFunc feval, PetscScalar x[])
2639: {
2640:   PetscScalar y[3] = {x[0], x[1], x[2]}; // Initial guess

2642:   for (PetscInt iter = 0; iter < 10; iter++) {
2643:     PetscScalar res[3], J[9];
2644:     PetscReal   resnorm;
2645:     TPSNearestPointResJac(feval, x, y, res, J);
2646:     resnorm = PetscSqrtReal(PetscSqr(PetscRealPart(res[0])) + PetscSqr(PetscRealPart(res[1])) + PetscSqr(PetscRealPart(res[2])));
2647:     if (0) { // Turn on this monitor if you need to confirm quadratic convergence
2648:       PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT "] res [%g %g %g]\n", iter, (double)PetscRealPart(res[0]), (double)PetscRealPart(res[1]), (double)PetscRealPart(res[2]));
2649:     }
2650:     if (resnorm < PETSC_SMALL) break;

2652:     // Take the Newton step
2653:     PetscKernel_A_gets_inverse_A_3(J, 0., PETSC_FALSE, NULL);
2654:     PetscKernel_v_gets_v_minus_A_times_w_3(y, J, res);
2655:   }
2656:   for (PetscInt i = 0; i < 3; i++) x[i] = y[i];
2657:   return 0;
2658: }

2660: const char *const DMPlexTPSTypes[] = {"SCHWARZ_P", "GYROID", "DMPlexTPSType", "DMPLEX_TPS_", NULL};

2662: static PetscErrorCode DMPlexCreateTPSMesh_Internal(DM dm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness)
2663: {
2664:   PetscMPIInt rank;
2665:   PetscInt    topoDim = 2, spaceDim = 3, numFaces = 0, numVertices = 0, numEdges = 0;
2666:   PetscInt(*edges)[2] = NULL, *edgeSets = NULL;
2667:   PetscInt            *cells_flat = NULL;
2668:   PetscReal           *vtxCoords  = NULL;
2669:   TPSEvaluateFunc      evalFunc   = NULL;
2670:   PetscSimplePointFunc normalFunc = NULL;
2671:   DMLabel              label;

2673:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
2675:   switch (tpstype) {
2676:   case DMPLEX_TPS_SCHWARZ_P:
2678:     if (rank == 0) {
2679:       PetscInt(*cells)[6][4][4] = NULL; // [junction, junction-face, cell, conn]
2680:       PetscInt  Njunctions = 0, Ncuts = 0, Npipes[3], vcount;
2681:       PetscReal L = 1;

2683:       Npipes[0]   = (extent[0] + 1) * extent[1] * extent[2];
2684:       Npipes[1]   = extent[0] * (extent[1] + 1) * extent[2];
2685:       Npipes[2]   = extent[0] * extent[1] * (extent[2] + 1);
2686:       Njunctions  = extent[0] * extent[1] * extent[2];
2687:       Ncuts       = 2 * (extent[0] * extent[1] + extent[1] * extent[2] + extent[2] * extent[0]);
2688:       numVertices = 4 * (Npipes[0] + Npipes[1] + Npipes[2]) + 8 * Njunctions;
2689:       PetscMalloc1(3 * numVertices, &vtxCoords);
2690:       PetscMalloc1(Njunctions, &cells);
2691:       PetscMalloc1(Ncuts * 4, &edges);
2692:       PetscMalloc1(Ncuts * 4, &edgeSets);
2693:       // x-normal pipes
2694:       vcount = 0;
2695:       for (PetscInt i = 0; i < extent[0] + 1; i++) {
2696:         for (PetscInt j = 0; j < extent[1]; j++) {
2697:           for (PetscInt k = 0; k < extent[2]; k++) {
2698:             for (PetscInt l = 0; l < 4; l++) {
2699:               vtxCoords[vcount++] = (2 * i - 1) * L;
2700:               vtxCoords[vcount++] = 2 * j * L + PetscCosReal((2 * l + 1) * PETSC_PI / 4) * L / 2;
2701:               vtxCoords[vcount++] = 2 * k * L + PetscSinReal((2 * l + 1) * PETSC_PI / 4) * L / 2;
2702:             }
2703:           }
2704:         }
2705:       }
2706:       // y-normal pipes
2707:       for (PetscInt i = 0; i < extent[0]; i++) {
2708:         for (PetscInt j = 0; j < extent[1] + 1; j++) {
2709:           for (PetscInt k = 0; k < extent[2]; k++) {
2710:             for (PetscInt l = 0; l < 4; l++) {
2711:               vtxCoords[vcount++] = 2 * i * L + PetscSinReal((2 * l + 1) * PETSC_PI / 4) * L / 2;
2712:               vtxCoords[vcount++] = (2 * j - 1) * L;
2713:               vtxCoords[vcount++] = 2 * k * L + PetscCosReal((2 * l + 1) * PETSC_PI / 4) * L / 2;
2714:             }
2715:           }
2716:         }
2717:       }
2718:       // z-normal pipes
2719:       for (PetscInt i = 0; i < extent[0]; i++) {
2720:         for (PetscInt j = 0; j < extent[1]; j++) {
2721:           for (PetscInt k = 0; k < extent[2] + 1; k++) {
2722:             for (PetscInt l = 0; l < 4; l++) {
2723:               vtxCoords[vcount++] = 2 * i * L + PetscCosReal((2 * l + 1) * PETSC_PI / 4) * L / 2;
2724:               vtxCoords[vcount++] = 2 * j * L + PetscSinReal((2 * l + 1) * PETSC_PI / 4) * L / 2;
2725:               vtxCoords[vcount++] = (2 * k - 1) * L;
2726:             }
2727:           }
2728:         }
2729:       }
2730:       // junctions
2731:       for (PetscInt i = 0; i < extent[0]; i++) {
2732:         for (PetscInt j = 0; j < extent[1]; j++) {
2733:           for (PetscInt k = 0; k < extent[2]; k++) {
2734:             const PetscInt J = (i * extent[1] + j) * extent[2] + k, Jvoff = (Npipes[0] + Npipes[1] + Npipes[2]) * 4 + J * 8;
2736:             for (PetscInt ii = 0; ii < 2; ii++) {
2737:               for (PetscInt jj = 0; jj < 2; jj++) {
2738:                 for (PetscInt kk = 0; kk < 2; kk++) {
2739:                   double Ls           = (1 - sqrt(2) / 4) * L;
2740:                   vtxCoords[vcount++] = 2 * i * L + (2 * ii - 1) * Ls;
2741:                   vtxCoords[vcount++] = 2 * j * L + (2 * jj - 1) * Ls;
2742:                   vtxCoords[vcount++] = 2 * k * L + (2 * kk - 1) * Ls;
2743:                 }
2744:               }
2745:             }
2746:             const PetscInt jfaces[3][2][4] = {
2747:               {{3, 1, 0, 2}, {7, 5, 4, 6}}, // x-aligned
2748:               {{5, 4, 0, 1}, {7, 6, 2, 3}}, // y-aligned
2749:               {{6, 2, 0, 4}, {7, 3, 1, 5}}  // z-aligned
2750:             };
2751:             const PetscInt pipe_lo[3] = {// vertex numbers of pipes
2752:                                          ((i * extent[1] + j) * extent[2] + k) * 4, ((i * (extent[1] + 1) + j) * extent[2] + k + Npipes[0]) * 4, ((i * extent[1] + j) * (extent[2] + 1) + k + Npipes[0] + Npipes[1]) * 4};
2753:             const PetscInt pipe_hi[3] = {// vertex numbers of pipes
2754:                                          (((i + 1) * extent[1] + j) * extent[2] + k) * 4, ((i * (extent[1] + 1) + j + 1) * extent[2] + k + Npipes[0]) * 4, ((i * extent[1] + j) * (extent[2] + 1) + k + 1 + Npipes[0] + Npipes[1]) * 4};
2755:             for (PetscInt dir = 0; dir < 3; dir++) { // x,y,z
2756:               const PetscInt ijk[3] = {i, j, k};
2757:               for (PetscInt l = 0; l < 4; l++) { // rotations
2758:                 cells[J][dir * 2 + 0][l][0] = pipe_lo[dir] + l;
2759:                 cells[J][dir * 2 + 0][l][1] = Jvoff + jfaces[dir][0][l];
2760:                 cells[J][dir * 2 + 0][l][2] = Jvoff + jfaces[dir][0][(l - 1 + 4) % 4];
2761:                 cells[J][dir * 2 + 0][l][3] = pipe_lo[dir] + (l - 1 + 4) % 4;
2762:                 cells[J][dir * 2 + 1][l][0] = Jvoff + jfaces[dir][1][l];
2763:                 cells[J][dir * 2 + 1][l][1] = pipe_hi[dir] + l;
2764:                 cells[J][dir * 2 + 1][l][2] = pipe_hi[dir] + (l - 1 + 4) % 4;
2765:                 cells[J][dir * 2 + 1][l][3] = Jvoff + jfaces[dir][1][(l - 1 + 4) % 4];
2766:                 if (ijk[dir] == 0) {
2767:                   edges[numEdges][0] = pipe_lo[dir] + l;
2768:                   edges[numEdges][1] = pipe_lo[dir] + (l + 1) % 4;
2769:                   edgeSets[numEdges] = dir * 2 + 1;
2770:                   numEdges++;
2771:                 }
2772:                 if (ijk[dir] + 1 == extent[dir]) {
2773:                   edges[numEdges][0] = pipe_hi[dir] + l;
2774:                   edges[numEdges][1] = pipe_hi[dir] + (l + 1) % 4;
2775:                   edgeSets[numEdges] = dir * 2 + 2;
2776:                   numEdges++;
2777:                 }
2778:               }
2779:             }
2780:           }
2781:         }
2782:       }
2784:       numFaces   = 24 * Njunctions;
2785:       cells_flat = cells[0][0][0];
2786:     }
2787:     evalFunc   = TPSEvaluate_SchwarzP;
2788:     normalFunc = TPSExtrudeNormalFunc_SchwarzP;
2789:     break;
2790:   case DMPLEX_TPS_GYROID:
2791:     if (rank == 0) {
2792:       // This is a coarse mesh approximation of the gyroid shifted to being the zero of the level set
2793:       //
2794:       //     sin(pi*x)*cos(pi*(y+1/2)) + sin(pi*(y+1/2))*cos(pi*(z+1/4)) + sin(pi*(z+1/4))*cos(x)
2795:       //
2796:       // on the cell [0,2]^3.
2797:       //
2798:       // Think about dividing that cell into four columns, and focus on the column [0,1]x[0,1]x[0,2].
2799:       // If you looked at the gyroid in that column at different slices of z you would see that it kind of spins
2800:       // like a boomerang:
2801:       //
2802:       //     z = 0          z = 1/4        z = 1/2        z = 3/4     //
2803:       //     -----          -------        -------        -------     //
2804:       //                                                              //
2805:       //     +       +      +       +      +       +      +   \   +   //
2806:       //      \                                   /            \      //
2807:       //       \            `-_   _-'            /              }     //
2808:       //        *-_            `-'            _-'              /      //
2809:       //     +     `-+      +       +      +-'     +      +   /   +   //
2810:       //                                                              //
2811:       //                                                              //
2812:       //     z = 1          z = 5/4        z = 3/2        z = 7/4     //
2813:       //     -----          -------        -------        -------     //
2814:       //                                                              //
2815:       //     +-_     +      +       +      +     _-+      +   /   +   //
2816:       //        `-_            _-_            _-`            /        //
2817:       //           \        _-'   `-_        /              {         //
2818:       //            \                       /                \        //
2819:       //     +       +      +       +      +       +      +   \   +   //
2820:       //
2821:       //
2822:       // This course mesh approximates each of these slices by two line segments,
2823:       // and then connects the segments in consecutive layers with quadrilateral faces.
2824:       // All of the end points of the segments are multiples of 1/4 except for the
2825:       // point * in the picture for z = 0 above and the similar points in other layers.
2826:       // That point is at (gamma, gamma, 0), where gamma is calculated below.
2827:       //
2828:       // The column  [1,2]x[1,2]x[0,2] looks the same as this column;
2829:       // The columns [1,2]x[0,1]x[0,2] and [0,1]x[1,2]x[0,2] are mirror images.
2830:       //
2831:       // As for how this method turned into the names given to the vertices:
2832:       // that was not systematic, it was just the way it worked out in my handwritten notes.

2834:       PetscInt facesPerBlock = 64;
2835:       PetscInt vertsPerBlock = 56;
2836:       PetscInt extentPlus[3];
2837:       PetscInt numBlocks, numBlocksPlus;
2838:       const PetscInt A = 0, B = 1, C = 2, D = 3, E = 4, F = 5, G = 6, H = 7, II = 8, J = 9, K = 10, L = 11, M = 12, N = 13, O = 14, P = 15, Q = 16, R = 17, S = 18, T = 19, U = 20, V = 21, W = 22, X = 23, Y = 24, Z = 25, Ap = 26, Bp = 27, Cp = 28, Dp = 29, Ep = 30, Fp = 31, Gp = 32, Hp = 33, Ip = 34, Jp = 35, Kp = 36, Lp = 37, Mp = 38, Np = 39, Op = 40, Pp = 41, Qp = 42, Rp = 43, Sp = 44, Tp = 45, Up = 46, Vp = 47, Wp = 48, Xp = 49, Yp = 50, Zp = 51, Aq = 52, Bq = 53, Cq = 54, Dq = 55;
2839:       const PetscInt pattern[64][4] = {
2840:   /* face to vertex within the coarse discretization of a single gyroid block */
2841:   /* layer 0 */
2842:         {A,           C,           K,           G          },
2843:         {C,           B,           II,          K          },
2844:         {D,           A,           H,           L          },
2845:         {B + 56 * 1,  D,           L,           J          },
2846:         {E,           B + 56 * 1,  J,           N          },
2847:         {A + 56 * 2,  E,           N,           H + 56 * 2 },
2848:         {F,           A + 56 * 2,  G + 56 * 2,  M          },
2849:         {B,           F,           M,           II         },
2850:  /* layer 1 */
2851:         {G,           K,           Q,           O          },
2852:         {K,           II,          P,           Q          },
2853:         {L,           H,           O + 56 * 1,  R          },
2854:         {J,           L,           R,           P          },
2855:         {N,           J,           P,           S          },
2856:         {H + 56 * 2,  N,           S,           O + 56 * 3 },
2857:         {M,           G + 56 * 2,  O + 56 * 2,  T          },
2858:         {II,          M,           T,           P          },
2859:  /* layer 2 */
2860:         {O,           Q,           Y,           U          },
2861:         {Q,           P,           W,           Y          },
2862:         {R,           O + 56 * 1,  U + 56 * 1,  Ap         },
2863:         {P,           R,           Ap,          W          },
2864:         {S,           P,           X,           Bp         },
2865:         {O + 56 * 3,  S,           Bp,          V + 56 * 1 },
2866:         {T,           O + 56 * 2,  V,           Z          },
2867:         {P,           T,           Z,           X          },
2868:  /* layer 3 */
2869:         {U,           Y,           Ep,          Dp         },
2870:         {Y,           W,           Cp,          Ep         },
2871:         {Ap,          U + 56 * 1,  Dp + 56 * 1, Gp         },
2872:         {W,           Ap,          Gp,          Cp         },
2873:         {Bp,          X,           Cp + 56 * 2, Fp         },
2874:         {V + 56 * 1,  Bp,          Fp,          Dp + 56 * 1},
2875:         {Z,           V,           Dp,          Hp         },
2876:         {X,           Z,           Hp,          Cp + 56 * 2},
2877:  /* layer 4 */
2878:         {Dp,          Ep,          Mp,          Kp         },
2879:         {Ep,          Cp,          Ip,          Mp         },
2880:         {Gp,          Dp + 56 * 1, Lp,          Np         },
2881:         {Cp,          Gp,          Np,          Jp         },
2882:         {Fp,          Cp + 56 * 2, Jp + 56 * 2, Pp         },
2883:         {Dp + 56 * 1, Fp,          Pp,          Lp         },
2884:         {Hp,          Dp,          Kp,          Op         },
2885:         {Cp + 56 * 2, Hp,          Op,          Ip + 56 * 2},
2886:  /* layer 5 */
2887:         {Kp,          Mp,          Sp,          Rp         },
2888:         {Mp,          Ip,          Qp,          Sp         },
2889:         {Np,          Lp,          Rp,          Tp         },
2890:         {Jp,          Np,          Tp,          Qp + 56 * 1},
2891:         {Pp,          Jp + 56 * 2, Qp + 56 * 3, Up         },
2892:         {Lp,          Pp,          Up,          Rp         },
2893:         {Op,          Kp,          Rp,          Vp         },
2894:         {Ip + 56 * 2, Op,          Vp,          Qp + 56 * 2},
2895:  /* layer 6 */
2896:         {Rp,          Sp,          Aq,          Yp         },
2897:         {Sp,          Qp,          Wp,          Aq         },
2898:         {Tp,          Rp,          Yp,          Cq         },
2899:         {Qp + 56 * 1, Tp,          Cq,          Wp + 56 * 1},
2900:         {Up,          Qp + 56 * 3, Xp + 56 * 1, Dq         },
2901:         {Rp,          Up,          Dq,          Zp         },
2902:         {Vp,          Rp,          Zp,          Bq         },
2903:         {Qp + 56 * 2, Vp,          Bq,          Xp         },
2904:  /* layer 7 (the top is the periodic image of the bottom of layer 0) */
2905:         {Yp,          Aq,          C + 56 * 4,  A + 56 * 4 },
2906:         {Aq,          Wp,          B + 56 * 4,  C + 56 * 4 },
2907:         {Cq,          Yp,          A + 56 * 4,  D + 56 * 4 },
2908:         {Wp + 56 * 1, Cq,          D + 56 * 4,  B + 56 * 5 },
2909:         {Dq,          Xp + 56 * 1, B + 56 * 5,  E + 56 * 4 },
2910:         {Zp,          Dq,          E + 56 * 4,  A + 56 * 6 },
2911:         {Bq,          Zp,          A + 56 * 6,  F + 56 * 4 },
2912:         {Xp,          Bq,          F + 56 * 4,  B + 56 * 4 }
2913:       };
2914:       const PetscReal gamma                = PetscAcosReal((PetscSqrtReal(3.) - 1.) / PetscSqrtReal(2.)) / PETSC_PI;
2915:       const PetscReal patternCoords[56][3] = {
2916:         {1.,        0.,        0.  }, /* A  */
2917:         {0.,        1.,        0.  }, /* B  */
2918:         {gamma,     gamma,     0.  }, /* C  */
2919:         {1 + gamma, 1 - gamma, 0.  }, /* D  */
2920:         {2 - gamma, 2 - gamma, 0.  }, /* E  */
2921:         {1 - gamma, 1 + gamma, 0.  }, /* F  */

2923:         {.5,        0,         .25 }, /* G  */
2924:         {1.5,       0.,        .25 }, /* H  */
2925:         {.5,        1.,        .25 }, /* II */
2926:         {1.5,       1.,        .25 }, /* J  */
2927:         {.25,       .5,        .25 }, /* K  */
2928:         {1.25,      .5,        .25 }, /* L  */
2929:         {.75,       1.5,       .25 }, /* M  */
2930:         {1.75,      1.5,       .25 }, /* N  */

2932:         {0.,        0.,        .5  }, /* O  */
2933:         {1.,        1.,        .5  }, /* P  */
2934:         {gamma,     1 - gamma, .5  }, /* Q  */
2935:         {1 + gamma, gamma,     .5  }, /* R  */
2936:         {2 - gamma, 1 + gamma, .5  }, /* S  */
2937:         {1 - gamma, 2 - gamma, .5  }, /* T  */

2939:         {0.,        .5,        .75 }, /* U  */
2940:         {0.,        1.5,       .75 }, /* V  */
2941:         {1.,        .5,        .75 }, /* W  */
2942:         {1.,        1.5,       .75 }, /* X  */
2943:         {.5,        .75,       .75 }, /* Y  */
2944:         {.5,        1.75,      .75 }, /* Z  */
2945:         {1.5,       .25,       .75 }, /* Ap */
2946:         {1.5,       1.25,      .75 }, /* Bp */

2948:         {1.,        0.,        1.  }, /* Cp */
2949:         {0.,        1.,        1.  }, /* Dp */
2950:         {1 - gamma, 1 - gamma, 1.  }, /* Ep */
2951:         {1 + gamma, 1 + gamma, 1.  }, /* Fp */
2952:         {2 - gamma, gamma,     1.  }, /* Gp */
2953:         {gamma,     2 - gamma, 1.  }, /* Hp */

2955:         {.5,        0.,        1.25}, /* Ip */
2956:         {1.5,       0.,        1.25}, /* Jp */
2957:         {.5,        1.,        1.25}, /* Kp */
2958:         {1.5,       1.,        1.25}, /* Lp */
2959:         {.75,       .5,        1.25}, /* Mp */
2960:         {1.75,      .5,        1.25}, /* Np */
2961:         {.25,       1.5,       1.25}, /* Op */
2962:         {1.25,      1.5,       1.25}, /* Pp */

2964:         {0.,        0.,        1.5 }, /* Qp */
2965:         {1.,        1.,        1.5 }, /* Rp */
2966:         {1 - gamma, gamma,     1.5 }, /* Sp */
2967:         {2 - gamma, 1 - gamma, 1.5 }, /* Tp */
2968:         {1 + gamma, 2 - gamma, 1.5 }, /* Up */
2969:         {gamma,     1 + gamma, 1.5 }, /* Vp */

2971:         {0.,        .5,        1.75}, /* Wp */
2972:         {0.,        1.5,       1.75}, /* Xp */
2973:         {1.,        .5,        1.75}, /* Yp */
2974:         {1.,        1.5,       1.75}, /* Zp */
2975:         {.5,        .25,       1.75}, /* Aq */
2976:         {.5,        1.25,      1.75}, /* Bq */
2977:         {1.5,       .75,       1.75}, /* Cq */
2978:         {1.5,       1.75,      1.75}, /* Dq */
2979:       };
2980:       PetscInt(*cells)[64][4] = NULL;
2981:       PetscBool *seen;
2982:       PetscInt  *vertToTrueVert;
2983:       PetscInt   count;

2985:       for (PetscInt i = 0; i < 3; i++) extentPlus[i] = extent[i] + 1;
2986:       numBlocks = 1;
2987:       for (PetscInt i = 0; i < 3; i++) numBlocks *= extent[i];
2988:       numBlocksPlus = 1;
2989:       for (PetscInt i = 0; i < 3; i++) numBlocksPlus *= extentPlus[i];
2990:       numFaces = numBlocks * facesPerBlock;
2991:       PetscMalloc1(numBlocks, &cells);
2992:       PetscCalloc1(numBlocksPlus * vertsPerBlock, &seen);
2993:       for (PetscInt k = 0; k < extent[2]; k++) {
2994:         for (PetscInt j = 0; j < extent[1]; j++) {
2995:           for (PetscInt i = 0; i < extent[0]; i++) {
2996:             for (PetscInt f = 0; f < facesPerBlock; f++) {
2997:               for (PetscInt v = 0; v < 4; v++) {
2998:                 PetscInt vertRaw     = pattern[f][v];
2999:                 PetscInt blockidx    = vertRaw / 56;
3000:                 PetscInt patternvert = vertRaw % 56;
3001:                 PetscInt xplus       = (blockidx & 1);
3002:                 PetscInt yplus       = (blockidx & 2) >> 1;
3003:                 PetscInt zplus       = (blockidx & 4) >> 2;
3004:                 PetscInt zcoord      = (periodic && periodic[2] == DM_BOUNDARY_PERIODIC) ? ((k + zplus) % extent[2]) : (k + zplus);
3005:                 PetscInt ycoord      = (periodic && periodic[1] == DM_BOUNDARY_PERIODIC) ? ((j + yplus) % extent[1]) : (j + yplus);
3006:                 PetscInt xcoord      = (periodic && periodic[0] == DM_BOUNDARY_PERIODIC) ? ((i + xplus) % extent[0]) : (i + xplus);
3007:                 PetscInt vert        = ((zcoord * extentPlus[1] + ycoord) * extentPlus[0] + xcoord) * 56 + patternvert;

3009:                 cells[(k * extent[1] + j) * extent[0] + i][f][v] = vert;
3010:                 seen[vert]                                       = PETSC_TRUE;
3011:               }
3012:             }
3013:           }
3014:         }
3015:       }
3016:       for (PetscInt i = 0; i < numBlocksPlus * vertsPerBlock; i++)
3017:         if (seen[i]) numVertices++;
3018:       count = 0;
3019:       PetscMalloc1(numBlocksPlus * vertsPerBlock, &vertToTrueVert);
3020:       PetscMalloc1(numVertices * 3, &vtxCoords);
3021:       for (PetscInt i = 0; i < numBlocksPlus * vertsPerBlock; i++) vertToTrueVert[i] = -1;
3022:       for (PetscInt k = 0; k < extentPlus[2]; k++) {
3023:         for (PetscInt j = 0; j < extentPlus[1]; j++) {
3024:           for (PetscInt i = 0; i < extentPlus[0]; i++) {
3025:             for (PetscInt v = 0; v < vertsPerBlock; v++) {
3026:               PetscInt vIdx = ((k * extentPlus[1] + j) * extentPlus[0] + i) * vertsPerBlock + v;

3028:               if (seen[vIdx]) {
3029:                 PetscInt thisVert;

3031:                 vertToTrueVert[vIdx] = thisVert = count++;

3033:                 for (PetscInt d = 0; d < 3; d++) vtxCoords[3 * thisVert + d] = patternCoords[v][d];
3034:                 vtxCoords[3 * thisVert + 0] += i * 2;
3035:                 vtxCoords[3 * thisVert + 1] += j * 2;
3036:                 vtxCoords[3 * thisVert + 2] += k * 2;
3037:               }
3038:             }
3039:           }
3040:         }
3041:       }
3042:       for (PetscInt i = 0; i < numBlocks; i++) {
3043:         for (PetscInt f = 0; f < facesPerBlock; f++) {
3044:           for (PetscInt v = 0; v < 4; v++) cells[i][f][v] = vertToTrueVert[cells[i][f][v]];
3045:         }
3046:       }
3047:       PetscFree(vertToTrueVert);
3048:       PetscFree(seen);
3049:       cells_flat = cells[0][0];
3050:       numEdges   = 0;
3051:       for (PetscInt i = 0; i < numFaces; i++) {
3052:         for (PetscInt e = 0; e < 4; e++) {
3053:           PetscInt         ev[]       = {cells_flat[i * 4 + e], cells_flat[i * 4 + ((e + 1) % 4)]};
3054:           const PetscReal *evCoords[] = {&vtxCoords[3 * ev[0]], &vtxCoords[3 * ev[1]]};

3056:           for (PetscInt d = 0; d < 3; d++) {
3057:             if (!periodic || periodic[0] != DM_BOUNDARY_PERIODIC) {
3058:               if (evCoords[0][d] == 0. && evCoords[1][d] == 0.) numEdges++;
3059:               if (evCoords[0][d] == 2. * extent[d] && evCoords[1][d] == 2. * extent[d]) numEdges++;
3060:             }
3061:           }
3062:         }
3063:       }
3064:       PetscMalloc1(numEdges, &edges);
3065:       PetscMalloc1(numEdges, &edgeSets);
3066:       for (PetscInt edge = 0, i = 0; i < numFaces; i++) {
3067:         for (PetscInt e = 0; e < 4; e++) {
3068:           PetscInt         ev[]       = {cells_flat[i * 4 + e], cells_flat[i * 4 + ((e + 1) % 4)]};
3069:           const PetscReal *evCoords[] = {&vtxCoords[3 * ev[0]], &vtxCoords[3 * ev[1]]};

3071:           for (PetscInt d = 0; d < 3; d++) {
3072:             if (!periodic || periodic[d] != DM_BOUNDARY_PERIODIC) {
3073:               if (evCoords[0][d] == 0. && evCoords[1][d] == 0.) {
3074:                 edges[edge][0]   = ev[0];
3075:                 edges[edge][1]   = ev[1];
3076:                 edgeSets[edge++] = 2 * d;
3077:               }
3078:               if (evCoords[0][d] == 2. * extent[d] && evCoords[1][d] == 2. * extent[d]) {
3079:                 edges[edge][0]   = ev[0];
3080:                 edges[edge][1]   = ev[1];
3081:                 edgeSets[edge++] = 2 * d + 1;
3082:               }
3083:             }
3084:           }
3085:         }
3086:       }
3087:     }
3088:     evalFunc   = TPSEvaluate_Gyroid;
3089:     normalFunc = TPSExtrudeNormalFunc_Gyroid;
3090:     break;
3091:   }

3093:   DMSetDimension(dm, topoDim);
3094:   if (rank == 0) DMPlexBuildFromCellList(dm, numFaces, numVertices, 4, cells_flat);
3095:   else DMPlexBuildFromCellList(dm, 0, 0, 0, NULL);
3096:   PetscFree(cells_flat);
3097:   {
3098:     DM idm;
3099:     DMPlexInterpolate(dm, &idm);
3100:     DMPlexReplace_Internal(dm, &idm);
3101:   }
3102:   if (rank == 0) DMPlexBuildCoordinatesFromCellList(dm, spaceDim, vtxCoords);
3103:   else DMPlexBuildCoordinatesFromCellList(dm, spaceDim, NULL);
3104:   PetscFree(vtxCoords);

3106:   DMCreateLabel(dm, "Face Sets");
3107:   DMGetLabel(dm, "Face Sets", &label);
3108:   for (PetscInt e = 0; e < numEdges; e++) {
3109:     PetscInt        njoin;
3110:     const PetscInt *join, verts[] = {numFaces + edges[e][0], numFaces + edges[e][1]};
3111:     DMPlexGetJoin(dm, 2, verts, &njoin, &join);
3113:     DMLabelSetValue(label, join[0], edgeSets[e]);
3114:     DMPlexRestoreJoin(dm, 2, verts, &njoin, &join);
3115:   }
3116:   PetscFree(edges);
3117:   PetscFree(edgeSets);
3118:   if (tps_distribute) {
3119:     DM               pdm = NULL;
3120:     PetscPartitioner part;

3122:     DMPlexGetPartitioner(dm, &part);
3123:     PetscPartitionerSetFromOptions(part);
3124:     DMPlexDistribute(dm, 0, NULL, &pdm);
3125:     if (pdm) DMPlexReplace_Internal(dm, &pdm);
3126:     // Do not auto-distribute again
3127:     DMPlexDistributeSetDefault(dm, PETSC_FALSE);
3128:   }

3130:   DMPlexSetRefinementUniform(dm, PETSC_TRUE);
3131:   for (PetscInt refine = 0; refine < refinements; refine++) {
3132:     PetscInt     m;
3133:     DM           dmf;
3134:     Vec          X;
3135:     PetscScalar *x;
3136:     DMRefine(dm, MPI_COMM_NULL, &dmf);
3137:     DMPlexReplace_Internal(dm, &dmf);

3139:     DMGetCoordinatesLocal(dm, &X);
3140:     VecGetLocalSize(X, &m);
3141:     VecGetArray(X, &x);
3142:     for (PetscInt i = 0; i < m; i += 3) TPSNearestPoint(evalFunc, &x[i]);
3143:     VecRestoreArray(X, &x);
3144:   }

3146:   // Face Sets has already been propagated to new vertices during refinement; this propagates to the initial vertices.
3147:   DMGetLabel(dm, "Face Sets", &label);
3148:   DMPlexLabelComplete(dm, label);

3150:   if (thickness > 0) {
3151:     DM              edm, cdm, ecdm;
3152:     DMPlexTransform tr;
3153:     const char     *prefix;
3154:     PetscOptions    options;
3155:     // Code from DMPlexExtrude
3156:     DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr);
3157:     DMPlexTransformSetDM(tr, dm);
3158:     DMPlexTransformSetType(tr, DMPLEXEXTRUDE);
3159:     PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
3160:     PetscObjectSetOptionsPrefix((PetscObject)tr, prefix);
3161:     PetscObjectGetOptions((PetscObject)dm, &options);
3162:     PetscObjectSetOptions((PetscObject)tr, options);
3163:     DMPlexTransformExtrudeSetLayers(tr, layers);
3164:     DMPlexTransformExtrudeSetThickness(tr, thickness);
3165:     DMPlexTransformExtrudeSetTensor(tr, PETSC_FALSE);
3166:     DMPlexTransformExtrudeSetSymmetric(tr, PETSC_TRUE);
3167:     DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc);
3168:     DMPlexTransformSetFromOptions(tr);
3169:     PetscObjectSetOptions((PetscObject)tr, NULL);
3170:     DMPlexTransformSetUp(tr);
3171:     PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_tps_transform_view");
3172:     DMPlexTransformApply(tr, dm, &edm);
3173:     DMCopyDisc(dm, edm);
3174:     DMGetCoordinateDM(dm, &cdm);
3175:     DMGetCoordinateDM(edm, &ecdm);
3176:     DMCopyDisc(cdm, ecdm);
3177:     DMPlexTransformCreateDiscLabels(tr, edm);
3178:     DMPlexTransformDestroy(&tr);
3179:     if (edm) {
3180:       ((DM_Plex *)edm->data)->printFEM    = ((DM_Plex *)dm->data)->printFEM;
3181:       ((DM_Plex *)edm->data)->printL2     = ((DM_Plex *)dm->data)->printL2;
3182:       ((DM_Plex *)edm->data)->printLocate = ((DM_Plex *)dm->data)->printLocate;
3183:     }
3184:     DMPlexReplace_Internal(dm, &edm);
3185:   }
3186:   return 0;
3187: }

3189: /*@
3190:   DMPlexCreateTPSMesh - Create a distributed, interpolated mesh of a triply-periodic surface

3192:   Collective

3194:   Input Parameters:
3195: + comm   - The communicator for the DM object
3196: . tpstype - Type of triply-periodic surface
3197: . extent - Array of length 3 containing number of periods in each direction
3198: . periodic - array of length 3 with periodicity, or NULL for non-periodic
3199: . tps_distribute - Distribute 2D manifold mesh prior to refinement and extrusion (more scalable)
3200: . refinements - Number of factor-of-2 refinements of 2D manifold mesh
3201: . layers - Number of cell layers extruded in normal direction
3202: - thickness - Thickness in normal direction

3204:   Output Parameter:
3205: . dm  - The DM object

3207:   Notes:
3208:   This meshes the surface of the Schwarz P or Gyroid surfaces.  Schwarz P is is the simplest member of the triply-periodic minimal surfaces.
3209:   https://en.wikipedia.org/wiki/Schwarz_minimal_surface#Schwarz_P_(%22Primitive%22) and can be cut with "clean" boundaries.
3210:   The Gyroid (https://en.wikipedia.org/wiki/Gyroid) is another triply-periodic minimal surface with applications in additive manufacturing; it is much more difficult to "cut" since there are no planes of symmetry.
3211:   Our implementation creates a very coarse mesh of the surface and refines (by 4-way splitting) as many times as requested.
3212:   On each refinement, all vertices are projected to their nearest point on the surface.
3213:   This projection could readily be extended to related surfaces.

3215:   The face (edge) sets for the Schwarz P surface are numbered 1(-x), 2(+x), 3(-y), 4(+y), 5(-z), 6(+z).
3216:   When the mesh is refined, "Face Sets" contain the new vertices (created during refinement).  Use DMPlexLabelComplete() to propagate to coarse-level vertices.

3218:   References:
3219: . * - Maskery et al, Insights into the mechanical properties of several triply periodic minimal surface lattice structures made by polymer additive manufacturing, 2017. https://doi.org/10.1016/j.polymer.2017.11.049

3221:   Developer Notes:
3222:   The Gyroid mesh does not currently mark boundary sets.

3224:   Level: beginner

3226: .seealso: `DMPlexCreateSphereMesh()`, `DMSetType()`, `DMCreate()`
3227: @*/
3228: PetscErrorCode DMPlexCreateTPSMesh(MPI_Comm comm, DMPlexTPSType tpstype, const PetscInt extent[], const DMBoundaryType periodic[], PetscBool tps_distribute, PetscInt refinements, PetscInt layers, PetscReal thickness, DM *dm)
3229: {
3230:   DMCreate(comm, dm);
3231:   DMSetType(*dm, DMPLEX);
3232:   DMPlexCreateTPSMesh_Internal(*dm, tpstype, extent, periodic, tps_distribute, refinements, layers, thickness);
3233:   return 0;
3234: }

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

3239:   Collective

3241:   Input Parameters:
3242: + comm    - The communicator for the DM object
3243: . dim     - The dimension
3244: . simplex - Use simplices, or tensor product cells
3245: - R       - The radius

3247:   Output Parameter:
3248: . dm  - The DM object

3250:   Level: beginner

3252: .seealso: `DMPlexCreateBallMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()`
3253: @*/
3254: PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm)
3255: {
3257:   DMCreate(comm, dm);
3258:   DMSetType(*dm, DMPLEX);
3259:   DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R);
3260:   return 0;
3261: }

3263: static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R)
3264: {
3265:   DM      sdm, vol;
3266:   DMLabel bdlabel;

3268:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
3269:   DMSetType(sdm, DMPLEX);
3270:   PetscObjectSetOptionsPrefix((PetscObject)sdm, "bd_");
3271:   DMPlexCreateSphereMesh_Internal(sdm, dim - 1, PETSC_TRUE, R);
3272:   DMSetFromOptions(sdm);
3273:   DMViewFromOptions(sdm, NULL, "-dm_view");
3274:   DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol);
3275:   DMDestroy(&sdm);
3276:   DMPlexReplace_Internal(dm, &vol);
3277:   DMCreateLabel(dm, "marker");
3278:   DMGetLabel(dm, "marker", &bdlabel);
3279:   DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel);
3280:   DMPlexLabelComplete(dm, bdlabel);
3281:   return 0;
3282: }

3284: /*@
3285:   DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d.

3287:   Collective

3289:   Input Parameters:
3290: + comm  - The communicator for the DM object
3291: . dim   - The dimension
3292: - R     - The radius

3294:   Output Parameter:
3295: . dm  - The DM object

3297:   Options Database Keys:
3298: - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry

3300:   Level: beginner

3302: .seealso: `DMPlexCreateSphereMesh()`, `DMPlexCreateBoxMesh()`, `DMSetType()`, `DMCreate()`
3303: @*/
3304: PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm)
3305: {
3306:   DMCreate(comm, dm);
3307:   DMSetType(*dm, DMPLEX);
3308:   DMPlexCreateBallMesh_Internal(*dm, dim, R);
3309:   return 0;
3310: }

3312: static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct)
3313: {
3314:   switch (ct) {
3315:   case DM_POLYTOPE_POINT: {
3316:     PetscInt    numPoints[1]        = {1};
3317:     PetscInt    coneSize[1]         = {0};
3318:     PetscInt    cones[1]            = {0};
3319:     PetscInt    coneOrientations[1] = {0};
3320:     PetscScalar vertexCoords[1]     = {0.0};

3322:     DMSetDimension(rdm, 0);
3323:     DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3324:   } break;
3325:   case DM_POLYTOPE_SEGMENT: {
3326:     PetscInt    numPoints[2]        = {2, 1};
3327:     PetscInt    coneSize[3]         = {2, 0, 0};
3328:     PetscInt    cones[2]            = {1, 2};
3329:     PetscInt    coneOrientations[2] = {0, 0};
3330:     PetscScalar vertexCoords[2]     = {-1.0, 1.0};

3332:     DMSetDimension(rdm, 1);
3333:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3334:   } break;
3335:   case DM_POLYTOPE_POINT_PRISM_TENSOR: {
3336:     PetscInt    numPoints[2]        = {2, 1};
3337:     PetscInt    coneSize[3]         = {2, 0, 0};
3338:     PetscInt    cones[2]            = {1, 2};
3339:     PetscInt    coneOrientations[2] = {0, 0};
3340:     PetscScalar vertexCoords[2]     = {-1.0, 1.0};

3342:     DMSetDimension(rdm, 1);
3343:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3344:   } break;
3345:   case DM_POLYTOPE_TRIANGLE: {
3346:     PetscInt    numPoints[2]        = {3, 1};
3347:     PetscInt    coneSize[4]         = {3, 0, 0, 0};
3348:     PetscInt    cones[3]            = {1, 2, 3};
3349:     PetscInt    coneOrientations[3] = {0, 0, 0};
3350:     PetscScalar vertexCoords[6]     = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0};

3352:     DMSetDimension(rdm, 2);
3353:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3354:   } break;
3355:   case DM_POLYTOPE_QUADRILATERAL: {
3356:     PetscInt    numPoints[2]        = {4, 1};
3357:     PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
3358:     PetscInt    cones[4]            = {1, 2, 3, 4};
3359:     PetscInt    coneOrientations[4] = {0, 0, 0, 0};
3360:     PetscScalar vertexCoords[8]     = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};

3362:     DMSetDimension(rdm, 2);
3363:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3364:   } break;
3365:   case DM_POLYTOPE_SEG_PRISM_TENSOR: {
3366:     PetscInt    numPoints[2]        = {4, 1};
3367:     PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
3368:     PetscInt    cones[4]            = {1, 2, 3, 4};
3369:     PetscInt    coneOrientations[4] = {0, 0, 0, 0};
3370:     PetscScalar vertexCoords[8]     = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0};

3372:     DMSetDimension(rdm, 2);
3373:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3374:   } break;
3375:   case DM_POLYTOPE_TETRAHEDRON: {
3376:     PetscInt    numPoints[2]        = {4, 1};
3377:     PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
3378:     PetscInt    cones[4]            = {1, 2, 3, 4};
3379:     PetscInt    coneOrientations[4] = {0, 0, 0, 0};
3380:     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};

3382:     DMSetDimension(rdm, 3);
3383:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3384:   } break;
3385:   case DM_POLYTOPE_HEXAHEDRON: {
3386:     PetscInt    numPoints[2]        = {8, 1};
3387:     PetscInt    coneSize[9]         = {8, 0, 0, 0, 0, 0, 0, 0, 0};
3388:     PetscInt    cones[8]            = {1, 2, 3, 4, 5, 6, 7, 8};
3389:     PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
3390:     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, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};

3392:     DMSetDimension(rdm, 3);
3393:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3394:   } break;
3395:   case DM_POLYTOPE_TRI_PRISM: {
3396:     PetscInt    numPoints[2]        = {6, 1};
3397:     PetscInt    coneSize[7]         = {6, 0, 0, 0, 0, 0, 0};
3398:     PetscInt    cones[6]            = {1, 2, 3, 4, 5, 6};
3399:     PetscInt    coneOrientations[6] = {0, 0, 0, 0, 0, 0};
3400:     PetscScalar vertexCoords[18]    = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0};

3402:     DMSetDimension(rdm, 3);
3403:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3404:   } break;
3405:   case DM_POLYTOPE_TRI_PRISM_TENSOR: {
3406:     PetscInt    numPoints[2]        = {6, 1};
3407:     PetscInt    coneSize[7]         = {6, 0, 0, 0, 0, 0, 0};
3408:     PetscInt    cones[6]            = {1, 2, 3, 4, 5, 6};
3409:     PetscInt    coneOrientations[6] = {0, 0, 0, 0, 0, 0};
3410:     PetscScalar vertexCoords[18]    = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0};

3412:     DMSetDimension(rdm, 3);
3413:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3414:   } break;
3415:   case DM_POLYTOPE_QUAD_PRISM_TENSOR: {
3416:     PetscInt    numPoints[2]        = {8, 1};
3417:     PetscInt    coneSize[9]         = {8, 0, 0, 0, 0, 0, 0, 0, 0};
3418:     PetscInt    cones[8]            = {1, 2, 3, 4, 5, 6, 7, 8};
3419:     PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
3420:     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, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};

3422:     DMSetDimension(rdm, 3);
3423:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3424:   } break;
3425:   case DM_POLYTOPE_PYRAMID: {
3426:     PetscInt    numPoints[2]        = {5, 1};
3427:     PetscInt    coneSize[6]         = {5, 0, 0, 0, 0, 0};
3428:     PetscInt    cones[5]            = {1, 2, 3, 4, 5};
3429:     PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
3430:     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, 0.0, 0.0, 1.0};

3432:     DMSetDimension(rdm, 3);
3433:     DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
3434:   } break;
3435:   default:
3436:     SETERRQ(PetscObjectComm((PetscObject)rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]);
3437:   }
3438:   {
3439:     PetscInt Nv, v;

3441:     /* Must create the celltype label here so that we do not automatically try to compute the types */
3442:     DMCreateLabel(rdm, "celltype");
3443:     DMPlexSetCellType(rdm, 0, ct);
3444:     DMPlexGetChart(rdm, NULL, &Nv);
3445:     for (v = 1; v < Nv; ++v) DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);
3446:   }
3447:   DMPlexInterpolateInPlace_Internal(rdm);
3448:   PetscObjectSetName((PetscObject)rdm, DMPolytopeTypes[ct]);
3449:   return 0;
3450: }

3452: /*@
3453:   DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

3455:   Collective

3457:   Input Parameters:
3458: + comm - The communicator
3459: - ct   - The cell type of the reference cell

3461:   Output Parameter:
3462: . refdm - The reference cell

3464:   Level: intermediate

3466: .seealso: `DMPlexCreateReferenceCell()`, `DMPlexCreateBoxMesh()`
3467: @*/
3468: PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm)
3469: {
3470:   DMCreate(comm, refdm);
3471:   DMSetType(*refdm, DMPLEX);
3472:   DMPlexCreateReferenceCell_Internal(*refdm, ct);
3473:   return 0;
3474: }

3476: static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[])
3477: {
3478:   DM        plex;
3479:   DMLabel   label;
3480:   PetscBool hasLabel;

3482:   DMHasLabel(dm, name, &hasLabel);
3483:   if (hasLabel) return 0;
3484:   DMCreateLabel(dm, name);
3485:   DMGetLabel(dm, name, &label);
3486:   DMConvert(dm, DMPLEX, &plex);
3487:   DMPlexMarkBoundaryFaces(plex, 1, label);
3488:   DMPlexLabelComplete(plex, label);
3489:   DMDestroy(&plex);
3490:   return 0;
3491: }

3493: /*
3494:   We use the last coordinate as the radius, the inner radius is lower[dim-1] and the outer radius is upper[dim-1]. Then we map the first coordinate around the circle.

3496:     (x, y) -> (r, theta) = (x[1], (x[0] - lower[0]) * 2\pi/(upper[0] - lower[0]))
3497: */
3498: static void boxToAnnulus(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
3499: {
3500:   const PetscReal low = PetscRealPart(constants[0]);
3501:   const PetscReal upp = PetscRealPart(constants[1]);
3502:   const PetscReal r   = PetscRealPart(u[1]);
3503:   const PetscReal th  = 2. * PETSC_PI * (PetscRealPart(u[0]) - low) / (upp - low);

3505:   f0[0] = r * PetscCosReal(th);
3506:   f0[1] = r * PetscSinReal(th);
3507: }

3509: const char *const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "schwarz_p", "gyroid", "doublet", "annulus", "unknown", "DMPlexShape", "DM_SHAPE_", NULL};

3511: static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm)
3512: {
3513:   DMPlexShape    shape   = DM_SHAPE_BOX;
3514:   DMPolytopeType cell    = DM_POLYTOPE_TRIANGLE;
3515:   PetscInt       dim     = 2;
3516:   PetscBool      simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE;
3517:   PetscBool      flg, flg2, fflg, bdfflg, nameflg;
3518:   MPI_Comm       comm;
3519:   char           filename[PETSC_MAX_PATH_LEN]   = "<unspecified>";
3520:   char           bdFilename[PETSC_MAX_PATH_LEN] = "<unspecified>";
3521:   char           plexname[PETSC_MAX_PATH_LEN]   = "";

3523:   PetscObjectGetComm((PetscObject)dm, &comm);
3524:   /* TODO Turn this into a registration interface */
3525:   PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg);
3526:   PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg);
3527:   PetscOptionsString("-dm_plex_name", "Name of the mesh in the file", "DMPlexCreateFromFile", plexname, plexname, sizeof(plexname), &nameflg);
3528:   PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum)cell, (PetscEnum *)&cell, NULL);
3529:   PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL);
3530:   PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum)shape, (PetscEnum *)&shape, &flg);
3531:   PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0);
3533:   PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg);
3534:   PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg);
3535:   PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg);
3536:   PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2);
3537:   if (flg || flg2) DMSetBasicAdjacency(dm, adjCone, adjClosure);

3539:   switch (cell) {
3540:   case DM_POLYTOPE_POINT:
3541:   case DM_POLYTOPE_SEGMENT:
3542:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
3543:   case DM_POLYTOPE_TRIANGLE:
3544:   case DM_POLYTOPE_QUADRILATERAL:
3545:   case DM_POLYTOPE_TETRAHEDRON:
3546:   case DM_POLYTOPE_HEXAHEDRON:
3547:     *useCoordSpace = PETSC_TRUE;
3548:     break;
3549:   default:
3550:     *useCoordSpace = PETSC_FALSE;
3551:     break;
3552:   }

3554:   if (fflg) {
3555:     DM dmnew;

3557:     DMPlexCreateFromFile(PetscObjectComm((PetscObject)dm), filename, plexname, interpolate, &dmnew);
3558:     DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew);
3559:     DMPlexReplace_Internal(dm, &dmnew);
3560:   } else if (refDomain) {
3561:     DMPlexCreateReferenceCell_Internal(dm, cell);
3562:   } else if (bdfflg) {
3563:     DM bdm, dmnew;

3565:     DMPlexCreateFromFile(PetscObjectComm((PetscObject)dm), bdFilename, plexname, interpolate, &bdm);
3566:     PetscObjectSetOptionsPrefix((PetscObject)bdm, "bd_");
3567:     DMSetFromOptions(bdm);
3568:     DMPlexGenerate(bdm, NULL, interpolate, &dmnew);
3569:     DMDestroy(&bdm);
3570:     DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew);
3571:     DMPlexReplace_Internal(dm, &dmnew);
3572:   } else {
3573:     PetscObjectSetName((PetscObject)dm, DMPlexShapes[shape]);
3574:     switch (shape) {
3575:     case DM_SHAPE_BOX:
3576:     case DM_SHAPE_ANNULUS: {
3577:       PetscInt       faces[3]  = {0, 0, 0};
3578:       PetscReal      lower[3]  = {0, 0, 0};
3579:       PetscReal      upper[3]  = {1, 1, 1};
3580:       DMBoundaryType bdt[3]    = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
3581:       PetscBool      isAnnular = shape == DM_SHAPE_ANNULUS ? PETSC_TRUE : PETSC_FALSE;
3582:       PetscInt       i, n;

3584:       n = dim;
3585:       for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4 - dim);
3586:       PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);
3587:       n = 3;
3588:       PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);
3590:       n = 3;
3591:       PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);
3593:       n = 3;
3594:       PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *)bdt, &n, &flg);

3598:       if (isAnnular)
3599:         for (i = 0; i < dim - 1; ++i) bdt[i] = DM_BOUNDARY_PERIODIC;

3601:       switch (cell) {
3602:       case DM_POLYTOPE_TRI_PRISM_TENSOR:
3603:         DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt);
3604:         if (!interpolate) {
3605:           DM udm;

3607:           DMPlexUninterpolate(dm, &udm);
3608:           DMPlexReplace_Internal(dm, &udm);
3609:         }
3610:         break;
3611:       default:
3612:         DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate);
3613:         break;
3614:       }
3615:       if (isAnnular) {
3616:         DM          cdm;
3617:         PetscDS     cds;
3618:         PetscScalar bounds[2] = {lower[0], upper[0]};

3620:         // Fix coordinates for annular region
3621:         DMSetPeriodicity(dm, NULL, NULL, NULL);
3622:         DMSetCellCoordinatesLocal(dm, NULL);
3623:         DMSetCellCoordinates(dm, NULL);
3624:         DMPlexCreateCoordinateSpace(dm, 1, NULL);
3625:         DMGetCoordinateDM(dm, &cdm);
3626:         DMGetDS(cdm, &cds);
3627:         PetscDSSetConstants(cds, 2, bounds);
3628:         DMPlexRemapGeometry(dm, 0.0, boxToAnnulus);
3629:       }
3630:     } break;
3631:     case DM_SHAPE_BOX_SURFACE: {
3632:       PetscInt  faces[3] = {0, 0, 0};
3633:       PetscReal lower[3] = {0, 0, 0};
3634:       PetscReal upper[3] = {1, 1, 1};
3635:       PetscInt  i, n;

3637:       n = dim + 1;
3638:       for (i = 0; i < dim + 1; ++i) faces[i] = (dim + 1 == 1 ? 1 : 4 - (dim + 1));
3639:       PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);
3640:       n = 3;
3641:       PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);
3643:       n = 3;
3644:       PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);
3646:       DMPlexCreateBoxSurfaceMesh_Internal(dm, dim + 1, faces, lower, upper, interpolate);
3647:     } break;
3648:     case DM_SHAPE_SPHERE: {
3649:       PetscReal R = 1.0;

3651:       PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg);
3652:       DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R);
3653:     } break;
3654:     case DM_SHAPE_BALL: {
3655:       PetscReal R = 1.0;

3657:       PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg);
3658:       DMPlexCreateBallMesh_Internal(dm, dim, R);
3659:     } break;
3660:     case DM_SHAPE_CYLINDER: {
3661:       DMBoundaryType bdt = DM_BOUNDARY_NONE;
3662:       PetscInt       Nw  = 6;

3664:       PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum)bdt, (PetscEnum *)&bdt, NULL);
3665:       PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL);
3666:       switch (cell) {
3667:       case DM_POLYTOPE_TRI_PRISM_TENSOR:
3668:         DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate);
3669:         break;
3670:       default:
3671:         DMPlexCreateHexCylinderMesh_Internal(dm, bdt);
3672:         break;
3673:       }
3674:     } break;
3675:     case DM_SHAPE_SCHWARZ_P: // fallthrough
3676:     case DM_SHAPE_GYROID: {
3677:       PetscInt       extent[3] = {1, 1, 1}, refine = 0, layers = 0, three;
3678:       PetscReal      thickness   = 0.;
3679:       DMBoundaryType periodic[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
3680:       DMPlexTPSType  tps_type    = shape == DM_SHAPE_SCHWARZ_P ? DMPLEX_TPS_SCHWARZ_P : DMPLEX_TPS_GYROID;
3681:       PetscBool      tps_distribute;
3682:       PetscOptionsIntArray("-dm_plex_tps_extent", "Number of replicas for each of three dimensions", NULL, extent, (three = 3, &three), NULL);
3683:       PetscOptionsInt("-dm_plex_tps_refine", "Number of refinements", NULL, refine, &refine, NULL);
3684:       PetscOptionsEnumArray("-dm_plex_tps_periodic", "Periodicity in each of three dimensions", NULL, DMBoundaryTypes, (PetscEnum *)periodic, (three = 3, &three), NULL);
3685:       PetscOptionsInt("-dm_plex_tps_layers", "Number of layers in volumetric extrusion (or zero to not extrude)", NULL, layers, &layers, NULL);
3686:       PetscOptionsReal("-dm_plex_tps_thickness", "Thickness of volumetric extrusion", NULL, thickness, &thickness, NULL);
3687:       DMPlexDistributeGetDefault(dm, &tps_distribute);
3688:       PetscOptionsBool("-dm_plex_tps_distribute", "Distribute the 2D mesh prior to refinement and extrusion", NULL, tps_distribute, &tps_distribute, NULL);
3689:       DMPlexCreateTPSMesh_Internal(dm, tps_type, extent, periodic, tps_distribute, refine, layers, thickness);
3690:     } break;
3691:     case DM_SHAPE_DOUBLET: {
3692:       DM        dmnew;
3693:       PetscReal rl = 0.0;

3695:       PetscOptionsReal("-dm_plex_doublet_refinementlimit", "Refinement limit", NULL, rl, &rl, NULL);
3696:       DMPlexCreateDoublet(PetscObjectComm((PetscObject)dm), dim, simplex, interpolate, rl, &dmnew);
3697:       DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_FALSE, dmnew);
3698:       DMPlexReplace_Internal(dm, &dmnew);
3699:     } break;
3700:     default:
3701:       SETERRQ(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]);
3702:     }
3703:   }
3704:   DMPlexSetRefinementUniform(dm, PETSC_TRUE);
3705:   if (!((PetscObject)dm)->name && nameflg) PetscObjectSetName((PetscObject)dm, plexname);
3706:   return 0;
3707: }

3709: PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM dm, PetscOptionItems *PetscOptionsObject)
3710: {
3711:   DM_Plex  *mesh = (DM_Plex *)dm->data;
3712:   PetscBool flg, flg2;
3713:   char      bdLabel[PETSC_MAX_PATH_LEN];

3715:   /* Handle viewing */
3716:   PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);
3717:   PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL, 0);
3718:   PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);
3719:   PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL, 0);
3720:   PetscOptionsBoundedInt("-dm_plex_print_locate", "Debug output level all point location computations", "DMLocatePoints", 0, &mesh->printLocate, NULL, 0);
3721:   DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);
3722:   if (flg) PetscLogDefaultBegin();
3723:   /* Labeling */
3724:   PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg);
3725:   if (flg) DMPlexCreateBoundaryLabel_Private(dm, bdLabel);
3726:   /* Point Location */
3727:   PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);
3728:   /* Partitioning and distribution */
3729:   PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);
3730:   /* Generation and remeshing */
3731:   PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);
3732:   /* Projection behavior */
3733:   PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL, 0);
3734:   PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);
3735:   /* Checking structure */
3736:   {
3737:     PetscBool all = PETSC_FALSE;

3739:     PetscOptionsBool("-dm_plex_check_all", "Perform all basic checks", "DMPlexCheck", PETSC_FALSE, &all, NULL);
3740:     if (all) {
3741:       DMPlexCheck(dm);
3742:     } else {
3743:       PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);
3744:       if (flg && flg2) DMPlexCheckSymmetry(dm);
3745:       PetscOptionsBool("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices (only for homogeneous simplex or tensor meshes)", "DMPlexCheckSkeleton", PETSC_FALSE, &flg, &flg2);
3746:       if (flg && flg2) DMPlexCheckSkeleton(dm, 0);
3747:       PetscOptionsBool("-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", PETSC_FALSE, &flg, &flg2);
3748:       if (flg && flg2) DMPlexCheckFaces(dm, 0);
3749:       PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);
3750:       if (flg && flg2) DMPlexCheckGeometry(dm);
3751:       PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);
3752:       if (flg && flg2) DMPlexCheckPointSF(dm, NULL, PETSC_FALSE);
3753:       PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2);
3754:       if (flg && flg2) DMPlexCheckInterfaceCones(dm);
3755:     }
3756:     PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);
3757:     if (flg && flg2) DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);
3758:   }
3759:   {
3760:     PetscReal scale = 1.0;

3762:     PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg);
3763:     if (flg) {
3764:       Vec coordinates, coordinatesLocal;

3766:       DMGetCoordinates(dm, &coordinates);
3767:       DMGetCoordinatesLocal(dm, &coordinatesLocal);
3768:       VecScale(coordinates, scale);
3769:       VecScale(coordinatesLocal, scale);
3770:     }
3771:   }
3772:   PetscPartitionerSetFromOptions(mesh->partitioner);
3773:   return 0;
3774: }

3776: PetscErrorCode DMSetFromOptions_Overlap_Plex(DM dm, PetscOptionItems *PetscOptionsObject, PetscInt *overlap)
3777: {
3778:   PetscInt  numOvLabels = 16, numOvExLabels = 16;
3779:   char     *ovLabelNames[16], *ovExLabelNames[16];
3780:   PetscInt  numOvValues = 16, numOvExValues = 16, l;
3781:   PetscBool flg;

3783:   PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMPlexDistribute", *overlap, overlap, NULL, 0);
3784:   PetscOptionsStringArray("-dm_distribute_overlap_labels", "List of overlap label names", "DMPlexDistribute", ovLabelNames, &numOvLabels, &flg);
3785:   if (!flg) numOvLabels = 0;
3786:   if (numOvLabels) {
3787:     ((DM_Plex *)dm->data)->numOvLabels = numOvLabels;
3788:     for (l = 0; l < numOvLabels; ++l) {
3789:       DMGetLabel(dm, ovLabelNames[l], &((DM_Plex *)dm->data)->ovLabels[l]);
3791:       PetscFree(ovLabelNames[l]);
3792:     }
3793:     PetscOptionsIntArray("-dm_distribute_overlap_values", "List of overlap label values", "DMPlexDistribute", ((DM_Plex *)dm->data)->ovValues, &numOvValues, &flg);
3794:     if (!flg) numOvValues = 0;

3797:     PetscOptionsStringArray("-dm_distribute_overlap_exclude_labels", "List of overlap exclude label names", "DMPlexDistribute", ovExLabelNames, &numOvExLabels, &flg);
3798:     if (!flg) numOvExLabels = 0;
3799:     ((DM_Plex *)dm->data)->numOvExLabels = numOvExLabels;
3800:     for (l = 0; l < numOvExLabels; ++l) {
3801:       DMGetLabel(dm, ovExLabelNames[l], &((DM_Plex *)dm->data)->ovExLabels[l]);
3803:       PetscFree(ovExLabelNames[l]);
3804:     }
3805:     PetscOptionsIntArray("-dm_distribute_overlap_exclude_values", "List of overlap exclude label values", "DMPlexDistribute", ((DM_Plex *)dm->data)->ovExValues, &numOvExValues, &flg);
3806:     if (!flg) numOvExValues = 0;
3808:   }
3809:   return 0;
3810: }

3812: static PetscErrorCode DMSetFromOptions_Plex(DM dm, PetscOptionItems *PetscOptionsObject)
3813: {
3814:   PetscFunctionList        ordlist;
3815:   char                     oname[256];
3816:   PetscReal                volume    = -1.0;
3817:   PetscInt                 prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim;
3818:   PetscBool                uniformOrig, created = PETSC_FALSE, uniform = PETSC_TRUE, distribute, interpolate = PETSC_TRUE, coordSpace = PETSC_TRUE, remap = PETSC_TRUE, ghostCells = PETSC_FALSE, isHierarchy, ignoreModel = PETSC_FALSE, flg;
3819:   DMPlexReorderDefaultFlag reorder;

3821:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Options");
3822:   /* Handle automatic creation */
3823:   DMGetDimension(dm, &dim);
3824:   if (dim < 0) {
3825:     DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm);
3826:     created = PETSC_TRUE;
3827:   }
3828:   DMGetDimension(dm, &dim);
3829:   /* Handle interpolation before distribution */
3830:   PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg);
3831:   if (flg) {
3832:     DMPlexInterpolatedFlag interpolated;

3834:     DMPlexIsInterpolated(dm, &interpolated);
3835:     if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) {
3836:       DM udm;

3838:       DMPlexUninterpolate(dm, &udm);
3839:       DMPlexReplace_Internal(dm, &udm);
3840:     } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) {
3841:       DM idm;

3843:       DMPlexInterpolate(dm, &idm);
3844:       DMPlexReplace_Internal(dm, &idm);
3845:     }
3846:   }
3847:   /* Handle DMPlex refinement before distribution */
3848:   PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg);
3849:   if (flg) ((DM_Plex *)dm->data)->ignoreModel = ignoreModel;
3850:   DMPlexGetRefinementUniform(dm, &uniformOrig);
3851:   PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL, 0);
3852:   PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);
3853:   PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);
3854:   if (flg) DMPlexSetRefinementUniform(dm, uniform);
3855:   PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);
3856:   if (flg) {
3857:     DMPlexSetRefinementUniform(dm, PETSC_FALSE);
3858:     DMPlexSetRefinementLimit(dm, volume);
3859:     prerefine = PetscMax(prerefine, 1);
3860:   }
3861:   for (r = 0; r < prerefine; ++r) {
3862:     DM             rdm;
3863:     PetscPointFunc coordFunc = ((DM_Plex *)dm->data)->coordFunc;

3865:     DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
3866:     DMRefine(dm, PetscObjectComm((PetscObject)dm), &rdm);
3867:     DMPlexReplace_Internal(dm, &rdm);
3868:     DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
3869:     if (coordFunc && remap) {
3870:       DMPlexRemapGeometry(dm, 0.0, coordFunc);
3871:       ((DM_Plex *)dm->data)->coordFunc = coordFunc;
3872:     }
3873:   }
3874:   DMPlexSetRefinementUniform(dm, uniformOrig);
3875:   /* Handle DMPlex extrusion before distribution */
3876:   PetscOptionsBoundedInt("-dm_extrude", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0);
3877:   if (extLayers) {
3878:     DM edm;

3880:     DMExtrude(dm, extLayers, &edm);
3881:     DMPlexReplace_Internal(dm, &edm);
3882:     ((DM_Plex *)dm->data)->coordFunc = NULL;
3883:     DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
3884:     extLayers = 0;
3885:   }
3886:   /* Handle DMPlex reordering before distribution */
3887:   DMPlexReorderGetDefault(dm, &reorder);
3888:   MatGetOrderingList(&ordlist);
3889:   PetscStrncpy(oname, MATORDERINGNATURAL, sizeof(oname));
3890:   PetscOptionsFList("-dm_plex_reorder", "Set mesh reordering type", "DMPlexGetOrdering", ordlist, MATORDERINGNATURAL, oname, sizeof(oname), &flg);
3891:   if (reorder == DMPLEX_REORDER_DEFAULT_TRUE || flg) {
3892:     DM pdm;
3893:     IS perm;

3895:     DMPlexGetOrdering(dm, oname, NULL, &perm);
3896:     DMPlexPermute(dm, perm, &pdm);
3897:     ISDestroy(&perm);
3898:     DMPlexReplace_Internal(dm, &pdm);
3899:     DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
3900:   }
3901:   /* Handle DMPlex distribution */
3902:   DMPlexDistributeGetDefault(dm, &distribute);
3903:   PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMPlexDistribute", distribute, &distribute, NULL);
3904:   DMSetFromOptions_Overlap_Plex(dm, PetscOptionsObject, &overlap);
3905:   if (distribute) {
3906:     DM               pdm = NULL;
3907:     PetscPartitioner part;

3909:     DMPlexGetPartitioner(dm, &part);
3910:     PetscPartitionerSetFromOptions(part);
3911:     DMPlexDistribute(dm, overlap, NULL, &pdm);
3912:     if (pdm) DMPlexReplace_Internal(dm, &pdm);
3913:   }
3914:   /* Create coordinate space */
3915:   if (created) {
3916:     DM_Plex  *mesh   = (DM_Plex *)dm->data;
3917:     PetscInt  degree = 1;
3918:     PetscBool flg;

3920:     PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg);
3921:     PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, &degree, NULL);
3922:     if (coordSpace) DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc);
3923:     if (flg && !coordSpace) {
3924:       DM           cdm;
3925:       PetscDS      cds;
3926:       PetscObject  obj;
3927:       PetscClassId id;

3929:       DMGetCoordinateDM(dm, &cdm);
3930:       DMGetDS(cdm, &cds);
3931:       PetscDSGetDiscretization(cds, 0, &obj);
3932:       PetscObjectGetClassId(obj, &id);
3933:       if (id == PETSCFE_CLASSID) {
3934:         PetscContainer dummy;

3936:         PetscContainerCreate(PETSC_COMM_SELF, &dummy);
3937:         PetscObjectSetName((PetscObject)dummy, "coordinates");
3938:         DMSetField(cdm, 0, NULL, (PetscObject)dummy);
3939:         PetscContainerDestroy(&dummy);
3940:         DMClearDS(cdm);
3941:       }
3942:       mesh->coordFunc = NULL;
3943:     }
3944:     PetscOptionsBool("-dm_sparse_localize", "Localize only necessary cells", "", dm->sparseLocalize, &dm->sparseLocalize, &flg);
3945:     DMLocalizeCoordinates(dm);
3946:   }
3947:   /* Handle DMPlex refinement */
3948:   remap = PETSC_TRUE;
3949:   PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL, 0);
3950:   PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);
3951:   PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy, 0);
3952:   if (refine) DMPlexSetRefinementUniform(dm, PETSC_TRUE);
3953:   if (refine && isHierarchy) {
3954:     DM *dms, coarseDM;

3956:     DMGetCoarseDM(dm, &coarseDM);
3957:     PetscObjectReference((PetscObject)coarseDM);
3958:     PetscMalloc1(refine, &dms);
3959:     DMRefineHierarchy(dm, refine, dms);
3960:     /* Total hack since we do not pass in a pointer */
3961:     DMPlexSwap_Static(dm, dms[refine - 1]);
3962:     if (refine == 1) {
3963:       DMSetCoarseDM(dm, dms[0]);
3964:       DMPlexSetRegularRefinement(dm, PETSC_TRUE);
3965:     } else {
3966:       DMSetCoarseDM(dm, dms[refine - 2]);
3967:       DMPlexSetRegularRefinement(dm, PETSC_TRUE);
3968:       DMSetCoarseDM(dms[0], dms[refine - 1]);
3969:       DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);
3970:     }
3971:     DMSetCoarseDM(dms[refine - 1], coarseDM);
3972:     PetscObjectDereference((PetscObject)coarseDM);
3973:     /* Free DMs */
3974:     for (r = 0; r < refine; ++r) {
3975:       DMSetFromOptions_NonRefinement_Plex(dms[r], PetscOptionsObject);
3976:       DMDestroy(&dms[r]);
3977:     }
3978:     PetscFree(dms);
3979:   } else {
3980:     for (r = 0; r < refine; ++r) {
3981:       DM             rdm;
3982:       PetscPointFunc coordFunc = ((DM_Plex *)dm->data)->coordFunc;

3984:       DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
3985:       DMRefine(dm, PetscObjectComm((PetscObject)dm), &rdm);
3986:       /* Total hack since we do not pass in a pointer */
3987:       DMPlexReplace_Internal(dm, &rdm);
3988:       DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
3989:       if (coordFunc && remap) {
3990:         DMPlexRemapGeometry(dm, 0.0, coordFunc);
3991:         ((DM_Plex *)dm->data)->coordFunc = coordFunc;
3992:       }
3993:     }
3994:   }
3995:   /* Handle DMPlex coarsening */
3996:   PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL, 0);
3997:   PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy, 0);
3998:   if (coarsen && isHierarchy) {
3999:     DM *dms;

4001:     PetscMalloc1(coarsen, &dms);
4002:     DMCoarsenHierarchy(dm, coarsen, dms);
4003:     /* Free DMs */
4004:     for (r = 0; r < coarsen; ++r) {
4005:       DMSetFromOptions_NonRefinement_Plex(dms[r], PetscOptionsObject);
4006:       DMDestroy(&dms[r]);
4007:     }
4008:     PetscFree(dms);
4009:   } else {
4010:     for (r = 0; r < coarsen; ++r) {
4011:       DM             cdm;
4012:       PetscPointFunc coordFunc = ((DM_Plex *)dm->data)->coordFunc;

4014:       DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
4015:       DMCoarsen(dm, PetscObjectComm((PetscObject)dm), &cdm);
4016:       /* Total hack since we do not pass in a pointer */
4017:       DMPlexReplace_Internal(dm, &cdm);
4018:       DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
4019:       if (coordFunc) {
4020:         DMPlexRemapGeometry(dm, 0.0, coordFunc);
4021:         ((DM_Plex *)dm->data)->coordFunc = coordFunc;
4022:       }
4023:     }
4024:   }
4025:   /* Handle ghost cells */
4026:   PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL);
4027:   if (ghostCells) {
4028:     DM   gdm;
4029:     char lname[PETSC_MAX_PATH_LEN];

4031:     lname[0] = '\0';
4032:     PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg);
4033:     DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm);
4034:     DMPlexReplace_Internal(dm, &gdm);
4035:   }
4036:   /* Handle 1D order */
4037:   if (reorder != DMPLEX_REORDER_DEFAULT_FALSE && dim == 1) {
4038:     DM           cdm, rdm;
4039:     PetscDS      cds;
4040:     PetscObject  obj;
4041:     PetscClassId id = PETSC_OBJECT_CLASSID;
4042:     IS           perm;
4043:     PetscInt     Nf;
4044:     PetscBool    distributed;

4046:     DMPlexIsDistributed(dm, &distributed);
4047:     DMGetCoordinateDM(dm, &cdm);
4048:     DMGetDS(cdm, &cds);
4049:     PetscDSGetNumFields(cds, &Nf);
4050:     if (Nf) {
4051:       PetscDSGetDiscretization(cds, 0, &obj);
4052:       PetscObjectGetClassId(obj, &id);
4053:     }
4054:     if (!distributed && id != PETSCFE_CLASSID) {
4055:       DMPlexGetOrdering1D(dm, &perm);
4056:       DMPlexPermute(dm, perm, &rdm);
4057:       DMPlexReplace_Internal(dm, &rdm);
4058:       ISDestroy(&perm);
4059:     }
4060:   }
4061:   /* Handle */
4062:   DMSetFromOptions_NonRefinement_Plex(dm, PetscOptionsObject);
4063:   PetscOptionsHeadEnd();
4064:   return 0;
4065: }

4067: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm, Vec *vec)
4068: {
4069:   DMCreateGlobalVector_Section_Private(dm, vec);
4070:   /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
4071:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);
4072:   VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void))VecView_Plex_Native);
4073:   VecSetOperation(*vec, VECOP_LOAD, (void (*)(void))VecLoad_Plex);
4074:   VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void))VecLoad_Plex_Native);
4075:   return 0;
4076: }

4078: static PetscErrorCode DMCreateLocalVector_Plex(DM dm, Vec *vec)
4079: {
4080:   DMCreateLocalVector_Section_Private(dm, vec);
4081:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex_Local);
4082:   VecSetOperation(*vec, VECOP_LOAD, (void (*)(void))VecLoad_Plex_Local);
4083:   return 0;
4084: }

4086: static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4087: {
4088:   PetscInt depth, d;

4090:   DMPlexGetDepth(dm, &depth);
4091:   if (depth == 1) {
4092:     DMGetDimension(dm, &d);
4093:     if (dim == 0) DMPlexGetDepthStratum(dm, dim, pStart, pEnd);
4094:     else if (dim == d) DMPlexGetDepthStratum(dm, 1, pStart, pEnd);
4095:     else {
4096:       *pStart = 0;
4097:       *pEnd   = 0;
4098:     }
4099:   } else {
4100:     DMPlexGetDepthStratum(dm, dim, pStart, pEnd);
4101:   }
4102:   return 0;
4103: }

4105: static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[])
4106: {
4107:   PetscSF            sf;
4108:   PetscInt           niranks, njranks, n;
4109:   const PetscMPIInt *iranks, *jranks;
4110:   DM_Plex           *data = (DM_Plex *)dm->data;

4112:   DMGetPointSF(dm, &sf);
4113:   if (!data->neighbors) {
4114:     PetscSFSetUp(sf);
4115:     PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);
4116:     PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);
4117:     PetscMalloc1(njranks + niranks + 1, &data->neighbors);
4118:     PetscArraycpy(data->neighbors + 1, jranks, njranks);
4119:     PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);
4120:     n = njranks + niranks;
4121:     PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);
4122:     /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */
4123:     PetscMPIIntCast(n, data->neighbors);
4124:   }
4125:   if (nranks) *nranks = data->neighbors[0];
4126:   if (ranks) {
4127:     if (data->neighbors[0]) *ranks = data->neighbors + 1;
4128:     else *ranks = NULL;
4129:   }
4130:   return 0;
4131: }

4133: PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec);

4135: static PetscErrorCode DMInitialize_Plex(DM dm)
4136: {
4137:   dm->ops->view                      = DMView_Plex;
4138:   dm->ops->load                      = DMLoad_Plex;
4139:   dm->ops->setfromoptions            = DMSetFromOptions_Plex;
4140:   dm->ops->clone                     = DMClone_Plex;
4141:   dm->ops->setup                     = DMSetUp_Plex;
4142:   dm->ops->createlocalsection        = DMCreateLocalSection_Plex;
4143:   dm->ops->createdefaultconstraints  = DMCreateDefaultConstraints_Plex;
4144:   dm->ops->createglobalvector        = DMCreateGlobalVector_Plex;
4145:   dm->ops->createlocalvector         = DMCreateLocalVector_Plex;
4146:   dm->ops->getlocaltoglobalmapping   = NULL;
4147:   dm->ops->createfieldis             = NULL;
4148:   dm->ops->createcoordinatedm        = DMCreateCoordinateDM_Plex;
4149:   dm->ops->createcoordinatefield     = DMCreateCoordinateField_Plex;
4150:   dm->ops->getcoloring               = NULL;
4151:   dm->ops->creatematrix              = DMCreateMatrix_Plex;
4152:   dm->ops->createinterpolation       = DMCreateInterpolation_Plex;
4153:   dm->ops->createmassmatrix          = DMCreateMassMatrix_Plex;
4154:   dm->ops->createmassmatrixlumped    = DMCreateMassMatrixLumped_Plex;
4155:   dm->ops->createinjection           = DMCreateInjection_Plex;
4156:   dm->ops->refine                    = DMRefine_Plex;
4157:   dm->ops->coarsen                   = DMCoarsen_Plex;
4158:   dm->ops->refinehierarchy           = DMRefineHierarchy_Plex;
4159:   dm->ops->coarsenhierarchy          = DMCoarsenHierarchy_Plex;
4160:   dm->ops->extrude                   = DMExtrude_Plex;
4161:   dm->ops->globaltolocalbegin        = NULL;
4162:   dm->ops->globaltolocalend          = NULL;
4163:   dm->ops->localtoglobalbegin        = NULL;
4164:   dm->ops->localtoglobalend          = NULL;
4165:   dm->ops->destroy                   = DMDestroy_Plex;
4166:   dm->ops->createsubdm               = DMCreateSubDM_Plex;
4167:   dm->ops->createsuperdm             = DMCreateSuperDM_Plex;
4168:   dm->ops->getdimpoints              = DMGetDimPoints_Plex;
4169:   dm->ops->locatepoints              = DMLocatePoints_Plex;
4170:   dm->ops->projectfunctionlocal      = DMProjectFunctionLocal_Plex;
4171:   dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex;
4172:   dm->ops->projectfieldlocal         = DMProjectFieldLocal_Plex;
4173:   dm->ops->projectfieldlabellocal    = DMProjectFieldLabelLocal_Plex;
4174:   dm->ops->projectbdfieldlabellocal  = DMProjectBdFieldLabelLocal_Plex;
4175:   dm->ops->computel2diff             = DMComputeL2Diff_Plex;
4176:   dm->ops->computel2gradientdiff     = DMComputeL2GradientDiff_Plex;
4177:   dm->ops->computel2fielddiff        = DMComputeL2FieldDiff_Plex;
4178:   dm->ops->getneighbors              = DMGetNeighbors_Plex;
4179:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_Plex);
4180:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertTimeDerviativeBoundaryValues_C", DMPlexInsertTimeDerivativeBoundaryValues_Plex);
4181:   PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Plex);
4182:   PetscObjectComposeFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", DMCreateNeumannOverlap_Plex);
4183:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", DMPlexGetOverlap_Plex);
4184:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexDistributeGetDefault_C", DMPlexDistributeGetDefault_Plex);
4185:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexDistributeSetDefault_C", DMPlexDistributeSetDefault_Plex);
4186:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexReorderGetDefault_C", DMPlexReorderGetDefault_Plex);
4187:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexReorderSetDefault_C", DMPlexReorderSetDefault_Plex);
4188:   PetscObjectComposeFunction((PetscObject)dm, "DMInterpolateSolution_C", DMInterpolateSolution_Plex);
4189:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", DMPlexGetOverlap_Plex);
4190:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexSetOverlap_C", DMPlexSetOverlap_Plex);
4191:   return 0;
4192: }

4194: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
4195: {
4196:   DM_Plex *mesh = (DM_Plex *)dm->data;

4198:   mesh->refct++;
4199:   (*newdm)->data = mesh;
4200:   PetscObjectChangeTypeName((PetscObject)*newdm, DMPLEX);
4201:   DMInitialize_Plex(*newdm);
4202:   return 0;
4203: }

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

4211:   Options Database Keys:
4212: + -dm_refine_pre                     - Refine mesh before distribution
4213: + -dm_refine_uniform_pre             - Choose uniform or generator-based refinement
4214: + -dm_refine_volume_limit_pre        - Cell volume limit after pre-refinement using generator
4215: . -dm_distribute                     - Distribute mesh across processes
4216: . -dm_distribute_overlap             - Number of cells to overlap for distribution
4217: . -dm_refine                         - Refine mesh after distribution
4218: . -dm_plex_hash_location             - Use grid hashing for point location
4219: . -dm_plex_hash_box_faces <n,m,p>    - The number of divisions in each direction of the grid hash
4220: . -dm_plex_partition_balance         - Attempt to evenly divide points on partition boundary between processes
4221: . -dm_plex_remesh_bd                 - Allow changes to the boundary on remeshing
4222: . -dm_plex_max_projection_height     - Maxmimum mesh point height used to project locally
4223: . -dm_plex_regular_refinement        - Use special nested projection algorithm for regular refinement
4224: . -dm_plex_check_all                 - Perform all shecks below
4225: . -dm_plex_check_symmetry            - Check that the adjacency information in the mesh is symmetric
4226: . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices
4227: . -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
4228: . -dm_plex_check_geometry            - Check that cells have positive volume
4229: . -dm_view :mesh.tex:ascii_latex     - View the mesh in LaTeX/TikZ
4230: . -dm_plex_view_scale <num>          - Scale the TikZ
4231: - -dm_plex_print_fem <num>           - View FEM assembly information, such as element vectors and matrices

4233:   Level: intermediate

4235: .seealso: `DMType`, `DMPlexCreate()`, `DMCreate()`, `DMSetType()`
4236: M*/

4238: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
4239: {
4240:   DM_Plex *mesh;
4241:   PetscInt unit;

4244:   PetscNew(&mesh);
4245:   dm->data = mesh;

4247:   mesh->refct = 1;
4248:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
4249:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
4250:   mesh->refinementUniform      = PETSC_TRUE;
4251:   mesh->refinementLimit        = -1.0;
4252:   mesh->distDefault            = PETSC_TRUE;
4253:   mesh->reorderDefault         = DMPLEX_REORDER_DEFAULT_NOTSET;
4254:   mesh->distributionName       = NULL;
4255:   mesh->interpolated           = DMPLEX_INTERPOLATED_INVALID;
4256:   mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID;

4258:   PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);
4259:   mesh->remeshBd = PETSC_FALSE;

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

4263:   mesh->depthState    = -1;
4264:   mesh->celltypeState = -1;
4265:   mesh->printTol      = 1.0e-10;

4267:   DMInitialize_Plex(dm);
4268:   return 0;
4269: }

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

4274:   Collective

4276:   Input Parameter:
4277: . comm - The communicator for the DMPlex object

4279:   Output Parameter:
4280: . mesh  - The DMPlex object

4282:   Level: beginner

4284: @*/
4285: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
4286: {
4288:   DMCreate(comm, mesh);
4289:   DMSetType(*mesh, DMPLEX);
4290:   return 0;
4291: }

4293: /*@C
4294:   DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output)

4296:   Input Parameters:
4297: + dm - The DM
4298: . numCells - The number of cells owned by this process
4299: . numVertices - The number of vertices to be owned by this process, or PETSC_DECIDE
4300: . NVertices - The global number of vertices, or PETSC_DETERMINE
4301: . numCorners - The number of vertices for each cell
4302: - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell

4304:   Output Parameters:
4305: + vertexSF - (Optional) SF describing complete vertex ownership
4306: - verticesAdjSaved - (Optional) vertex adjacency array

4308:   Notes:
4309:   Two triangles sharing a face
4310: $
4311: $        2
4312: $      / | \
4313: $     /  |  \
4314: $    /   |   \
4315: $   0  0 | 1  3
4316: $    \   |   /
4317: $     \  |  /
4318: $      \ | /
4319: $        1
4320: would have input
4321: $  numCells = 2, numVertices = 4
4322: $  cells = [0 1 2  1 3 2]
4323: $
4324: which would result in the DMPlex
4325: $
4326: $        4
4327: $      / | \
4328: $     /  |  \
4329: $    /   |   \
4330: $   2  0 | 1  5
4331: $    \   |   /
4332: $     \  |  /
4333: $      \ | /
4334: $        3

4336:   Vertices are implicitly numbered consecutively 0,...,NVertices.
4337:   Each rank owns a chunk of numVertices consecutive vertices.
4338:   If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout.
4339:   If NVertices is PETSC_DETERMINE and numVertices is PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1.
4340:   If only NVertices is PETSC_DETERMINE, it is computed as the sum of numVertices over all ranks.

4342:   The cell distribution is arbitrary non-overlapping, independent of the vertex distribution.

4344:   Not currently supported in Fortran.

4346:   Level: advanced

4348: .seealso: `DMPlexBuildFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildCoordinatesFromCellListParallel()`
4349: @*/
4350: PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF, PetscInt **verticesAdjSaved)
4351: {
4352:   PetscSF     sfPoint;
4353:   PetscLayout layout;
4354:   PetscInt    numVerticesAdj, *verticesAdj, *cones, c, p;

4357:   PetscLogEventBegin(DMPLEX_BuildFromCellList, dm, 0, 0, 0);
4358:   /* Get/check global number of vertices */
4359:   {
4360:     PetscInt       NVerticesInCells, i;
4361:     const PetscInt len = numCells * numCorners;

4363:     /* NVerticesInCells = max(cells) + 1 */
4364:     NVerticesInCells = PETSC_MIN_INT;
4365:     for (i = 0; i < len; i++)
4366:       if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i];
4367:     ++NVerticesInCells;
4368:     MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));

4370:     if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells;
4371:     else
4373:   }
4374:   /* Count locally unique vertices */
4375:   {
4376:     PetscHSetI vhash;
4377:     PetscInt   off = 0;

4379:     PetscHSetICreate(&vhash);
4380:     for (c = 0; c < numCells; ++c) {
4381:       for (p = 0; p < numCorners; ++p) PetscHSetIAdd(vhash, cells[c * numCorners + p]);
4382:     }
4383:     PetscHSetIGetSize(vhash, &numVerticesAdj);
4384:     if (!verticesAdjSaved) PetscMalloc1(numVerticesAdj, &verticesAdj);
4385:     else verticesAdj = *verticesAdjSaved;
4386:     PetscHSetIGetElems(vhash, &off, verticesAdj);
4387:     PetscHSetIDestroy(&vhash);
4389:   }
4390:   PetscSortInt(numVerticesAdj, verticesAdj);
4391:   /* Create cones */
4392:   DMPlexSetChart(dm, 0, numCells + numVerticesAdj);
4393:   for (c = 0; c < numCells; ++c) DMPlexSetConeSize(dm, c, numCorners);
4394:   DMSetUp(dm);
4395:   DMPlexGetCones(dm, &cones);
4396:   for (c = 0; c < numCells; ++c) {
4397:     for (p = 0; p < numCorners; ++p) {
4398:       const PetscInt gv = cells[c * numCorners + p];
4399:       PetscInt       lv;

4401:       /* Positions within verticesAdj form 0-based local vertex numbering;
4402:          we need to shift it by numCells to get correct DAG points (cells go first) */
4403:       PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);
4405:       cones[c * numCorners + p] = lv + numCells;
4406:     }
4407:   }
4408:   /* Build point sf */
4409:   PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout);
4410:   PetscLayoutSetSize(layout, NVertices);
4411:   PetscLayoutSetLocalSize(layout, numVertices);
4412:   PetscLayoutSetBlockSize(layout, 1);
4413:   PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint);
4414:   PetscLayoutDestroy(&layout);
4415:   if (!verticesAdjSaved) PetscFree(verticesAdj);
4416:   PetscObjectSetName((PetscObject)sfPoint, "point SF");
4417:   if (dm->sf) {
4418:     const char *prefix;

4420:     PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix);
4421:     PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix);
4422:   }
4423:   DMSetPointSF(dm, sfPoint);
4424:   PetscSFDestroy(&sfPoint);
4425:   if (vertexSF) PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF");
4426:   /* Fill in the rest of the topology structure */
4427:   DMPlexSymmetrize(dm);
4428:   DMPlexStratify(dm);
4429:   PetscLogEventEnd(DMPLEX_BuildFromCellList, dm, 0, 0, 0);
4430:   return 0;
4431: }

4433: /*@C
4434:   DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output)

4436:   Input Parameters:
4437: + dm - The DM
4438: . spaceDim - The spatial dimension used for coordinates
4439: . sfVert - SF describing complete vertex ownership
4440: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex

4442:   Level: advanced

4444:   Notes:
4445:   Not currently supported in Fortran.

4447: .seealso: `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellListParallel()`
4448: @*/
4449: PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[])
4450: {
4451:   PetscSection coordSection;
4452:   Vec          coordinates;
4453:   PetscScalar *coords;
4454:   PetscInt     numVertices, numVerticesAdj, coordSize, v, vStart, vEnd;

4456:   PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList, dm, 0, 0, 0);
4457:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
4459:   DMSetCoordinateDim(dm, spaceDim);
4460:   PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);
4462:   DMGetCoordinateSection(dm, &coordSection);
4463:   PetscSectionSetNumFields(coordSection, 1);
4464:   PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
4465:   PetscSectionSetChart(coordSection, vStart, vEnd);
4466:   for (v = vStart; v < vEnd; ++v) {
4467:     PetscSectionSetDof(coordSection, v, spaceDim);
4468:     PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
4469:   }
4470:   PetscSectionSetUp(coordSection);
4471:   PetscSectionGetStorageSize(coordSection, &coordSize);
4472:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
4473:   VecSetBlockSize(coordinates, spaceDim);
4474:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
4475:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
4476:   VecSetType(coordinates, VECSTANDARD);
4477:   VecGetArray(coordinates, &coords);
4478:   {
4479:     MPI_Datatype coordtype;

4481:     /* Need a temp buffer for coords if we have complex/single */
4482:     MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);
4483:     MPI_Type_commit(&coordtype);
4484: #if defined(PETSC_USE_COMPLEX)
4485:     {
4486:       PetscScalar *svertexCoords;
4487:       PetscInt     i;
4488:       PetscMalloc1(numVertices * spaceDim, &svertexCoords);
4489:       for (i = 0; i < numVertices * spaceDim; i++) svertexCoords[i] = vertexCoords[i];
4490:       PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords, MPI_REPLACE);
4491:       PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords, MPI_REPLACE);
4492:       PetscFree(svertexCoords);
4493:     }
4494: #else
4495:     PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords, MPI_REPLACE);
4496:     PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords, MPI_REPLACE);
4497: #endif
4498:     MPI_Type_free(&coordtype);
4499:   }
4500:   VecRestoreArray(coordinates, &coords);
4501:   DMSetCoordinatesLocal(dm, coordinates);
4502:   VecDestroy(&coordinates);
4503:   PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList, dm, 0, 0, 0);
4504:   return 0;
4505: }

4507: /*@
4508:   DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output)

4510:   Input Parameters:
4511: + comm - The communicator
4512: . dim - The topological dimension of the mesh
4513: . numCells - The number of cells owned by this process
4514: . numVertices - The number of vertices owned by this process, or PETSC_DECIDE
4515: . NVertices - The global number of vertices, or PETSC_DECIDE
4516: . numCorners - The number of vertices for each cell
4517: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
4518: . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell
4519: . spaceDim - The spatial dimension used for coordinates
4520: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex

4522:   Output Parameters:
4523: + dm - The DM
4524: . vertexSF - (Optional) SF describing complete vertex ownership
4525: - verticesAdjSaved - (Optional) vertex adjacency array

4527:   Notes:
4528:   This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(),
4529:   DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel()

4531:   See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters.
4532:   See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters.

4534:   Level: intermediate

4536: .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()`
4537: @*/
4538: PetscErrorCode DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, PetscInt **verticesAdj, DM *dm)
4539: {
4540:   PetscSF sfVert;

4542:   DMCreate(comm, dm);
4543:   DMSetType(*dm, DMPLEX);
4546:   DMSetDimension(*dm, dim);
4547:   DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert, verticesAdj);
4548:   if (interpolate) {
4549:     DM idm;

4551:     DMPlexInterpolate(*dm, &idm);
4552:     DMDestroy(dm);
4553:     *dm = idm;
4554:   }
4555:   DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);
4556:   if (vertexSF) *vertexSF = sfVert;
4557:   else PetscSFDestroy(&sfVert);
4558:   return 0;
4559: }

4561: /*@C
4562:   DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output)

4564:   Input Parameters:
4565: + dm - The DM
4566: . numCells - The number of cells owned by this process
4567: . numVertices - The number of vertices owned by this process, or PETSC_DETERMINE
4568: . numCorners - The number of vertices for each cell
4569: - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell

4571:   Level: advanced

4573:   Notes:
4574:   Two triangles sharing a face
4575: $
4576: $        2
4577: $      / | \
4578: $     /  |  \
4579: $    /   |   \
4580: $   0  0 | 1  3
4581: $    \   |   /
4582: $     \  |  /
4583: $      \ | /
4584: $        1
4585: would have input
4586: $  numCells = 2, numVertices = 4
4587: $  cells = [0 1 2  1 3 2]
4588: $
4589: which would result in the DMPlex
4590: $
4591: $        4
4592: $      / | \
4593: $     /  |  \
4594: $    /   |   \
4595: $   2  0 | 1  5
4596: $    \   |   /
4597: $     \  |  /
4598: $      \ | /
4599: $        3

4601:   If numVertices is PETSC_DETERMINE, it is computed by PETSc as the maximum vertex index in cells + 1.

4603:   Not currently supported in Fortran.

4605: .seealso: `DMPlexBuildFromCellListParallel()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromCellListPetsc()`
4606: @*/
4607: PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[])
4608: {
4609:   PetscInt *cones, c, p, dim;

4611:   PetscLogEventBegin(DMPLEX_BuildFromCellList, dm, 0, 0, 0);
4612:   DMGetDimension(dm, &dim);
4613:   /* Get/check global number of vertices */
4614:   {
4615:     PetscInt       NVerticesInCells, i;
4616:     const PetscInt len = numCells * numCorners;

4618:     /* NVerticesInCells = max(cells) + 1 */
4619:     NVerticesInCells = PETSC_MIN_INT;
4620:     for (i = 0; i < len; i++)
4621:       if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i];
4622:     ++NVerticesInCells;

4624:     if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells;
4625:     else
4627:   }
4628:   DMPlexSetChart(dm, 0, numCells + numVertices);
4629:   for (c = 0; c < numCells; ++c) DMPlexSetConeSize(dm, c, numCorners);
4630:   DMSetUp(dm);
4631:   DMPlexGetCones(dm, &cones);
4632:   for (c = 0; c < numCells; ++c) {
4633:     for (p = 0; p < numCorners; ++p) cones[c * numCorners + p] = cells[c * numCorners + p] + numCells;
4634:   }
4635:   DMPlexSymmetrize(dm);
4636:   DMPlexStratify(dm);
4637:   PetscLogEventEnd(DMPLEX_BuildFromCellList, dm, 0, 0, 0);
4638:   return 0;
4639: }

4641: /*@C
4642:   DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output)

4644:   Input Parameters:
4645: + dm - The DM
4646: . spaceDim - The spatial dimension used for coordinates
4647: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex

4649:   Level: advanced

4651:   Notes:
4652:   Not currently supported in Fortran.

4654: .seealso: `DMPlexBuildCoordinatesFromCellListParallel()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexBuildFromCellList()`
4655: @*/
4656: PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[])
4657: {
4658:   PetscSection coordSection;
4659:   Vec          coordinates;
4660:   DM           cdm;
4661:   PetscScalar *coords;
4662:   PetscInt     v, vStart, vEnd, d;

4664:   PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList, dm, 0, 0, 0);
4665:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
4667:   DMSetCoordinateDim(dm, spaceDim);
4668:   DMGetCoordinateSection(dm, &coordSection);
4669:   PetscSectionSetNumFields(coordSection, 1);
4670:   PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
4671:   PetscSectionSetChart(coordSection, vStart, vEnd);
4672:   for (v = vStart; v < vEnd; ++v) {
4673:     PetscSectionSetDof(coordSection, v, spaceDim);
4674:     PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
4675:   }
4676:   PetscSectionSetUp(coordSection);

4678:   DMGetCoordinateDM(dm, &cdm);
4679:   DMCreateLocalVector(cdm, &coordinates);
4680:   VecSetBlockSize(coordinates, spaceDim);
4681:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
4682:   VecGetArrayWrite(coordinates, &coords);
4683:   for (v = 0; v < vEnd - vStart; ++v) {
4684:     for (d = 0; d < spaceDim; ++d) coords[v * spaceDim + d] = vertexCoords[v * spaceDim + d];
4685:   }
4686:   VecRestoreArrayWrite(coordinates, &coords);
4687:   DMSetCoordinatesLocal(dm, coordinates);
4688:   VecDestroy(&coordinates);
4689:   PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList, dm, 0, 0, 0);
4690:   return 0;
4691: }

4693: /*@
4694:   DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input

4696:   Collective on comm

4698:   Input Parameters:
4699: + comm - The communicator
4700: . dim - The topological dimension of the mesh
4701: . numCells - The number of cells, only on process 0
4702: . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0
4703: . numCorners - The number of vertices for each cell, only on process 0
4704: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
4705: . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0
4706: . spaceDim - The spatial dimension used for coordinates
4707: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0

4709:   Output Parameter:
4710: . dm - The DM, which only has points on process 0

4712:   Notes:
4713:   This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(),
4714:   DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList()

4716:   See DMPlexBuildFromCellList() for an example and details about the topology-related parameters.
4717:   See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters.
4718:   See DMPlexCreateFromCellListParallelPetsc() for parallel input

4720:   Level: intermediate

4722: .seealso: `DMPlexCreateFromCellListParallelPetsc()`, `DMPlexBuildFromCellList()`, `DMPlexBuildCoordinatesFromCellList()`, `DMPlexCreateFromDAG()`, `DMPlexCreate()`
4723: @*/
4724: PetscErrorCode DMPlexCreateFromCellListPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], DM *dm)
4725: {
4726:   PetscMPIInt rank;

4729:   MPI_Comm_rank(comm, &rank);
4730:   DMCreate(comm, dm);
4731:   DMSetType(*dm, DMPLEX);
4732:   DMSetDimension(*dm, dim);
4733:   if (rank == 0) DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);
4734:   else DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL);
4735:   if (interpolate) {
4736:     DM idm;

4738:     DMPlexInterpolate(*dm, &idm);
4739:     DMDestroy(dm);
4740:     *dm = idm;
4741:   }
4742:   if (rank == 0) DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);
4743:   else DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL);
4744:   return 0;
4745: }

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

4750:   Input Parameters:
4751: + dm - The empty DM object, usually from DMCreate() and DMSetDimension()
4752: . depth - The depth of the DAG
4753: . numPoints - Array of size depth + 1 containing the number of points at each depth
4754: . coneSize - The cone size of each point
4755: . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
4756: . coneOrientations - The orientation of each cone point
4757: - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim()

4759:   Output Parameter:
4760: . dm - The DM

4762:   Note: Two triangles sharing a face would have input
4763: $  depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
4764: $  cones = [2 3 4  3 5 4], coneOrientations = [0 0 0  0 0 0]
4765: $ vertexCoords = [-1.0 0.0  0.0 -1.0  0.0 1.0  1.0 0.0]
4766: $
4767: which would result in the DMPlex
4768: $
4769: $        4
4770: $      / | \
4771: $     /  |  \
4772: $    /   |   \
4773: $   2  0 | 1  5
4774: $    \   |   /
4775: $     \  |  /
4776: $      \ | /
4777: $        3
4778: $
4779: $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc()

4781:   Level: advanced

4783: .seealso: `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()`
4784: @*/
4785: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
4786: {
4787:   Vec          coordinates;
4788:   PetscSection coordSection;
4789:   PetscScalar *coords;
4790:   PetscInt     coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off;

4792:   DMGetDimension(dm, &dim);
4793:   DMGetCoordinateDim(dm, &dimEmbed);
4795:   for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
4796:   DMPlexSetChart(dm, pStart, pEnd);
4797:   for (p = pStart; p < pEnd; ++p) {
4798:     DMPlexSetConeSize(dm, p, coneSize[p - pStart]);
4799:     if (firstVertex < 0 && !coneSize[p - pStart]) firstVertex = p - pStart;
4800:   }
4802:   DMSetUp(dm); /* Allocate space for cones */
4803:   for (p = pStart, off = 0; p < pEnd; off += coneSize[p - pStart], ++p) {
4804:     DMPlexSetCone(dm, p, &cones[off]);
4805:     DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
4806:   }
4807:   DMPlexSymmetrize(dm);
4808:   DMPlexStratify(dm);
4809:   /* Build coordinates */
4810:   DMGetCoordinateSection(dm, &coordSection);
4811:   PetscSectionSetNumFields(coordSection, 1);
4812:   PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
4813:   PetscSectionSetChart(coordSection, firstVertex, firstVertex + numPoints[0]);
4814:   for (v = firstVertex; v < firstVertex + numPoints[0]; ++v) {
4815:     PetscSectionSetDof(coordSection, v, dimEmbed);
4816:     PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
4817:   }
4818:   PetscSectionSetUp(coordSection);
4819:   PetscSectionGetStorageSize(coordSection, &coordSize);
4820:   VecCreate(PETSC_COMM_SELF, &coordinates);
4821:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
4822:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
4823:   VecSetBlockSize(coordinates, dimEmbed);
4824:   VecSetType(coordinates, VECSTANDARD);
4825:   if (vertexCoords) {
4826:     VecGetArray(coordinates, &coords);
4827:     for (v = 0; v < numPoints[0]; ++v) {
4828:       PetscInt off;

4830:       PetscSectionGetOffset(coordSection, v + firstVertex, &off);
4831:       for (d = 0; d < dimEmbed; ++d) coords[off + d] = vertexCoords[v * dimEmbed + d];
4832:     }
4833:   }
4834:   VecRestoreArray(coordinates, &coords);
4835:   DMSetCoordinatesLocal(dm, coordinates);
4836:   VecDestroy(&coordinates);
4837:   return 0;
4838: }

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

4843: + comm        - The MPI communicator
4844: . filename    - Name of the .dat file
4845: - interpolate - Create faces and edges in the mesh

4847:   Output Parameter:
4848: . dm  - The DM object representing the mesh

4850:   Note: The format is the simplest possible:
4851: $ Ne
4852: $ v0 v1 ... vk
4853: $ Nv
4854: $ x y z marker

4856:   Level: beginner

4858: .seealso: `DMPlexCreateFromFile()`, `DMPlexCreateMedFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
4859: @*/
4860: PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
4861: {
4862:   DMLabel      marker;
4863:   PetscViewer  viewer;
4864:   Vec          coordinates;
4865:   PetscSection coordSection;
4866:   PetscScalar *coords;
4867:   char         line[PETSC_MAX_PATH_LEN];
4868:   PetscInt     dim = 3, cdim = 3, coordSize, v, c, d;
4869:   PetscMPIInt  rank;
4870:   int          snum, Nv, Nc, Ncn, Nl;

4872:   MPI_Comm_rank(comm, &rank);
4873:   PetscViewerCreate(comm, &viewer);
4874:   PetscViewerSetType(viewer, PETSCVIEWERASCII);
4875:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
4876:   PetscViewerFileSetName(viewer, filename);
4877:   if (rank == 0) {
4878:     PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);
4879:     snum = sscanf(line, "%d %d %d %d", &Nc, &Nv, &Ncn, &Nl);
4881:   } else {
4882:     Nc = Nv = Ncn = Nl = 0;
4883:   }
4884:   DMCreate(comm, dm);
4885:   DMSetType(*dm, DMPLEX);
4886:   DMPlexSetChart(*dm, 0, Nc + Nv);
4887:   DMSetDimension(*dm, dim);
4888:   DMSetCoordinateDim(*dm, cdim);
4889:   /* Read topology */
4890:   if (rank == 0) {
4891:     char     format[PETSC_MAX_PATH_LEN];
4892:     PetscInt cone[8];
4893:     int      vbuf[8], v;

4895:     for (c = 0; c < Ncn; ++c) {
4896:       format[c * 3 + 0] = '%';
4897:       format[c * 3 + 1] = 'd';
4898:       format[c * 3 + 2] = ' ';
4899:     }
4900:     format[Ncn * 3 - 1] = '\0';
4901:     for (c = 0; c < Nc; ++c) DMPlexSetConeSize(*dm, c, Ncn);
4902:     DMSetUp(*dm);
4903:     for (c = 0; c < Nc; ++c) {
4904:       PetscViewerRead(viewer, line, Ncn, NULL, PETSC_STRING);
4905:       switch (Ncn) {
4906:       case 2:
4907:         snum = sscanf(line, format, &vbuf[0], &vbuf[1]);
4908:         break;
4909:       case 3:
4910:         snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2]);
4911:         break;
4912:       case 4:
4913:         snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3]);
4914:         break;
4915:       case 6:
4916:         snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5]);
4917:         break;
4918:       case 8:
4919:         snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);
4920:         break;
4921:       default:
4922:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell shape with %d vertices", Ncn);
4923:       }
4925:       for (v = 0; v < Ncn; ++v) cone[v] = vbuf[v] + Nc;
4926:       /* Hexahedra are inverted */
4927:       if (Ncn == 8) {
4928:         PetscInt tmp = cone[1];
4929:         cone[1]      = cone[3];
4930:         cone[3]      = tmp;
4931:       }
4932:       DMPlexSetCone(*dm, c, cone);
4933:     }
4934:   }
4935:   DMPlexSymmetrize(*dm);
4936:   DMPlexStratify(*dm);
4937:   /* Read coordinates */
4938:   DMGetCoordinateSection(*dm, &coordSection);
4939:   PetscSectionSetNumFields(coordSection, 1);
4940:   PetscSectionSetFieldComponents(coordSection, 0, cdim);
4941:   PetscSectionSetChart(coordSection, Nc, Nc + Nv);
4942:   for (v = Nc; v < Nc + Nv; ++v) {
4943:     PetscSectionSetDof(coordSection, v, cdim);
4944:     PetscSectionSetFieldDof(coordSection, v, 0, cdim);
4945:   }
4946:   PetscSectionSetUp(coordSection);
4947:   PetscSectionGetStorageSize(coordSection, &coordSize);
4948:   VecCreate(PETSC_COMM_SELF, &coordinates);
4949:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
4950:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
4951:   VecSetBlockSize(coordinates, cdim);
4952:   VecSetType(coordinates, VECSTANDARD);
4953:   VecGetArray(coordinates, &coords);
4954:   if (rank == 0) {
4955:     char   format[PETSC_MAX_PATH_LEN];
4956:     double x[3];
4957:     int    l, val[3];

4959:     if (Nl) {
4960:       for (l = 0; l < Nl; ++l) {
4961:         format[l * 3 + 0] = '%';
4962:         format[l * 3 + 1] = 'd';
4963:         format[l * 3 + 2] = ' ';
4964:       }
4965:       format[Nl * 3 - 1] = '\0';
4966:       DMCreateLabel(*dm, "marker");
4967:       DMGetLabel(*dm, "marker", &marker);
4968:     }
4969:     for (v = 0; v < Nv; ++v) {
4970:       PetscViewerRead(viewer, line, 3 + Nl, NULL, PETSC_STRING);
4971:       snum = sscanf(line, "%lg %lg %lg", &x[0], &x[1], &x[2]);
4973:       switch (Nl) {
4974:       case 0:
4975:         snum = 0;
4976:         break;
4977:       case 1:
4978:         snum = sscanf(line, format, &val[0]);
4979:         break;
4980:       case 2:
4981:         snum = sscanf(line, format, &val[0], &val[1]);
4982:         break;
4983:       case 3:
4984:         snum = sscanf(line, format, &val[0], &val[1], &val[2]);
4985:         break;
4986:       default:
4987:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Request support for %d labels", Nl);
4988:       }
4990:       for (d = 0; d < cdim; ++d) coords[v * cdim + d] = x[d];
4991:       for (l = 0; l < Nl; ++l) DMLabelSetValue(marker, v + Nc, val[l]);
4992:     }
4993:   }
4994:   VecRestoreArray(coordinates, &coords);
4995:   DMSetCoordinatesLocal(*dm, coordinates);
4996:   VecDestroy(&coordinates);
4997:   PetscViewerDestroy(&viewer);
4998:   if (interpolate) {
4999:     DM      idm;
5000:     DMLabel bdlabel;

5002:     DMPlexInterpolate(*dm, &idm);
5003:     DMDestroy(dm);
5004:     *dm = idm;

5006:     if (!Nl) {
5007:       DMCreateLabel(*dm, "marker");
5008:       DMGetLabel(*dm, "marker", &bdlabel);
5009:       DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);
5010:       DMPlexLabelComplete(*dm, bdlabel);
5011:     }
5012:   }
5013:   return 0;
5014: }

5016: /*@C
5017:   DMPlexCreateFromFile - This takes a filename and produces a DM

5019:   Input Parameters:
5020: + comm - The communicator
5021: . filename - A file name
5022: . plexname - The object name of the resulting DM, also used for intra-datafile lookup by some formats
5023: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

5025:   Output Parameter:
5026: . dm - The DM

5028:   Options Database Keys:
5029: . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5

5031:   Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g.
5032: $ -dm_plex_create_viewer_hdf5_collective

5034:   Notes:
5035:   Using PETSCVIEWERHDF5 type with PETSC_VIEWER_HDF5_PETSC format, one can save multiple DMPlex
5036:   meshes in a single HDF5 file. This in turn requires one to name the DMPlex object with PetscObjectSetName()
5037:   before saving it with DMView() and before loading it with DMLoad() for identification of the mesh object.
5038:   The input parameter name is thus used to name the DMPlex object when DMPlexCreateFromFile() internally
5039:   calls DMLoad(). Currently, name is ignored for other viewer types and/or formats.

5041:   Level: beginner

5043: .seealso: `DMPlexCreateFromDAG()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCreate()`, `PetscObjectSetName()`, `DMView()`, `DMLoad()`
5044: @*/
5045: PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], const char plexname[], PetscBool interpolate, DM *dm)
5046: {
5047:   const char  extGmsh[]      = ".msh";
5048:   const char  extGmsh2[]     = ".msh2";
5049:   const char  extGmsh4[]     = ".msh4";
5050:   const char  extCGNS[]      = ".cgns";
5051:   const char  extExodus[]    = ".exo";
5052:   const char  extExodus_e[]  = ".e";
5053:   const char  extGenesis[]   = ".gen";
5054:   const char  extFluent[]    = ".cas";
5055:   const char  extHDF5[]      = ".h5";
5056:   const char  extMed[]       = ".med";
5057:   const char  extPLY[]       = ".ply";
5058:   const char  extEGADSLite[] = ".egadslite";
5059:   const char  extEGADS[]     = ".egads";
5060:   const char  extIGES[]      = ".igs";
5061:   const char  extSTEP[]      = ".stp";
5062:   const char  extCV[]        = ".dat";
5063:   size_t      len;
5064:   PetscBool   isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV;
5065:   PetscMPIInt rank;

5070:   DMInitializePackage();
5071:   PetscLogEventBegin(DMPLEX_CreateFromFile, 0, 0, 0, 0);
5072:   MPI_Comm_rank(comm, &rank);
5073:   PetscStrlen(filename, &len);

5076: #define CheckExtension(extension__, is_extension__) \
5077:   do { \
5078:     PetscAssert(sizeof(extension__), comm, PETSC_ERR_PLIB, "Zero-size extension: %s", extension__); \
5079:     /* don't count the null-terminator at the end */ \
5080:     const size_t ext_len = sizeof(extension__) - 1; \
5081:     if (len < ext_len) { \
5082:       is_extension__ = PETSC_FALSE; \
5083:     } else { \
5084:       PetscStrncmp(filename + len - ext_len, extension__, ext_len, &is_extension__); \
5085:     } \
5086:   } while (0)

5088:   CheckExtension(extGmsh, isGmsh);
5089:   CheckExtension(extGmsh2, isGmsh2);
5090:   CheckExtension(extGmsh4, isGmsh4);
5091:   CheckExtension(extCGNS, isCGNS);
5092:   CheckExtension(extExodus, isExodus);
5093:   if (!isExodus) CheckExtension(extExodus_e, isExodus);
5094:   CheckExtension(extGenesis, isGenesis);
5095:   CheckExtension(extFluent, isFluent);
5096:   CheckExtension(extHDF5, isHDF5);
5097:   CheckExtension(extMed, isMed);
5098:   CheckExtension(extPLY, isPLY);
5099:   CheckExtension(extEGADSLite, isEGADSLite);
5100:   CheckExtension(extEGADS, isEGADS);
5101:   CheckExtension(extIGES, isIGES);
5102:   CheckExtension(extSTEP, isSTEP);
5103:   CheckExtension(extCV, isCV);

5105: #undef CheckExtension

5107:   if (isGmsh || isGmsh2 || isGmsh4) {
5108:     DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);
5109:   } else if (isCGNS) {
5110:     DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);
5111:   } else if (isExodus || isGenesis) {
5112:     DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);
5113:   } else if (isFluent) {
5114:     DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);
5115:   } else if (isHDF5) {
5116:     PetscBool   load_hdf5_xdmf = PETSC_FALSE;
5117:     PetscViewer viewer;

5119:     /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */
5120:     PetscStrncmp(&filename[PetscMax(0, len - 8)], ".xdmf", 5, &load_hdf5_xdmf);
5121:     PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);
5122:     PetscViewerCreate(comm, &viewer);
5123:     PetscViewerSetType(viewer, PETSCVIEWERHDF5);
5124:     PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");
5125:     PetscViewerSetFromOptions(viewer);
5126:     PetscViewerFileSetMode(viewer, FILE_MODE_READ);
5127:     PetscViewerFileSetName(viewer, filename);

5129:     DMCreate(comm, dm);
5130:     PetscObjectSetName((PetscObject)(*dm), plexname);
5131:     DMSetType(*dm, DMPLEX);
5132:     if (load_hdf5_xdmf) PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);
5133:     DMLoad(*dm, viewer);
5134:     if (load_hdf5_xdmf) PetscViewerPopFormat(viewer);
5135:     PetscViewerDestroy(&viewer);

5137:     if (interpolate) {
5138:       DM idm;

5140:       DMPlexInterpolate(*dm, &idm);
5141:       DMDestroy(dm);
5142:       *dm = idm;
5143:     }
5144:   } else if (isMed) {
5145:     DMPlexCreateMedFromFile(comm, filename, interpolate, dm);
5146:   } else if (isPLY) {
5147:     DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);
5148:   } else if (isEGADSLite || isEGADS || isIGES || isSTEP) {
5149:     if (isEGADSLite) DMPlexCreateEGADSLiteFromFile(comm, filename, dm);
5150:     else DMPlexCreateEGADSFromFile(comm, filename, dm);
5151:     if (!interpolate) {
5152:       DM udm;

5154:       DMPlexUninterpolate(*dm, &udm);
5155:       DMDestroy(dm);
5156:       *dm = udm;
5157:     }
5158:   } else if (isCV) {
5159:     DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);
5160:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename);
5161:   PetscStrlen(plexname, &len);
5162:   if (len) PetscObjectSetName((PetscObject)(*dm), plexname);
5163:   PetscLogEventEnd(DMPLEX_CreateFromFile, 0, 0, 0, 0);
5164:   return 0;
5165: }