Actual source code: meshcreate.c

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

  7: PetscErrorCode  DMSetFromOptions_Mesh(DM dm)
  8: {
  9:   PetscBool      flg;

 14:   PetscOptionsHead("DMMesh Options");
 15:   /* Handle DMMesh refinement */
 16:   /* Handle associated vectors */
 17:   /* Handle viewing */
 18:   PetscOptionsBool("-dm_mesh_view_vtk", "Output mesh in VTK format", "DMView", PETSC_FALSE, &flg, NULL);
 19:   if (flg) {
 20:     PetscViewer viewer;

 22:     PetscViewerCreate(((PetscObject) dm)->comm, &viewer);
 23:     PetscViewerSetType(viewer, PETSCVIEWERASCII);
 24:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
 25:     PetscViewerFileSetName(viewer, "mesh.vtk");
 26:     DMView(dm, viewer);
 27:     PetscViewerDestroy(&viewer);
 28:   }
 29:   PetscOptionsBool("-dm_mesh_view", "Exhaustive mesh description", "DMView", PETSC_FALSE, &flg, NULL);
 30:   if (flg) {
 31:     PetscViewer viewer;

 33:     PetscViewerCreate(((PetscObject) dm)->comm, &viewer);
 34:     PetscViewerSetType(viewer, PETSCVIEWERASCII);
 35:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_INFO_DETAIL);
 36:     DMView(dm, viewer);
 37:     PetscViewerDestroy(&viewer);
 38:   }
 39:   PetscOptionsTail();
 40:   return(0);
 41: }

 43: #include <sieve/DMBuilder.hh>

 47: /*
 48:  Simple square boundary:

 50:  18--5-17--4--16
 51:   |     |     |
 52:   6    10     3
 53:   |     |     |
 54:  19-11-20--9--15
 55:   |     |     |
 56:   7     8     2
 57:   |     |     |
 58:  12--0-13--1--14
 59: */
 60: PetscErrorCode DMMeshCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
 61: {
 62:   DM_Mesh        *mesh       = (DM_Mesh*) dm->data;
 63:   PetscInt       numVertices = (edges[0]+1)*(edges[1]+1);
 64:   PetscInt       numEdges    = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
 65:   PetscScalar    *coords;
 66:   PetscInt       coordSize;
 67:   PetscMPIInt    rank;
 68:   PetscInt       v, vx, vy;

 72:   MPI_Comm_rank(((PetscObject) dm)->comm, &rank);
 73:   if (!rank) {
 74:     PetscInt e, ex, ey;

 76:     DMMeshSetChart(dm, 0, numEdges+numVertices);
 77:     for (e = 0; e < numEdges; ++e) {
 78:       DMMeshSetConeSize(dm, e, 2);
 79:     }
 80:     DMMeshSetUp(dm); /* Allocate space for cones */
 81:     for (vy = 0; vy <= edges[1]; vy++) {
 82:       for (ex = 0; ex < edges[0]; ex++) {
 83:         PetscInt edge    = vy*edges[0]     + ex;
 84:         PetscInt vertex  = vy*(edges[0]+1) + ex + numEdges;
 85:         PetscInt cone[2] = {vertex, vertex+1};

 87:         DMMeshSetCone(dm, edge, cone);
 88:         if ((vy == 0) || (vy == edges[1])) {
 89:           DMMeshSetLabelValue(dm, "marker", edge,    1);
 90:           DMMeshSetLabelValue(dm, "marker", cone[0], 1);
 91:           if (ex == edges[0]-1) {
 92:             DMMeshSetLabelValue(dm, "marker", cone[1], 1);
 93:           }
 94:         }
 95:       }
 96:     }
 97:     for (vx = 0; vx <= edges[0]; vx++) {
 98:       for (ey = 0; ey < edges[1]; ey++) {
 99:         PetscInt edge    = vx*edges[1] + ey + edges[0]*(edges[1]+1);
100:         PetscInt vertex  = ey*(edges[0]+1) + vx + numEdges;
101:         PetscInt cone[2] = {vertex, vertex+edges[0]+1};

103:         DMMeshSetCone(dm, edge, cone);
104:         if ((vx == 0) || (vx == edges[0])) {
105:           DMMeshSetLabelValue(dm, "marker", edge,    1);
106:           DMMeshSetLabelValue(dm, "marker", cone[0], 1);
107:           if (ey == edges[1]-1) {
108:             DMMeshSetLabelValue(dm, "marker", cone[1], 1);
109:           }
110:         }
111:       }
112:     }
113:   }
114:   DMMeshSymmetrize(dm);
115:   DMMeshStratify(dm);
116:   /* Build coordinates */
117:   PetscSectionSetChart(mesh->coordSection, numEdges, numEdges + numVertices);
118:   for (v = numEdges; v < numEdges+numVertices; ++v) {
119:     PetscSectionSetDof(mesh->coordSection, v, 2);
120:   }
121:   PetscSectionSetUp(mesh->coordSection);
122:   PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
123:   VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
124:   VecSetFromOptions(mesh->coordinates);
125:   VecGetArray(mesh->coordinates, &coords);
126:   for (vy = 0; vy <= edges[1]; ++vy) {
127:     for (vx = 0; vx <= edges[0]; ++vx) {
128:       coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
129:       coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
130:     }
131:   }
132:   VecRestoreArray(mesh->coordinates, &coords);
133:   return(0);
134: }

