Actual source code: plexcreate.c
petsc-3.5.4 2015-05-23
1: #define PETSCDM_DLL
2: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
3: #include <petscdmda.h>
4: #include <petscsf.h>
8: /*@
9: DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement.
11: Collective on MPI_Comm
13: Input Parameters:
14: + comm - The communicator for the DM object
15: . dim - The spatial dimension
16: . simplex - Flag for simplicial cells, otherwise they are tensor product cells
17: . interpolate - Flag to create intermediate mesh pieces (edges, faces)
18: . refinementUniform - Flag for uniform parallel refinement
19: - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell
21: Output Parameter:
22: . dm - The DM object
24: Level: beginner
26: .keywords: DM, create
27: .seealso: DMSetType(), DMCreate()
28: @*/
29: PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm)
30: {
31: DM dm;
32: PetscInt p;
33: PetscMPIInt rank;
37: DMCreate(comm, &dm);
38: DMSetType(dm, DMPLEX);
39: DMPlexSetDimension(dm, dim);
40: MPI_Comm_rank(comm, &rank);
41: switch (dim) {
42: case 2:
43: if (simplex) {PetscObjectSetName((PetscObject) dm, "triangular");}
44: else {PetscObjectSetName((PetscObject) dm, "quadrilateral");}
45: break;
46: case 3:
47: if (simplex) {PetscObjectSetName((PetscObject) dm, "tetrahedral");}
48: else {PetscObjectSetName((PetscObject) dm, "hexahedral");}
49: break;
50: default:
51: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
52: }
53: if (rank) {
54: PetscInt numPoints[2] = {0, 0};
55: DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);
56: } else {
57: switch (dim) {
58: case 2:
59: if (simplex) {
60: PetscInt numPoints[2] = {4, 2};
61: PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0};
62: PetscInt cones[6] = {2, 3, 4, 5, 4, 3};
63: PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
64: PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
65: PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1};
67: DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
68: for (p = 0; p < 4; ++p) {DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);}
69: } else {
70: PetscInt numPoints[2] = {6, 2};
71: PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0};
72: PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4};
73: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
74: PetscScalar vertexCoords[12] = {-1.0, -0.5, 0.0, -0.5, 0.0, 0.5, -1.0, 0.5, 1.0, -0.5, 1.0, 0.5};
76: DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
77: }
78: break;
79: case 3:
80: if (simplex) {
81: PetscInt numPoints[2] = {5, 2};
82: PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0};
83: PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6};
84: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
85: PetscScalar vertexCoords[15] = {-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0};
86: PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
88: DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
89: for (p = 0; p < 5; ++p) {DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);}
90: } else {
91: PetscInt numPoints[2] = {12, 2};
92: PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
93: PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8};
94: PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
95: PetscScalar vertexCoords[36] = {-1.0, -0.5, -0.5, -1.0, 0.5, -0.5, 0.0, 0.5, -0.5, 0.0, -0.5, -0.5,
96: -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5,
97: 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5};
99: DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
100: }
101: break;
102: default:
103: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
104: }
105: }
106: *newdm = dm;
107: if (refinementLimit > 0.0) {
108: DM rdm;
109: const char *name;
111: DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);
112: DMPlexSetRefinementLimit(*newdm, refinementLimit);
113: DMRefine(*newdm, comm, &rdm);
114: PetscObjectGetName((PetscObject) *newdm, &name);
115: PetscObjectSetName((PetscObject) rdm, name);
116: DMDestroy(newdm);
117: *newdm = rdm;
118: }
119: if (interpolate) {
120: DM idm = NULL;
121: const char *name;
123: DMPlexInterpolate(*newdm, &idm);
124: PetscObjectGetName((PetscObject) *newdm, &name);
125: PetscObjectSetName((PetscObject) idm, name);
126: DMPlexCopyCoordinates(*newdm, idm);
127: DMPlexCopyLabels(*newdm, idm);
128: DMDestroy(newdm);
129: *newdm = idm;
130: }
131: {
132: DM refinedMesh = NULL;
133: DM distributedMesh = NULL;
135: /* Distribute mesh over processes */
136: DMPlexDistribute(*newdm, NULL, 0, NULL, &distributedMesh);
137: if (distributedMesh) {
138: DMDestroy(newdm);
139: *newdm = distributedMesh;
140: }
141: if (refinementUniform) {
142: DMPlexSetRefinementUniform(*newdm, refinementUniform);
143: DMRefine(*newdm, comm, &refinedMesh);
144: if (refinedMesh) {
145: DMDestroy(newdm);
146: *newdm = refinedMesh;
147: }
148: }
149: }
150: return(0);
151: }
155: /*@
156: DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice.
158: Collective on MPI_Comm
160: Input Parameters:
161: + comm - The communicator for the DM object
162: . lower - The lower left corner coordinates
163: . upper - The upper right corner coordinates
164: - edges - The number of cells in each direction
166: Output Parameter:
167: . dm - The DM object
169: Note: Here is the numbering returned for 2 cells in each direction:
170: $ 18--5-17--4--16
171: $ | | |
172: $ 6 10 3
173: $ | | |
174: $ 19-11-20--9--15
175: $ | | |
176: $ 7 8 2
177: $ | | |
178: $ 12--0-13--1--14
180: Level: beginner
182: .keywords: DM, create
183: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate()
184: @*/
185: PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
186: {
187: PetscInt numVertices = (edges[0]+1)*(edges[1]+1);
188: PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
189: PetscInt markerTop = 1;
190: PetscInt markerBottom = 1;
191: PetscInt markerRight = 1;
192: PetscInt markerLeft = 1;
193: PetscBool markerSeparate = PETSC_FALSE;
194: Vec coordinates;
195: PetscSection coordSection;
196: PetscScalar *coords;
197: PetscInt coordSize;
198: PetscMPIInt rank;
199: PetscInt v, vx, vy;
203: PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
204: if (markerSeparate) {
205: markerTop = 1;
206: markerBottom = 0;
207: markerRight = 0;
208: markerLeft = 0;
209: }
210: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
211: if (!rank) {
212: PetscInt e, ex, ey;
214: DMPlexSetChart(dm, 0, numEdges+numVertices);
215: for (e = 0; e < numEdges; ++e) {
216: DMPlexSetConeSize(dm, e, 2);
217: }
218: DMSetUp(dm); /* Allocate space for cones */
219: for (vx = 0; vx <= edges[0]; vx++) {
220: for (ey = 0; ey < edges[1]; ey++) {
221: PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1);
222: PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
223: PetscInt cone[2];
225: cone[0] = vertex; cone[1] = vertex+edges[0]+1;
226: DMPlexSetCone(dm, edge, cone);
227: if (vx == edges[0]) {
228: DMPlexSetLabelValue(dm, "marker", edge, markerRight);
229: DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
230: if (ey == edges[1]-1) {
231: DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
232: }
233: } else if (vx == 0) {
234: DMPlexSetLabelValue(dm, "marker", edge, markerLeft);
235: DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
236: if (ey == edges[1]-1) {
237: DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
238: }
239: }
240: }
241: }
242: for (vy = 0; vy <= edges[1]; vy++) {
243: for (ex = 0; ex < edges[0]; ex++) {
244: PetscInt edge = vy*edges[0] + ex;
245: PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
246: PetscInt cone[2];
248: cone[0] = vertex; cone[1] = vertex+1;
249: DMPlexSetCone(dm, edge, cone);
250: if (vy == edges[1]) {
251: DMPlexSetLabelValue(dm, "marker", edge, markerTop);
252: DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
253: if (ex == edges[0]-1) {
254: DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
255: }
256: } else if (vy == 0) {
257: DMPlexSetLabelValue(dm, "marker", edge, markerBottom);
258: DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
259: if (ex == edges[0]-1) {
260: DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
261: }
262: }
263: }
264: }
265: }
266: DMPlexSymmetrize(dm);
267: DMPlexStratify(dm);
268: /* Build coordinates */
269: DMGetCoordinateSection(dm, &coordSection);
270: PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);
271: for (v = numEdges; v < numEdges+numVertices; ++v) {
272: PetscSectionSetDof(coordSection, v, 2);
273: }
274: PetscSectionSetUp(coordSection);
275: PetscSectionGetStorageSize(coordSection, &coordSize);
276: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
277: VecSetBlockSize(coordinates, 2);
278: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
279: VecSetType(coordinates,VECSTANDARD);
280: VecGetArray(coordinates, &coords);
281: for (vy = 0; vy <= edges[1]; ++vy) {
282: for (vx = 0; vx <= edges[0]; ++vx) {
283: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
284: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
285: }
286: }
287: VecRestoreArray(coordinates, &coords);
288: DMSetCoordinatesLocal(dm, coordinates);
289: VecDestroy(&coordinates);
290: return(0);
291: }
295: /*@
296: DMPlexCreateCubeBoundary - Creates a 2D mesh the is the boundary of a cubic lattice.
298: Collective on MPI_Comm
300: Input Parameters:
301: + comm - The communicator for the DM object
302: . lower - The lower left front corner coordinates
303: . upper - The upper right back corner coordinates
304: - edges - The number of cells in each direction
306: Output Parameter:
307: . dm - The DM object
309: Level: beginner
311: .keywords: DM, create
312: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate()
313: @*/
314: PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
315: {
316: PetscInt vertices[3], numVertices;
317: PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2];
318: Vec coordinates;
319: PetscSection coordSection;
320: PetscScalar *coords;
321: PetscInt coordSize;
322: PetscMPIInt rank;
323: PetscInt v, vx, vy, vz;
324: PetscInt voffset, iface=0, cone[4];
328: 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");
329: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
330: vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1;
331: numVertices = vertices[0]*vertices[1]*vertices[2];
332: if (!rank) {
333: PetscInt f;
335: DMPlexSetChart(dm, 0, numFaces+numVertices);
336: for (f = 0; f < numFaces; ++f) {
337: DMPlexSetConeSize(dm, f, 4);
338: }
339: DMSetUp(dm); /* Allocate space for cones */
340: for (v = 0; v < numFaces+numVertices; ++v) {
341: DMPlexSetLabelValue(dm, "marker", v, 1);
342: }
344: /* Side 0 (Top) */
345: for (vy = 0; vy < faces[1]; vy++) {
346: for (vx = 0; vx < faces[0]; vx++) {
347: voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx;
348: cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0];
349: DMPlexSetCone(dm, iface, cone);
350: iface++;
351: }
352: }
354: /* Side 1 (Bottom) */
355: for (vy = 0; vy < faces[1]; vy++) {
356: for (vx = 0; vx < faces[0]; vx++) {
357: voffset = numFaces + vy*(faces[0]+1) + vx;
358: cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1;
359: DMPlexSetCone(dm, iface, cone);
360: iface++;
361: }
362: }
364: /* Side 2 (Front) */
365: for (vz = 0; vz < faces[2]; vz++) {
366: for (vx = 0; vx < faces[0]; vx++) {
367: voffset = numFaces + vz*vertices[0]*vertices[1] + vx;
368: cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1];
369: DMPlexSetCone(dm, iface, cone);
370: iface++;
371: }
372: }
374: /* Side 3 (Back) */
375: for (vz = 0; vz < faces[2]; vz++) {
376: for (vx = 0; vx < faces[0]; vx++) {
377: voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx;
378: cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1;
379: cone[2] = voffset+1; cone[3] = voffset;
380: DMPlexSetCone(dm, iface, cone);
381: iface++;
382: }
383: }
385: /* Side 4 (Left) */
386: for (vz = 0; vz < faces[2]; vz++) {
387: for (vy = 0; vy < faces[1]; vy++) {
388: voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0];
389: cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1];
390: cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0];
391: DMPlexSetCone(dm, iface, cone);
392: iface++;
393: }
394: }
396: /* Side 5 (Right) */
397: for (vz = 0; vz < faces[2]; vz++) {
398: for (vy = 0; vy < faces[1]; vy++) {
399: voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + vx;
400: cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset;
401: cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0];
402: DMPlexSetCone(dm, iface, cone);
403: iface++;
404: }
405: }
406: }
407: DMPlexSymmetrize(dm);
408: DMPlexStratify(dm);
409: /* Build coordinates */
410: DMGetCoordinateSection(dm, &coordSection);
411: PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);
412: for (v = numFaces; v < numFaces+numVertices; ++v) {
413: PetscSectionSetDof(coordSection, v, 3);
414: }
415: PetscSectionSetUp(coordSection);
416: PetscSectionGetStorageSize(coordSection, &coordSize);
417: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
418: VecSetBlockSize(coordinates, 3);
419: PetscObjectSetName((PetscObject) coordinates, "coordinates");
420: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
421: VecSetType(coordinates,VECSTANDARD);
422: VecGetArray(coordinates, &coords);
423: for (vz = 0; vz <= faces[2]; ++vz) {
424: for (vy = 0; vy <= faces[1]; ++vy) {
425: for (vx = 0; vx <= faces[0]; ++vx) {
426: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
427: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
428: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
429: }
430: }
431: }
432: VecRestoreArray(coordinates, &coords);
433: DMSetCoordinatesLocal(dm, coordinates);
434: VecDestroy(&coordinates);
435: return(0);
436: }
440: /*@
441: DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice.
443: Collective on MPI_Comm
445: Input Parameters:
446: + comm - The communicator for the DM object
447: . lower - The lower left corner coordinates
448: . upper - The upper right corner coordinates
449: . edges - The number of cells in each direction
450: . bdX - The boundary type for the X direction
451: - bdY - The boundary type for the Y direction
453: Output Parameter:
454: . dm - The DM object
456: Note: Here is the numbering returned for 2 cells in each direction:
457: $ 22--8-23--9--24
458: $ | | |
459: $ 13 2 14 3 15
460: $ | | |
461: $ 19--6-20--7--21
462: $ | | |
463: $ 10 0 11 1 12
464: $ | | |
465: $ 16--4-17--5--18
467: Level: beginner
469: .keywords: DM, create
470: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate()
471: @*/
472: PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY)
473: {
474: PetscInt markerTop = 1;
475: PetscInt markerBottom = 1;
476: PetscInt markerRight = 1;
477: PetscInt markerLeft = 1;
478: PetscBool markerSeparate = PETSC_FALSE;
479: PetscMPIInt rank;
483: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
484: PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
485: if (markerSeparate) {
486: markerTop = 3;
487: markerBottom = 1;
488: markerRight = 2;
489: markerLeft = 4;
490: }
491: {
492: const PetscInt numXEdges = !rank ? edges[0] : 0;
493: const PetscInt numYEdges = !rank ? edges[1] : 0;
494: const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0;
495: const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0;
496: const PetscInt numTotXEdges = numXEdges*numYVertices;
497: const PetscInt numTotYEdges = numYEdges*numXVertices;
498: const PetscInt numVertices = numXVertices*numYVertices;
499: const PetscInt numEdges = numTotXEdges + numTotYEdges;
500: const PetscInt numFaces = numXEdges*numYEdges;
501: const PetscInt firstVertex = numFaces;
502: const PetscInt firstXEdge = numFaces + numVertices;
503: const PetscInt firstYEdge = numFaces + numVertices + numTotXEdges;
504: Vec coordinates;
505: PetscSection coordSection;
506: PetscScalar *coords;
507: PetscInt coordSize;
508: PetscInt v, vx, vy;
509: PetscInt f, fx, fy, e, ex, ey;
511: DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);
512: for (f = 0; f < numFaces; ++f) {
513: DMPlexSetConeSize(dm, f, 4);
514: }
515: for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
516: DMPlexSetConeSize(dm, e, 2);
517: }
518: DMSetUp(dm); /* Allocate space for cones */
519: /* Build faces */
520: for (fy = 0; fy < numYEdges; ++fy) {
521: for (fx = 0; fx < numXEdges; ++fx) {
522: PetscInt face = fy*numXEdges + fx;
523: PetscInt edgeL = firstYEdge + fx *numYEdges + fy;
524: PetscInt edgeR = firstYEdge + ((fx+1)%numXVertices)*numYEdges + fy;
525: PetscInt edgeB = firstXEdge + fy *numXEdges + fx;
526: PetscInt edgeT = firstXEdge + ((fy+1)%numYVertices)*numXEdges + fx;
527: PetscInt ornt[4] = {0, 0, -2, -2};
528: PetscInt cone[4];
530: if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;}
531: if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;}
532: cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
533: DMPlexSetCone(dm, face, cone);
534: DMPlexSetConeOrientation(dm, face, ornt);
535: }
536: }
537: /* Build Y edges*/
538: for (vx = 0; vx < numXVertices; vx++) {
539: for (ey = 0; ey < numYEdges; ey++) {
540: const PetscInt nextv = (bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : ((ey+1)%numYVertices)*numXVertices + vx;
541: const PetscInt edge = firstYEdge + vx*numYEdges + ey;
542: const PetscInt vertexB = firstVertex + ey*numXVertices + vx;
543: const PetscInt vertexT = firstVertex + nextv;
544: PetscInt cone[2];
546: cone[0] = vertexB; cone[1] = vertexT;
547: DMPlexSetCone(dm, edge, cone);
548: if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) {
549: if (vx == numXVertices-1) {
550: DMPlexSetLabelValue(dm, "marker", edge, markerRight);
551: DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
552: if (ey == numYEdges-1) {
553: DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
554: }
555: } else if (vx == 0) {
556: DMPlexSetLabelValue(dm, "marker", edge, markerLeft);
557: DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
558: if (ey == numYEdges-1) {
559: DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
560: }
561: }
562: }
563: }
564: }
565: /* Build X edges*/
566: for (vy = 0; vy < numYVertices; vy++) {
567: for (ex = 0; ex < numXEdges; ex++) {
568: const PetscInt nextv = (bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : vy*numXVertices + (ex+1)%numXVertices;
569: const PetscInt edge = firstXEdge + vy*numXEdges + ex;
570: const PetscInt vertexL = firstVertex + vy*numXVertices + ex;
571: const PetscInt vertexR = firstVertex + nextv;
572: PetscInt cone[2];
574: cone[0] = vertexL; cone[1] = vertexR;
575: DMPlexSetCone(dm, edge, cone);
576: if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) {
577: if (vy == numYVertices-1) {
578: DMPlexSetLabelValue(dm, "marker", edge, markerTop);
579: DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
580: if (ex == numXEdges-1) {
581: DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
582: }
583: } else if (vy == 0) {
584: DMPlexSetLabelValue(dm, "marker", edge, markerBottom);
585: DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
586: if (ex == numXEdges-1) {
587: DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
588: }
589: }
590: }
591: }
592: }
593: DMPlexSymmetrize(dm);
594: DMPlexStratify(dm);
595: /* Build coordinates */
596: DMGetCoordinateSection(dm, &coordSection);
597: PetscSectionSetNumFields(coordSection, 1);
598: PetscSectionSetFieldComponents(coordSection, 0, 2);
599: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);
600: for (v = firstVertex; v < firstVertex+numVertices; ++v) {
601: PetscSectionSetDof(coordSection, v, 2);
602: PetscSectionSetFieldDof(coordSection, v, 0, 2);
603: }
604: PetscSectionSetUp(coordSection);
605: PetscSectionGetStorageSize(coordSection, &coordSize);
606: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
607: VecSetBlockSize(coordinates, 2);
608: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
609: VecSetType(coordinates,VECSTANDARD);
610: VecGetArray(coordinates, &coords);
611: for (vy = 0; vy < numYVertices; ++vy) {
612: for (vx = 0; vx < numXVertices; ++vx) {
613: coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
614: coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
615: }
616: }
617: VecRestoreArray(coordinates, &coords);
618: DMSetCoordinatesLocal(dm, coordinates);
619: VecDestroy(&coordinates);
620: }
621: return(0);
622: }
626: /*@
627: DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices.
629: Collective on MPI_Comm
631: Input Parameters:
632: + comm - The communicator for the DM object
633: . dim - The spatial dimension
634: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
636: Output Parameter:
637: . dm - The DM object
639: Level: beginner
641: .keywords: DM, create
642: .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate()
643: @*/
644: PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm)
645: {
646: DM boundary;
651: DMCreate(comm, &boundary);
653: DMSetType(boundary, DMPLEX);
654: DMPlexSetDimension(boundary, dim-1);
655: switch (dim) {
656: case 2:
657: {
658: PetscReal lower[2] = {0.0, 0.0};
659: PetscReal upper[2] = {1.0, 1.0};
660: PetscInt edges[2] = {2, 2};
662: DMPlexCreateSquareBoundary(boundary, lower, upper, edges);
663: break;
664: }
665: case 3:
666: {
667: PetscReal lower[3] = {0.0, 0.0, 0.0};
668: PetscReal upper[3] = {1.0, 1.0, 1.0};
669: PetscInt faces[3] = {1, 1, 1};
671: DMPlexCreateCubeBoundary(boundary, lower, upper, faces);
672: break;
673: }
674: default:
675: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
676: }
677: DMPlexGenerate(boundary, NULL, interpolate, dm);
678: DMDestroy(&boundary);
679: return(0);
680: }
684: /*@
685: DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra.
687: Collective on MPI_Comm
689: Input Parameters:
690: + comm - The communicator for the DM object
691: . dim - The spatial dimension
692: . periodicX - The boundary type for the X direction
693: . periodicY - The boundary type for the Y direction
694: . periodicZ - The boundary type for the Z direction
695: - cells - The number of cells in each direction
697: Output Parameter:
698: . dm - The DM object
700: Level: beginner
702: .keywords: DM, create
703: .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
704: @*/
705: PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm)
706: {
711: DMCreate(comm, dm);
713: DMSetType(*dm, DMPLEX);
714: DMPlexSetDimension(*dm, dim);
715: switch (dim) {
716: case 2:
717: {
718: PetscReal lower[2] = {0.0, 0.0};
719: PetscReal upper[2] = {1.0, 1.0};
721: DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);
722: break;
723: }
724: #if 0
725: case 3:
726: {
727: PetscReal lower[3] = {0.0, 0.0, 0.0};
728: PetscReal upper[3] = {1.0, 1.0, 1.0};
730: DMPlexCreateCubeMesh(boundary, lower, upper, cells, periodicX, periodicY, periodicZ);
731: break;
732: }
733: #endif
734: default:
735: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
736: }
737: return(0);
738: }
740: /* External function declarations here */
741: extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
742: extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, VecScatter *ctx);
743: extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm);
744: extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
745: extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
746: extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined);
747: extern PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened);
748: extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]);
749: extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
750: extern PetscErrorCode DMSetUp_Plex(DM dm);
751: extern PetscErrorCode DMDestroy_Plex(DM dm);
752: extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
753: extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
754: extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
755: extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS);
759: /* Replace dm with the contents of dmNew
760: - Share the DM_Plex structure
761: - Share the coordinates
762: - Share the SF
763: */
764: static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew)
765: {
766: PetscSF sf;
767: PetscSection coordSection;
768: Vec coords;
772: DMGetPointSF(dmNew, &sf);
773: DMSetPointSF(dm, sf);
774: DMGetCoordinateSection(dmNew, &coordSection);
775: DMGetCoordinatesLocal(dmNew, &coords);
776: DMSetCoordinateSection(dm, coordSection);
777: DMSetCoordinatesLocal(dm, coords);
778: DMDestroy_Plex(dm);
779: dm->data = dmNew->data;
780: ((DM_Plex *) dmNew->data)->refct++;
781: return(0);
782: }
786: /* Swap dm with the contents of dmNew
787: - Swap the DM_Plex structure
788: - Swap the coordinates
789: */
790: static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB)
791: {
792: DM coordDMA, coordDMB;
793: Vec coordsA, coordsB;
794: void *tmp;
798: DMGetCoordinateDM(dmA, &coordDMA);
799: DMGetCoordinateDM(dmB, &coordDMB);
800: PetscObjectReference((PetscObject) coordDMA);
801: DMSetCoordinateDM(dmA, coordDMB);
802: DMSetCoordinateDM(dmB, coordDMA);
803: PetscObjectDereference((PetscObject) coordDMA);
805: DMGetCoordinatesLocal(dmA, &coordsA);
806: DMGetCoordinatesLocal(dmB, &coordsB);
807: PetscObjectReference((PetscObject) coordsA);
808: DMSetCoordinatesLocal(dmA, coordsB);
809: DMSetCoordinatesLocal(dmB, coordsA);
810: PetscObjectDereference((PetscObject) coordsA);
811: tmp = dmA->data;
812: dmA->data = dmB->data;
813: dmB->data = tmp;
814: return(0);
815: }
819: PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM dm)
820: {
821: DM_Plex *mesh = (DM_Plex*) dm->data;
822: DMBoundary b;
826: /* Handle boundary conditions */
827: PetscOptionsBegin(PetscObjectComm((PetscObject) dm), NULL, "Boundary condition options", "");
828: for (b = mesh->boundary; b; b = b->next) {
829: char optname[1024];
830: PetscInt ids[1024], len = 1024, i;
831: PetscBool flg;
833: PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
834: PetscMemzero(ids, sizeof(ids));
835: PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
836: if (flg) {
837: DMLabel label;
839: DMPlexGetLabel(dm, b->labelname, &label);
840: for (i = 0; i < len; ++i) {
841: PetscBool has;
843: DMLabelHasValue(label, ids[i], &has);
844: if (!has) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Boundary id %D is not present in the label %s", ids[i], b->name);
845: }
846: b->numids = len;
847: PetscFree(b->ids);
848: PetscMalloc1(len, &b->ids);
849: PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));
850: }
851: }
852: PetscOptionsEnd();
853: /* Handle viewing */
854: PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);
855: PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);
856: PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", PETSC_FALSE, &mesh->printTol, NULL);
857: return(0);
858: }
862: PetscErrorCode DMSetFromOptions_Plex(DM dm)
863: {
864: PetscInt refine = 0, r;
865: PetscBool isHierarchy;
870: PetscOptionsHead("DMPlex Options");
871: /* Handle DMPlex refinement */
872: PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);
873: PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);
874: if (refine) {DMPlexSetRefinementUniform(dm, PETSC_TRUE);}
875: if (refine && isHierarchy) {
876: DM *dms;
878: PetscMalloc1(refine,&dms);
879: DMRefineHierarchy(dm, refine, dms);
880: /* Total hack since we do not pass in a pointer */
881: DMPlexSwap_Static(dm, dms[refine-1]);
882: if (refine == 1) {
883: DMPlexSetCoarseDM(dm, dms[0]);
884: } else {
885: DMPlexSetCoarseDM(dm, dms[refine-2]);
886: DMPlexSetCoarseDM(dms[0], dms[refine-1]);
887: }
888: /* Free DMs */
889: for (r = 0; r < refine; ++r) {
890: DMSetFromOptions_NonRefinement_Plex(dm);
891: DMDestroy(&dms[r]);
892: }
893: PetscFree(dms);
894: } else {
895: for (r = 0; r < refine; ++r) {
896: DM refinedMesh;
898: DMSetFromOptions_NonRefinement_Plex(dm);
899: DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);
900: /* Total hack since we do not pass in a pointer */
901: DMPlexReplace_Static(dm, refinedMesh);
902: DMDestroy(&refinedMesh);
903: }
904: }
905: DMSetFromOptions_NonRefinement_Plex(dm);
906: PetscOptionsTail();
907: return(0);
908: }
912: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
913: {
917: DMCreateGlobalVector_Section_Private(dm,vec);
918: /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
919: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);
920: VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);
921: return(0);
922: }
926: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
927: {
931: DMCreateLocalVector_Section_Private(dm,vec);
932: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);
933: VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);
934: return(0);
935: }
939: PetscErrorCode DMInitialize_Plex(DM dm)
940: {
942: dm->ops->view = DMView_Plex;
943: dm->ops->load = DMLoad_Plex;
944: dm->ops->setfromoptions = DMSetFromOptions_Plex;
945: dm->ops->clone = DMClone_Plex;
946: dm->ops->setup = DMSetUp_Plex;
947: dm->ops->createdefaultsection = DMCreateDefaultSection_Plex;
948: dm->ops->createglobalvector = DMCreateGlobalVector_Plex;
949: dm->ops->createlocalvector = DMCreateLocalVector_Plex;
950: dm->ops->getlocaltoglobalmapping = NULL;
951: dm->ops->createfieldis = NULL;
952: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex;
953: dm->ops->getcoloring = NULL;
954: dm->ops->creatematrix = DMCreateMatrix_Plex;
955: dm->ops->createinterpolation = DMCreateInterpolation_Plex;
956: dm->ops->getaggregates = NULL;
957: dm->ops->getinjection = DMCreateInjection_Plex;
958: dm->ops->refine = DMRefine_Plex;
959: dm->ops->coarsen = DMCoarsen_Plex;
960: dm->ops->refinehierarchy = DMRefineHierarchy_Plex;
961: dm->ops->coarsenhierarchy = NULL;
962: dm->ops->globaltolocalbegin = NULL;
963: dm->ops->globaltolocalend = NULL;
964: dm->ops->localtoglobalbegin = NULL;
965: dm->ops->localtoglobalend = NULL;
966: dm->ops->destroy = DMDestroy_Plex;
967: dm->ops->createsubdm = DMCreateSubDM_Plex;
968: dm->ops->locatepoints = DMLocatePoints_Plex;
969: return(0);
970: }
974: PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
975: {
976: DM_Plex *mesh = (DM_Plex *) dm->data;
980: mesh->refct++;
981: (*newdm)->data = mesh;
982: PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
983: DMInitialize_Plex(*newdm);
984: return(0);
985: }
987: /*MC
988: DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram.
989: In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is
990: specified by a PetscSection object. Ownership in the global representation is determined by
991: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
993: Level: intermediate
995: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
996: M*/
1000: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
1001: {
1002: DM_Plex *mesh;
1003: PetscInt unit, d;
1008: PetscNewLog(dm,&mesh);
1009: dm->data = mesh;
1011: mesh->refct = 1;
1012: mesh->dim = 0;
1013: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
1014: mesh->maxConeSize = 0;
1015: mesh->cones = NULL;
1016: mesh->coneOrientations = NULL;
1017: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
1018: mesh->maxSupportSize = 0;
1019: mesh->supports = NULL;
1020: mesh->refinementUniform = PETSC_TRUE;
1021: mesh->refinementLimit = -1.0;
1023: mesh->facesTmp = NULL;
1025: mesh->subpointMap = NULL;
1027: for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0;
1029: mesh->labels = NULL;
1030: mesh->depthLabel = NULL;
1031: mesh->globalVertexNumbers = NULL;
1032: mesh->globalCellNumbers = NULL;
1033: for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
1034: mesh->vtkCellHeight = 0;
1035: mesh->useCone = PETSC_FALSE;
1036: mesh->useClosure = PETSC_TRUE;
1038: mesh->printSetValues = PETSC_FALSE;
1039: mesh->printFEM = 0;
1040: mesh->printTol = 1.0e-10;
1042: DMInitialize_Plex(dm);
1043: return(0);
1044: }
1048: /*@
1049: DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.
1051: Collective on MPI_Comm
1053: Input Parameter:
1054: . comm - The communicator for the DMPlex object
1056: Output Parameter:
1057: . mesh - The DMPlex object
1059: Level: beginner
1061: .keywords: DMPlex, create
1062: @*/
1063: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
1064: {
1069: DMCreate(comm, mesh);
1070: DMSetType(*mesh, DMPLEX);
1071: return(0);
1072: }
1076: /*
1077: This takes as input the common mesh generator output, a list of the vertices for each cell
1078: */
1079: PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[])
1080: {
1081: PetscInt *cone, c, p;
1085: DMPlexSetChart(dm, 0, numCells+numVertices);
1086: for (c = 0; c < numCells; ++c) {
1087: DMPlexSetConeSize(dm, c, numCorners);
1088: }
1089: DMSetUp(dm);
1090: DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);
1091: for (c = 0; c < numCells; ++c) {
1092: for (p = 0; p < numCorners; ++p) {
1093: cone[p] = cells[c*numCorners+p]+numCells;
1094: }
1095: DMPlexSetCone(dm, c, cone);
1096: }
1097: DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);
1098: DMPlexSymmetrize(dm);
1099: DMPlexStratify(dm);
1100: return(0);
1101: }
1105: /*
1106: This takes as input the coordinates for each vertex
1107: */
1108: PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
1109: {
1110: PetscSection coordSection;
1111: Vec coordinates;
1112: PetscScalar *coords;
1113: PetscInt coordSize, v, d;
1117: DMGetCoordinateSection(dm, &coordSection);
1118: PetscSectionSetNumFields(coordSection, 1);
1119: PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
1120: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1121: for (v = numCells; v < numCells+numVertices; ++v) {
1122: PetscSectionSetDof(coordSection, v, spaceDim);
1123: PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
1124: }
1125: PetscSectionSetUp(coordSection);
1126: PetscSectionGetStorageSize(coordSection, &coordSize);
1127: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1128: VecSetBlockSize(coordinates, spaceDim);
1129: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1130: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1131: VecSetType(coordinates,VECSTANDARD);
1132: VecGetArray(coordinates, &coords);
1133: for (v = 0; v < numVertices; ++v) {
1134: for (d = 0; d < spaceDim; ++d) {
1135: coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
1136: }
1137: }
1138: VecRestoreArray(coordinates, &coords);
1139: DMSetCoordinatesLocal(dm, coordinates);
1140: VecDestroy(&coordinates);
1141: return(0);
1142: }
1146: /*@C
1147: DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM
1149: Input Parameters:
1150: + comm - The communicator
1151: . dim - The topological dimension of the mesh
1152: . numCells - The number of cells
1153: . numVertices - The number of vertices
1154: . numCorners - The number of vertices for each cell
1155: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
1156: . cells - An array of numCells*numCorners numbers, the vertices for each cell
1157: . spaceDim - The spatial dimension used for coordinates
1158: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
1160: Output Parameter:
1161: . dm - The DM
1163: Note: Two triangles sharing a face
1164: $
1165: $ 2
1166: $ / | \
1167: $ / | \
1168: $ / | \
1169: $ 0 0 | 1 3
1170: $ \ | /
1171: $ \ | /
1172: $ \ | /
1173: $ 1
1174: would have input
1175: $ numCells = 2, numVertices = 4
1176: $ cells = [0 1 2 1 3 2]
1177: $
1178: which would result in the DMPlex
1179: $
1180: $ 4
1181: $ / | \
1182: $ / | \
1183: $ / | \
1184: $ 2 0 | 1 5
1185: $ \ | /
1186: $ \ | /
1187: $ \ | /
1188: $ 3
1190: Level: beginner
1192: .seealso: DMPlexCreateFromDAG(), DMPlexCreate()
1193: @*/
1194: 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)
1195: {
1199: DMCreate(comm, dm);
1200: DMSetType(*dm, DMPLEX);
1201: DMPlexSetDimension(*dm, dim);
1202: DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);
1203: if (interpolate) {
1204: DM idm = NULL;
1206: DMPlexInterpolate(*dm, &idm);
1207: DMDestroy(dm);
1208: *dm = idm;
1209: }
1210: DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);
1211: return(0);
1212: }
1216: /*@
1217: DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM
1219: Input Parameters:
1220: + dm - The empty DM object, usually from DMCreate() and DMPlexSetDimension()
1221: . depth - The depth of the DAG
1222: . numPoints - The number of points at each depth
1223: . coneSize - The cone size of each point
1224: . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
1225: . coneOrientations - The orientation of each cone point
1226: - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex
1228: Output Parameter:
1229: . dm - The DM
1231: Note: Two triangles sharing a face would have input
1232: $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
1233: $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0]
1234: $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0]
1235: $
1236: which would result in the DMPlex
1237: $
1238: $ 4
1239: $ / | \
1240: $ / | \
1241: $ / | \
1242: $ 2 0 | 1 5
1243: $ \ | /
1244: $ \ | /
1245: $ \ | /
1246: $ 3
1247: $
1248: $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList()
1250: Level: advanced
1252: .seealso: DMPlexCreateFromCellList(), DMPlexCreate()
1253: @*/
1254: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
1255: {
1256: Vec coordinates;
1257: PetscSection coordSection;
1258: PetscScalar *coords;
1259: PetscInt coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off;
1263: DMPlexGetDimension(dm, &dim);
1264: for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
1265: DMPlexSetChart(dm, pStart, pEnd);
1266: for (p = pStart; p < pEnd; ++p) {
1267: DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
1268: }
1269: DMSetUp(dm); /* Allocate space for cones */
1270: for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
1271: DMPlexSetCone(dm, p, &cones[off]);
1272: DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
1273: }
1274: DMPlexSymmetrize(dm);
1275: DMPlexStratify(dm);
1276: /* Build coordinates */
1277: DMGetCoordinateSection(dm, &coordSection);
1278: PetscSectionSetNumFields(coordSection, 1);
1279: PetscSectionSetFieldComponents(coordSection, 0, dim);
1280: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
1281: for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
1282: PetscSectionSetDof(coordSection, v, dim);
1283: PetscSectionSetFieldDof(coordSection, v, 0, dim);
1284: }
1285: PetscSectionSetUp(coordSection);
1286: PetscSectionGetStorageSize(coordSection, &coordSize);
1287: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1288: VecSetBlockSize(coordinates, dim);
1289: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1290: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1291: VecSetType(coordinates,VECSTANDARD);
1292: VecGetArray(coordinates, &coords);
1293: for (v = 0; v < numPoints[0]; ++v) {
1294: PetscInt off;
1296: PetscSectionGetOffset(coordSection, v+firstVertex, &off);
1297: for (d = 0; d < dim; ++d) {
1298: coords[off+d] = vertexCoords[v*dim+d];
1299: }
1300: }
1301: VecRestoreArray(coordinates, &coords);
1302: DMSetCoordinatesLocal(dm, coordinates);
1303: VecDestroy(&coordinates);
1304: return(0);
1305: }