Actual source code: plexcreate.c

petsc-3.4.4 2014-03-13
  1: #define PETSCDM_DLL
  2: #include <petsc-private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
  3: #include <petscdmda.h>
  4: #include <petscsf.h>

  8: PetscErrorCode  DMSetFromOptions_Plex(DM dm)
  9: {
 10:   DM_Plex        *mesh = (DM_Plex*) dm->data;

 15:   PetscOptionsHead("DMPlex Options");
 16:   /* Handle DMPlex refinement */
 17:   /* Handle associated vectors */
 18:   /* Handle viewing */
 19:   PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);
 20:   PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);
 21:   PetscOptionsTail();
 22:   return(0);
 23: }

 27: /*
 28:  Simple square boundary:

 30:  18--5-17--4--16
 31:   |     |     |
 32:   6    10     3
 33:   |     |     |
 34:  19-11-20--9--15
 35:   |     |     |
 36:   7     8     2
 37:   |     |     |
 38:  12--0-13--1--14
 39: */
 40: PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
 41: {
 42:   PetscInt       numVertices    = (edges[0]+1)*(edges[1]+1);
 43:   PetscInt       numEdges       = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
 44:   PetscInt       markerTop      = 1;
 45:   PetscInt       markerBottom   = 1;
 46:   PetscInt       markerRight    = 1;
 47:   PetscInt       markerLeft     = 1;
 48:   PetscBool      markerSeparate = PETSC_FALSE;
 49:   Vec            coordinates;
 50:   PetscSection   coordSection;
 51:   PetscScalar    *coords;
 52:   PetscInt       coordSize;
 53:   PetscMPIInt    rank;
 54:   PetscInt       v, vx, vy;

 58:   PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
 59:   if (markerSeparate) {
 60:     markerTop    = 1;
 61:     markerBottom = 0;
 62:     markerRight  = 0;
 63:     markerLeft   = 0;
 64:   }
 65:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
 66:   if (!rank) {
 67:     PetscInt e, ex, ey;

 69:     DMPlexSetChart(dm, 0, numEdges+numVertices);
 70:     for (e = 0; e < numEdges; ++e) {
 71:       DMPlexSetConeSize(dm, e, 2);
 72:     }
 73:     DMSetUp(dm); /* Allocate space for cones */
 74:     for (vx = 0; vx <= edges[0]; vx++) {
 75:       for (ey = 0; ey < edges[1]; ey++) {
 76:         PetscInt edge   = vx*edges[1] + ey + edges[0]*(edges[1]+1);
 77:         PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
 78:         PetscInt cone[2];

 80:         cone[0] = vertex; cone[1] = vertex+edges[0]+1;
 81:         DMPlexSetCone(dm, edge, cone);
 82:         if (vx == edges[0]) {
 83:           DMPlexSetLabelValue(dm, "marker", edge,    markerRight);
 84:           DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
 85:           if (ey == edges[1]-1) {
 86:             DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
 87:           }
 88:         } else if (vx == 0) {
 89:           DMPlexSetLabelValue(dm, "marker", edge,    markerLeft);
 90:           DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
 91:           if (ey == edges[1]-1) {
 92:             DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
 93:           }
 94:         }
 95:       }
 96:     }
 97:     for (vy = 0; vy <= edges[1]; vy++) {
 98:       for (ex = 0; ex < edges[0]; ex++) {
 99:         PetscInt edge   = vy*edges[0]     + ex;
100:         PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
101:         PetscInt cone[2];

103:         cone[0] = vertex; cone[1] = vertex+1;
104:         DMPlexSetCone(dm, edge, cone);
105:         if (vy == edges[1]) {
106:           DMPlexSetLabelValue(dm, "marker", edge,    markerTop);
107:           DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
108:           if (ex == edges[0]-1) {
109:             DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
110:           }
111:         } else if (vy == 0) {
112:           DMPlexSetLabelValue(dm, "marker", edge,    markerBottom);
113:           DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
114:           if (ex == edges[0]-1) {
115:             DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
116:           }
117:         }
118:       }
119:     }
120:   }
121:   DMPlexSymmetrize(dm);
122:   DMPlexStratify(dm);
123:   /* Build coordinates */
124:   DMPlexGetCoordinateSection(dm, &coordSection);
125:   PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);
126:   for (v = numEdges; v < numEdges+numVertices; ++v) {
127:     PetscSectionSetDof(coordSection, v, 2);
128:   }
129:   PetscSectionSetUp(coordSection);
130:   PetscSectionGetStorageSize(coordSection, &coordSize);
131:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
132:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
133:   VecSetFromOptions(coordinates);
134:   VecGetArray(coordinates, &coords);
135:   for (vy = 0; vy <= edges[1]; ++vy) {
136:     for (vx = 0; vx <= edges[0]; ++vx) {
137:       coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
138:       coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
139:     }
140:   }
141:   VecRestoreArray(coordinates, &coords);
142:   DMSetCoordinatesLocal(dm, coordinates);
143:   VecDestroy(&coordinates);
144:   return(0);
145: }