138: /*
139:  Simple cubic boundary:

141:      2-------3
142:     /|      /|
143:    6-------7 |
144:    | |     | |
145:    | 0-----|-1
146:    |/      |/
147:    4-------5
148: */
149: PetscErrorCode DMMeshCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
150: {
151:   DM_Mesh        *mesh       = (DM_Mesh*) dm->data;
152:   PetscInt       numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
153:   PetscInt       numFaces    = 6;
154:   PetscScalar    *coords;
155:   PetscInt       coordSize;
156:   PetscMPIInt    rank;
157:   PetscInt       v, vx, vy, vz;

161:   if ((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1)) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Must have at least 1 face per side");
162:   if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Currently can't handle more than 1 face per side");
163:   PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);
164:   MPI_Comm_rank(((PetscObject) dm)->comm, &rank);
165:   if (!rank) {
166:     PetscInt f;

168:     DMMeshSetChart(dm, 0, numFaces+numVertices);
169:     for (f = 0; f < numFaces; ++f) {
170:       DMMeshSetConeSize(dm, f, 4);
171:     }
172:     DMMeshSetUp(dm); /* Allocate space for cones */
173:     for (v = 0; v < numFaces+numVertices; ++v) {
174:       DMMeshSetLabelValue(dm, "marker", v, 1);
175:     }
176:     { // Side 0 (Front)
177:       PetscInt cone[4] = {numFaces+4, numFaces+5, numFaces+7, numFaces+6};
178:       DMMeshSetCone(dm, 0, cone);
179:     }
180:     { // Side 1 (Back)
181:       PetscInt cone[4] = {numFaces+1, numFaces+0, numFaces+2, numFaces+3};
182:       DMMeshSetCone(dm, 0, cone);
183:     }
184:     { // Side 0 (Bottom)
185:       PetscInt cone[4] = {numFaces+0, numFaces+1, numFaces+5, numFaces+4};
186:       DMMeshSetCone(dm, 0, cone);
187:     }
188:     { // Side 0 (Top)
189:       PetscInt cone[4] = {numFaces+6, numFaces+7, numFaces+3, numFaces+2};
190:       DMMeshSetCone(dm, 0, cone);
191:     }
192:     { // Side 0 (Left)
193:       PetscInt cone[4] = {numFaces+0, numFaces+4, numFaces+6, numFaces+2};
194:       DMMeshSetCone(dm, 0, cone);
195:     }
196:     { // Side 0 (Right)
197:       PetscInt cone[4] = {numFaces+5, numFaces+1, numFaces+3, numFaces+7};
198:       DMMeshSetCone(dm, 0, cone);
199:     }
200:   }
201:   DMMeshSymmetrize(dm);
202:   DMMeshStratify(dm);
203:   /* Build coordinates */
204:   PetscSectionSetChart(mesh->coordSection, numFaces, numFaces + numVertices);
205:   for (v = numFaces; v < numFaces+numVertices; ++v) {
206:     PetscSectionSetDof(mesh->coordSection, v, 3);
207:   }
208:   PetscSectionSetUp(mesh->coordSection);
209:   PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
210:   VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
211:   VecGetArray(mesh->coordinates, &coords);
212:   for (vz = 0; vz <= faces[2]; ++vz) {
213:     for (vy = 0; vy <= faces[1]; ++vy) {
214:       for (vx = 0; vx <= faces[0]; ++vx) {
215:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
216:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
217:         coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
218:       }
219:     }
220:   }
221:   VecRestoreArray(mesh->coordinates, &coords);
222:   return(0);
223: }

