Actual source code: plexcreate.c
petsc-3.7.5 2017-01-01
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: DMSetDimension(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) {DMSetLabelValue(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) {DMSetLabelValue(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: DMCopyLabels(*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, 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: const PetscInt numVertices = (edges[0]+1)*(edges[1]+1);
188: const 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)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
204: if (markerSeparate) {
205: markerTop = 3;
206: markerBottom = 1;
207: markerRight = 2;
208: markerLeft = 4;
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: DMSetLabelValue(dm, "marker", edge, markerRight);
229: DMSetLabelValue(dm, "marker", cone[0], markerRight);
230: if (ey == edges[1]-1) {
231: DMSetLabelValue(dm, "marker", cone[1], markerRight);
232: }
233: } else if (vx == 0) {
234: DMSetLabelValue(dm, "marker", edge, markerLeft);
235: DMSetLabelValue(dm, "marker", cone[0], markerLeft);
236: if (ey == edges[1]-1) {
237: DMSetLabelValue(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: DMSetLabelValue(dm, "marker", edge, markerTop);
252: DMSetLabelValue(dm, "marker", cone[0], markerTop);
253: if (ex == edges[0]-1) {
254: DMSetLabelValue(dm, "marker", cone[1], markerTop);
255: }
256: } else if (vy == 0) {
257: DMSetLabelValue(dm, "marker", edge, markerBottom);
258: DMSetLabelValue(dm, "marker", cone[0], markerBottom);
259: if (ex == edges[0]-1) {
260: DMSetLabelValue(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(PETSC_COMM_SELF, &coordinates);
277: PetscObjectSetName((PetscObject) coordinates, "coordinates");
278: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
279: VecSetBlockSize(coordinates, 2);
280: VecSetType(coordinates,VECSTANDARD);
281: VecGetArray(coordinates, &coords);
282: for (vy = 0; vy <= edges[1]; ++vy) {
283: for (vx = 0; vx <= edges[0]; ++vx) {
284: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
285: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
286: }
287: }
288: VecRestoreArray(coordinates, &coords);
289: DMSetCoordinatesLocal(dm, coordinates);
290: VecDestroy(&coordinates);
291: return(0);
292: }
296: /*@
297: DMPlexCreateCubeBoundary - Creates a 2D mesh the is the boundary of a cubic lattice.
299: Collective on MPI_Comm
301: Input Parameters:
302: + comm - The communicator for the DM object
303: . lower - The lower left front corner coordinates
304: . upper - The upper right back corner coordinates
305: - edges - The number of cells in each direction
307: Output Parameter:
308: . dm - The DM object
310: Level: beginner
312: .keywords: DM, create
313: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate()
314: @*/
315: PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
316: {
317: PetscInt vertices[3], numVertices;
318: PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2];
319: Vec coordinates;
320: PetscSection coordSection;
321: PetscScalar *coords;
322: PetscInt coordSize;
323: PetscMPIInt rank;
324: PetscInt v, vx, vy, vz;
325: PetscInt voffset, iface=0, cone[4];
329: 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");
330: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
331: vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1;
332: numVertices = vertices[0]*vertices[1]*vertices[2];
333: if (!rank) {
334: PetscInt f;
336: DMPlexSetChart(dm, 0, numFaces+numVertices);
337: for (f = 0; f < numFaces; ++f) {
338: DMPlexSetConeSize(dm, f, 4);
339: }
340: DMSetUp(dm); /* Allocate space for cones */
341: for (v = 0; v < numFaces+numVertices; ++v) {
342: DMSetLabelValue(dm, "marker", v, 1);
343: }
345: /* Side 0 (Top) */
346: for (vy = 0; vy < faces[1]; vy++) {
347: for (vx = 0; vx < faces[0]; vx++) {
348: voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx;
349: cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0];
350: DMPlexSetCone(dm, iface, cone);
351: iface++;
352: }
353: }
355: /* Side 1 (Bottom) */
356: for (vy = 0; vy < faces[1]; vy++) {
357: for (vx = 0; vx < faces[0]; vx++) {
358: voffset = numFaces + vy*(faces[0]+1) + vx;
359: cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1;
360: DMPlexSetCone(dm, iface, cone);
361: iface++;
362: }
363: }
365: /* Side 2 (Front) */
366: for (vz = 0; vz < faces[2]; vz++) {
367: for (vx = 0; vx < faces[0]; vx++) {
368: voffset = numFaces + vz*vertices[0]*vertices[1] + vx;
369: cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1];
370: DMPlexSetCone(dm, iface, cone);
371: iface++;
372: }
373: }
375: /* Side 3 (Back) */
376: for (vz = 0; vz < faces[2]; vz++) {
377: for (vx = 0; vx < faces[0]; vx++) {
378: voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx;
379: cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1;
380: cone[2] = voffset+1; cone[3] = voffset;
381: DMPlexSetCone(dm, iface, cone);
382: iface++;
383: }
384: }
386: /* Side 4 (Left) */
387: for (vz = 0; vz < faces[2]; vz++) {
388: for (vy = 0; vy < faces[1]; vy++) {
389: voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0];
390: cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1];
391: cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0];
392: DMPlexSetCone(dm, iface, cone);
393: iface++;
394: }
395: }
397: /* Side 5 (Right) */
398: for (vz = 0; vz < faces[2]; vz++) {
399: for (vy = 0; vy < faces[1]; vy++) {
400: voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + vx;
401: cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset;
402: cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0];
403: DMPlexSetCone(dm, iface, cone);
404: iface++;
405: }
406: }
407: }
408: DMPlexSymmetrize(dm);
409: DMPlexStratify(dm);
410: /* Build coordinates */
411: DMGetCoordinateSection(dm, &coordSection);
412: PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);
413: for (v = numFaces; v < numFaces+numVertices; ++v) {
414: PetscSectionSetDof(coordSection, v, 3);
415: }
416: PetscSectionSetUp(coordSection);
417: PetscSectionGetStorageSize(coordSection, &coordSize);
418: VecCreate(PETSC_COMM_SELF, &coordinates);
419: PetscObjectSetName((PetscObject) coordinates, "coordinates");
420: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
421: VecSetBlockSize(coordinates, 3);
422: VecSetType(coordinates,VECSTANDARD);
423: VecGetArray(coordinates, &coords);
424: for (vz = 0; vz <= faces[2]; ++vz) {
425: for (vy = 0; vy <= faces[1]; ++vy) {
426: for (vx = 0; vx <= faces[0]; ++vx) {
427: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
428: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
429: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
430: }
431: }
432: }
433: VecRestoreArray(coordinates, &coords);
434: DMSetCoordinatesLocal(dm, coordinates);
435: VecDestroy(&coordinates);
436: return(0);
437: }
441: static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ)
442: {
443: PetscInt markerTop = 1, faceMarkerTop = 1;
444: PetscInt markerBottom = 1, faceMarkerBottom = 1;
445: PetscInt markerFront = 1, faceMarkerFront = 1;
446: PetscInt markerBack = 1, faceMarkerBack = 1;
447: PetscInt markerRight = 1, faceMarkerRight = 1;
448: PetscInt markerLeft = 1, faceMarkerLeft = 1;
449: PetscInt dim;
450: PetscBool markerSeparate = PETSC_FALSE;
451: PetscMPIInt rank;
455: DMGetDimension(dm,&dim);
456: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
457: switch (dim) {
458: case 2:
459: faceMarkerTop = 3;
460: faceMarkerBottom = 1;
461: faceMarkerRight = 2;
462: faceMarkerLeft = 4;
463: break;
464: case 3:
465: faceMarkerBottom = 1;
466: faceMarkerTop = 2;
467: faceMarkerFront = 3;
468: faceMarkerBack = 4;
469: faceMarkerRight = 5;
470: faceMarkerLeft = 6;
471: break;
472: default:
473: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim);
474: break;
475: }
476: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
477: if (markerSeparate) {
478: markerBottom = faceMarkerBottom;
479: markerTop = faceMarkerTop;
480: markerFront = faceMarkerFront;
481: markerBack = faceMarkerBack;
482: markerRight = faceMarkerRight;
483: markerLeft = faceMarkerLeft;
484: }
485: {
486: const PetscInt numXEdges = !rank ? edges[0] : 0;
487: const PetscInt numYEdges = !rank ? edges[1] : 0;
488: const PetscInt numZEdges = !rank ? edges[2] : 0;
489: const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0;
490: const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0;
491: const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0;
492: const PetscInt numCells = numXEdges*numYEdges*numZEdges;
493: const PetscInt numXFaces = numYEdges*numZEdges;
494: const PetscInt numYFaces = numXEdges*numZEdges;
495: const PetscInt numZFaces = numXEdges*numYEdges;
496: const PetscInt numTotXFaces = numXVertices*numXFaces;
497: const PetscInt numTotYFaces = numYVertices*numYFaces;
498: const PetscInt numTotZFaces = numZVertices*numZFaces;
499: const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces;
500: const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices;
501: const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices;
502: const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices;
503: const PetscInt numVertices = numXVertices*numYVertices*numZVertices;
504: const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges;
505: const PetscInt firstVertex = (dim == 2) ? numFaces : numCells;
506: const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices;
507: const PetscInt firstYFace = firstXFace + numTotXFaces;
508: const PetscInt firstZFace = firstYFace + numTotYFaces;
509: const PetscInt firstXEdge = numCells + numFaces + numVertices;
510: const PetscInt firstYEdge = firstXEdge + numTotXEdges;
511: const PetscInt firstZEdge = firstYEdge + numTotYEdges;
512: Vec coordinates;
513: PetscSection coordSection;
514: PetscScalar *coords;
515: PetscInt coordSize;
516: PetscInt v, vx, vy, vz;
517: PetscInt c, f, fx, fy, fz, e, ex, ey, ez;
519: DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);
520: for (c = 0; c < numCells; c++) {
521: DMPlexSetConeSize(dm, c, 6);
522: }
523: for (f = firstXFace; f < firstXFace+numFaces; ++f) {
524: DMPlexSetConeSize(dm, f, 4);
525: }
526: for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
527: DMPlexSetConeSize(dm, e, 2);
528: }
529: DMSetUp(dm); /* Allocate space for cones */
530: /* Build cells */
531: for (fz = 0; fz < numZEdges; ++fz) {
532: for (fy = 0; fy < numYEdges; ++fy) {
533: for (fx = 0; fx < numXEdges; ++fx) {
534: PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx;
535: PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz;
536: PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices);
537: PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy;
538: PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices);
539: PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx;
540: PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices);
541: /* B, T, F, K, R, L */
542: PetscInt ornt[6] = {-4, 0, 0, -1, 0, -4}; /* ??? */
543: PetscInt cone[6];
545: /* no boundary twisting in 3D */
546: cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL;
547: DMPlexSetCone(dm, cell, cone);
548: DMPlexSetConeOrientation(dm, cell, ornt);
549: }
550: }
551: }
552: /* Build x faces */
553: for (fz = 0; fz < numZEdges; ++fz) {
554: for (fy = 0; fy < numYEdges; ++fy) {
555: for (fx = 0; fx < numXVertices; ++fx) {
556: PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx;
557: PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz;
558: PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz;
559: PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy;
560: PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy;
561: PetscInt ornt[4] = {0, 0, -2, -2};
562: PetscInt cone[4];
564: if (dim == 3) {
565: /* markers */
566: if (bdX != DM_BOUNDARY_PERIODIC) {
567: if (fx == numXVertices-1) {
568: DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);
569: DMSetLabelValue(dm, "marker", face, markerRight);
570: }
571: else if (fx == 0) {
572: DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);
573: DMSetLabelValue(dm, "marker", face, markerLeft);
574: }
575: }
576: }
577: cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
578: DMPlexSetCone(dm, face, cone);
579: DMPlexSetConeOrientation(dm, face, ornt);
580: }
581: }
582: }
583: /* Build y faces */
584: for (fz = 0; fz < numZEdges; ++fz) {
585: for (fx = 0; fx < numXEdges; ++fx) {
586: for (fy = 0; fy < numYVertices; ++fy) {
587: PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy;
588: PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz;
589: PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz;
590: PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx;
591: PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx;
592: PetscInt ornt[4] = {0, 0, -2, -2};
593: PetscInt cone[4];
595: if (dim == 3) {
596: /* markers */
597: if (bdY != DM_BOUNDARY_PERIODIC) {
598: if (fy == numYVertices-1) {
599: DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);
600: DMSetLabelValue(dm, "marker", face, markerBack);
601: }
602: else if (fy == 0) {
603: DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);
604: DMSetLabelValue(dm, "marker", face, markerFront);
605: }
606: }
607: }
608: cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
609: DMPlexSetCone(dm, face, cone);
610: DMPlexSetConeOrientation(dm, face, ornt);
611: }
612: }
613: }
614: /* Build z faces */
615: for (fy = 0; fy < numYEdges; ++fy) {
616: for (fx = 0; fx < numXEdges; ++fx) {
617: for (fz = 0; fz < numZVertices; fz++) {
618: PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz;
619: PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy;
620: PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy;
621: PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx;
622: PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx;
623: PetscInt ornt[4] = {0, 0, -2, -2};
624: PetscInt cone[4];
626: if (dim == 2) {
627: if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;}
628: if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;}
629: }
630: else {
631: /* markers */
632: if (bdZ != DM_BOUNDARY_PERIODIC) {
633: if (fz == numZVertices-1) {
634: DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);
635: DMSetLabelValue(dm, "marker", face, markerTop);
636: }
637: else if (fz == 0) {
638: DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);
639: DMSetLabelValue(dm, "marker", face, markerBottom);
640: }
641: }
642: }
643: cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
644: DMPlexSetCone(dm, face, cone);
645: DMPlexSetConeOrientation(dm, face, ornt);
646: }
647: }
648: }
649: /* Build Z edges*/
650: for (vy = 0; vy < numYVertices; vy++) {
651: for (vx = 0; vx < numXVertices; vx++) {
652: for (ez = 0; ez < numZEdges; ez++) {
653: const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez;
654: const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx;
655: const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx;
656: PetscInt cone[2];
658: if (dim == 3) {
659: if (bdX != DM_BOUNDARY_PERIODIC) {
660: if (vx == numXVertices-1) {
661: DMSetLabelValue(dm, "marker", edge, markerRight);
662: }
663: else if (vx == 0) {
664: DMSetLabelValue(dm, "marker", edge, markerLeft);
665: }
666: }
667: if (bdY != DM_BOUNDARY_PERIODIC) {
668: if (vy == numYVertices-1) {
669: DMSetLabelValue(dm, "marker", edge, markerBack);
670: }
671: else if (vy == 0) {
672: DMSetLabelValue(dm, "marker", edge, markerFront);
673: }
674: }
675: }
676: cone[0] = vertexB; cone[1] = vertexT;
677: DMPlexSetCone(dm, edge, cone);
678: }
679: }
680: }
681: /* Build Y edges*/
682: for (vz = 0; vz < numZVertices; vz++) {
683: for (vx = 0; vx < numXVertices; vx++) {
684: for (ey = 0; ey < numYEdges; ey++) {
685: const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx;
686: const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey;
687: const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx;
688: const PetscInt vertexK = firstVertex + nextv;
689: PetscInt cone[2];
691: cone[0] = vertexF; cone[1] = vertexK;
692: DMPlexSetCone(dm, edge, cone);
693: if (dim == 2) {
694: if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) {
695: if (vx == numXVertices-1) {
696: DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);
697: DMSetLabelValue(dm, "marker", edge, markerRight);
698: DMSetLabelValue(dm, "marker", cone[0], markerRight);
699: if (ey == numYEdges-1) {
700: DMSetLabelValue(dm, "marker", cone[1], markerRight);
701: }
702: }
703: else if (vx == 0) {
704: DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);
705: DMSetLabelValue(dm, "marker", edge, markerLeft);
706: DMSetLabelValue(dm, "marker", cone[0], markerLeft);
707: if (ey == numYEdges-1) {
708: DMSetLabelValue(dm, "marker", cone[1], markerLeft);
709: }
710: }
711: }
712: }
713: else {
714: if (bdX != DM_BOUNDARY_PERIODIC) {
715: if (vx == numXVertices-1) {
716: DMSetLabelValue(dm, "marker", edge, markerRight);
717: }
718: else if (vx == 0) {
719: DMSetLabelValue(dm, "marker", edge, markerLeft);
720: }
721: }
722: if (bdZ != DM_BOUNDARY_PERIODIC) {
723: if (vz == numZVertices-1) {
724: DMSetLabelValue(dm, "marker", edge, markerTop);
725: }
726: else if (vz == 0) {
727: DMSetLabelValue(dm, "marker", edge, markerBottom);
728: }
729: }
730: }
731: }
732: }
733: }
734: /* Build X edges*/
735: for (vz = 0; vz < numZVertices; vz++) {
736: for (vy = 0; vy < numYVertices; vy++) {
737: for (ex = 0; ex < numXEdges; ex++) {
738: const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices;
739: const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex;
740: const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex;
741: const PetscInt vertexR = firstVertex + nextv;
742: PetscInt cone[2];
744: cone[0] = vertexL; cone[1] = vertexR;
745: DMPlexSetCone(dm, edge, cone);
746: if (dim == 2) {
747: if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) {
748: if (vy == numYVertices-1) {
749: DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);
750: DMSetLabelValue(dm, "marker", edge, markerTop);
751: DMSetLabelValue(dm, "marker", cone[0], markerTop);
752: if (ex == numXEdges-1) {
753: DMSetLabelValue(dm, "marker", cone[1], markerTop);
754: }
755: }
756: else if (vy == 0) {
757: DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);
758: DMSetLabelValue(dm, "marker", edge, markerBottom);
759: DMSetLabelValue(dm, "marker", cone[0], markerBottom);
760: if (ex == numXEdges-1) {
761: DMSetLabelValue(dm, "marker", cone[1], markerBottom);
762: }
763: }
764: }
765: }
766: else {
767: if (bdY != DM_BOUNDARY_PERIODIC) {
768: if (vy == numYVertices-1) {
769: DMSetLabelValue(dm, "marker", edge, markerBack);
770: }
771: else if (vy == 0) {
772: DMSetLabelValue(dm, "marker", edge, markerFront);
773: }
774: }
775: if (bdZ != DM_BOUNDARY_PERIODIC) {
776: if (vz == numZVertices-1) {
777: DMSetLabelValue(dm, "marker", edge, markerTop);
778: }
779: else if (vz == 0) {
780: DMSetLabelValue(dm, "marker", edge, markerBottom);
781: }
782: }
783: }
784: }
785: }
786: }
787: DMPlexSymmetrize(dm);
788: DMPlexStratify(dm);
789: /* Build coordinates */
790: DMGetCoordinateSection(dm, &coordSection);
791: PetscSectionSetNumFields(coordSection, 1);
792: PetscSectionSetFieldComponents(coordSection, 0, dim);
793: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);
794: for (v = firstVertex; v < firstVertex+numVertices; ++v) {
795: PetscSectionSetDof(coordSection, v, dim);
796: PetscSectionSetFieldDof(coordSection, v, 0, dim);
797: }
798: PetscSectionSetUp(coordSection);
799: PetscSectionGetStorageSize(coordSection, &coordSize);
800: VecCreate(PETSC_COMM_SELF, &coordinates);
801: PetscObjectSetName((PetscObject) coordinates, "coordinates");
802: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
803: VecSetBlockSize(coordinates, dim);
804: VecSetType(coordinates,VECSTANDARD);
805: VecGetArray(coordinates, &coords);
806: for (vz = 0; vz < numZVertices; ++vz) {
807: for (vy = 0; vy < numYVertices; ++vy) {
808: for (vx = 0; vx < numXVertices; ++vx) {
809: coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
810: coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
811: if (dim == 3) {
812: coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz;
813: }
814: }
815: }
816: }
817: VecRestoreArray(coordinates, &coords);
818: DMSetCoordinatesLocal(dm, coordinates);
819: VecDestroy(&coordinates);
820: }
821: return(0);
822: }
826: /*@
827: DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice.
829: Collective on MPI_Comm
831: Input Parameters:
832: + comm - The communicator for the DM object
833: . lower - The lower left corner coordinates
834: . upper - The upper right corner coordinates
835: . edges - The number of cells in each direction
836: . bdX - The boundary type for the X direction
837: - bdY - The boundary type for the Y direction
839: Output Parameter:
840: . dm - The DM object
842: Note: Here is the numbering returned for 2 cells in each direction:
843: $ 22--8-23--9--24
844: $ | | |
845: $ 13 2 14 3 15
846: $ | | |
847: $ 19--6-20--7--21
848: $ | | |
849: $ 10 0 11 1 12
850: $ | | |
851: $ 16--4-17--5--18
853: Level: beginner
855: .keywords: DM, create
856: .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate()
857: @*/
858: PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY)
859: {
860: PetscReal lower3[3], upper3[3];
861: PetscInt edges3[3];
865: lower3[0] = lower[0]; lower3[1] = lower[1]; lower3[2] = 0.;
866: upper3[0] = upper[0]; upper3[1] = upper[1]; upper3[2] = 0.;
867: edges3[0] = edges[0]; edges3[1] = edges[1]; edges3[2] = 0;
868: DMPlexCreateCubeMesh_Internal(dm, lower3, upper3, edges3, bdX, bdY, DM_BOUNDARY_NONE);
869: return(0);
870: }
874: /*@
875: DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices.
877: Collective on MPI_Comm
879: Input Parameters:
880: + comm - The communicator for the DM object
881: . dim - The spatial dimension
882: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
884: Output Parameter:
885: . dm - The DM object
887: Level: beginner
889: .keywords: DM, create
890: .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate()
891: @*/
892: PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm)
893: {
894: DM boundary;
899: DMCreate(comm, &boundary);
901: DMSetType(boundary, DMPLEX);
902: DMSetDimension(boundary, dim-1);
903: switch (dim) {
904: case 2:
905: {
906: PetscReal lower[2] = {0.0, 0.0};
907: PetscReal upper[2] = {1.0, 1.0};
908: PetscInt edges[2] = {2, 2};
910: DMPlexCreateSquareBoundary(boundary, lower, upper, edges);
911: break;
912: }
913: case 3:
914: {
915: PetscReal lower[3] = {0.0, 0.0, 0.0};
916: PetscReal upper[3] = {1.0, 1.0, 1.0};
917: PetscInt faces[3] = {1, 1, 1};
919: DMPlexCreateCubeBoundary(boundary, lower, upper, faces);
920: break;
921: }
922: default:
923: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
924: }
925: DMPlexGenerate(boundary, NULL, interpolate, dm);
926: DMDestroy(&boundary);
927: return(0);
928: }
932: /*@
933: DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra.
935: Collective on MPI_Comm
937: Input Parameters:
938: + comm - The communicator for the DM object
939: . dim - The spatial dimension
940: . periodicX - The boundary type for the X direction
941: . periodicY - The boundary type for the Y direction
942: . periodicZ - The boundary type for the Z direction
943: - cells - The number of cells in each direction
945: Output Parameter:
946: . dm - The DM object
948: Level: beginner
950: .keywords: DM, create
951: .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
952: @*/
953: PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm)
954: {
959: DMCreate(comm, dm);
961: DMSetType(*dm, DMPLEX);
962: DMSetDimension(*dm, dim);
963: switch (dim) {
964: case 2:
965: {
966: PetscReal lower[2] = {0.0, 0.0};
967: PetscReal upper[2] = {1.0, 1.0};
969: DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);
970: break;
971: }
972: case 3:
973: {
974: PetscReal lower[3] = {0.0, 0.0, 0.0};
975: PetscReal upper[3] = {1.0, 1.0, 1.0};
977: DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);
978: break;
979: }
980: default:
981: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
982: }
983: return(0);
984: }
986: /* External function declarations here */
987: extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
988: extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
989: extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm);
990: extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
991: extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
992: extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
993: extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined);
994: extern PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened);
995: extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]);
996: extern PetscErrorCode DMCoarsenHierarchy_Plex(DM dm, PetscInt nlevels, DM dmCoarsened[]);
997: extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
998: extern PetscErrorCode DMSetUp_Plex(DM dm);
999: extern PetscErrorCode DMDestroy_Plex(DM dm);
1000: extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
1001: extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
1002: extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
1003: extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, PetscSF cellSF);
1004: extern PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
1005: extern PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
1006: extern PetscErrorCode DMProjectFieldLocal_Plex(DM,Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscScalar[]),InsertMode,Vec);
1007: extern PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
1008: extern PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
1009: extern PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
1013: /* Replace dm with the contents of dmNew
1014: - Share the DM_Plex structure
1015: - Share the coordinates
1016: - Share the SF
1017: */
1018: static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew)
1019: {
1020: PetscSF sf;
1021: DM coordDM, coarseDM;
1022: Vec coords;
1023: const PetscReal *maxCell, *L;
1024: const DMBoundaryType *bd;
1025: PetscErrorCode ierr;
1028: DMGetPointSF(dmNew, &sf);
1029: DMSetPointSF(dm, sf);
1030: DMGetCoordinateDM(dmNew, &coordDM);
1031: DMGetCoordinatesLocal(dmNew, &coords);
1032: DMSetCoordinateDM(dm, coordDM);
1033: DMSetCoordinatesLocal(dm, coords);
1034: DMGetPeriodicity(dm, &maxCell, &L, &bd);
1035: if (L) {DMSetPeriodicity(dmNew, maxCell, L, bd);}
1036: DMDestroy_Plex(dm);
1037: dm->data = dmNew->data;
1038: ((DM_Plex *) dmNew->data)->refct++;
1039: dmNew->labels->refct++;
1040: if (!--(dm->labels->refct)) {
1041: DMLabelLink next = dm->labels->next;
1043: /* destroy the labels */
1044: while (next) {
1045: DMLabelLink tmp = next->next;
1047: DMLabelDestroy(&next->label);
1048: PetscFree(next);
1049: next = tmp;
1050: }
1051: PetscFree(dm->labels);
1052: }
1053: dm->labels = dmNew->labels;
1054: dm->depthLabel = dmNew->depthLabel;
1055: DMGetCoarseDM(dmNew,&coarseDM);
1056: DMSetCoarseDM(dm,coarseDM);
1057: return(0);
1058: }
1062: /* Swap dm with the contents of dmNew
1063: - Swap the DM_Plex structure
1064: - Swap the coordinates
1065: - Swap the point PetscSF
1066: */
1067: static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB)
1068: {
1069: DM coordDMA, coordDMB;
1070: Vec coordsA, coordsB;
1071: PetscSF sfA, sfB;
1072: void *tmp;
1073: DMLabelLinkList listTmp;
1074: DMLabel depthTmp;
1075: PetscErrorCode ierr;
1078: DMGetPointSF(dmA, &sfA);
1079: DMGetPointSF(dmB, &sfB);
1080: PetscObjectReference((PetscObject) sfA);
1081: DMSetPointSF(dmA, sfB);
1082: DMSetPointSF(dmB, sfA);
1083: PetscObjectDereference((PetscObject) sfA);
1085: DMGetCoordinateDM(dmA, &coordDMA);
1086: DMGetCoordinateDM(dmB, &coordDMB);
1087: PetscObjectReference((PetscObject) coordDMA);
1088: DMSetCoordinateDM(dmA, coordDMB);
1089: DMSetCoordinateDM(dmB, coordDMA);
1090: PetscObjectDereference((PetscObject) coordDMA);
1092: DMGetCoordinatesLocal(dmA, &coordsA);
1093: DMGetCoordinatesLocal(dmB, &coordsB);
1094: PetscObjectReference((PetscObject) coordsA);
1095: DMSetCoordinatesLocal(dmA, coordsB);
1096: DMSetCoordinatesLocal(dmB, coordsA);
1097: PetscObjectDereference((PetscObject) coordsA);
1099: tmp = dmA->data;
1100: dmA->data = dmB->data;
1101: dmB->data = tmp;
1102: listTmp = dmA->labels;
1103: dmA->labels = dmB->labels;
1104: dmB->labels = listTmp;
1105: depthTmp = dmA->depthLabel;
1106: dmA->depthLabel = dmB->depthLabel;
1107: dmB->depthLabel = depthTmp;
1108: return(0);
1109: }
1113: PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm)
1114: {
1115: DM_Plex *mesh = (DM_Plex*) dm->data;
1116: DMBoundary b;
1120: /* Handle boundary conditions */
1121: PetscOptionsBegin(PetscObjectComm((PetscObject) dm), NULL, "Boundary condition options", "");
1122: for (b = dm->boundary->next; b; b = b->next) {
1123: char optname[1024];
1124: PetscInt ids[1024], len = 1024, i;
1125: PetscBool flg;
1127: PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
1128: PetscMemzero(ids, sizeof(ids));
1129: PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
1130: if (flg) {
1131: DMLabel label;
1133: DMGetLabel(dm, b->labelname, &label);
1134: for (i = 0; i < len; ++i) {
1135: PetscBool has;
1137: DMLabelHasValue(label, ids[i], &has);
1138: if (!has) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Boundary id %D is not present in the label %s", ids[i], b->name);
1139: }
1140: b->numids = len;
1141: PetscFree(b->ids);
1142: PetscMalloc1(len, &b->ids);
1143: PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));
1144: }
1145: PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
1146: PetscMemzero(ids, sizeof(ids));
1147: PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
1148: if (flg) {
1149: b->numcomps = len;
1150: PetscFree(b->comps);
1151: PetscMalloc1(len, &b->comps);
1152: PetscMemcpy(b->comps, ids, len*sizeof(PetscInt));
1153: }
1154: }
1155: PetscOptionsEnd();
1156: /* Handle viewing */
1157: PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);
1158: PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);
1159: PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);
1160: /* Point Location */
1161: PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);
1162: /* Projection behavior */
1163: PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);
1164: PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);
1165: return(0);
1166: }
1170: PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm)
1171: {
1172: PetscInt refine = 0, coarsen = 0, r;
1173: PetscBool isHierarchy;
1178: PetscOptionsHead(PetscOptionsObject,"DMPlex Options");
1179: /* Handle DMPlex refinement */
1180: PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);
1181: PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);
1182: if (refine) {DMPlexSetRefinementUniform(dm, PETSC_TRUE);}
1183: if (refine && isHierarchy) {
1184: DM *dms, coarseDM;
1186: DMGetCoarseDM(dm, &coarseDM);
1187: PetscObjectReference((PetscObject)coarseDM);
1188: PetscMalloc1(refine,&dms);
1189: DMRefineHierarchy(dm, refine, dms);
1190: /* Total hack since we do not pass in a pointer */
1191: DMPlexSwap_Static(dm, dms[refine-1]);
1192: if (refine == 1) {
1193: DMSetCoarseDM(dm, dms[0]);
1194: DMPlexSetRegularRefinement(dm, PETSC_TRUE);
1195: } else {
1196: DMSetCoarseDM(dm, dms[refine-2]);
1197: DMPlexSetRegularRefinement(dm, PETSC_TRUE);
1198: DMSetCoarseDM(dms[0], dms[refine-1]);
1199: DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);
1200: }
1201: DMSetCoarseDM(dms[refine-1], coarseDM);
1202: PetscObjectDereference((PetscObject)coarseDM);
1203: /* Free DMs */
1204: for (r = 0; r < refine; ++r) {
1205: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);
1206: DMDestroy(&dms[r]);
1207: }
1208: PetscFree(dms);
1209: } else {
1210: for (r = 0; r < refine; ++r) {
1211: DM refinedMesh;
1213: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
1214: DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);
1215: /* Total hack since we do not pass in a pointer */
1216: DMPlexReplace_Static(dm, refinedMesh);
1217: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
1218: DMDestroy(&refinedMesh);
1219: }
1220: }
1221: /* Handle DMPlex coarsening */
1222: PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);
1223: PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);
1224: if (coarsen && isHierarchy) {
1225: DM *dms;
1227: PetscMalloc1(coarsen, &dms);
1228: DMCoarsenHierarchy(dm, coarsen, dms);
1229: /* Free DMs */
1230: for (r = 0; r < coarsen; ++r) {
1231: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);
1232: DMDestroy(&dms[r]);
1233: }
1234: PetscFree(dms);
1235: } else {
1236: for (r = 0; r < coarsen; ++r) {
1237: DM coarseMesh;
1239: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
1240: DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);
1241: /* Total hack since we do not pass in a pointer */
1242: DMPlexReplace_Static(dm, coarseMesh);
1243: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
1244: DMDestroy(&coarseMesh);
1245: }
1246: }
1247: /* Handle */
1248: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
1249: PetscOptionsTail();
1250: return(0);
1251: }
1255: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
1256: {
1260: DMCreateGlobalVector_Section_Private(dm,vec);
1261: /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
1262: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);
1263: VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);
1264: VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);
1265: VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);
1266: return(0);
1267: }
1271: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
1272: {
1276: DMCreateLocalVector_Section_Private(dm,vec);
1277: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);
1278: VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);
1279: return(0);
1280: }
1284: static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
1285: {
1286: PetscInt depth, d;
1290: DMPlexGetDepth(dm, &depth);
1291: if (depth == 1) {
1292: DMGetDimension(dm, &d);
1293: if (dim == 0) {DMPlexGetDepthStratum(dm, dim, pStart, pEnd);}
1294: else if (dim == d) {DMPlexGetDepthStratum(dm, 1, pStart, pEnd);}
1295: else {*pStart = 0; *pEnd = 0;}
1296: } else {
1297: DMPlexGetDepthStratum(dm, dim, pStart, pEnd);
1298: }
1299: return(0);
1300: }
1304: PetscErrorCode DMInitialize_Plex(DM dm)
1305: {
1307: dm->ops->view = DMView_Plex;
1308: dm->ops->load = DMLoad_Plex;
1309: dm->ops->setfromoptions = DMSetFromOptions_Plex;
1310: dm->ops->clone = DMClone_Plex;
1311: dm->ops->setup = DMSetUp_Plex;
1312: dm->ops->createdefaultsection = DMCreateDefaultSection_Plex;
1313: dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex;
1314: dm->ops->createglobalvector = DMCreateGlobalVector_Plex;
1315: dm->ops->createlocalvector = DMCreateLocalVector_Plex;
1316: dm->ops->getlocaltoglobalmapping = NULL;
1317: dm->ops->createfieldis = NULL;
1318: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex;
1319: dm->ops->getcoloring = NULL;
1320: dm->ops->creatematrix = DMCreateMatrix_Plex;
1321: dm->ops->createinterpolation = DMCreateInterpolation_Plex;
1322: dm->ops->getaggregates = NULL;
1323: dm->ops->getinjection = DMCreateInjection_Plex;
1324: dm->ops->refine = DMRefine_Plex;
1325: dm->ops->coarsen = DMCoarsen_Plex;
1326: dm->ops->refinehierarchy = DMRefineHierarchy_Plex;
1327: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex;
1328: dm->ops->globaltolocalbegin = NULL;
1329: dm->ops->globaltolocalend = NULL;
1330: dm->ops->localtoglobalbegin = NULL;
1331: dm->ops->localtoglobalend = NULL;
1332: dm->ops->destroy = DMDestroy_Plex;
1333: dm->ops->createsubdm = DMCreateSubDM_Plex;
1334: dm->ops->getdimpoints = DMGetDimPoints_Plex;
1335: dm->ops->locatepoints = DMLocatePoints_Plex;
1336: dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex;
1337: dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex;
1338: dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex;
1339: dm->ops->computel2diff = DMComputeL2Diff_Plex;
1340: dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex;
1341: dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex;
1342: return(0);
1343: }
1347: PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
1348: {
1349: DM_Plex *mesh = (DM_Plex *) dm->data;
1353: mesh->refct++;
1354: (*newdm)->data = mesh;
1355: PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
1356: DMInitialize_Plex(*newdm);
1357: return(0);
1358: }
1360: /*MC
1361: DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram.
1362: In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is
1363: specified by a PetscSection object. Ownership in the global representation is determined by
1364: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
1366: Level: intermediate
1368: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
1369: M*/
1373: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
1374: {
1375: DM_Plex *mesh;
1376: PetscInt unit, d;
1381: PetscNewLog(dm,&mesh);
1382: dm->dim = 0;
1383: dm->data = mesh;
1385: mesh->refct = 1;
1386: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
1387: mesh->maxConeSize = 0;
1388: mesh->cones = NULL;
1389: mesh->coneOrientations = NULL;
1390: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
1391: mesh->maxSupportSize = 0;
1392: mesh->supports = NULL;
1393: mesh->refinementUniform = PETSC_TRUE;
1394: mesh->refinementLimit = -1.0;
1396: mesh->facesTmp = NULL;
1398: mesh->tetgenOpts = NULL;
1399: mesh->triangleOpts = NULL;
1400: PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);
1401: PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);
1403: mesh->subpointMap = NULL;
1405: for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0;
1407: mesh->regularRefinement = PETSC_FALSE;
1408: mesh->depthState = -1;
1409: mesh->globalVertexNumbers = NULL;
1410: mesh->globalCellNumbers = NULL;
1411: mesh->anchorSection = NULL;
1412: mesh->anchorIS = NULL;
1413: mesh->createanchors = NULL;
1414: mesh->computeanchormatrix = NULL;
1415: mesh->parentSection = NULL;
1416: mesh->parents = NULL;
1417: mesh->childIDs = NULL;
1418: mesh->childSection = NULL;
1419: mesh->children = NULL;
1420: mesh->referenceTree = NULL;
1421: mesh->getchildsymmetry = NULL;
1422: for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
1423: mesh->vtkCellHeight = 0;
1424: mesh->useCone = PETSC_FALSE;
1425: mesh->useClosure = PETSC_TRUE;
1426: mesh->useAnchors = PETSC_FALSE;
1428: mesh->maxProjectionHeight = 0;
1430: mesh->printSetValues = PETSC_FALSE;
1431: mesh->printFEM = 0;
1432: mesh->printTol = 1.0e-10;
1434: DMInitialize_Plex(dm);
1435: return(0);
1436: }
1440: /*@
1441: DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.
1443: Collective on MPI_Comm
1445: Input Parameter:
1446: . comm - The communicator for the DMPlex object
1448: Output Parameter:
1449: . mesh - The DMPlex object
1451: Level: beginner
1453: .keywords: DMPlex, create
1454: @*/
1455: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
1456: {
1461: DMCreate(comm, mesh);
1462: DMSetType(*mesh, DMPLEX);
1463: return(0);
1464: }
1468: /*
1469: This takes as input the common mesh generator output, a list of the vertices for each cell
1470: */
1471: PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[])
1472: {
1473: PetscInt *cone, c, p;
1477: DMPlexSetChart(dm, 0, numCells+numVertices);
1478: for (c = 0; c < numCells; ++c) {
1479: DMPlexSetConeSize(dm, c, numCorners);
1480: }
1481: DMSetUp(dm);
1482: DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);
1483: for (c = 0; c < numCells; ++c) {
1484: for (p = 0; p < numCorners; ++p) {
1485: cone[p] = cells[c*numCorners+p]+numCells;
1486: }
1487: DMPlexSetCone(dm, c, cone);
1488: }
1489: DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);
1490: DMPlexSymmetrize(dm);
1491: DMPlexStratify(dm);
1492: return(0);
1493: }
1497: /*
1498: This takes as input the coordinates for each vertex
1499: */
1500: PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
1501: {
1502: PetscSection coordSection;
1503: Vec coordinates;
1504: DM cdm;
1505: PetscScalar *coords;
1506: PetscInt v, d;
1510: DMGetCoordinateSection(dm, &coordSection);
1511: PetscSectionSetNumFields(coordSection, 1);
1512: PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
1513: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1514: for (v = numCells; v < numCells+numVertices; ++v) {
1515: PetscSectionSetDof(coordSection, v, spaceDim);
1516: PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
1517: }
1518: PetscSectionSetUp(coordSection);
1520: DMGetCoordinateDM(dm, &cdm);
1521: DMCreateLocalVector(cdm, &coordinates);
1522: VecSetBlockSize(coordinates, spaceDim);
1523: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1524: VecGetArray(coordinates, &coords);
1525: for (v = 0; v < numVertices; ++v) {
1526: for (d = 0; d < spaceDim; ++d) {
1527: coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
1528: }
1529: }
1530: VecRestoreArray(coordinates, &coords);
1531: DMSetCoordinatesLocal(dm, coordinates);
1532: VecDestroy(&coordinates);
1533: return(0);
1534: }
1538: /*@C
1539: DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM
1541: Input Parameters:
1542: + comm - The communicator
1543: . dim - The topological dimension of the mesh
1544: . numCells - The number of cells
1545: . numVertices - The number of vertices
1546: . numCorners - The number of vertices for each cell
1547: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
1548: . cells - An array of numCells*numCorners numbers, the vertices for each cell
1549: . spaceDim - The spatial dimension used for coordinates
1550: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
1552: Output Parameter:
1553: . dm - The DM
1555: Note: Two triangles sharing a face
1556: $
1557: $ 2
1558: $ / | \
1559: $ / | \
1560: $ / | \
1561: $ 0 0 | 1 3
1562: $ \ | /
1563: $ \ | /
1564: $ \ | /
1565: $ 1
1566: would have input
1567: $ numCells = 2, numVertices = 4
1568: $ cells = [0 1 2 1 3 2]
1569: $
1570: which would result in the DMPlex
1571: $
1572: $ 4
1573: $ / | \
1574: $ / | \
1575: $ / | \
1576: $ 2 0 | 1 5
1577: $ \ | /
1578: $ \ | /
1579: $ \ | /
1580: $ 3
1582: Level: beginner
1584: .seealso: DMPlexCreateFromDAG(), DMPlexCreate()
1585: @*/
1586: 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)
1587: {
1591: DMCreate(comm, dm);
1592: DMSetType(*dm, DMPLEX);
1593: DMSetDimension(*dm, dim);
1594: DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);
1595: if (interpolate) {
1596: DM idm = NULL;
1598: DMPlexInterpolate(*dm, &idm);
1599: DMDestroy(dm);
1600: *dm = idm;
1601: }
1602: DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);
1603: return(0);
1604: }
1608: /*@
1609: DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM
1611: Input Parameters:
1612: + dm - The empty DM object, usually from DMCreate() and DMSetDimension()
1613: . depth - The depth of the DAG
1614: . numPoints - The number of points at each depth
1615: . coneSize - The cone size of each point
1616: . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
1617: . coneOrientations - The orientation of each cone point
1618: - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex
1620: Output Parameter:
1621: . dm - The DM
1623: Note: Two triangles sharing a face would have input
1624: $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
1625: $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0]
1626: $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0]
1627: $
1628: which would result in the DMPlex
1629: $
1630: $ 4
1631: $ / | \
1632: $ / | \
1633: $ / | \
1634: $ 2 0 | 1 5
1635: $ \ | /
1636: $ \ | /
1637: $ \ | /
1638: $ 3
1639: $
1640: $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList()
1642: Level: advanced
1644: .seealso: DMPlexCreateFromCellList(), DMPlexCreate()
1645: @*/
1646: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
1647: {
1648: Vec coordinates;
1649: PetscSection coordSection;
1650: PetscScalar *coords;
1651: PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off;
1655: DMGetDimension(dm, &dim);
1656: DMGetCoordinateDim(dm, &dimEmbed);
1657: if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim);
1658: for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
1659: DMPlexSetChart(dm, pStart, pEnd);
1660: for (p = pStart; p < pEnd; ++p) {
1661: DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
1662: if (firstVertex < 0 && !coneSize[p - pStart]) {
1663: firstVertex = p - pStart;
1664: }
1665: }
1666: if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]);
1667: DMSetUp(dm); /* Allocate space for cones */
1668: for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
1669: DMPlexSetCone(dm, p, &cones[off]);
1670: DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
1671: }
1672: DMPlexSymmetrize(dm);
1673: DMPlexStratify(dm);
1674: /* Build coordinates */
1675: DMGetCoordinateSection(dm, &coordSection);
1676: PetscSectionSetNumFields(coordSection, 1);
1677: PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
1678: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
1679: for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
1680: PetscSectionSetDof(coordSection, v, dimEmbed);
1681: PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
1682: }
1683: PetscSectionSetUp(coordSection);
1684: PetscSectionGetStorageSize(coordSection, &coordSize);
1685: VecCreate(PETSC_COMM_SELF, &coordinates);
1686: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1687: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1688: VecSetBlockSize(coordinates, dimEmbed);
1689: VecSetType(coordinates,VECSTANDARD);
1690: VecGetArray(coordinates, &coords);
1691: for (v = 0; v < numPoints[0]; ++v) {
1692: PetscInt off;
1694: PetscSectionGetOffset(coordSection, v+firstVertex, &off);
1695: for (d = 0; d < dimEmbed; ++d) {
1696: coords[off+d] = vertexCoords[v*dimEmbed+d];
1697: }
1698: }
1699: VecRestoreArray(coordinates, &coords);
1700: DMSetCoordinatesLocal(dm, coordinates);
1701: VecDestroy(&coordinates);
1702: return(0);
1703: }
1707: /*@C
1708: DMPlexCreateFromFile - This takes a filename and produces a DM
1710: Input Parameters:
1711: + comm - The communicator
1712: . filename - A file name
1713: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
1715: Output Parameter:
1716: . dm - The DM
1718: Level: beginner
1720: .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate()
1721: @*/
1722: PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1723: {
1724: const char *extGmsh = ".msh";
1725: const char *extCGNS = ".cgns";
1726: const char *extExodus = ".exo";
1727: const char *extFluent = ".cas";
1728: const char *extHDF5 = ".h5";
1729: size_t len;
1730: PetscBool isGmsh, isCGNS, isExodus, isFluent, isHDF5;
1731: PetscMPIInt rank;
1737: MPI_Comm_rank(comm, &rank);
1738: PetscStrlen(filename, &len);
1739: if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path");
1740: PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);
1741: PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);
1742: PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);
1743: PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);
1744: PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);
1745: if (isGmsh) {
1746: DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);
1747: } else if (isCGNS) {
1748: DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);
1749: } else if (isExodus) {
1750: DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);
1751: } else if (isFluent) {
1752: DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);
1753: } else if (isHDF5) {
1754: PetscViewer viewer;
1756: PetscViewerCreate(comm, &viewer);
1757: PetscViewerSetType(viewer, PETSCVIEWERHDF5);
1758: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1759: PetscViewerFileSetName(viewer, filename);
1760: DMCreate(comm, dm);
1761: DMSetType(*dm, DMPLEX);
1762: DMLoad(*dm, viewer);
1763: PetscViewerDestroy(&viewer);
1764: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename);
1765: return(0);
1766: }
1770: /*@
1771: DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
1773: Collective on comm
1775: Input Parameters:
1776: + comm - The communicator
1777: . dim - The spatial dimension
1778: - simplex - Flag for simplex, otherwise use a tensor-product cell
1780: Output Parameter:
1781: . refdm - The reference cell
1783: Level: intermediate
1785: .keywords: reference cell
1786: .seealso:
1787: @*/
1788: PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm)
1789: {
1790: DM rdm;
1791: Vec coords;
1795: DMCreate(comm, &rdm);
1796: DMSetType(rdm, DMPLEX);
1797: DMSetDimension(rdm, dim);
1798: switch (dim) {
1799: case 0:
1800: {
1801: PetscInt numPoints[1] = {1};
1802: PetscInt coneSize[1] = {0};
1803: PetscInt cones[1] = {0};
1804: PetscInt coneOrientations[1] = {0};
1805: PetscScalar vertexCoords[1] = {0.0};
1807: DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1808: }
1809: break;
1810: case 1:
1811: {
1812: PetscInt numPoints[2] = {2, 1};
1813: PetscInt coneSize[3] = {2, 0, 0};
1814: PetscInt cones[2] = {1, 2};
1815: PetscInt coneOrientations[2] = {0, 0};
1816: PetscScalar vertexCoords[2] = {-1.0, 1.0};
1818: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1819: }
1820: break;
1821: case 2:
1822: if (simplex) {
1823: PetscInt numPoints[2] = {3, 1};
1824: PetscInt coneSize[4] = {3, 0, 0, 0};
1825: PetscInt cones[3] = {1, 2, 3};
1826: PetscInt coneOrientations[3] = {0, 0, 0};
1827: PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0};
1829: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1830: } else {
1831: PetscInt numPoints[2] = {4, 1};
1832: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
1833: PetscInt cones[4] = {1, 2, 3, 4};
1834: PetscInt coneOrientations[4] = {0, 0, 0, 0};
1835: PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
1837: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1838: }
1839: break;
1840: case 3:
1841: if (simplex) {
1842: PetscInt numPoints[2] = {4, 1};
1843: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
1844: PetscInt cones[4] = {1, 3, 2, 4};
1845: PetscInt coneOrientations[4] = {0, 0, 0, 0};
1846: PetscScalar vertexCoords[12] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0};
1848: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1849: } else {
1850: PetscInt numPoints[2] = {8, 1};
1851: PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0};
1852: PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8};
1853: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
1854: PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0,
1855: -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
1857: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1858: }
1859: break;
1860: default:
1861: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim);
1862: }
1863: *refdm = NULL;
1864: DMPlexInterpolate(rdm, refdm);
1865: if (rdm->coordinateDM) {
1866: DM ncdm;
1867: PetscSection cs;
1868: PetscInt pEnd = -1;
1870: DMGetDefaultSection(rdm->coordinateDM, &cs);
1871: if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
1872: if (pEnd >= 0) {
1873: DMClone(rdm->coordinateDM, &ncdm);
1874: DMSetDefaultSection(ncdm, cs);
1875: DMSetCoordinateDM(*refdm, ncdm);
1876: DMDestroy(&ncdm);
1877: }
1878: }
1879: DMGetCoordinatesLocal(rdm, &coords);
1880: if (coords) {
1881: DMSetCoordinatesLocal(*refdm, coords);
1882: } else {
1883: DMGetCoordinates(rdm, &coords);
1884: if (coords) {DMSetCoordinates(*refdm, coords);}
1885: }
1886: DMDestroy(&rdm);
1887: return(0);
1888: }