Actual source code: plexcreate.c
petsc-3.4.5 2014-06-29
1: #define PETSCDM_DLL
2: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
3: #include <petscdmda.h>
4: #include <petscsf.h>
8: PetscErrorCode DMSetFromOptions_Plex(DM dm)
9: {
10: DM_Plex *mesh = (DM_Plex*) dm->data;
15: PetscOptionsHead("DMPlex Options");
16: /* Handle DMPlex refinement */
17: /* Handle associated vectors */
18: /* Handle viewing */
19: PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);
20: PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);
21: PetscOptionsTail();
22: return(0);
23: }
27: /*
28: Simple square boundary:
30: 18--5-17--4--16
31: | | |
32: 6 10 3
33: | | |
34: 19-11-20--9--15
35: | | |
36: 7 8 2
37: | | |
38: 12--0-13--1--14
39: */
40: PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
41: {
42: PetscInt numVertices = (edges[0]+1)*(edges[1]+1);
43: PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
44: PetscInt markerTop = 1;
45: PetscInt markerBottom = 1;
46: PetscInt markerRight = 1;
47: PetscInt markerLeft = 1;
48: PetscBool markerSeparate = PETSC_FALSE;
49: Vec coordinates;
50: PetscSection coordSection;
51: PetscScalar *coords;
52: PetscInt coordSize;
53: PetscMPIInt rank;
54: PetscInt v, vx, vy;
58: PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
59: if (markerSeparate) {
60: markerTop = 1;
61: markerBottom = 0;
62: markerRight = 0;
63: markerLeft = 0;
64: }
65: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
66: if (!rank) {
67: PetscInt e, ex, ey;
69: DMPlexSetChart(dm, 0, numEdges+numVertices);
70: for (e = 0; e < numEdges; ++e) {
71: DMPlexSetConeSize(dm, e, 2);
72: }
73: DMSetUp(dm); /* Allocate space for cones */
74: for (vx = 0; vx <= edges[0]; vx++) {
75: for (ey = 0; ey < edges[1]; ey++) {
76: PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1);
77: PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
78: PetscInt cone[2];
80: cone[0] = vertex; cone[1] = vertex+edges[0]+1;
81: DMPlexSetCone(dm, edge, cone);
82: if (vx == edges[0]) {
83: DMPlexSetLabelValue(dm, "marker", edge, markerRight);
84: DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
85: if (ey == edges[1]-1) {
86: DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
87: }
88: } else if (vx == 0) {
89: DMPlexSetLabelValue(dm, "marker", edge, markerLeft);
90: DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
91: if (ey == edges[1]-1) {
92: DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
93: }
94: }
95: }
96: }
97: for (vy = 0; vy <= edges[1]; vy++) {
98: for (ex = 0; ex < edges[0]; ex++) {
99: PetscInt edge = vy*edges[0] + ex;
100: PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
101: PetscInt cone[2];
103: cone[0] = vertex; cone[1] = vertex+1;
104: DMPlexSetCone(dm, edge, cone);
105: if (vy == edges[1]) {
106: DMPlexSetLabelValue(dm, "marker", edge, markerTop);
107: DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
108: if (ex == edges[0]-1) {
109: DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
110: }
111: } else if (vy == 0) {
112: DMPlexSetLabelValue(dm, "marker", edge, markerBottom);
113: DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
114: if (ex == edges[0]-1) {
115: DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
116: }
117: }
118: }
119: }
120: }
121: DMPlexSymmetrize(dm);
122: DMPlexStratify(dm);
123: /* Build coordinates */
124: DMPlexGetCoordinateSection(dm, &coordSection);
125: PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);
126: for (v = numEdges; v < numEdges+numVertices; ++v) {
127: PetscSectionSetDof(coordSection, v, 2);
128: }
129: PetscSectionSetUp(coordSection);
130: PetscSectionGetStorageSize(coordSection, &coordSize);
131: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
132: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
133: VecSetFromOptions(coordinates);
134: VecGetArray(coordinates, &coords);
135: for (vy = 0; vy <= edges[1]; ++vy) {
136: for (vx = 0; vx <= edges[0]; ++vx) {
137: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
138: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
139: }
140: }
141: VecRestoreArray(coordinates, &coords);
142: DMSetCoordinatesLocal(dm, coordinates);
143: VecDestroy(&coordinates);
144: return(0);
145: }
149: /*
150: Simple cubic boundary:
152: 2-------3
153: /| /|
154: 6-------7 |
155: | | | |
156: | 0-----|-1
157: |/ |/
158: 4-------5
159: */
160: PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
161: {
162: PetscInt numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
163: PetscInt numFaces = 6;
164: Vec coordinates;
165: PetscSection coordSection;
166: PetscScalar *coords;
167: PetscInt coordSize;
168: PetscMPIInt rank;
169: PetscInt v, vx, vy, vz;
173: 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");
174: if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Currently can't handle more than 1 face per side");
175: PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);
176: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
177: if (!rank) {
178: PetscInt f;
180: DMPlexSetChart(dm, 0, numFaces+numVertices);
181: for (f = 0; f < numFaces; ++f) {
182: DMPlexSetConeSize(dm, f, 4);
183: }
184: DMSetUp(dm); /* Allocate space for cones */
185: for (v = 0; v < numFaces+numVertices; ++v) {
186: DMPlexSetLabelValue(dm, "marker", v, 1);
187: }
188: { /* Side 0 (Front) */
189: PetscInt cone[4];
190: cone[0] = numFaces+4; cone[1] = numFaces+5; cone[2] = numFaces+7; cone[3] = numFaces+6;
191: DMPlexSetCone(dm, 0, cone);
192: }
193: { /* Side 1 (Back) */
194: PetscInt cone[4];
195: cone[0] = numFaces+1; cone[1] = numFaces+0; cone[2] = numFaces+2; cone[3] = numFaces+3;
196: DMPlexSetCone(dm, 1, cone);
197: }
198: { /* Side 2 (Bottom) */
199: PetscInt cone[4];
200: cone[0] = numFaces+0; cone[1] = numFaces+1; cone[2] = numFaces+5; cone[3] = numFaces+4;
201: DMPlexSetCone(dm, 2, cone);
202: }
203: { /* Side 3 (Top) */
204: PetscInt cone[4];
205: cone[0] = numFaces+6; cone[1] = numFaces+7; cone[2] = numFaces+3; cone[3] = numFaces+2;
206: DMPlexSetCone(dm, 3, cone);
207: }
208: { /* Side 4 (Left) */
209: PetscInt cone[4];
210: cone[0] = numFaces+0; cone[1] = numFaces+4; cone[2] = numFaces+6; cone[3] = numFaces+2;
211: DMPlexSetCone(dm, 4, cone);
212: }
213: { /* Side 5 (Right) */
214: PetscInt cone[4];
215: cone[0] = numFaces+5; cone[1] = numFaces+1; cone[2] = numFaces+3; cone[3] = numFaces+7;
216: DMPlexSetCone(dm, 5, cone);
217: }
218: }
219: DMPlexSymmetrize(dm);
220: DMPlexStratify(dm);
221: /* Build coordinates */
222: DMPlexGetCoordinateSection(dm, &coordSection);
223: PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);
224: for (v = numFaces; v < numFaces+numVertices; ++v) {
225: PetscSectionSetDof(coordSection, v, 3);
226: }
227: PetscSectionSetUp(coordSection);
228: PetscSectionGetStorageSize(coordSection, &coordSize);
229: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
230: PetscObjectSetName((PetscObject) coordinates, "coordinates");
231: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
232: VecSetFromOptions(coordinates);
233: VecGetArray(coordinates, &coords);
234: for (vz = 0; vz <= faces[2]; ++vz) {
235: for (vy = 0; vy <= faces[1]; ++vy) {
236: for (vx = 0; vx <= faces[0]; ++vx) {
237: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
238: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
239: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
240: }
241: }
242: }
243: VecRestoreArray(coordinates, &coords);
244: DMSetCoordinatesLocal(dm, coordinates);
245: VecDestroy(&coordinates);
246: return(0);
247: }
251: /*
252: Simple square mesh:
254: 22--8-23--9--24
255: | | |
256: 13 2 14 3 15
257: | | |
258: 19--6-20--7--21
259: | | |
260: 10 0 11 1 12
261: | | |
262: 16--4-17--5--18
263: */
264: PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
265: {
266: PetscInt markerTop = 1;
267: PetscInt markerBottom = 1;
268: PetscInt markerRight = 1;
269: PetscInt markerLeft = 1;
270: PetscBool markerSeparate = PETSC_FALSE;
271: PetscMPIInt rank;
275: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
276: PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
277: if (markerSeparate) {
278: markerTop = 3;
279: markerBottom = 1;
280: markerRight = 2;
281: markerLeft = 4;
282: }
283: {
284: const PetscInt numXEdges = !rank ? edges[0] : 0;
285: const PetscInt numYEdges = !rank ? edges[1] : 0;
286: const PetscInt numXVertices = !rank ? edges[0]+1 : 0;
287: const PetscInt numYVertices = !rank ? edges[1]+1 : 0;
288: const PetscInt numTotXEdges = numXEdges*numYVertices;
289: const PetscInt numTotYEdges = numYEdges*numXVertices;
290: const PetscInt numVertices = numXVertices*numYVertices;
291: const PetscInt numEdges = numTotXEdges + numTotYEdges;
292: const PetscInt numFaces = numXEdges*numYEdges;
293: const PetscInt firstVertex = numFaces;
294: const PetscInt firstXEdge = numFaces + numVertices;
295: const PetscInt firstYEdge = numFaces + numVertices + numTotXEdges;
296: Vec coordinates;
297: PetscSection coordSection;
298: PetscScalar *coords;
299: PetscInt coordSize;
300: PetscInt v, vx, vy;
301: PetscInt f, fx, fy, e, ex, ey;
303: DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);
304: for (f = 0; f < numFaces; ++f) {
305: DMPlexSetConeSize(dm, f, 4);
306: }
307: for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
308: DMPlexSetConeSize(dm, e, 2);
309: }
310: DMSetUp(dm); /* Allocate space for cones */
311: /* Build faces */
312: for (fy = 0; fy < numYEdges; fy++) {
313: for (fx = 0; fx < numXEdges; fx++) {
314: const PetscInt face = fy*numXEdges + fx;
315: const PetscInt edgeL = firstYEdge + fx*numYEdges + fy;
316: const PetscInt edgeB = firstXEdge + fy*numXEdges + fx;
317: const PetscInt ornt[4] = {0, 0, -2, -2};
318: PetscInt cone[4];
320: cone[0] = edgeB; cone[1] = edgeL+numYEdges; cone[2] = edgeB+numXEdges; cone[3] = edgeL;
321: DMPlexSetCone(dm, face, cone);
322: DMPlexSetConeOrientation(dm, face, ornt);
323: }
324: }
325: /* Build Y edges*/
326: for (vx = 0; vx < numXVertices; vx++) {
327: for (ey = 0; ey < numYEdges; ey++) {
328: const PetscInt edge = firstYEdge + vx*numYEdges + ey;
329: const PetscInt vertex = firstVertex + ey*numXVertices + vx;
330: PetscInt cone[2];
332: cone[0] = vertex; cone[1] = vertex+numXVertices;
333: DMPlexSetCone(dm, edge, cone);
334: if (vx == numXVertices-1) {
335: DMPlexSetLabelValue(dm, "marker", edge, markerRight);
336: DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);
337: if (ey == numYEdges-1) {
338: DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);
339: }
340: } else if (vx == 0) {
341: DMPlexSetLabelValue(dm, "marker", edge, markerLeft);
342: DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);
343: if (ey == numYEdges-1) {
344: DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);
345: }
346: }
347: }
348: }
349: /* Build X edges*/
350: for (vy = 0; vy < numYVertices; vy++) {
351: for (ex = 0; ex < numXEdges; ex++) {
352: const PetscInt edge = firstXEdge + vy*numXEdges + ex;
353: const PetscInt vertex = firstVertex + vy*numXVertices + ex;
354: PetscInt cone[2];
356: cone[0] = vertex; cone[1] = vertex+1;
357: DMPlexSetCone(dm, edge, cone);
358: if (vy == numYVertices-1) {
359: DMPlexSetLabelValue(dm, "marker", edge, markerTop);
360: DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);
361: if (ex == numXEdges-1) {
362: DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);
363: }
364: } else if (vy == 0) {
365: DMPlexSetLabelValue(dm, "marker", edge, markerBottom);
366: DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);
367: if (ex == numXEdges-1) {
368: DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);
369: }
370: }
371: }
372: }
373: DMPlexSymmetrize(dm);
374: DMPlexStratify(dm);
375: /* Build coordinates */
376: DMPlexGetCoordinateSection(dm, &coordSection);
377: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);
378: for (v = firstVertex; v < firstVertex+numVertices; ++v) {
379: PetscSectionSetDof(coordSection, v, 2);
380: }
381: PetscSectionSetUp(coordSection);
382: PetscSectionGetStorageSize(coordSection, &coordSize);
383: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
384: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
385: VecSetFromOptions(coordinates);
386: VecGetArray(coordinates, &coords);
387: for (vy = 0; vy < numYVertices; ++vy) {
388: for (vx = 0; vx < numXVertices; ++vx) {
389: coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
390: coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
391: }
392: }
393: VecRestoreArray(coordinates, &coords);
394: DMSetCoordinatesLocal(dm, coordinates);
395: VecDestroy(&coordinates);
396: }
397: return(0);
398: }
402: /*@
403: DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box).
405: Collective on MPI_Comm
407: Input Parameters:
408: + comm - The communicator for the DM object
409: . dim - The spatial dimension
410: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
412: Output Parameter:
413: . dm - The DM object
415: Level: beginner
417: .keywords: DM, create
418: .seealso: DMSetType(), DMCreate()
419: @*/
420: PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm)
421: {
422: DM boundary;
427: DMCreate(comm, &boundary);
429: DMSetType(boundary, DMPLEX);
430: DMPlexSetDimension(boundary, dim-1);
431: switch (dim) {
432: case 2:
433: {
434: PetscReal lower[2] = {0.0, 0.0};
435: PetscReal upper[2] = {1.0, 1.0};
436: PetscInt edges[2] = {2, 2};
438: DMPlexCreateSquareBoundary(boundary, lower, upper, edges);
439: break;
440: }
441: case 3:
442: {
443: PetscReal lower[3] = {0.0, 0.0, 0.0};
444: PetscReal upper[3] = {1.0, 1.0, 1.0};
445: PetscInt faces[3] = {1, 1, 1};
447: DMPlexCreateCubeBoundary(boundary, lower, upper, faces);
448: break;
449: }
450: default:
451: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
452: }
453: DMPlexGenerate(boundary, NULL, interpolate, dm);
454: DMDestroy(&boundary);
455: return(0);
456: }
460: PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm)
461: {
466: DMCreate(comm, dm);
468: DMSetType(*dm, DMPLEX);
469: DMPlexSetDimension(*dm, dim);
470: switch (dim) {
471: case 2:
472: {
473: PetscReal lower[2] = {0.0, 0.0};
474: PetscReal upper[2] = {1.0, 1.0};
476: DMPlexCreateSquareMesh(*dm, lower, upper, cells);
477: break;
478: }
479: #if 0
480: case 3:
481: {
482: PetscReal lower[3] = {0.0, 0.0, 0.0};
483: PetscReal upper[3] = {1.0, 1.0, 1.0};
485: DMPlexCreateCubeMesh(boundary, lower, upper, cells);
486: break;
487: }
488: #endif
489: default:
490: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
491: }
492: return(0);
493: }
495: /* External function declarations here */
496: extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
497: extern PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J);
498: extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
499: extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined);
500: extern PetscErrorCode DMSetUp_Plex(DM dm);
501: extern PetscErrorCode DMDestroy_Plex(DM dm);
502: extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
503: extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
504: extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS);
508: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
509: {
513: DMCreateGlobalVector_Section_Private(dm,vec);
514: /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
515: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);
516: return(0);
517: }
521: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
522: {
526: DMCreateLocalVector_Section_Private(dm,vec);
527: VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);
528: return(0);
529: }
534: PetscErrorCode DMInitialize_Plex(DM dm)
535: {
539: PetscStrallocpy(VECSTANDARD, (char**)&dm->vectype);
541: dm->ops->view = DMView_Plex;
542: dm->ops->setfromoptions = DMSetFromOptions_Plex;
543: dm->ops->setup = DMSetUp_Plex;
544: dm->ops->createglobalvector = DMCreateGlobalVector_Plex;
545: dm->ops->createlocalvector = DMCreateLocalVector_Plex;
546: dm->ops->getlocaltoglobalmapping = NULL;
547: dm->ops->getlocaltoglobalmappingblock = NULL;
548: dm->ops->createfieldis = NULL;
549: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex;
550: dm->ops->getcoloring = 0;
551: dm->ops->creatematrix = DMCreateMatrix_Plex;
552: dm->ops->createinterpolation = 0;
553: dm->ops->getaggregates = 0;
554: dm->ops->getinjection = 0;
555: dm->ops->refine = DMRefine_Plex;
556: dm->ops->coarsen = 0;
557: dm->ops->refinehierarchy = 0;
558: dm->ops->coarsenhierarchy = 0;
559: dm->ops->globaltolocalbegin = NULL;
560: dm->ops->globaltolocalend = NULL;
561: dm->ops->localtoglobalbegin = NULL;
562: dm->ops->localtoglobalend = NULL;
563: dm->ops->destroy = DMDestroy_Plex;
564: dm->ops->createsubdm = DMCreateSubDM_Plex;
565: dm->ops->locatepoints = DMLocatePoints_Plex;
566: return(0);
567: }
569: /*MC
570: DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram.
571: In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is
572: specified by a PetscSection object. Ownership in the global representation is determined by
573: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
575: Level: intermediate
577: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
578: M*/
582: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
583: {
584: DM_Plex *mesh;
585: PetscInt unit, d;
590: PetscNewLog(dm, DM_Plex, &mesh);
591: dm->data = mesh;
593: mesh->refct = 1;
594: mesh->dim = 0;
595: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
596: mesh->maxConeSize = 0;
597: mesh->cones = NULL;
598: mesh->coneOrientations = NULL;
599: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
600: mesh->maxSupportSize = 0;
601: mesh->supports = NULL;
602: mesh->refinementUniform = PETSC_TRUE;
603: mesh->refinementLimit = -1.0;
605: mesh->facesTmp = NULL;
607: mesh->subpointMap = NULL;
609: for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0;
611: mesh->labels = NULL;
612: mesh->globalVertexNumbers = NULL;
613: mesh->globalCellNumbers = NULL;
614: for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE;
615: mesh->vtkCellHeight = 0;
616: mesh->preallocCenterDim = -1;
618: mesh->integrateResidualFEM = NULL;
619: mesh->integrateJacobianActionFEM = NULL;
620: mesh->integrateJacobianFEM = NULL;
622: mesh->printSetValues = PETSC_FALSE;
623: mesh->printFEM = 0;
625: DMInitialize_Plex(dm);
626: return(0);
627: }
631: /*@
632: DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.
634: Collective on MPI_Comm
636: Input Parameter:
637: . comm - The communicator for the DMPlex object
639: Output Parameter:
640: . mesh - The DMPlex object
642: Level: beginner
644: .keywords: DMPlex, create
645: @*/
646: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
647: {
652: DMCreate(comm, mesh);
653: DMSetType(*mesh, DMPLEX);
654: return(0);
655: }
659: /*@
660: DMPlexClone - Creates a DMPlex object with the same mesh as the original.
662: Collective on MPI_Comm
664: Input Parameter:
665: . dm - The original DMPlex object
667: Output Parameter:
668: . newdm - The new DMPlex object
670: Level: beginner
672: .keywords: DMPlex, create
673: @*/
674: PetscErrorCode DMPlexClone(DM dm, DM *newdm)
675: {
676: DM_Plex *mesh;
677: Vec coords;
678: void *ctx;
684: DMCreate(PetscObjectComm((PetscObject)dm), newdm);
685: PetscSFDestroy(&(*newdm)->sf);
686: PetscObjectReference((PetscObject) dm->sf);
687: (*newdm)->sf = dm->sf;
688: mesh = (DM_Plex*) dm->data;
689: mesh->refct++;
690: (*newdm)->data = mesh;
691: PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
692: DMInitialize_Plex(*newdm);
693: DMGetApplicationContext(dm, &ctx);
694: DMSetApplicationContext(*newdm, ctx);
695: DMGetCoordinatesLocal(dm, &coords);
696: if (coords) {
697: DMSetCoordinatesLocal(*newdm, coords);
698: } else {
699: DMGetCoordinates(dm, &coords);
700: if (coords) {DMSetCoordinates(*newdm, coords);}
701: }
702: return(0);
703: }
707: /*
708: This takes as input the common mesh generator output, a list of the vertices for each cell
709: */
710: PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[])
711: {
712: PetscInt *cone, c, p;
716: DMPlexSetChart(dm, 0, numCells+numVertices);
717: for (c = 0; c < numCells; ++c) {
718: DMPlexSetConeSize(dm, c, numCorners);
719: }
720: DMSetUp(dm);
721: DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);
722: for (c = 0; c < numCells; ++c) {
723: for (p = 0; p < numCorners; ++p) {
724: cone[p] = cells[c*numCorners+p]+numCells;
725: }
726: DMPlexSetCone(dm, c, cone);
727: }
728: DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);
729: DMPlexSymmetrize(dm);
730: DMPlexStratify(dm);
731: return(0);
732: }
736: /*
737: This takes as input the coordinates for each vertex
738: */
739: PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[])
740: {
741: PetscSection coordSection;
742: Vec coordinates;
743: PetscScalar *coords;
744: PetscInt coordSize, v, d;
748: DMPlexGetCoordinateSection(dm, &coordSection);
749: PetscSectionSetNumFields(coordSection, 1);
750: PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
751: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
752: for (v = numCells; v < numCells+numVertices; ++v) {
753: PetscSectionSetDof(coordSection, v, spaceDim);
754: PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
755: }
756: PetscSectionSetUp(coordSection);
757: PetscSectionGetStorageSize(coordSection, &coordSize);
758: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
759: PetscObjectSetName((PetscObject) coordinates, "coordinates");
760: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
761: VecSetFromOptions(coordinates);
762: VecGetArray(coordinates, &coords);
763: for (v = 0; v < numVertices; ++v) {
764: for (d = 0; d < spaceDim; ++d) {
765: coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
766: }
767: }
768: VecRestoreArray(coordinates, &coords);
769: DMSetCoordinatesLocal(dm, coordinates);
770: VecDestroy(&coordinates);
771: return(0);
772: }
776: /*@C
777: DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM
779: Input Parameters:
780: + comm - The communicator
781: . dim - The topological dimension of the mesh
782: . numCells - The number of cells
783: . numVertices - The number of vertices
784: . numCorners - The number of vertices for each cell
785: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
786: . cells - An array of numCells*numCorners numbers, the vertices for each cell
787: . spaceDim - The spatial dimension used for coordinates
788: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
790: Output Parameter:
791: . dm - The DM
793: Note: Two triangles sharing a face
794: $
795: $ 2
796: $ / | \
797: $ / | \
798: $ / | \
799: $ 0 0 | 1 3
800: $ \ | /
801: $ \ | /
802: $ \ | /
803: $ 1
804: would have input
805: $ numCells = 2, numVertices = 4
806: $ cells = [0 1 2 1 3 2]
807: $
808: which would result in the DMPlex
809: $
810: $ 4
811: $ / | \
812: $ / | \
813: $ / | \
814: $ 2 0 | 1 5
815: $ \ | /
816: $ \ | /
817: $ \ | /
818: $ 3
820: Level: beginner
822: .seealso: DMPlexCreate()
823: @*/
824: 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)
825: {
829: DMCreate(comm, dm);
830: DMSetType(*dm, DMPLEX);
831: DMPlexSetDimension(*dm, dim);
832: DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);
833: if (interpolate) {
834: DM idm;
836: DMPlexInterpolate(*dm, &idm);
837: DMDestroy(dm);
838: *dm = idm;
839: }
840: DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);
841: return(0);
842: }
846: /*
847: This takes as input the raw Hasse Diagram data
848: */
849: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
850: {
851: Vec coordinates;
852: PetscSection coordSection;
853: PetscScalar *coords;
854: PetscInt coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off;
858: DMPlexGetDimension(dm, &dim);
859: for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
860: DMPlexSetChart(dm, pStart, pEnd);
861: for (p = pStart; p < pEnd; ++p) {
862: DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
863: }
864: DMSetUp(dm); /* Allocate space for cones */
865: for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
866: DMPlexSetCone(dm, p, &cones[off]);
867: DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
868: }
869: DMPlexSymmetrize(dm);
870: DMPlexStratify(dm);
871: /* Build coordinates */
872: DMPlexGetCoordinateSection(dm, &coordSection);
873: PetscSectionSetNumFields(coordSection, 1);
874: PetscSectionSetFieldComponents(coordSection, 0, dim);
875: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
876: for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
877: PetscSectionSetDof(coordSection, v, dim);
878: PetscSectionSetFieldDof(coordSection, v, 0, dim);
879: }
880: PetscSectionSetUp(coordSection);
881: PetscSectionGetStorageSize(coordSection, &coordSize);
882: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
883: PetscObjectSetName((PetscObject) coordinates, "coordinates");
884: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
885: VecSetFromOptions(coordinates);
886: VecGetArray(coordinates, &coords);
887: for (v = 0; v < numPoints[0]; ++v) {
888: PetscInt off;
890: PetscSectionGetOffset(coordSection, v+firstVertex, &off);
891: for (d = 0; d < dim; ++d) {
892: coords[off+d] = vertexCoords[v*dim+d];
893: }
894: }
895: VecRestoreArray(coordinates, &coords);
896: DMSetCoordinatesLocal(dm, coordinates);
897: VecDestroy(&coordinates);
898: return(0);
899: }