149: /*
150:  Simple cubic boundary:

152:      2-------3
153:     /|      /|
154:    6-------7 |
155:    | |     | |
156:    | 0-----|-1
157:    |/      |/
158:    4-------5
159: */
160: PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
161: {
162:   PetscInt       numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
163:   PetscInt       numFaces    = 6;
164:   Vec            coordinates;
165:   PetscSection   coordSection;
166:   PetscScalar    *coords;
167:   PetscInt       coordSize;
168:   PetscMPIInt    rank;
169:   PetscInt       v, vx, vy, vz;

173:   if ((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Must have at least 1 face per side");
174:   if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Currently can't handle more than 1 face per side");
175:   PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);
176:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
177:   if (!rank) {
178:     PetscInt f;

180:     DMPlexSetChart(dm, 0, numFaces+numVertices);
181:     for (f = 0; f < numFaces; ++f) {
182:       DMPlexSetConeSize(dm, f, 4);
183:     }
184:     DMSetUp(dm); /* Allocate space for cones */
185:     for (v = 0; v < numFaces+numVertices; ++v) {
186:       DMPlexSetLabelValue(dm, "marker", v, 1);
187:     }
188:     { /* Side 0 (Front) */
189:       PetscInt cone[4];
190:       cone[0] = numFaces+4; cone[1] = numFaces+5; cone[2] = numFaces+7; cone[3] = numFaces+6;
191:       DMPlexSetCone(dm, 0, cone);
192:     }
193:     { /* Side 1 (Back) */
194:       PetscInt cone[4];
195:       cone[0] = numFaces+1; cone[1] = numFaces+0; cone[2] = numFaces+2; cone[3] = numFaces+3;
196:       DMPlexSetCone(dm, 1, cone);
197:     }
198:     { /* Side 2 (Bottom) */
199:       PetscInt cone[4];
200:       cone[0] = numFaces+0; cone[1] = numFaces+1; cone[2] = numFaces+5; cone[3] = numFaces+4;
201:       DMPlexSetCone(dm, 2, cone);
202:     }
203:     { /* Side 3 (Top) */
204:       PetscInt cone[4];
205:       cone[0] = numFaces+6; cone[1] = numFaces+7; cone[2] = numFaces+3; cone[3] = numFaces+2;
206:       DMPlexSetCone(dm, 3, cone);
207:     }
208:     { /* Side 4 (Left) */
209:       PetscInt cone[4];
210:       cone[0] = numFaces+0; cone[1] = numFaces+4; cone[2] = numFaces+6; cone[3] = numFaces+2;
211:       DMPlexSetCone(dm, 4, cone);
212:     }
213:     { /* Side 5 (Right) */
214:       PetscInt cone[4];
215:       cone[0] = numFaces+5; cone[1] = numFaces+1; cone[2] = numFaces+3; cone[3] = numFaces+7;
216:       DMPlexSetCone(dm, 5, cone);
217:     }
218:   }
219:   DMPlexSymmetrize(dm);
220:   DMPlexStratify(dm);
221:   /* Build coordinates */
222:   DMPlexGetCoordinateSection(dm, &coordSection);
223:   PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);
224:   for (v = numFaces; v < numFaces+numVertices; ++v) {
225:     PetscSectionSetDof(coordSection, v, 3);
226:   }
227:   PetscSectionSetUp(coordSection);
228:   PetscSectionGetStorageSize(coordSection, &coordSize);
229:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
230:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
231:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
232:   VecSetFromOptions(coordinates);
233:   VecGetArray(coordinates, &coords);
234:   for (vz = 0; vz <= faces[2]; ++vz) {
235:     for (vy = 0; vy <= faces[1]; ++vy) {
236:       for (vx = 0; vx <= faces[0]; ++vx) {
237:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
238:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
239:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
240:       }
241:     }
242:   }
243:   VecRestoreArray(coordinates, &coords);
244:   DMSetCoordinatesLocal(dm, coordinates);
245:   VecDestroy(&coordinates);
246:   return(0);
247: }

