Actual source code: meshcreate.c
petsc-3.3-p7 2013-05-11
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, PETSC_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, PETSC_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: PetscBool flg;
233: PetscOptionsBool("-dm_mesh_new_impl", "Use the new C unstructured mesh implementation", "DMView", PETSC_FALSE, &flg, PETSC_NULL);
234: if (flg) {
235: DM boundary;
237: if (interpolate) {SETERRQ(comm, PETSC_ERR_SUP, "Interpolation (creation of faces and edges) is not yet supported.");}
238: DMCreate(comm, &boundary);
240: DMSetType(boundary, DMMESH);
241: DMMeshSetDimension(boundary, dim-1);
242: switch(dim) {
243: case 2:
244: {
245: PetscReal lower[2] = {0.0, 0.0};
246: PetscReal upper[2] = {1.0, 1.0};
247: PetscInt edges[2] = {2, 2};
249: DMMeshCreateSquareBoundary(boundary, lower, upper, edges);
250: break;
251: }
252: case 3:
253: {
254: PetscReal lower[3] = {0.0, 0.0, 0.0};
255: PetscReal upper[3] = {1.0, 1.0, 1.0};
256: PetscInt faces[3] = {1, 1, 1};
258: DMMeshCreateCubeBoundary(boundary, lower, upper, faces);
259: break;
260: }
261: default:
262: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
263: }
264: DMMeshGenerate(boundary, interpolate, dm);
265: DMDestroy(&boundary);
266: } else {
267: PetscInt debug = 0;
268: try {
269: ALE::DMBuilder::createBoxMesh(comm, dim, false, interpolate, debug, dm);
270: } catch(ALE::Exception e) {
271: SETERRQ1(comm, PETSC_ERR_PLIB, "Unable to create mesh: %s", e.message());
272: }
273: }
274: return(0);
275: }
279: /*@
280: DMMeshCreateMeshFromAdjacency - Create an unstrctured mesh from a list of the vertices for each cell, and the coordinates for each vertex.
282: Collective on comm
284: Input Parameters:
285: + comm - An MPI communicator
286: . dim - The dimension of the cells, e.g. triangles have dimension 2
287: . numCells - The number of cells in the mesh
288: . numCorners - The number of vertices in each cell
289: . cellVertices - An array of the vertices for each cell, numbered 0 to numVertices-1
290: . spatialDim - The dimension for coordinates, e.g. for a triangle in 3D this would be 3
291: . numVertices - The number of mesh vertices
292: . coordinates - An array of the coordinates for each vertex
293: - interpolate - Flag to create faces and edges
295: Output Parameter:
296: . dm - The DMMesh object
298: Level: beginner
300: .seealso DMMESH, DMMeshCreateMeshFromAdjacencyHybrid(), DMMeshCreateBoxMesh()
301: @*/
302: PetscErrorCode DMMeshCreateMeshFromAdjacency(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numCorners, PetscInt cellVertices[], PetscInt spatialDim, PetscInt numVertices, const PetscReal coordinates[], PetscBool interpolate, DM *dm) {
303: PetscInt *cone;
304: PetscInt *coneO;
305: PetscInt debug = 0;
312: if (interpolate) {SETERRQ(comm, PETSC_ERR_SUP, "Interpolation (creation of faces and edges) is not yet supported.");}
313: PetscOptionsGetInt(PETSC_NULL, "-dm_mesh_debug", &debug, PETSC_NULL);
314: Obj<PETSC_MESH_TYPE> mesh = new PETSC_MESH_TYPE(comm, dim, debug);
315: Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, 0, numCells+numVertices, debug);
317: mesh->setSieve(sieve);
318: for(PetscInt c = 0; c < numCells; ++c) {
319: sieve->setConeSize(c, numCorners);
320: }
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) {
325: coneO[v] = 1;
326: }
327: for(PetscInt c = 0; c < numCells; ++c) {
328: for(PetscInt v = 0; v < numCorners; ++v) {
329: cone[v] = cellVertices[c*numCorners+v]+numCells;
330: }
331: sieve->setCone(cone, c);
332: sieve->setConeOrientation(coneO, c);
333: }
334: PetscFree2(cone,coneO);
335: sieve->symmetrize();
336: mesh->stratify();
337: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, spatialDim, coordinates, numCells);
338: DMCreate(comm, dm);
339: DMSetType(*dm, DMMESH);
340: DMMeshSetMesh(*dm, mesh);
341: return(0);
342: }
346: PetscErrorCode DMMeshCreateMeshFromAdjacencyHybrid(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numCorners[], PetscInt cellVertices[], PetscInt spatialDim, PetscInt numVertices, const PetscReal coordinates[], PetscBool interpolate, DM *dm) {
347: PetscInt debug = 0;
355: if (interpolate) {SETERRQ(comm, PETSC_ERR_SUP, "Interpolation (creation of faces and edges) is not yet supported.");}
356: PetscOptionsGetInt(PETSC_NULL, "-dmmesh_debug", &debug, PETSC_NULL);
357: Obj<PETSC_MESH_TYPE> mesh = new PETSC_MESH_TYPE(comm, dim, debug);
358: Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, 0, numCells+numVertices, debug);
360: mesh->setSieve(sieve);
361: for(PetscInt c = 0; c < numCells; ++c) {
362: sieve->setConeSize(c, numCorners[c]);
363: }
364: //sieve->symmetrizeSizes(numCells, numCorners, cellVertices);
365: sieve->allocate();
366: for(PetscInt c = 0, offset = 0; c < numCells; offset += numCorners[c], ++c) {
367: sieve->setCone(&cellVertices[offset], c);
368: }
369: sieve->symmetrize();
370: return(0);
371: }
373: /* External function declarations here */
374: extern PetscErrorCode DMGlobalToLocalBegin_Mesh(DM dm, Vec g, InsertMode mode, Vec l);
375: extern PetscErrorCode DMGlobalToLocalEnd_Mesh(DM dm, Vec g, InsertMode mode, Vec l);
376: extern PetscErrorCode DMLocalToGlobalBegin_Mesh(DM dm, Vec l, InsertMode mode, Vec g);
377: extern PetscErrorCode DMLocalToGlobalEnd_Mesh(DM dm, Vec l, InsertMode mode, Vec g);
378: extern PetscErrorCode DMCreateGlobalVector_Mesh(DM dm, Vec *gvec);
379: extern PetscErrorCode DMCreateLocalVector_Mesh(DM dm, Vec *lvec);
380: extern PetscErrorCode DMCreateLocalToGlobalMapping_Mesh(DM dm);
381: extern PetscErrorCode DMCreateInterpolation_Mesh(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
382: extern PetscErrorCode DMCreateMatrix_Mesh(DM dm, const MatType mtype, Mat *J);
383: extern PetscErrorCode DMRefine_Mesh(DM dm, MPI_Comm comm, DM *dmRefined);
384: extern PetscErrorCode DMCoarsenHierarchy_Mesh(DM dm, int numLevels, DM *coarseHierarchy);
385: extern PetscErrorCode DMDestroy_Mesh(DM dm);
386: extern PetscErrorCode DMView_Mesh(DM dm, PetscViewer viewer);
388: EXTERN_C_BEGIN
391: PetscErrorCode DMConvert_DA_Mesh(DM dm, const DMType newtype, DM *dmNew)
392: {
393: PetscSection section;
394: DM cda;
395: DMDALocalInfo info;
396: Vec coordinates;
397: PetscInt *cone, *coneO;
398: PetscInt dim, M, N, P, numCells, numGlobalCells, numCorners, numVertices, c = 0, v = 0;
399: PetscInt ye, ze;
400: PetscInt debug = 0;
404: PetscOptionsGetInt(PETSC_NULL, "-dm_mesh_debug", &debug, PETSC_NULL);
405: DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);
406: DMDAGetLocalInfo(dm, &info);
407: if (info.sw > 1) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Currently, only DMDAs with unti stencil width can be converted to DMMeshes.");
408: /* In order to get a partition of cells, rather than vertices, we give each process the cells between vertices it owns
409: and also higher numbered ghost vertices (vertices to the right and up) */
410: numCorners = 1 << dim;
411: numCells = ((info.gxm+info.gxs - info.xs) - 1);
412: if (dim > 1) {numCells *= ((info.gym+info.gys - info.ys) - 1);}
413: if (dim > 2) {numCells *= ((info.gzm+info.gzs - info.zs) - 1);}
414: numVertices = (info.gxm+info.gxs - info.xs);
415: if (dim > 1) {numVertices *= (info.gym+info.gys - info.ys);}
416: if (dim > 2) {numVertices *= (info.gzm+info.gzs - info.zs);}
417: numGlobalCells = M-1;
418: if (dim > 1) {numGlobalCells *= N-1;}
419: if (dim > 2) {numGlobalCells *= P-1;}
421: ALE::Obj<PETSC_MESH_TYPE> mesh = new PETSC_MESH_TYPE(((PetscObject) dm)->comm, info.dim, debug);
422: ALE::Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(((PetscObject) dm)->comm, 0, numCells+numVertices, debug);
423: PETSC_MESH_TYPE::renumbering_type renumbering;
425: mesh->setSieve(sieve);
426: /* Number each cell for the vertex in the lower left corner */
427: if (dim < 3) {ze = 1; P = 1;} else {ze = info.gzs+info.gzm-1;}
428: if (dim < 2) {ye = 1; N = 1;} else {ye = info.gys+info.gym-1;}
429: for(PetscInt k = info.zs; k < ze; ++k) {
430: for(PetscInt j = info.ys; j < ye; ++j) {
431: for(PetscInt i = info.xs; i < info.gxs+info.gxm-1; ++i, ++c) {
432: PetscInt globalC = (k*(N-1) + j)*(M-1) + i;
434: renumbering[globalC] = c;
435: sieve->setConeSize(c, numCorners);
436: }
437: }
438: }
439: if (c != numCells) {SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Error in generated cell numbering, %d should be %d", c, numCells);}
440: /* Get vertex renumbering */
441: for(PetscInt k = info.zs; k < info.gzs+info.gzm; ++k) {
442: for(PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
443: for(PetscInt i = info.xs; i < info.gxs+info.gxm; ++i, ++v) {
444: PetscInt globalV = (k*N + j)*M + i + numGlobalCells;
446: renumbering[globalV] = v+numCells;
447: }
448: }
449: }
450: if (v != numVertices) {SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Error in generated vertex numbering, %d should be %d", v, numVertices);}
451: /* Calculate support sizes */
452: for(PetscInt k = info.zs; k < ze; ++k, ++c) {
453: for(PetscInt j = info.ys; j < ye; ++j) {
454: for(PetscInt i = info.xs; i < info.gxs+info.gxm-1; ++i) {
455: for(PetscInt kp = k; kp <= k+(dim>2); ++kp) {
456: for(PetscInt jp = j; jp <= j+(dim>1); ++jp) {
457: for(PetscInt ip = i; ip <= i+1; ++ip) {
458: PetscInt globalV = (kp*N + jp)*M + ip + numGlobalCells;
460: sieve->addSupportSize(renumbering[globalV], 1);
461: }
462: }
463: }
464: }
465: }
466: }
467: sieve->allocate();
468: PetscMalloc2(numCorners,PetscInt,&cone,numCorners,PetscInt,&coneO);
469: for(PetscInt v = 0; v < numCorners; ++v) {
470: coneO[v] = 1;
471: }
472: for(PetscInt k = info.zs; k < ze; ++k) {
473: for(PetscInt j = info.ys; j < ye; ++j) {
474: for(PetscInt i = info.xs; i < info.gxs+info.gxm-1; ++i) {
475: PetscInt globalC = (k*(N-1) + j)*(M-1) + i;
476: PetscInt v = 0;
478: cone[v++] = renumbering[(k*N + j)*M + i+0 + numGlobalCells];
479: cone[v++] = renumbering[(k*N + j)*M + i+1 + numGlobalCells];
480: if (dim > 1) {
481: cone[v++] = renumbering[(k*N + j+1)*M + i+0 + numGlobalCells];
482: cone[v++] = renumbering[(k*N + j+1)*M + i+1 + numGlobalCells];
483: }
484: if (dim > 2) {
485: cone[v++] = renumbering[((k+1)*N + j+0)*M + i+0 + numGlobalCells];
486: cone[v++] = renumbering[((k+1)*N + j+0)*M + i+1 + numGlobalCells];
487: cone[v++] = renumbering[((k+1)*N + j+1)*M + i+0 + numGlobalCells];
488: cone[v++] = renumbering[((k+1)*N + j+1)*M + i+1 + numGlobalCells];
489: }
490: sieve->setCone(cone, renumbering[globalC]);
491: sieve->setConeOrientation(coneO, renumbering[globalC]);
492: }
493: }
494: }
495: PetscFree2(cone,coneO);
496: sieve->symmetrize();
497: mesh->stratify();
498: /* Create boundary marker */
499: {
500: const Obj<PETSC_MESH_TYPE::label_type>& boundary = mesh->createLabel("marker");
502: for(PetscInt k = info.zs; k < info.gzs+info.gzm; ++k) {
503: for(PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
504: if (info.xs == 0) {
505: PetscInt globalV = (k*N + j)*M + info.xs + numGlobalCells;
507: mesh->setValue(boundary, renumbering[globalV], 1);
508: }
509: if (info.gxs+info.gxm-1 == M-1) {
510: PetscInt globalV = (k*N + j)*M + info.gxs+info.gxm-1 + numGlobalCells;
512: mesh->setValue(boundary, renumbering[globalV], 1);
513: }
514: }
515: }
516: if (dim > 1) {
517: for(PetscInt k = info.zs; k < info.gzs+info.gzm; ++k) {
518: for(PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
519: if (info.ys == 0) {
520: PetscInt globalV = (k*N + info.ys)*M + i + numGlobalCells;
522: mesh->setValue(boundary, renumbering[globalV], 1);
523: }
524: if (info.gys+info.gym-1 == N-1) {
525: PetscInt globalV = (k*N + info.gys+info.gym-1)*M + i + numGlobalCells;
527: mesh->setValue(boundary, renumbering[globalV], 1);
528: }
529: }
530: }
531: }
532: if (dim > 2) {
533: for(PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
534: for(PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
535: if (info.zs == 0) {
536: PetscInt globalV = (info.zs*N + j)*M + i + numGlobalCells;
538: mesh->setValue(boundary, renumbering[globalV], 1);
539: }
540: if (info.gzs+info.gzm-1 == P-1) {
541: PetscInt globalV = ((info.gzs+info.gzm-1)*N + j)*M + i + numGlobalCells;
543: mesh->setValue(boundary, renumbering[globalV], 1);
544: }
545: }
546: }
547: }
548: }
549: /* Create new DM */
550: DMMeshCreate(((PetscObject) dm)->comm, dmNew);
551: DMMeshSetMesh(*dmNew, mesh);
552: /* Set coordinates */
553: PetscSectionCreate(((PetscObject) dm)->comm, §ion);
554: PetscSectionSetChart(section, numCells, numCells+numVertices);
555: for(PetscInt v = numCells; v < numCells+numVertices; ++v) {
556: PetscSectionSetDof(section, v, dim);
557: }
558: PetscSectionSetUp(section);
559: DMMeshSetCoordinateSection(*dmNew, section);
560: DMDAGetCoordinateDA(dm, &cda);
561: DMDAGetGhostedCoordinates(dm, &coordinates);
562: {
563: Obj<PETSC_MESH_TYPE::real_section_type> coordSection = mesh->getRealSection("coordinates");
565: switch(dim) {
566: case 1:
567: {
568: PetscScalar **coords;
570: DMDAVecGetArrayDOF(cda, coordinates, &coords);
571: for(PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
572: PetscInt globalV = i + numGlobalCells;
574: coordSection->updatePoint(renumbering[globalV], coords[i]);
575: }
576: DMDAVecRestoreArrayDOF(cda, coordinates, &coords);
577: break;
578: }
579: case 2:
580: {
581: PetscScalar ***coords;
583: DMDAVecGetArrayDOF(cda, coordinates, &coords);
584: for(PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
585: for(PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
586: PetscInt globalV = j*M + i + numGlobalCells;
588: coordSection->updatePoint(renumbering[globalV], coords[j][i]);
589: }
590: }
591: DMDAVecRestoreArrayDOF(cda, coordinates, &coords);
592: break;
593: }
594: case 3:
595: {
596: PetscScalar ****coords;
598: DMDAVecGetArrayDOF(cda, coordinates, &coords);
599: for(PetscInt k = info.zs; k < info.gzs+info.gzm; ++k, ++v) {
600: for(PetscInt j = info.ys; j < info.gys+info.gym; ++j) {
601: for(PetscInt i = info.xs; i < info.gxs+info.gxm; ++i) {
602: PetscInt globalV = (k*N + j)*M + i + numGlobalCells;
604: coordSection->updatePoint(renumbering[globalV], coords[k][j][i]);
605: }
606: }
607: }
608: DMDAVecRestoreArrayDOF(cda, coordinates, &coords);
609: break;
610: }
611: default:
612: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid DMDA dimension %d", dim);
613: }
614: }
615: /* Get overlap for interdomain communication */
616: {
617: typedef PETSC_MESH_TYPE::point_type point_type;
618: PETSc::Log::Event("CreateOverlap").begin();
619: ALE::Obj<PETSC_MESH_TYPE::send_overlap_type> sendParallelMeshOverlap = mesh->getSendOverlap();
620: ALE::Obj<PETSC_MESH_TYPE::recv_overlap_type> recvParallelMeshOverlap = mesh->getRecvOverlap();
621: // Can I figure this out in a nicer way?
622: ALE::SetFromMap<std::map<point_type,point_type> > globalPoints(renumbering);
624: ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering, sendParallelMeshOverlap, recvParallelMeshOverlap);
625: if (debug) {
626: sendParallelMeshOverlap->view("Send Overlap");
627: recvParallelMeshOverlap->view("Recieve Overlap");
628: }
629: mesh->setCalculatedOverlap(true);
630: PETSc::Log::Event("CreateOverlap").end();
631: }
632: return(0);
633: }
637: PetscErrorCode DMCreate_Mesh(DM dm)
638: {
639: DM_Mesh *mesh;
644: PetscNewLog(dm, DM_Mesh, &mesh);
645: dm->data = mesh;
647: new(&mesh->m) ALE::Obj<PETSC_MESH_TYPE>(PETSC_NULL);
649: mesh->globalScatter = PETSC_NULL;
650: mesh->defaultSection = PETSC_NULL;
651: mesh->lf = PETSC_NULL;
652: mesh->lj = PETSC_NULL;
654: mesh->useNewImpl = PETSC_FALSE;
655: mesh->dim = 0;
656: mesh->sf = PETSC_NULL;
657: PetscSectionCreate(((PetscObject) dm)->comm, &mesh->coneSection);
658: mesh->maxConeSize = 0;
659: mesh->cones = PETSC_NULL;
660: PetscSectionCreate(((PetscObject) dm)->comm, &mesh->supportSection);
661: mesh->maxSupportSize = 0;
662: mesh->supports = PETSC_NULL;
663: PetscSectionCreate(((PetscObject) dm)->comm, &mesh->coordSection);
664: VecCreate(((PetscObject) dm)->comm, &mesh->coordinates);
665: PetscObjectSetName((PetscObject) mesh->coordinates, "coordinates");
667: mesh->meetTmpA = PETSC_NULL;
668: mesh->meetTmpB = PETSC_NULL;
669: mesh->joinTmpA = PETSC_NULL;
670: mesh->joinTmpB = PETSC_NULL;
671: mesh->closureTmpA = PETSC_NULL;
672: mesh->closureTmpB = PETSC_NULL;
674: PetscStrallocpy(VECSTANDARD, &dm->vectype);
675: dm->ops->view = DMView_Mesh;
676: dm->ops->setfromoptions = DMSetFromOptions_Mesh;
677: dm->ops->setup = 0;
678: dm->ops->createglobalvector = DMCreateGlobalVector_Mesh;
679: dm->ops->createlocalvector = DMCreateLocalVector_Mesh;
680: dm->ops->createlocaltoglobalmapping = DMCreateLocalToGlobalMapping_Mesh;
681: dm->ops->createlocaltoglobalmappingblock = 0;
683: dm->ops->getcoloring = 0;
684: dm->ops->creatematrix = DMCreateMatrix_Mesh;
685: dm->ops->createinterpolation = DMCreateInterpolation_Mesh;
686: dm->ops->getaggregates = 0;
687: dm->ops->getinjection = 0;
689: dm->ops->refine = DMRefine_Mesh;
690: dm->ops->coarsen = 0;
691: dm->ops->refinehierarchy = 0;
692: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Mesh;
694: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Mesh;
695: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Mesh;
696: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Mesh;
697: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Mesh;
699: dm->ops->destroy = DMDestroy_Mesh;
701: PetscObjectComposeFunction((PetscObject) dm, "DMConvert_da_mesh_C", "DMConvert_DA_Mesh", (void (*)(void)) DMConvert_DA_Mesh);
703: /* NEW_MESH_IMPL */
704: PetscOptionsBool("-dm_mesh_new_impl", "Use the new C unstructured mesh implementation", "DMCreate", PETSC_FALSE, &mesh->useNewImpl, PETSC_NULL);
705: return(0);
706: }
707: EXTERN_C_END
711: /*@
712: DMMeshCreate - Creates a DMMesh object.
714: Collective on MPI_Comm
716: Input Parameter:
717: . comm - The communicator for the DMMesh object
719: Output Parameter:
720: . mesh - The DMMesh object
722: Level: beginner
724: .keywords: DMMesh, create
725: @*/
726: PetscErrorCode DMMeshCreate(MPI_Comm comm, DM *mesh)
727: {
732: DMCreate(comm, mesh);
733: DMSetType(*mesh, DMMESH);
734: return(0);
735: }