Actual source code: plexcreate.c
petsc-3.6.1 2015-08-06
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) {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, 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)->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: 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: PetscObjectSetName((PetscObject) coordinates, "coordinates");
279: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
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: DMPlexSetLabelValue(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(PetscObjectComm((PetscObject)dm), &coordinates);
419: VecSetBlockSize(coordinates, 3);
420: PetscObjectSetName((PetscObject) coordinates, "coordinates");
421: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
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)->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 || bdY == 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[8] = {-4, 0, 0, -1, 0, -4}; /* ??? */
543: PetscInt cone[8];
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: DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerRight);
569: DMPlexSetLabelValue(dm, "marker", face, markerRight);
570: }
571: else if (fx == 0) {
572: DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);
573: DMPlexSetLabelValue(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 < numYEdges; ++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: DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerBack);
600: DMPlexSetLabelValue(dm, "marker", face, markerBack);
601: }
602: else if (fy == 0) {
603: DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerFront);
604: DMPlexSetLabelValue(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: DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerTop);
635: DMPlexSetLabelValue(dm, "marker", face, markerTop);
636: }
637: else if (fz == 0) {
638: DMPlexSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);
639: DMPlexSetLabelValue(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: DMPlexSetLabelValue(dm, "marker", edge, markerRight);
662: }
663: else if (vx == 0) {
664: DMPlexSetLabelValue(dm, "marker", edge, markerLeft);
665: }
666: }
667: if (bdY != DM_BOUNDARY_PERIODIC) {
668: if (vy == numYVertices-1) {
669: DMPlexSetLabelValue(dm, "marker", edge, markerBack);
670: }
671: else if (vy == 0) {
672: DMPlexSetLabelValue(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: DMPlexSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);
697: DMPlexSetLabelValue(dm, "marker", edge, markerRight);
698: DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
699: if (ey == numYEdges-1) {
700: DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
701: }
702: }
703: else if (vx == 0) {
704: DMPlexSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);
705: DMPlexSetLabelValue(dm, "marker", edge, markerLeft);
706: DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
707: if (ey == numYEdges-1) {
708: DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
709: }
710: }
711: }
712: }
713: else {
714: if (bdX != DM_BOUNDARY_PERIODIC) {
715: if (vx == numXVertices-1) {
716: DMPlexSetLabelValue(dm, "marker", edge, markerRight);
717: }
718: else if (vx == 0) {
719: DMPlexSetLabelValue(dm, "marker", edge, markerLeft);
720: }
721: }
722: if (bdZ != DM_BOUNDARY_PERIODIC) {
723: if (vz == numZVertices-1) {
724: DMPlexSetLabelValue(dm, "marker", edge, markerTop);
725: }
726: else if (vz == 0) {
727: DMPlexSetLabelValue(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: DMPlexSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);
750: DMPlexSetLabelValue(dm, "marker", edge, markerTop);
751: DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
752: if (ex == numXEdges-1) {
753: DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
754: }
755: }
756: else if (vy == 0) {
757: DMPlexSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);
758: DMPlexSetLabelValue(dm, "marker", edge, markerBottom);
759: DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
760: if (ex == numXEdges-1) {
761: DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
762: }
763: }
764: }
765: }
766: else {
767: if (bdY != DM_BOUNDARY_PERIODIC) {
768: if (vy == numYVertices-1) {
769: DMPlexSetLabelValue(dm, "marker", edge, markerBack);
770: }
771: else if (vy == 0) {
772: DMPlexSetLabelValue(dm, "marker", edge, markerFront);
773: }
774: }
775: if (bdZ != DM_BOUNDARY_PERIODIC) {
776: if (vz == numZVertices-1) {
777: DMPlexSetLabelValue(dm, "marker", edge, markerTop);
778: }
779: else if (vz == 0) {
780: DMPlexSetLabelValue(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(PetscObjectComm((PetscObject)dm), &coordinates);
801: VecSetBlockSize(coordinates, dim);
802: PetscObjectSetName((PetscObject) coordinates, "coordinates");
803: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
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 DMClone_Plex(DM dm, DM *newdm);
997: extern PetscErrorCode DMSetUp_Plex(DM dm);
998: extern PetscErrorCode DMDestroy_Plex(DM dm);
999: extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
1000: extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
1001: extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
1002: extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS);
1006: /* Replace dm with the contents of dmNew
1007: - Share the DM_Plex structure
1008: - Share the coordinates
1009: - Share the SF
1010: */
1011: static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew)
1012: {
1013: PetscSF sf;
1014: DM coordDM;
1015: Vec coords;
1016: const PetscReal *maxCell, *L;
1017: const DMBoundaryType *bd;
1018: PetscErrorCode ierr;
1021: DMGetPointSF(dmNew, &sf);
1022: DMSetPointSF(dm, sf);
1023: DMGetCoordinateDM(dmNew, &coordDM);
1024: DMGetCoordinatesLocal(dmNew, &coords);
1025: DMSetCoordinateDM(dm, coordDM);
1026: DMSetCoordinatesLocal(dm, coords);
1027: DMGetPeriodicity(dm, &maxCell, &L, &bd);
1028: if (L) {DMSetPeriodicity(dmNew, maxCell, L, bd);}
1029: DMDestroy_Plex(dm);
1030: dm->data = dmNew->data;
1031: ((DM_Plex *) dmNew->data)->refct++;
1032: return(0);
1033: }
1037: /* Swap dm with the contents of dmNew
1038: - Swap the DM_Plex structure
1039: - Swap the coordinates
1040: */
1041: static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB)
1042: {
1043: DM coordDMA, coordDMB;
1044: Vec coordsA, coordsB;
1045: void *tmp;
1049: DMGetCoordinateDM(dmA, &coordDMA);
1050: DMGetCoordinateDM(dmB, &coordDMB);
1051: PetscObjectReference((PetscObject) coordDMA);
1052: DMSetCoordinateDM(dmA, coordDMB);
1053: DMSetCoordinateDM(dmB, coordDMA);
1054: PetscObjectDereference((PetscObject) coordDMA);
1056: DMGetCoordinatesLocal(dmA, &coordsA);
1057: DMGetCoordinatesLocal(dmB, &coordsB);
1058: PetscObjectReference((PetscObject) coordsA);
1059: DMSetCoordinatesLocal(dmA, coordsB);
1060: DMSetCoordinatesLocal(dmB, coordsA);
1061: PetscObjectDereference((PetscObject) coordsA);
1062: tmp = dmA->data;
1063: dmA->data = dmB->data;
1064: dmB->data = tmp;
1065: return(0);
1066: }
1070: PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptions *PetscOptionsObject,DM dm)
1071: {
1072: DM_Plex *mesh = (DM_Plex*) dm->data;
1073: DMBoundary b;
1077: /* Handle boundary conditions */
1078: PetscOptionsBegin(PetscObjectComm((PetscObject) dm), NULL, "Boundary condition options", "");
1079: for (b = mesh->boundary; b; b = b->next) {
1080: char optname[1024];
1081: PetscInt ids[1024], len = 1024, i;
1082: PetscBool flg;
1084: PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
1085: PetscMemzero(ids, sizeof(ids));
1086: PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
1087: if (flg) {
1088: DMLabel label;
1090: DMPlexGetLabel(dm, b->labelname, &label);
1091: for (i = 0; i < len; ++i) {
1092: PetscBool has;
1094: DMLabelHasValue(label, ids[i], &has);
1095: if (!has) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Boundary id %D is not present in the label %s", ids[i], b->name);
1096: }
1097: b->numids = len;
1098: PetscFree(b->ids);
1099: PetscMalloc1(len, &b->ids);
1100: PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));
1101: }
1102: PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
1103: PetscMemzero(ids, sizeof(ids));
1104: PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
1105: if (flg) {
1106: b->numcomps = len;
1107: PetscFree(b->comps);
1108: PetscMalloc1(len, &b->comps);
1109: PetscMemcpy(b->comps, ids, len*sizeof(PetscInt));
1110: }
1111: }
1112: PetscOptionsEnd();
1113: /* Handle viewing */
1114: PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);
1115: PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);
1116: PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);
1117: /* Projection behavior */
1118: PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);
1119: return(0);
1120: }
1124: PetscErrorCode DMSetFromOptions_Plex(PetscOptions *PetscOptionsObject,DM dm)
1125: {
1126: PetscInt refine = 0, r;
1127: PetscBool isHierarchy;
1132: PetscOptionsHead(PetscOptionsObject,"DMPlex Options");
1133: /* Handle DMPlex refinement */
1134: PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);
1135: PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);
1136: if (refine) {DMPlexSetRefinementUniform(dm, PETSC_TRUE);}
1137: if (refine && isHierarchy) {
1138: DM *dms;
1140: PetscMalloc1(refine,&dms);
1141: DMRefineHierarchy(dm, refine, dms);
1142: /* Total hack since we do not pass in a pointer */
1143: DMPlexSwap_Static(dm, dms[refine-1]);
1144: if (refine == 1) {
1145: DMPlexSetCoarseDM(dm, dms[0]);
1146: } else {
1147: DMPlexSetCoarseDM(dm, dms[refine-2]);
1148: DMPlexSetCoarseDM(dms[0], dms[refine-1]);
1149: }
1150: /* Free DMs */
1151: for (r = 0; r < refine; ++r) {
1152: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject,dm);
1153: DMDestroy(&dms[r]);
1154: }
1155: PetscFree(dms);
1156: } else {
1157: for (r = 0; r < refine; ++r) {
1158: DM refinedMesh;
1160: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject,dm);
1161: DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);
1162: /* Total hack since we do not pass in a pointer */
1163: DMPlexReplace_Static(dm, refinedMesh);
1164: DMDestroy(&refinedMesh);
1165: }
1166: }
1167: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject,dm);
1168: PetscOptionsTail();
1169: return(0);
1170: }
1174: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
1175: {
1179: DMCreateGlobalVector_Section_Private(dm,vec);
1180: /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
1181: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);
1182: VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);
1183: return(0);
1184: }
1188: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
1189: {
1193: DMCreateLocalVector_Section_Private(dm,vec);
1194: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);
1195: VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);
1196: return(0);
1197: }
1201: static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
1202: {
1203: PetscInt depth, d;
1207: DMPlexGetDepth(dm, &depth);
1208: if (depth == 1) {
1209: DMGetDimension(dm, &d);
1210: if (dim == 0) {DMPlexGetDepthStratum(dm, dim, pStart, pEnd);}
1211: else if (dim == d) {DMPlexGetDepthStratum(dm, 1, pStart, pEnd);}
1212: else {*pStart = 0; *pEnd = 0;}
1213: } else {
1214: DMPlexGetDepthStratum(dm, dim, pStart, pEnd);
1215: }
1216: return(0);
1217: }
1221: PetscErrorCode DMInitialize_Plex(DM dm)
1222: {
1224: dm->ops->view = DMView_Plex;
1225: dm->ops->load = DMLoad_Plex;
1226: dm->ops->setfromoptions = DMSetFromOptions_Plex;
1227: dm->ops->clone = DMClone_Plex;
1228: dm->ops->setup = DMSetUp_Plex;
1229: dm->ops->createdefaultsection = DMCreateDefaultSection_Plex;
1230: dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex;
1231: dm->ops->createglobalvector = DMCreateGlobalVector_Plex;
1232: dm->ops->createlocalvector = DMCreateLocalVector_Plex;
1233: dm->ops->getlocaltoglobalmapping = NULL;
1234: dm->ops->createfieldis = NULL;
1235: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex;
1236: dm->ops->getcoloring = NULL;
1237: dm->ops->creatematrix = DMCreateMatrix_Plex;
1238: dm->ops->createinterpolation = DMCreateInterpolation_Plex;
1239: dm->ops->getaggregates = NULL;
1240: dm->ops->getinjection = DMCreateInjection_Plex;
1241: dm->ops->refine = DMRefine_Plex;
1242: dm->ops->coarsen = DMCoarsen_Plex;
1243: dm->ops->refinehierarchy = DMRefineHierarchy_Plex;
1244: dm->ops->coarsenhierarchy = NULL;
1245: dm->ops->globaltolocalbegin = NULL;
1246: dm->ops->globaltolocalend = NULL;
1247: dm->ops->localtoglobalbegin = NULL;
1248: dm->ops->localtoglobalend = NULL;
1249: dm->ops->destroy = DMDestroy_Plex;
1250: dm->ops->createsubdm = DMCreateSubDM_Plex;
1251: dm->ops->getdimpoints = DMGetDimPoints_Plex;
1252: dm->ops->locatepoints = DMLocatePoints_Plex;
1253: return(0);
1254: }
1258: PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
1259: {
1260: DM_Plex *mesh = (DM_Plex *) dm->data;
1264: mesh->refct++;
1265: (*newdm)->data = mesh;
1266: PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
1267: DMInitialize_Plex(*newdm);
1268: return(0);
1269: }
1271: /*MC
1272: DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram.
1273: In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is
1274: specified by a PetscSection object. Ownership in the global representation is determined by
1275: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
1277: Level: intermediate
1279: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
1280: M*/
1284: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
1285: {
1286: DM_Plex *mesh;
1287: PetscInt unit, d;
1292: PetscNewLog(dm,&mesh);
1293: dm->dim = 0;
1294: dm->data = mesh;
1296: mesh->refct = 1;
1297: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
1298: mesh->maxConeSize = 0;
1299: mesh->cones = NULL;
1300: mesh->coneOrientations = NULL;
1301: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
1302: mesh->maxSupportSize = 0;
1303: mesh->supports = NULL;
1304: mesh->refinementUniform = PETSC_TRUE;
1305: mesh->refinementLimit = -1.0;
1307: mesh->facesTmp = NULL;
1309: mesh->tetgenOpts = NULL;
1310: mesh->triangleOpts = NULL;
1311: PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);
1312: PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);
1314: mesh->subpointMap = NULL;
1316: for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0;
1318: mesh->labels = NULL;
1319: mesh->depthLabel = NULL;
1320: mesh->globalVertexNumbers = NULL;
1321: mesh->globalCellNumbers = NULL;
1322: mesh->anchorSection = NULL;
1323: mesh->anchorIS = NULL;
1324: mesh->createanchors = NULL;
1325: mesh->computeanchormatrix = NULL;
1326: mesh->parentSection = NULL;
1327: mesh->parents = NULL;
1328: mesh->childIDs = NULL;
1329: mesh->childSection = NULL;
1330: mesh->children = NULL;
1331: mesh->referenceTree = NULL;
1332: mesh->getchildsymmetry = NULL;
1333: for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
1334: mesh->vtkCellHeight = 0;
1335: mesh->useCone = PETSC_FALSE;
1336: mesh->useClosure = PETSC_TRUE;
1337: mesh->useAnchors = PETSC_FALSE;
1339: mesh->maxProjectionHeight = 0;
1341: mesh->printSetValues = PETSC_FALSE;
1342: mesh->printFEM = 0;
1343: mesh->printTol = 1.0e-10;
1345: DMInitialize_Plex(dm);
1346: return(0);
1347: }
1351: /*@
1352: DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.
1354: Collective on MPI_Comm
1356: Input Parameter:
1357: . comm - The communicator for the DMPlex object
1359: Output Parameter:
1360: . mesh - The DMPlex object
1362: Level: beginner
1364: .keywords: DMPlex, create
1365: @*/
1366: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
1367: {
1372: DMCreate(comm, mesh);
1373: DMSetType(*mesh, DMPLEX);
1374: return(0);
1375: }
1379: /*
1380: This takes as input the common mesh generator output, a list of the vertices for each cell
1381: */
1382: PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[])
1383: {
1384: PetscInt *cone, c, p;
1388: DMPlexSetChart(dm, 0, numCells+numVertices);
1389: for (c = 0; c < numCells; ++c) {
1390: DMPlexSetConeSize(dm, c, numCorners);
1391: }
1392: DMSetUp(dm);
1393: DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);
1394: for (c = 0; c < numCells; ++c) {
1395: for (p = 0; p < numCorners; ++p) {
1396: cone[p] = cells[c*numCorners+p]+numCells;
1397: }
1398: DMPlexSetCone(dm, c, cone);
1399: }
1400: DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);
1401: DMPlexSymmetrize(dm);
1402: DMPlexStratify(dm);
1403: return(0);
1404: }
1408: /*
1409: This takes as input the coordinates for each vertex
1410: */
1411: PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
1412: {
1413: PetscSection coordSection;
1414: Vec coordinates;
1415: PetscScalar *coords;
1416: PetscInt coordSize, v, d;
1420: DMGetCoordinateSection(dm, &coordSection);
1421: PetscSectionSetNumFields(coordSection, 1);
1422: PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
1423: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
1424: for (v = numCells; v < numCells+numVertices; ++v) {
1425: PetscSectionSetDof(coordSection, v, spaceDim);
1426: PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
1427: }
1428: PetscSectionSetUp(coordSection);
1429: PetscSectionGetStorageSize(coordSection, &coordSize);
1430: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1431: VecSetBlockSize(coordinates, spaceDim);
1432: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1433: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1434: VecSetType(coordinates,VECSTANDARD);
1435: VecGetArray(coordinates, &coords);
1436: for (v = 0; v < numVertices; ++v) {
1437: for (d = 0; d < spaceDim; ++d) {
1438: coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
1439: }
1440: }
1441: VecRestoreArray(coordinates, &coords);
1442: DMSetCoordinatesLocal(dm, coordinates);
1443: VecDestroy(&coordinates);
1444: return(0);
1445: }
1449: /*@C
1450: DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM
1452: Input Parameters:
1453: + comm - The communicator
1454: . dim - The topological dimension of the mesh
1455: . numCells - The number of cells
1456: . numVertices - The number of vertices
1457: . numCorners - The number of vertices for each cell
1458: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
1459: . cells - An array of numCells*numCorners numbers, the vertices for each cell
1460: . spaceDim - The spatial dimension used for coordinates
1461: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
1463: Output Parameter:
1464: . dm - The DM
1466: Note: Two triangles sharing a face
1467: $
1468: $ 2
1469: $ / | \
1470: $ / | \
1471: $ / | \
1472: $ 0 0 | 1 3
1473: $ \ | /
1474: $ \ | /
1475: $ \ | /
1476: $ 1
1477: would have input
1478: $ numCells = 2, numVertices = 4
1479: $ cells = [0 1 2 1 3 2]
1480: $
1481: which would result in the DMPlex
1482: $
1483: $ 4
1484: $ / | \
1485: $ / | \
1486: $ / | \
1487: $ 2 0 | 1 5
1488: $ \ | /
1489: $ \ | /
1490: $ \ | /
1491: $ 3
1493: Level: beginner
1495: .seealso: DMPlexCreateFromDAG(), DMPlexCreate()
1496: @*/
1497: 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)
1498: {
1502: DMCreate(comm, dm);
1503: DMSetType(*dm, DMPLEX);
1504: DMSetDimension(*dm, dim);
1505: DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);
1506: if (interpolate) {
1507: DM idm = NULL;
1509: DMPlexInterpolate(*dm, &idm);
1510: DMDestroy(dm);
1511: *dm = idm;
1512: }
1513: DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);
1514: return(0);
1515: }
1519: /*@
1520: DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM
1522: Input Parameters:
1523: + dm - The empty DM object, usually from DMCreate() and DMSetDimension()
1524: . depth - The depth of the DAG
1525: . numPoints - The number of points at each depth
1526: . coneSize - The cone size of each point
1527: . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
1528: . coneOrientations - The orientation of each cone point
1529: - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex
1531: Output Parameter:
1532: . dm - The DM
1534: Note: Two triangles sharing a face would have input
1535: $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
1536: $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0]
1537: $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0]
1538: $
1539: which would result in the DMPlex
1540: $
1541: $ 4
1542: $ / | \
1543: $ / | \
1544: $ / | \
1545: $ 2 0 | 1 5
1546: $ \ | /
1547: $ \ | /
1548: $ \ | /
1549: $ 3
1550: $
1551: $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList()
1553: Level: advanced
1555: .seealso: DMPlexCreateFromCellList(), DMPlexCreate()
1556: @*/
1557: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
1558: {
1559: Vec coordinates;
1560: PetscSection coordSection;
1561: PetscScalar *coords;
1562: PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off;
1566: DMGetDimension(dm, &dim);
1567: DMGetCoordinateDim(dm, &dimEmbed);
1568: if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim);
1569: for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
1570: DMPlexSetChart(dm, pStart, pEnd);
1571: for (p = pStart; p < pEnd; ++p) {
1572: DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
1573: if (firstVertex < 0 && !coneSize[p - pStart]) {
1574: firstVertex = p - pStart;
1575: }
1576: }
1577: if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]);
1578: DMSetUp(dm); /* Allocate space for cones */
1579: for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
1580: DMPlexSetCone(dm, p, &cones[off]);
1581: DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
1582: }
1583: DMPlexSymmetrize(dm);
1584: DMPlexStratify(dm);
1585: /* Build coordinates */
1586: DMGetCoordinateSection(dm, &coordSection);
1587: PetscSectionSetNumFields(coordSection, 1);
1588: PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
1589: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
1590: for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
1591: PetscSectionSetDof(coordSection, v, dimEmbed);
1592: PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
1593: }
1594: PetscSectionSetUp(coordSection);
1595: PetscSectionGetStorageSize(coordSection, &coordSize);
1596: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1597: VecSetBlockSize(coordinates, dimEmbed);
1598: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1599: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1600: VecSetType(coordinates,VECSTANDARD);
1601: VecGetArray(coordinates, &coords);
1602: for (v = 0; v < numPoints[0]; ++v) {
1603: PetscInt off;
1605: PetscSectionGetOffset(coordSection, v+firstVertex, &off);
1606: for (d = 0; d < dimEmbed; ++d) {
1607: coords[off+d] = vertexCoords[v*dimEmbed+d];
1608: }
1609: }
1610: VecRestoreArray(coordinates, &coords);
1611: DMSetCoordinatesLocal(dm, coordinates);
1612: VecDestroy(&coordinates);
1613: return(0);
1614: }
1618: /*@C
1619: DMPlexCreateFromFile - This takes a filename and produces a DM
1621: Input Parameters:
1622: + comm - The communicator
1623: . filename - A file name
1624: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
1626: Output Parameter:
1627: . dm - The DM
1629: Level: beginner
1631: .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate()
1632: @*/
1633: PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1634: {
1635: const char *extGmsh = ".msh";
1636: const char *extCGNS = ".cgns";
1637: const char *extExodus = ".exo";
1638: const char *extFluent = ".cas";
1639: size_t len;
1640: PetscBool isGmsh, isCGNS, isExodus, isFluent;
1641: PetscMPIInt rank;
1647: MPI_Comm_rank(comm, &rank);
1648: PetscStrlen(filename, &len);
1649: if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path");
1650: PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);
1651: PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);
1652: PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);
1653: PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);
1654: if (isGmsh) {
1655: DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);
1656: } else if (isCGNS) {
1657: DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);
1658: } else if (isExodus) {
1659: DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);
1660: } else if (isFluent) {
1661: DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);
1662: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename);
1663: return(0);
1664: }
1668: /*@
1669: DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
1671: Collective on comm
1673: Input Parameters:
1674: + comm - The communicator
1675: . dim - The spatial dimension
1676: - simplex - Flag for simplex, otherwise use a tensor-product cell
1678: Output Parameter:
1679: . refdm - The reference cell
1681: Level: intermediate
1683: .keywords: reference cell
1684: .seealso:
1685: @*/
1686: PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm)
1687: {
1688: DM rdm;
1692: DMCreate(comm, &rdm);
1693: DMSetType(rdm, DMPLEX);
1694: DMSetDimension(rdm, dim);
1695: switch (dim) {
1696: case 0:
1697: {
1698: PetscInt numPoints[1] = {1};
1699: PetscInt coneSize[1] = {0};
1700: PetscInt cones[1] = {0};
1701: PetscInt coneOrientations[1] = {0};
1702: PetscScalar vertexCoords[1] = {0.0};
1704: DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1705: }
1706: break;
1707: case 1:
1708: {
1709: PetscInt numPoints[2] = {2, 1};
1710: PetscInt coneSize[3] = {2, 0, 0};
1711: PetscInt cones[2] = {1, 2};
1712: PetscInt coneOrientations[2] = {0, 0};
1713: PetscScalar vertexCoords[2] = {-1.0, 1.0};
1715: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1716: }
1717: break;
1718: case 2:
1719: if (simplex) {
1720: PetscInt numPoints[2] = {3, 1};
1721: PetscInt coneSize[4] = {3, 0, 0, 0};
1722: PetscInt cones[3] = {1, 2, 3};
1723: PetscInt coneOrientations[3] = {0, 0, 0};
1724: PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0};
1726: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1727: } else {
1728: PetscInt numPoints[2] = {4, 1};
1729: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
1730: PetscInt cones[4] = {1, 2, 3, 4};
1731: PetscInt coneOrientations[4] = {0, 0, 0, 0};
1732: PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
1734: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1735: }
1736: break;
1737: case 3:
1738: if (simplex) {
1739: PetscInt numPoints[2] = {4, 1};
1740: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
1741: PetscInt cones[4] = {1, 3, 2, 4};
1742: PetscInt coneOrientations[4] = {0, 0, 0, 0};
1743: 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};
1745: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1746: } else {
1747: PetscInt numPoints[2] = {8, 1};
1748: PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0};
1749: PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8};
1750: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
1751: 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,
1752: -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
1754: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1755: }
1756: break;
1757: default:
1758: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim);
1759: }
1760: *refdm = NULL;
1761: DMPlexInterpolate(rdm, refdm);
1762: DMPlexCopyCoordinates(rdm, *refdm);
1763: DMDestroy(&rdm);
1764: return(0);
1765: }