251: /*
252:  Simple square mesh:

254:  22--8-23--9--24
255:   |     |     |
256:  13  2 14  3  15
257:   |     |     |
258:  19--6-20--7--21
259:   |     |     |
260:  10  0 11  1 12
261:   |     |     |
262:  16--4-17--5--18
263: */
264: PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
265: {
266:   PetscInt       markerTop      = 1;
267:   PetscInt       markerBottom   = 1;
268:   PetscInt       markerRight    = 1;
269:   PetscInt       markerLeft     = 1;
270:   PetscBool      markerSeparate = PETSC_FALSE;
271:   PetscMPIInt    rank;

275:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
276:   PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
277:   if (markerSeparate) {
278:     markerTop    = 3;
279:     markerBottom = 1;
280:     markerRight  = 2;
281:     markerLeft   = 4;
282:   }
283:   {
284:     const PetscInt numXEdges    = !rank ? edges[0]   : 0;
285:     const PetscInt numYEdges    = !rank ? edges[1]   : 0;
286:     const PetscInt numXVertices = !rank ? edges[0]+1 : 0;
287:     const PetscInt numYVertices = !rank ? edges[1]+1 : 0;
288:     const PetscInt numTotXEdges = numXEdges*numYVertices;
289:     const PetscInt numTotYEdges = numYEdges*numXVertices;
290:     const PetscInt numVertices  = numXVertices*numYVertices;
291:     const PetscInt numEdges     = numTotXEdges + numTotYEdges;
292:     const PetscInt numFaces     = numXEdges*numYEdges;
293:     const PetscInt firstVertex  = numFaces;
294:     const PetscInt firstXEdge   = numFaces + numVertices;
295:     const PetscInt firstYEdge   = numFaces + numVertices + numTotXEdges;
296:     Vec            coordinates;
297:     PetscSection   coordSection;
298:     PetscScalar    *coords;
299:     PetscInt       coordSize;
300:     PetscInt       v, vx, vy;
301:     PetscInt       f, fx, fy, e, ex, ey;

303:     DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);
304:     for (f = 0; f < numFaces; ++f) {
305:       DMPlexSetConeSize(dm, f, 4);
306:     }
307:     for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
308:       DMPlexSetConeSize(dm, e, 2);
309:     }
310:     DMSetUp(dm); /* Allocate space for cones */
311:     /* Build faces */
312:     for (fy = 0; fy < numYEdges; fy++) {
313:       for (fx = 0; fx < numXEdges; fx++) {
314:         const PetscInt face    = fy*numXEdges + fx;
315:         const PetscInt edgeL   = firstYEdge + fx*numYEdges + fy;
316:         const PetscInt edgeB   = firstXEdge + fy*numXEdges + fx;
317:         const PetscInt ornt[4] = {0,     0,               -2,              -2};
318:         PetscInt       cone[4];

320:         cone[0] = edgeB; cone[1] = edgeL+numYEdges; cone[2] = edgeB+numXEdges; cone[3] = edgeL;
321:         DMPlexSetCone(dm, face, cone);
322:         DMPlexSetConeOrientation(dm, face, ornt);
323:       }
324:     }
325:     /* Build Y edges*/
326:     for (vx = 0; vx < numXVertices; vx++) {
327:       for (ey = 0; ey < numYEdges; ey++) {
328:         const PetscInt edge   = firstYEdge  + vx*numYEdges + ey;
329:         const PetscInt vertex = firstVertex + ey*numXVertices + vx;
330:         PetscInt       cone[2];

332:         cone[0] = vertex; cone[1] = vertex+numXVertices;
333:         DMPlexSetCone(dm, edge, cone);
334:         if (vx == numXVertices-1) {
335:           DMPlexSetLabelValue(dm, "marker", edge,    markerRight);
336:           DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
337:           if (ey == numYEdges-1) {
338:             DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
339:           }
340:         } else if (vx == 0) {
341:           DMPlexSetLabelValue(dm, "marker", edge,    markerLeft);
342:           DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
343:           if (ey == numYEdges-1) {
344:             DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
345:           }
346:         }
347:       }
348:     }
349:     /* Build X edges*/
350:     for (vy = 0; vy < numYVertices; vy++) {
351:       for (ex = 0; ex < numXEdges; ex++) {
352:         const PetscInt edge   = firstXEdge  + vy*numXEdges + ex;
353:         const PetscInt vertex = firstVertex + vy*numXVertices + ex;
354:         PetscInt       cone[2];

356:         cone[0] = vertex; cone[1] = vertex+1;
357:         DMPlexSetCone(dm, edge, cone);
358:         if (vy == numYVertices-1) {
359:           DMPlexSetLabelValue(dm, "marker", edge,    markerTop);
360:           DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
361:           if (ex == numXEdges-1) {
362:             DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
363:           }
364:         } else if (vy == 0) {
365:           DMPlexSetLabelValue(dm, "marker", edge,    markerBottom);
366:           DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
367:           if (ex == numXEdges-1) {
368:             DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
369:           }
370:         }
371:       }
372:     }
373:     DMPlexSymmetrize(dm);
374:     DMPlexStratify(dm);
375:     /* Build coordinates */
376:     DMPlexGetCoordinateSection(dm, &coordSection);
377:     PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);
378:     for (v = firstVertex; v < firstVertex+numVertices; ++v) {
379:       PetscSectionSetDof(coordSection, v, 2);
380:     }
381:     PetscSectionSetUp(coordSection);
382:     PetscSectionGetStorageSize(coordSection, &coordSize);
383:     VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
384:     VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
385:     VecSetFromOptions(coordinates);
386:     VecGetArray(coordinates, &coords);
387:     for (vy = 0; vy < numYVertices; ++vy) {
388:       for (vx = 0; vx < numXVertices; ++vx) {
389:         coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
390:         coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
391:       }
392:     }
393:     VecRestoreArray(coordinates, &coords);
394:     DMSetCoordinatesLocal(dm, coordinates);
395:     VecDestroy(&coordinates);
396:   }
397:   return(0);
398: }