227: PetscErrorCode DMMeshCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm)
228: {
229:   PetscBool      flg;

234:   PetscOptionsBool("-dm_mesh_new_impl", "Use the new C unstructured mesh implementation", "DMView", PETSC_FALSE, &flg, NULL);
235:   if (flg) {
236:     DM boundary;

238:     if (interpolate) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation (creation of faces and edges) is not yet supported.");
239:     DMCreate(comm, &boundary);
241:     DMSetType(boundary, DMMESH);
242:     DMMeshSetDimension(boundary, dim-1);
243:     switch (dim) {
244:     case 2:
245:     {
246:       PetscReal lower[2] = {0.0, 0.0};
247:       PetscReal upper[2] = {1.0, 1.0};
248:       PetscInt  edges[2] = {2, 2};

250:       DMMeshCreateSquareBoundary(boundary, lower, upper, edges);
251:       break;
252:     }
253:     case 3:
254:     {
255:       PetscReal lower[3] = {0.0, 0.0, 0.0};
256:       PetscReal upper[3] = {1.0, 1.0, 1.0};
257:       PetscInt  faces[3] = {1, 1, 1};

259:       DMMeshCreateCubeBoundary(boundary, lower, upper, faces);
260:       break;
261:     }
262:     default:
263:       SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
264:     }
265:     DMMeshGenerate(boundary, interpolate, dm);
266:     DMDestroy(&boundary);
267:   } else {
268:     PetscInt debug = 0;
269:     try {
270:       ALE::DMBuilder::createBoxMesh(comm, dim, false, interpolate, debug, dm);
271:     } catch(ALE::Exception e) {
272:       SETERRQ1(comm, PETSC_ERR_PLIB, "Unable to create mesh: %s", e.message());
273:     }
274:   }
275:   return(0);
276: }

280: /*@
281:   DMMeshCreateMeshFromAdjacency - Create an unstrctured mesh from a list of the vertices for each cell, and the coordinates for each vertex.

283:  Collective on comm

285:   Input Parameters:
286: + comm - An MPI communicator
287: . dim - The dimension of the cells, e.g. triangles have dimension 2
288: . numCells - The number of cells in the mesh
289: . numCorners - The number of vertices in each cell
290: . cellVertices - An array of the vertices for each cell, numbered 0 to numVertices-1
291: . spatialDim - The dimension for coordinates, e.g. for a triangle in 3D this would be 3
292: . numVertices - The number of mesh vertices
293: . coordinates - An array of the coordinates for each vertex
294: - interpolate - Flag to create faces and edges

296:   Output Parameter:
297: . dm - The DMMesh object

299:   Level: beginner

301: .seealso DMMESH, DMMeshCreateMeshFromAdjacencyHybrid(), DMMeshCreateBoxMesh()
302: @*/
303: PetscErrorCode DMMeshCreateMeshFromAdjacency(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numCorners, PetscInt cellVertices[], PetscInt spatialDim, PetscInt numVertices, const PetscReal coordinates[], PetscBool interpolate, DM *dm)
304: {
305:   PetscInt       *cone;
306:   PetscInt       *coneO;
307:   PetscInt       debug = 0;

314:   if (interpolate) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation (creation of faces and edges) is not yet supported.");
315:   PetscOptionsGetInt(NULL, "-dm_mesh_debug", &debug, NULL);
316:   Obj<PETSC_MESH_TYPE>             mesh  = new PETSC_MESH_TYPE(comm, dim, debug);
317:   Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, 0, numCells+numVertices, debug);

