Actual source code: meshcreate.c
petsc-3.4.5 2014-06-29
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, §ion);
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: }