402: /*@
403:   DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box).

405:   Collective on MPI_Comm

407:   Input Parameters:
408: + comm - The communicator for the DM object
409: . dim - The spatial dimension
410: - interpolate - Flag to create intermediate mesh pieces (edges, faces)

412:   Output Parameter:
413: . dm  - The DM object

415:   Level: beginner

417: .keywords: DM, create
418: .seealso: DMSetType(), DMCreate()
419: @*/
420: PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm)
421: {
422:   DM             boundary;

427:   DMCreate(comm, &boundary);
429:   DMSetType(boundary, DMPLEX);
430:   DMPlexSetDimension(boundary, dim-1);
431:   switch (dim) {
432:   case 2:
433:   {
434:     PetscReal lower[2] = {0.0, 0.0};
435:     PetscReal upper[2] = {1.0, 1.0};
436:     PetscInt  edges[2] = {2, 2};

438:     DMPlexCreateSquareBoundary(boundary, lower, upper, edges);
439:     break;
440:   }
441:   case 3:
442:   {
443:     PetscReal lower[3] = {0.0, 0.0, 0.0};
444:     PetscReal upper[3] = {1.0, 1.0, 1.0};
445:     PetscInt  faces[3] = {1, 1, 1};

447:     DMPlexCreateCubeBoundary(boundary, lower, upper, faces);
448:     break;
449:   }
450:   default:
451:     SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
452:   }
453:   DMPlexGenerate(boundary, NULL, interpolate, dm);
454:   DMDestroy(&boundary);
455:   return(0);
456: }