319:   mesh->setSieve(sieve);
320:   for (PetscInt c = 0; c < numCells; ++c) sieve->setConeSize(c, numCorners);
321:   sieve->symmetrizeSizes(numCells, numCorners, cellVertices, numCells);
322:   sieve->allocate();
323:   PetscMalloc2(numCorners,PetscInt,&cone,numCorners,PetscInt,&coneO);
324:   for (PetscInt v = 0; v < numCorners; ++v) coneO[v] = 1;
325:   for (PetscInt c = 0; c < numCells; ++c) {
326:     for (PetscInt v = 0; v < numCorners; ++v) {
327:       cone[v] = cellVertices[c*numCorners+v]+numCells;
328:     }
329:     sieve->setCone(cone, c);
330:     sieve->setConeOrientation(coneO, c);
331:   }
332:   PetscFree2(cone,coneO);
333:   sieve->symmetrize();
334:   mesh->stratify();
335:   ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, spatialDim, coordinates, numCells);
336:   DMCreate(comm, dm);
337:   DMSetType(*dm, DMMESH);
338:   DMMeshSetMesh(*dm, mesh);
339:   return(0);
340: }

344: PetscErrorCode DMMeshCreateMeshFromAdjacencyHybrid(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numCorners[], PetscInt cellVertices[], PetscInt spatialDim, PetscInt numVertices, const PetscReal coordinates[], PetscBool interpolate, DM *dm)
345: {
346:   PetscInt       debug = 0;

354:   if (interpolate) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation (creation of faces and edges) is not yet supported.");
355:   PetscOptionsGetInt(NULL, "-dmmesh_debug", &debug, NULL);
356:   Obj<PETSC_MESH_TYPE>             mesh  = new PETSC_MESH_TYPE(comm, dim, debug);
357:   Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, 0, numCells+numVertices, debug);

359:   mesh->setSieve(sieve);
360:   for (PetscInt c = 0; c < numCells; ++c) {
361:     sieve->setConeSize(c, numCorners[c]);
362:   }
363:   //sieve->symmetrizeSizes(numCells, numCorners, cellVertices);
364:   sieve->allocate();
365:   for (PetscInt c = 0, offset = 0; c < numCells; offset += numCorners[c], ++c) {
366:     sieve->setCone(&cellVertices[offset], c);
367:   }
368:   sieve->symmetrize();
369:   return(0);
370: }