460: PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm)
461: {

466:   DMCreate(comm, dm);
468:   DMSetType(*dm, DMPLEX);
469:   DMPlexSetDimension(*dm, dim);
470:   switch (dim) {
471:   case 2:
472:   {
473:     PetscReal lower[2] = {0.0, 0.0};
474:     PetscReal upper[2] = {1.0, 1.0};

476:     DMPlexCreateSquareMesh(*dm, lower, upper, cells);
477:     break;
478:   }
479: #if 0
480:   case 3:
481:   {
482:     PetscReal lower[3] = {0.0, 0.0, 0.0};
483:     PetscReal upper[3] = {1.0, 1.0, 1.0};

485:     DMPlexCreateCubeMesh(boundary, lower, upper, cells);
486:     break;
487:   }
488: #endif
489:   default:
490:     SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
491:   }
492:   return(0);
493: }

495: /* External function declarations here */
496: extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
497: extern PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J);
498: extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
499: extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined);
500: extern PetscErrorCode DMSetUp_Plex(DM dm);
501: extern PetscErrorCode DMDestroy_Plex(DM dm);
502: extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
503: extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
504: extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS);

508: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
509: {

513:   DMCreateGlobalVector_Section_Private(dm,vec);
514:   /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
515:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);
516:   return(0);
517: }

521: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
522: {

526:   DMCreateLocalVector_Section_Private(dm,vec);
527:   VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);
528:   return(0);
529: }


534: PetscErrorCode DMInitialize_Plex(DM dm)
535: {

539:   PetscStrallocpy(VECSTANDARD, (char**)&dm->vectype);

541:   dm->ops->view                            = DMView_Plex;
542:   dm->ops->setfromoptions                  = DMSetFromOptions_Plex;
543:   dm->ops->setup                           = DMSetUp_Plex;
544:   dm->ops->createglobalvector              = DMCreateGlobalVector_Plex;
545:   dm->ops->createlocalvector               = DMCreateLocalVector_Plex;
546:   dm->ops->getlocaltoglobalmapping         = NULL;
547:   dm->ops->getlocaltoglobalmappingblock    = NULL;
548:   dm->ops->createfieldis                   = NULL;
549:   dm->ops->createcoordinatedm              = DMCreateCoordinateDM_Plex;
550:   dm->ops->getcoloring                     = 0;
551:   dm->ops->creatematrix                    = DMCreateMatrix_Plex;
552:   dm->ops->createinterpolation             = 0;
553:   dm->ops->getaggregates                   = 0;
554:   dm->ops->getinjection                    = 0;
555:   dm->ops->refine                          = DMRefine_Plex;
556:   dm->ops->coarsen                         = 0;
557:   dm->ops->refinehierarchy                 = 0;
558:   dm->ops->coarsenhierarchy                = 0;
559:   dm->ops->globaltolocalbegin              = NULL;
560:   dm->ops->globaltolocalend                = NULL;
561:   dm->ops->localtoglobalbegin              = NULL;
562:   dm->ops->localtoglobalend                = NULL;
563:   dm->ops->destroy                         = DMDestroy_Plex;
564:   dm->ops->createsubdm                     = DMCreateSubDM_Plex;
565:   dm->ops->locatepoints                    = DMLocatePoints_Plex;
566:   return(0);
567: }

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

575:   Level: intermediate

577: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
578: M*/

582: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
583: {
584:   DM_Plex        *mesh;
585:   PetscInt       unit, d;

590:   PetscNewLog(dm, DM_Plex, &mesh);
591:   dm->data = mesh;

593:   mesh->refct             = 1;
594:   mesh->dim               = 0;
595:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
596:   mesh->maxConeSize       = 0;
597:   mesh->cones             = NULL;
598:   mesh->coneOrientations  = NULL;
599:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
600:   mesh->maxSupportSize    = 0;
601:   mesh->supports          = NULL;
602:   mesh->refinementUniform = PETSC_TRUE;
603:   mesh->refinementLimit   = -1.0;

605:   mesh->facesTmp = NULL;

607:   mesh->subpointMap = NULL;

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

611:   mesh->labels              = NULL;
612:   mesh->globalVertexNumbers = NULL;
613:   mesh->globalCellNumbers   = NULL;
614:   for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
615:   mesh->vtkCellHeight     = 0;
616:   mesh->preallocCenterDim = -1;

618:   mesh->integrateResidualFEM       = NULL;
619:   mesh->integrateJacobianActionFEM = NULL;
620:   mesh->integrateJacobianFEM       = NULL;

622:   mesh->printSetValues = PETSC_FALSE;
623:   mesh->printFEM       = 0;

625:   DMInitialize_Plex(dm);
626:   return(0);
627: }

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

634:   Collective on MPI_Comm

636:   Input Parameter:
637: . comm - The communicator for the DMPlex object

639:   Output Parameter:
640: . mesh  - The DMPlex object

642:   Level: beginner

644: .keywords: DMPlex, create
645: @*/
646: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
647: {

652:   DMCreate(comm, mesh);
653:   DMSetType(*mesh, DMPLEX);
654:   return(0);
655: }

659: /*@
660:   DMPlexClone - Creates a DMPlex object with the same mesh as the original.

662:   Collective on MPI_Comm

664:   Input Parameter:
665: . dm - The original DMPlex object

667:   Output Parameter:
668: . newdm  - The new DMPlex object

670:   Level: beginner

672: .keywords: DMPlex, create
673: @*/
674: PetscErrorCode DMPlexClone(DM dm, DM *newdm)
675: {
676:   DM_Plex        *mesh;
677:   Vec             coords;
678:   void           *ctx;

684:   DMCreate(PetscObjectComm((PetscObject)dm), newdm);
685:   PetscSFDestroy(&(*newdm)->sf);
686:   PetscObjectReference((PetscObject) dm->sf);
687:   (*newdm)->sf = dm->sf;
688:   mesh         = (DM_Plex*) dm->data;
689:   mesh->refct++;
690:   (*newdm)->data = mesh;
691:   PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
692:   DMInitialize_Plex(*newdm);
693:   DMGetApplicationContext(dm, &ctx);
694:   DMSetApplicationContext(*newdm, ctx);
695:   DMGetCoordinatesLocal(dm, &coords);
696:   if (coords) {
697:     DMSetCoordinatesLocal(*newdm, coords);
698:   } else {
699:     DMGetCoordinates(dm, &coords);
700:     if (coords) {DMSetCoordinates(*newdm, coords);}
701:   }
702:   return(0);
703: }

707: /*
708:   This takes as input the common mesh generator output, a list of the vertices for each cell
709: */
710: PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[])
711: {
712:   PetscInt      *cone, c, p;

716:   DMPlexSetChart(dm, 0, numCells+numVertices);
717:   for (c = 0; c < numCells; ++c) {
718:     DMPlexSetConeSize(dm, c, numCorners);
719:   }
720:   DMSetUp(dm);
721:   DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);
722:   for (c = 0; c < numCells; ++c) {
723:     for (p = 0; p < numCorners; ++p) {
724:       cone[p] = cells[c*numCorners+p]+numCells;
725:     }
726:     DMPlexSetCone(dm, c, cone);
727:   }
728:   DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);
729:   DMPlexSymmetrize(dm);
730:   DMPlexStratify(dm);
731:   return(0);
732: }