372: /* External function declarations here */
373: extern PetscErrorCode DMGlobalToLocalBegin_Mesh(DM dm, Vec g, InsertMode mode, Vec l);
374: extern PetscErrorCode DMGlobalToLocalEnd_Mesh(DM dm, Vec g, InsertMode mode, Vec l);
375: extern PetscErrorCode DMLocalToGlobalBegin_Mesh(DM dm, Vec l, InsertMode mode, Vec g);
376: extern PetscErrorCode DMLocalToGlobalEnd_Mesh(DM dm, Vec l, InsertMode mode, Vec g);
377: extern PetscErrorCode DMCreateGlobalVector_Mesh(DM dm, Vec *gvec);
378: extern PetscErrorCode DMCreateLocalVector_Mesh(DM dm, Vec *lvec);
379: extern PetscErrorCode DMGetLocalToGlobalMapping_Mesh(DM dm);
380: extern PetscErrorCode DMCreateInterpolation_Mesh(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
381: extern PetscErrorCode DMCreateMatrix_Mesh(DM dm, MatType mtype, Mat *J);
382: extern PetscErrorCode DMRefine_Mesh(DM dm, MPI_Comm comm, DM *dmRefined);
383: extern PetscErrorCode DMCoarsenHierarchy_Mesh(DM dm, int numLevels, DM *coarseHierarchy);
384: extern PetscErrorCode DMDestroy_Mesh(DM dm);
385: extern PetscErrorCode DMView_Mesh(DM dm, PetscViewer viewer);

387: EXTERN_C_BEGIN
390: PetscErrorCode DMConvert_DA_Mesh(DM dm, DMType newtype, DM *dmNew)
391: {
392:   PetscSection   section;
393:   DM             cda;
394:   DMDALocalInfo  info;
395:   Vec            coordinates;
396:   PetscInt       *cone, *coneO;
397:   PetscInt       dim, M, N, P, numCells, numGlobalCells, numCorners, numVertices, c = 0, v = 0;
398:   PetscInt       ye, ze;
399:   PetscInt       debug = 0;

403:   PetscOptionsGetInt(NULL, "-dm_mesh_debug", &debug, NULL);
404:   DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);
405:   DMDAGetLocalInfo(dm, &info);
406:   if (info.sw  > 1) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Currently, only DMDAs with unti stencil width can be converted to DMMeshes.");
407:   /* In order to get a partition of cells, rather than vertices, we give each process the cells between vertices it owns
408:      and also higher numbered ghost vertices (vertices to the right and up) */
409:   numCorners  = 1 << dim;
410:   numCells    = ((info.gxm+info.gxs - info.xs) - 1);
411:   if (dim > 1) numCells *= ((info.gym+info.gys - info.ys) - 1);
412:   if (dim > 2) numCells *= ((info.gzm+info.gzs - info.zs) - 1);
413:   numVertices = (info.gxm+info.gxs - info.xs);
414:   if (dim > 1) numVertices *= (info.gym+info.gys - info.ys);
415:   if (dim > 2) numVertices *= (info.gzm+info.gzs - info.zs);
416:   numGlobalCells = M-1;
417:   if (dim > 1) numGlobalCells *= N-1;
418:   if (dim > 2) numGlobalCells *= P-1;

420:   ALE::Obj<PETSC_MESH_TYPE>             mesh  = new PETSC_MESH_TYPE(((PetscObject) dm)->comm, info.dim, debug);
421:   ALE::Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(((PetscObject) dm)->comm, 0, numCells+numVertices, debug);
422:   PETSC_MESH_TYPE::renumbering_type     renumbering;

424:   mesh->setSieve(sieve);
425:   /* Number each cell for the vertex in the lower left corner */
426:   if (dim < 3) {ze = 1; P = 1;} else ze = info.gzs+info.gzm-1;
427:   if (dim < 2) {ye = 1; N = 1;} else ye = info.gys+info.gym-1;
428:   for (PetscInt k = info.zs; k < ze; ++k) {
429:     for (PetscInt j = info.ys; j < ye; ++j) {
430:       for (PetscInt i = info.xs; i < info.gxs+info.gxm-1; ++i, ++c) {
431:         PetscInt globalC = (k*(N-1) + j)*(M-1) + i;

433:         renumbering[globalC] = c;
434:         sieve->setConeSize(c, numCorners);
435:       }
436:     }
437:   }
438:   if (c != numCells) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Error in generated cell numbering, %d should be %d", c, numCells);
439:   /* Get vertex renumbering */
440:   for (PetscInt k = info.zs; k < info.gzs+info.gzm; ++k) {
441:     for (PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
442:       for (PetscInt i = info.xs; i < info.gxs+info.gxm; ++i, ++v) {
443:         PetscInt globalV = (k*N + j)*M + i + numGlobalCells;

445:         renumbering[globalV] = v+numCells;
446:       }
447:     }
448:   }
449:   if (v != numVertices) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Error in generated vertex numbering, %d should be %d", v, numVertices);
450:   /* Calculate support sizes */
451:   for (PetscInt k = info.zs; k < ze; ++k, ++c) {
452:     for (PetscInt j = info.ys; j < ye; ++j) {
453:       for (PetscInt i = info.xs; i < info.gxs+info.gxm-1; ++i) {
454:         for (PetscInt kp = k; kp <= k+(dim>2); ++kp) {
455:           for (PetscInt jp = j; jp <= j+(dim>1); ++jp) {
456:             for (PetscInt ip = i; ip <= i+1; ++ip) {
457:               PetscInt globalV = (kp*N + jp)*M + ip + numGlobalCells;

459:               sieve->addSupportSize(renumbering[globalV], 1);
460:             }
461:           }
462:         }
463:       }
464:     }
465:   }
466:   sieve->allocate();
467:   PetscMalloc2(numCorners,PetscInt,&cone,numCorners,PetscInt,&coneO);
468:   for (PetscInt v = 0; v < numCorners; ++v) coneO[v] = 1;
469:   for (PetscInt k = info.zs; k < ze; ++k) {
470:     for (PetscInt j = info.ys; j < ye; ++j) {
471:       for (PetscInt i = info.xs; i < info.gxs+info.gxm-1; ++i) {
472:         PetscInt globalC = (k*(N-1) + j)*(M-1) + i;
473:         PetscInt v       = 0;

475:         cone[v++] = renumbering[(k*N + j)*M + i+0 + numGlobalCells];
476:         cone[v++] = renumbering[(k*N + j)*M + i+1 + numGlobalCells];
477:         if (dim > 1) {
478:           cone[v++] = renumbering[(k*N + j+1)*M + i+0 + numGlobalCells];
479:           cone[v++] = renumbering[(k*N + j+1)*M + i+1 + numGlobalCells];
480:         }
481:         if (dim > 2) {
482:           cone[v++] = renumbering[((k+1)*N + j+0)*M + i+0 + numGlobalCells];
483:           cone[v++] = renumbering[((k+1)*N + j+0)*M + i+1 + numGlobalCells];
484:           cone[v++] = renumbering[((k+1)*N + j+1)*M + i+0 + numGlobalCells];
485:           cone[v++] = renumbering[((k+1)*N + j+1)*M + i+1 + numGlobalCells];
486:         }
487:         sieve->setCone(cone, renumbering[globalC]);
488:         sieve->setConeOrientation(coneO, renumbering[globalC]);
489:       }
490:     }
491:   }
492:   PetscFree2(cone,coneO);
493:   sieve->symmetrize();
494:   mesh->stratify();
495:   /* Create boundary marker */
496:   {
497:     const Obj<PETSC_MESH_TYPE::label_type>& boundary = mesh->createLabel("marker");

499:     for (PetscInt k = info.zs; k < info.gzs+info.gzm; ++k) {
500:       for (PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
501:         if (info.xs == 0) {
502:           PetscInt globalV = (k*N + j)*M + info.xs + numGlobalCells;

504:           mesh->setValue(boundary, renumbering[globalV], 1);
505:         }
506:         if (info.gxs+info.gxm-1 == M-1) {
507:           PetscInt globalV = (k*N + j)*M + info.gxs+info.gxm-1 + numGlobalCells;

509:           mesh->setValue(boundary, renumbering[globalV], 1);
510:         }
511:       }
512:     }
513:     if (dim > 1) {
514:       for (PetscInt k = info.zs; k < info.gzs+info.gzm; ++k) {
515:         for (PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
516:           if (info.ys == 0) {
517:             PetscInt globalV = (k*N + info.ys)*M + i + numGlobalCells;

519:             mesh->setValue(boundary, renumbering[globalV], 1);
520:           }
521:           if (info.gys+info.gym-1 == N-1) {
522:             PetscInt globalV = (k*N + info.gys+info.gym-1)*M + i + numGlobalCells;

524:             mesh->setValue(boundary, renumbering[globalV], 1);
525:           }
526:         }
527:       }
528:     }
529:     if (dim > 2) {
530:       for (PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
531:         for (PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
532:           if (info.zs == 0) {
533:             PetscInt globalV = (info.zs*N + j)*M + i + numGlobalCells;

535:             mesh->setValue(boundary, renumbering[globalV], 1);
536:           }
537:           if (info.gzs+info.gzm-1 == P-1) {
538:             PetscInt globalV = ((info.gzs+info.gzm-1)*N + j)*M + i + numGlobalCells;

540:             mesh->setValue(boundary, renumbering[globalV], 1);
541:           }
542:         }
543:       }
544:     }
545:   }
546:   /* Create new DM */
547:   DMMeshCreate(((PetscObject) dm)->comm, dmNew);
548:   DMMeshSetMesh(*dmNew, mesh);
549:   /* Set coordinates */
550:   PetscSectionCreate(((PetscObject) dm)->comm, &section);
551:   PetscSectionSetChart(section, numCells, numCells+numVertices);
552:   for (PetscInt v = numCells; v < numCells+numVertices; ++v) {
553:     PetscSectionSetDof(section, v, dim);
554:   }
555:   PetscSectionSetUp(section);
556:   DMMeshSetCoordinateSection(*dmNew, section);
557:   DMGetCoordinateDM(dm, &cda);
558:   DMGetCoordinatesLocal(dm, &coordinates);
559:   {
560:     Obj<PETSC_MESH_TYPE::real_section_type> coordSection = mesh->getRealSection("coordinates");

562:     switch (dim) {
563:     case 1:
564:     {
565:       PetscScalar **coords;

567:       DMDAVecGetArrayDOF(cda, coordinates, &coords);
568:       for (PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
569:         PetscInt globalV = i + numGlobalCells;

571:         coordSection->updatePoint(renumbering[globalV], coords[i]);
572:       }
573:       DMDAVecRestoreArrayDOF(cda, coordinates, &coords);
574:       break;
575:     }
576:     case 2:
577:     {
578:       PetscScalar ***coords;

580:       DMDAVecGetArrayDOF(cda, coordinates, &coords);
581:       for (PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
582:         for (PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
583:           PetscInt globalV = j*M + i + numGlobalCells;

585:           coordSection->updatePoint(renumbering[globalV], coords[j][i]);
586:         }
587:       }
588:       DMDAVecRestoreArrayDOF(cda, coordinates, &coords);
589:       break;
590:     }
591:     case 3:
592:     {
593:       PetscScalar ****coords;

595:       DMDAVecGetArrayDOF(cda, coordinates, &coords);
596:       for (PetscInt k = info.zs; k < info.gzs+info.gzm; ++k, ++v) {
597:         for (PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
598:           for (PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
599:             PetscInt globalV = (k*N + j)*M + i + numGlobalCells;

601:             coordSection->updatePoint(renumbering[globalV], coords[k][j][i]);
602:           }
603:         }
604:       }
605:       DMDAVecRestoreArrayDOF(cda, coordinates, &coords);
606:       break;
607:     }
608:     default:
609:       SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid DMDA dimension %d", dim);
610:     }
611:   }
612:   /* Get overlap for interdomain communication */
613:   {
614:     typedef PETSC_MESH_TYPE::point_type point_type;
615:     PETSc::Log::Event("CreateOverlap").begin();
616:     ALE::Obj<PETSC_MESH_TYPE::send_overlap_type> sendParallelMeshOverlap = mesh->getSendOverlap();
617:     ALE::Obj<PETSC_MESH_TYPE::recv_overlap_type> recvParallelMeshOverlap = mesh->getRecvOverlap();
618:     //   Can I figure this out in a nicer way?
619:     ALE::SetFromMap<std::map<point_type,point_type> > globalPoints(renumbering);

621:     ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering, sendParallelMeshOverlap, recvParallelMeshOverlap);
622:     if (debug) {
623:       sendParallelMeshOverlap->view("Send Overlap");
624:       recvParallelMeshOverlap->view("Recieve Overlap");
625:     }
626:     mesh->setCalculatedOverlap(true);
627:     PETSc::Log::Event("CreateOverlap").end();
628:   }
629:   return(0);
630: }