736: /*
737:   This takes as input the coordinates for each vertex
738: */
739: PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
740: {
741:   PetscSection   coordSection;
742:   Vec            coordinates;
743:   PetscScalar   *coords;
744:   PetscInt       coordSize, v, d;

748:   DMPlexGetCoordinateSection(dm, &coordSection);
749:   PetscSectionSetNumFields(coordSection, 1);
750:   PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
751:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
752:   for (v = numCells; v < numCells+numVertices; ++v) {
753:     PetscSectionSetDof(coordSection, v, spaceDim);
754:     PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
755:   }
756:   PetscSectionSetUp(coordSection);
757:   PetscSectionGetStorageSize(coordSection, &coordSize);
758:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
759:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
760:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
761:   VecSetFromOptions(coordinates);
762:   VecGetArray(coordinates, &coords);
763:   for (v = 0; v < numVertices; ++v) {
764:     for (d = 0; d < spaceDim; ++d) {
765:       coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
766:     }
767:   }
768:   VecRestoreArray(coordinates, &coords);
769:   DMSetCoordinatesLocal(dm, coordinates);
770:   VecDestroy(&coordinates);
771:   return(0);
772: }

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

779:   Input Parameters:
780: + comm - The communicator
781: . dim - The topological dimension of the mesh
782: . numCells - The number of cells
783: . numVertices - The number of vertices
784: . numCorners - The number of vertices for each cell
785: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
786: . cells - An array of numCells*numCorners numbers, the vertices for each cell
787: . spaceDim - The spatial dimension used for coordinates
788: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex

790:   Output Parameter:
791: . dm - The DM

793:   Note: Two triangles sharing a face
794: $
795: $        2
796: $      / | \
797: $     /  |  \
798: $    /   |   \
799: $   0  0 | 1  3
800: $    \   |   /
801: $     \  |  /
802: $      \ | /
803: $        1
804: would have input
805: $  numCells = 2, numVertices = 4
806: $  cells = [0 1 2  1 3 2]
807: $
808: which would result in the DMPlex
809: $
810: $        4
811: $      / | \
812: $     /  |  \
813: $    /   |   \
814: $   2  0 | 1  5
815: $    \   |   /
816: $     \  |  /
817: $      \ | /
818: $        3

820:   Level: beginner

822: .seealso: DMPlexCreate()
823: @*/
824: PetscErrorCode DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm)
825: {

829:   DMCreate(comm, dm);
830:   DMSetType(*dm, DMPLEX);
831:   DMPlexSetDimension(*dm, dim);
832:   DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);
833:   if (interpolate) {
834:     DM idm;

836:     DMPlexInterpolate(*dm, &idm);
837:     DMDestroy(dm);
838:     *dm  = idm;
839:   }
840:   DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);
841:   return(0);
842: }

846: /*
847:   This takes as input the raw Hasse Diagram data
848: */
849: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
850: {
851:   Vec            coordinates;
852:   PetscSection   coordSection;
853:   PetscScalar    *coords;
854:   PetscInt       coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off;

858:   DMPlexGetDimension(dm, &dim);
859:   for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
860:   DMPlexSetChart(dm, pStart, pEnd);
861:   for (p = pStart; p < pEnd; ++p) {
862:     DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
863:   }
864:   DMSetUp(dm); /* Allocate space for cones */
865:   for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
866:     DMPlexSetCone(dm, p, &cones[off]);
867:     DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
868:   }
869:   DMPlexSymmetrize(dm);
870:   DMPlexStratify(dm);
871:   /* Build coordinates */
872:   DMPlexGetCoordinateSection(dm, &coordSection);
873:   PetscSectionSetNumFields(coordSection, 1);
874:   PetscSectionSetFieldComponents(coordSection, 0, dim);
875:   PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
876:   for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
877:     PetscSectionSetDof(coordSection, v, dim);
878:     PetscSectionSetFieldDof(coordSection, v, 0, dim);
879:   }
880:   PetscSectionSetUp(coordSection);
881:   PetscSectionGetStorageSize(coordSection, &coordSize);
882:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
883:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
884:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
885:   VecSetFromOptions(coordinates);
886:   VecGetArray(coordinates, &coords);
887:   for (v = 0; v < numPoints[0]; ++v) {
888:     PetscInt off;

890:     PetscSectionGetOffset(coordSection, v+firstVertex, &off);
891:     for (d = 0; d < dim; ++d) {
892:       coords[off+d] = vertexCoords[v*dim+d];
893:     }
894:   }
895:   VecRestoreArray(coordinates, &coords);
896:   DMSetCoordinatesLocal(dm, coordinates);
897:   VecDestroy(&coordinates);
898:   return(0);
899: }