634: PetscErrorCode DMCreate_Mesh(DM dm)
635: {
636:   DM_Mesh        *mesh;

641:   PetscNewLog(dm, DM_Mesh, &mesh);
642:   dm->data = mesh;

644:   new(&mesh->m) ALE::Obj<PETSC_MESH_TYPE>(NULL);

646:   mesh->globalScatter  = NULL;
647:   mesh->defaultSection = NULL;
648:   mesh->lf             = NULL;
649:   mesh->lj             = NULL;

651:   mesh->useNewImpl     = PETSC_FALSE;
652:   mesh->dim            = 0;
653:   mesh->sf             = NULL;

655:   PetscSectionCreate(((PetscObject) dm)->comm, &mesh->coneSection);

657:   mesh->maxConeSize    = 0;
658:   mesh->cones          = NULL;

660:   PetscSectionCreate(((PetscObject) dm)->comm, &mesh->supportSection);

662:   mesh->maxSupportSize = 0;
663:   mesh->supports       = NULL;

665:   PetscSectionCreate(((PetscObject) dm)->comm, &mesh->coordSection);
666:   VecCreate(((PetscObject) dm)->comm, &mesh->coordinates);
667:   PetscObjectSetName((PetscObject) mesh->coordinates, "coordinates");

669:   mesh->meetTmpA    = NULL;
670:   mesh->meetTmpB    = NULL;
671:   mesh->joinTmpA    = NULL;
672:   mesh->joinTmpB    = NULL;
673:   mesh->closureTmpA = NULL;
674:   mesh->closureTmpB = NULL;

676:   DMSetVecType(dm,VECSTANDARD);

678:   dm->ops->view                            = DMView_Mesh;
679:   dm->ops->setfromoptions                  = DMSetFromOptions_Mesh;
680:   dm->ops->setup                           = 0;
681:   dm->ops->createglobalvector              = DMCreateGlobalVector_Mesh;
682:   dm->ops->createlocalvector               = DMCreateLocalVector_Mesh;
683:   dm->ops->getlocaltoglobalmapping         = DMGetLocalToGlobalMapping_Mesh;
684:   dm->ops->getlocaltoglobalmappingblock    = 0;

686:   dm->ops->getcoloring         = 0;
687:   dm->ops->creatematrix        = DMCreateMatrix_Mesh;
688:   dm->ops->createinterpolation = DMCreateInterpolation_Mesh;
689:   dm->ops->getaggregates       = 0;
690:   dm->ops->getinjection        = 0;

692:   dm->ops->refine           = DMRefine_Mesh;
693:   dm->ops->coarsen          = 0;
694:   dm->ops->refinehierarchy  = 0;
695:   dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Mesh;

697:   dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Mesh;
698:   dm->ops->globaltolocalend   = DMGlobalToLocalEnd_Mesh;
699:   dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Mesh;
700:   dm->ops->localtoglobalend   = DMLocalToGlobalEnd_Mesh;

702:   dm->ops->destroy = DMDestroy_Mesh;

704:   PetscObjectComposeFunction((PetscObject) dm, "DMConvert_da_mesh_C", DMConvert_DA_Mesh);

706:   /* NEW_MESH_IMPL */
707:   PetscOptionsBool("-dm_mesh_new_impl", "Use the new C unstructured mesh implementation", "DMCreate", PETSC_FALSE, &mesh->useNewImpl, NULL);
708:   return(0);
709: }
710: EXTERN_C_END

714: /*@
715:   DMMeshCreate - Creates a DMMesh object.

717:   Collective on MPI_Comm

719:   Input Parameter:
720: . comm - The communicator for the DMMesh object

722:   Output Parameter:
723: . mesh  - The DMMesh object

725:   Level: beginner

727: .keywords: DMMesh, create
728: @*/
729: PetscErrorCode  DMMeshCreate(MPI_Comm comm, DM *mesh)
730: {

735:   DMCreate(comm, mesh);
736:   DMSetType(*mesh, DMMESH);
737:   return(0);
738: }