Actual source code: plexcreate.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h>
3: #include <petsc/private/hashseti.h>
4: #include <petscsf.h>
6: PetscLogEvent DMPLEX_CreateFromFile, DMPLEX_BuildFromCellList, DMPLEX_BuildCoordinatesFromCellList;
8: /* External function declarations here */
9: static PetscErrorCode DMInitialize_Plex(DM dm);
11: /* Replace dm with the contents of ndm, and then destroy ndm
12: - Share the DM_Plex structure
13: - Share the coordinates
14: - Share the SF
15: */
16: static PetscErrorCode DMPlexReplace_Static(DM dm, DM *ndm)
17: {
18: PetscSF sf;
19: DM dmNew = *ndm, coordDM, coarseDM;
20: Vec coords;
21: PetscBool isper;
22: const PetscReal *maxCell, *L;
23: const DMBoundaryType *bd;
24: PetscInt dim, cdim;
25: PetscErrorCode ierr;
28: if (dm == dmNew) {
29: DMDestroy(ndm);
30: return(0);
31: }
32: dm->setupcalled = dmNew->setupcalled;
33: DMGetDimension(dmNew, &dim);
34: DMSetDimension(dm, dim);
35: DMGetCoordinateDim(dmNew, &cdim);
36: DMSetCoordinateDim(dm, cdim);
37: DMGetPointSF(dmNew, &sf);
38: DMSetPointSF(dm, sf);
39: DMGetCoordinateDM(dmNew, &coordDM);
40: DMGetCoordinatesLocal(dmNew, &coords);
41: DMSetCoordinateDM(dm, coordDM);
42: DMSetCoordinatesLocal(dm, coords);
43: /* Do not want to create the coordinate field if it does not already exist, so do not call DMGetCoordinateField() */
44: DMFieldDestroy(&dm->coordinateField);
45: dm->coordinateField = dmNew->coordinateField;
46: ((DM_Plex *) dmNew->data)->coordFunc = ((DM_Plex *) dm->data)->coordFunc;
47: DMGetPeriodicity(dmNew, &isper, &maxCell, &L, &bd);
48: DMSetPeriodicity(dm, isper, maxCell, L, bd);
49: DMDestroy_Plex(dm);
50: DMInitialize_Plex(dm);
51: dm->data = dmNew->data;
52: ((DM_Plex *) dmNew->data)->refct++;
53: DMDestroyLabelLinkList_Internal(dm);
54: DMCopyLabels(dmNew, dm, PETSC_OWN_POINTER, PETSC_TRUE);
55: DMGetCoarseDM(dmNew,&coarseDM);
56: DMSetCoarseDM(dm,coarseDM);
57: DMDestroy(ndm);
58: return(0);
59: }
61: /* Swap dm with the contents of dmNew
62: - Swap the DM_Plex structure
63: - Swap the coordinates
64: - Swap the point PetscSF
65: */
66: static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB)
67: {
68: DM coordDMA, coordDMB;
69: Vec coordsA, coordsB;
70: PetscSF sfA, sfB;
71: DMField fieldTmp;
72: void *tmp;
73: DMLabelLink listTmp;
74: DMLabel depthTmp;
75: PetscInt tmpI;
76: PetscErrorCode ierr;
79: if (dmA == dmB) return(0);
80: DMGetPointSF(dmA, &sfA);
81: DMGetPointSF(dmB, &sfB);
82: PetscObjectReference((PetscObject) sfA);
83: DMSetPointSF(dmA, sfB);
84: DMSetPointSF(dmB, sfA);
85: PetscObjectDereference((PetscObject) sfA);
87: DMGetCoordinateDM(dmA, &coordDMA);
88: DMGetCoordinateDM(dmB, &coordDMB);
89: PetscObjectReference((PetscObject) coordDMA);
90: DMSetCoordinateDM(dmA, coordDMB);
91: DMSetCoordinateDM(dmB, coordDMA);
92: PetscObjectDereference((PetscObject) coordDMA);
94: DMGetCoordinatesLocal(dmA, &coordsA);
95: DMGetCoordinatesLocal(dmB, &coordsB);
96: PetscObjectReference((PetscObject) coordsA);
97: DMSetCoordinatesLocal(dmA, coordsB);
98: DMSetCoordinatesLocal(dmB, coordsA);
99: PetscObjectDereference((PetscObject) coordsA);
101: fieldTmp = dmA->coordinateField;
102: dmA->coordinateField = dmB->coordinateField;
103: dmB->coordinateField = fieldTmp;
104: tmp = dmA->data;
105: dmA->data = dmB->data;
106: dmB->data = tmp;
107: listTmp = dmA->labels;
108: dmA->labels = dmB->labels;
109: dmB->labels = listTmp;
110: depthTmp = dmA->depthLabel;
111: dmA->depthLabel = dmB->depthLabel;
112: dmB->depthLabel = depthTmp;
113: depthTmp = dmA->celltypeLabel;
114: dmA->celltypeLabel = dmB->celltypeLabel;
115: dmB->celltypeLabel = depthTmp;
116: tmpI = dmA->levelup;
117: dmA->levelup = dmB->levelup;
118: dmB->levelup = tmpI;
119: return(0);
120: }
122: static PetscErrorCode DMPlexInterpolateInPlace_Internal(DM dm)
123: {
124: DM idm;
128: DMPlexInterpolate(dm, &idm);
129: DMPlexCopyCoordinates(dm, idm);
130: DMPlexReplace_Static(dm, &idm);
131: return(0);
132: }
134: /*@C
135: DMPlexCreateCoordinateSpace - Creates a finite element space for the coordinates
137: Collective
139: Input Parameters:
140: + DM - The DM
141: . degree - The degree of the finite element
142: - coordFunc - An optional function to map new points from refinement to the surface
144: Level: advanced
146: .seealso: PetscFECreateLagrange(), DMGetCoordinateDM()
147: @*/
148: PetscErrorCode DMPlexCreateCoordinateSpace(DM dm, PetscInt degree, PetscPointFunc coordFunc)
149: {
150: DM_Plex *mesh = (DM_Plex *) dm->data;
151: DM cdm;
152: PetscDS cds;
153: PetscFE fe;
154: PetscClassId id;
158: DMGetCoordinateDM(dm, &cdm);
159: DMGetDS(cdm, &cds);
160: PetscDSGetDiscretization(cds, 0, (PetscObject *) &fe);
161: PetscObjectGetClassId((PetscObject) fe, &id);
162: if (id != PETSCFE_CLASSID) {
163: PetscBool simplex;
164: PetscInt dim, dE, qorder;
166: DMGetDimension(dm, &dim);
167: DMGetCoordinateDim(dm, &dE);
168: DMPlexIsSimplex(dm, &simplex);
169: qorder = degree;
170: PetscObjectOptionsBegin((PetscObject) cdm);
171: PetscOptionsBoundedInt("-coord_dm_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "DMPlexCreateCoordinateSpace", qorder, &qorder, NULL, 0);
172: PetscOptionsEnd();
173: PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, degree, qorder, &fe);
174: DMSetField(cdm, 0, NULL, (PetscObject) fe);
175: DMCreateDS(cdm);
176: DMProjectCoordinates(dm, fe);
177: PetscFEDestroy(&fe);
178: }
179: mesh->coordFunc = coordFunc;
180: return(0);
181: }
183: /*@
184: DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement.
186: Collective
188: Input Parameters:
189: + comm - The communicator for the DM object
190: . dim - The spatial dimension
191: . simplex - Flag for simplicial cells, otherwise they are tensor product cells
192: . interpolate - Flag to create intermediate mesh pieces (edges, faces)
193: - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell
195: Output Parameter:
196: . dm - The DM object
198: Level: beginner
200: .seealso: DMSetType(), DMCreate()
201: @*/
202: PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscReal refinementLimit, DM *newdm)
203: {
204: DM dm;
205: PetscMPIInt rank;
209: DMCreate(comm, &dm);
210: DMSetType(dm, DMPLEX);
211: DMSetDimension(dm, dim);
212: MPI_Comm_rank(comm, &rank);
213: switch (dim) {
214: case 2:
215: if (simplex) {PetscObjectSetName((PetscObject) dm, "triangular");}
216: else {PetscObjectSetName((PetscObject) dm, "quadrilateral");}
217: break;
218: case 3:
219: if (simplex) {PetscObjectSetName((PetscObject) dm, "tetrahedral");}
220: else {PetscObjectSetName((PetscObject) dm, "hexahedral");}
221: break;
222: default:
223: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", dim);
224: }
225: if (rank) {
226: PetscInt numPoints[2] = {0, 0};
227: DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);
228: } else {
229: switch (dim) {
230: case 2:
231: if (simplex) {
232: PetscInt numPoints[2] = {4, 2};
233: PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0};
234: PetscInt cones[6] = {2, 3, 4, 5, 4, 3};
235: PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
236: PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
238: DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
239: } else {
240: PetscInt numPoints[2] = {6, 2};
241: PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0};
242: PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4};
243: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
244: 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};
246: DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
247: }
248: break;
249: case 3:
250: if (simplex) {
251: PetscInt numPoints[2] = {5, 2};
252: PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0};
253: PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6};
254: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
255: 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};
257: DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
258: } else {
259: PetscInt numPoints[2] = {12, 2};
260: PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
261: PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8};
262: PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
263: 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,
264: -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5,
265: 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5};
267: DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
268: }
269: break;
270: default:
271: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", dim);
272: }
273: }
274: *newdm = dm;
275: if (refinementLimit > 0.0) {
276: DM rdm;
277: const char *name;
279: DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);
280: DMPlexSetRefinementLimit(*newdm, refinementLimit);
281: DMRefine(*newdm, comm, &rdm);
282: PetscObjectGetName((PetscObject) *newdm, &name);
283: PetscObjectSetName((PetscObject) rdm, name);
284: DMDestroy(newdm);
285: *newdm = rdm;
286: }
287: if (interpolate) {
288: DM idm;
290: DMPlexInterpolate(*newdm, &idm);
291: DMDestroy(newdm);
292: *newdm = idm;
293: }
294: return(0);
295: }
297: static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
298: {
299: const PetscInt numVertices = 2;
300: PetscInt markerRight = 1;
301: PetscInt markerLeft = 1;
302: PetscBool markerSeparate = PETSC_FALSE;
303: Vec coordinates;
304: PetscSection coordSection;
305: PetscScalar *coords;
306: PetscInt coordSize;
307: PetscMPIInt rank;
308: PetscInt cdim = 1, v;
312: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
313: if (markerSeparate) {
314: markerRight = 2;
315: markerLeft = 1;
316: }
317: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
318: if (!rank) {
319: DMPlexSetChart(dm, 0, numVertices);
320: DMSetUp(dm); /* Allocate space for cones */
321: DMSetLabelValue(dm, "marker", 0, markerLeft);
322: DMSetLabelValue(dm, "marker", 1, markerRight);
323: }
324: DMPlexSymmetrize(dm);
325: DMPlexStratify(dm);
326: /* Build coordinates */
327: DMSetCoordinateDim(dm, cdim);
328: DMGetCoordinateSection(dm, &coordSection);
329: PetscSectionSetNumFields(coordSection, 1);
330: PetscSectionSetChart(coordSection, 0, numVertices);
331: PetscSectionSetFieldComponents(coordSection, 0, cdim);
332: for (v = 0; v < numVertices; ++v) {
333: PetscSectionSetDof(coordSection, v, cdim);
334: PetscSectionSetFieldDof(coordSection, v, 0, cdim);
335: }
336: PetscSectionSetUp(coordSection);
337: PetscSectionGetStorageSize(coordSection, &coordSize);
338: VecCreate(PETSC_COMM_SELF, &coordinates);
339: PetscObjectSetName((PetscObject) coordinates, "coordinates");
340: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
341: VecSetBlockSize(coordinates, cdim);
342: VecSetType(coordinates,VECSTANDARD);
343: VecGetArray(coordinates, &coords);
344: coords[0] = lower[0];
345: coords[1] = upper[0];
346: VecRestoreArray(coordinates, &coords);
347: DMSetCoordinatesLocal(dm, coordinates);
348: VecDestroy(&coordinates);
349: return(0);
350: }
352: static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
353: {
354: const PetscInt numVertices = (edges[0]+1)*(edges[1]+1);
355: const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
356: PetscInt markerTop = 1;
357: PetscInt markerBottom = 1;
358: PetscInt markerRight = 1;
359: PetscInt markerLeft = 1;
360: PetscBool markerSeparate = PETSC_FALSE;
361: Vec coordinates;
362: PetscSection coordSection;
363: PetscScalar *coords;
364: PetscInt coordSize;
365: PetscMPIInt rank;
366: PetscInt v, vx, vy;
370: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
371: if (markerSeparate) {
372: markerTop = 3;
373: markerBottom = 1;
374: markerRight = 2;
375: markerLeft = 4;
376: }
377: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
378: if (rank == 0) {
379: PetscInt e, ex, ey;
381: DMPlexSetChart(dm, 0, numEdges+numVertices);
382: for (e = 0; e < numEdges; ++e) {
383: DMPlexSetConeSize(dm, e, 2);
384: }
385: DMSetUp(dm); /* Allocate space for cones */
386: for (vx = 0; vx <= edges[0]; vx++) {
387: for (ey = 0; ey < edges[1]; ey++) {
388: PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1);
389: PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
390: PetscInt cone[2];
392: cone[0] = vertex; cone[1] = vertex+edges[0]+1;
393: DMPlexSetCone(dm, edge, cone);
394: if (vx == edges[0]) {
395: DMSetLabelValue(dm, "marker", edge, markerRight);
396: DMSetLabelValue(dm, "marker", cone[0], markerRight);
397: if (ey == edges[1]-1) {
398: DMSetLabelValue(dm, "marker", cone[1], markerRight);
399: DMSetLabelValue(dm, "Face Sets", cone[1], markerRight);
400: }
401: } else if (vx == 0) {
402: DMSetLabelValue(dm, "marker", edge, markerLeft);
403: DMSetLabelValue(dm, "marker", cone[0], markerLeft);
404: if (ey == edges[1]-1) {
405: DMSetLabelValue(dm, "marker", cone[1], markerLeft);
406: DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft);
407: }
408: }
409: }
410: }
411: for (vy = 0; vy <= edges[1]; vy++) {
412: for (ex = 0; ex < edges[0]; ex++) {
413: PetscInt edge = vy*edges[0] + ex;
414: PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
415: PetscInt cone[2];
417: cone[0] = vertex; cone[1] = vertex+1;
418: DMPlexSetCone(dm, edge, cone);
419: if (vy == edges[1]) {
420: DMSetLabelValue(dm, "marker", edge, markerTop);
421: DMSetLabelValue(dm, "marker", cone[0], markerTop);
422: if (ex == edges[0]-1) {
423: DMSetLabelValue(dm, "marker", cone[1], markerTop);
424: DMSetLabelValue(dm, "Face Sets", cone[1], markerTop);
425: }
426: } else if (vy == 0) {
427: DMSetLabelValue(dm, "marker", edge, markerBottom);
428: DMSetLabelValue(dm, "marker", cone[0], markerBottom);
429: if (ex == edges[0]-1) {
430: DMSetLabelValue(dm, "marker", cone[1], markerBottom);
431: DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom);
432: }
433: }
434: }
435: }
436: }
437: DMPlexSymmetrize(dm);
438: DMPlexStratify(dm);
439: /* Build coordinates */
440: DMSetCoordinateDim(dm, 2);
441: DMGetCoordinateSection(dm, &coordSection);
442: PetscSectionSetNumFields(coordSection, 1);
443: PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);
444: PetscSectionSetFieldComponents(coordSection, 0, 2);
445: for (v = numEdges; v < numEdges+numVertices; ++v) {
446: PetscSectionSetDof(coordSection, v, 2);
447: PetscSectionSetFieldDof(coordSection, v, 0, 2);
448: }
449: PetscSectionSetUp(coordSection);
450: PetscSectionGetStorageSize(coordSection, &coordSize);
451: VecCreate(PETSC_COMM_SELF, &coordinates);
452: PetscObjectSetName((PetscObject) coordinates, "coordinates");
453: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
454: VecSetBlockSize(coordinates, 2);
455: VecSetType(coordinates,VECSTANDARD);
456: VecGetArray(coordinates, &coords);
457: for (vy = 0; vy <= edges[1]; ++vy) {
458: for (vx = 0; vx <= edges[0]; ++vx) {
459: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
460: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
461: }
462: }
463: VecRestoreArray(coordinates, &coords);
464: DMSetCoordinatesLocal(dm, coordinates);
465: VecDestroy(&coordinates);
466: return(0);
467: }
469: static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
470: {
471: PetscInt vertices[3], numVertices;
472: PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2];
473: Vec coordinates;
474: PetscSection coordSection;
475: PetscScalar *coords;
476: PetscInt coordSize;
477: PetscMPIInt rank;
478: PetscInt v, vx, vy, vz;
479: PetscInt voffset, iface=0, cone[4];
483: 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");
484: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
485: vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1;
486: numVertices = vertices[0]*vertices[1]*vertices[2];
487: if (rank == 0) {
488: PetscInt f;
490: DMPlexSetChart(dm, 0, numFaces+numVertices);
491: for (f = 0; f < numFaces; ++f) {
492: DMPlexSetConeSize(dm, f, 4);
493: }
494: DMSetUp(dm); /* Allocate space for cones */
496: /* Side 0 (Top) */
497: for (vy = 0; vy < faces[1]; vy++) {
498: for (vx = 0; vx < faces[0]; vx++) {
499: voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx;
500: cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0];
501: DMPlexSetCone(dm, iface, cone);
502: DMSetLabelValue(dm, "marker", iface, 1);
503: DMSetLabelValue(dm, "marker", voffset+0, 1);
504: DMSetLabelValue(dm, "marker", voffset+1, 1);
505: DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);
506: DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);
507: iface++;
508: }
509: }
511: /* Side 1 (Bottom) */
512: for (vy = 0; vy < faces[1]; vy++) {
513: for (vx = 0; vx < faces[0]; vx++) {
514: voffset = numFaces + vy*(faces[0]+1) + vx;
515: cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1;
516: DMPlexSetCone(dm, iface, cone);
517: DMSetLabelValue(dm, "marker", iface, 1);
518: DMSetLabelValue(dm, "marker", voffset+0, 1);
519: DMSetLabelValue(dm, "marker", voffset+1, 1);
520: DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);
521: DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);
522: iface++;
523: }
524: }
526: /* Side 2 (Front) */
527: for (vz = 0; vz < faces[2]; vz++) {
528: for (vx = 0; vx < faces[0]; vx++) {
529: voffset = numFaces + vz*vertices[0]*vertices[1] + vx;
530: cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1];
531: DMPlexSetCone(dm, iface, cone);
532: DMSetLabelValue(dm, "marker", iface, 1);
533: DMSetLabelValue(dm, "marker", voffset+0, 1);
534: DMSetLabelValue(dm, "marker", voffset+1, 1);
535: DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);
536: DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);
537: iface++;
538: }
539: }
541: /* Side 3 (Back) */
542: for (vz = 0; vz < faces[2]; vz++) {
543: for (vx = 0; vx < faces[0]; vx++) {
544: voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx;
545: cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1;
546: cone[2] = voffset+1; cone[3] = voffset;
547: DMPlexSetCone(dm, iface, cone);
548: DMSetLabelValue(dm, "marker", iface, 1);
549: DMSetLabelValue(dm, "marker", voffset+0, 1);
550: DMSetLabelValue(dm, "marker", voffset+1, 1);
551: DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);
552: DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);
553: iface++;
554: }
555: }
557: /* Side 4 (Left) */
558: for (vz = 0; vz < faces[2]; vz++) {
559: for (vy = 0; vy < faces[1]; vy++) {
560: voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0];
561: cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1];
562: cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0];
563: DMPlexSetCone(dm, iface, cone);
564: DMSetLabelValue(dm, "marker", iface, 1);
565: DMSetLabelValue(dm, "marker", voffset+0, 1);
566: DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);
567: DMSetLabelValue(dm, "marker", voffset+vertices[1]+0, 1);
568: DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);
569: iface++;
570: }
571: }
573: /* Side 5 (Right) */
574: for (vz = 0; vz < faces[2]; vz++) {
575: for (vy = 0; vy < faces[1]; vy++) {
576: voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + faces[0];
577: cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset;
578: cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0];
579: DMPlexSetCone(dm, iface, cone);
580: DMSetLabelValue(dm, "marker", iface, 1);
581: DMSetLabelValue(dm, "marker", voffset+0, 1);
582: DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);
583: DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);
584: DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);
585: iface++;
586: }
587: }
588: }
589: DMPlexSymmetrize(dm);
590: DMPlexStratify(dm);
591: /* Build coordinates */
592: DMSetCoordinateDim(dm, 3);
593: DMGetCoordinateSection(dm, &coordSection);
594: PetscSectionSetNumFields(coordSection, 1);
595: PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);
596: PetscSectionSetFieldComponents(coordSection, 0, 3);
597: for (v = numFaces; v < numFaces+numVertices; ++v) {
598: PetscSectionSetDof(coordSection, v, 3);
599: PetscSectionSetFieldDof(coordSection, v, 0, 3);
600: }
601: PetscSectionSetUp(coordSection);
602: PetscSectionGetStorageSize(coordSection, &coordSize);
603: VecCreate(PETSC_COMM_SELF, &coordinates);
604: PetscObjectSetName((PetscObject) coordinates, "coordinates");
605: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
606: VecSetBlockSize(coordinates, 3);
607: VecSetType(coordinates,VECSTANDARD);
608: VecGetArray(coordinates, &coords);
609: for (vz = 0; vz <= faces[2]; ++vz) {
610: for (vy = 0; vy <= faces[1]; ++vy) {
611: for (vx = 0; vx <= faces[0]; ++vx) {
612: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
613: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
614: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
615: }
616: }
617: }
618: VecRestoreArray(coordinates, &coords);
619: DMSetCoordinatesLocal(dm, coordinates);
620: VecDestroy(&coordinates);
621: return(0);
622: }
624: static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate)
625: {
630: DMSetDimension(dm, dim-1);
631: DMSetCoordinateDim(dm, dim);
632: switch (dim) {
633: case 1: DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(dm, lower, upper, faces);break;
634: case 2: DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(dm, lower, upper, faces);break;
635: case 3: DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(dm, lower, upper, faces);break;
636: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension not supported: %D", dim);
637: }
638: if (interpolate) {DMPlexInterpolateInPlace_Internal(dm);}
639: return(0);
640: }
642: /*@C
643: DMPlexCreateBoxSurfaceMesh - Creates a mesh on the surface of the tensor product of unit intervals (box) using tensor cells (hexahedra).
645: Collective
647: Input Parameters:
648: + comm - The communicator for the DM object
649: . dim - The spatial dimension of the box, so the resulting mesh is has dimension dim-1
650: . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D
651: . lower - The lower left corner, or NULL for (0, 0, 0)
652: . upper - The upper right corner, or NULL for (1, 1, 1)
653: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
655: Output Parameter:
656: . dm - The DM object
658: Level: beginner
660: .seealso: DMSetFromOptions(), DMPlexCreateBoxMesh(), DMPlexCreateFromFile(), DMSetType(), DMCreate()
661: @*/
662: PetscErrorCode DMPlexCreateBoxSurfaceMesh(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate, DM *dm)
663: {
664: PetscInt fac[3] = {1, 1, 1};
665: PetscReal low[3] = {0, 0, 0};
666: PetscReal upp[3] = {1, 1, 1};
670: DMCreate(comm,dm);
671: DMSetType(*dm,DMPLEX);
672: DMPlexCreateBoxSurfaceMesh_Internal(*dm, dim, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, interpolate);
673: return(0);
674: }
676: static PetscErrorCode DMPlexCreateLineMesh_Internal(DM dm,PetscInt segments,PetscReal lower,PetscReal upper,DMBoundaryType bd)
677: {
678: PetscInt i,fStart,fEnd,numCells = 0,numVerts = 0;
679: PetscInt numPoints[2],*coneSize,*cones,*coneOrientations;
680: PetscScalar *vertexCoords;
681: PetscReal L,maxCell;
682: PetscBool markerSeparate = PETSC_FALSE;
683: PetscInt markerLeft = 1, faceMarkerLeft = 1;
684: PetscInt markerRight = 1, faceMarkerRight = 2;
685: PetscBool wrap = (bd == DM_BOUNDARY_PERIODIC || bd == DM_BOUNDARY_TWIST) ? PETSC_TRUE : PETSC_FALSE;
686: PetscMPIInt rank;
692: DMSetDimension(dm,1);
693: DMCreateLabel(dm,"marker");
694: DMCreateLabel(dm,"Face Sets");
696: MPI_Comm_rank(PetscObjectComm((PetscObject) dm),&rank);
697: if (rank == 0) numCells = segments;
698: if (rank == 0) numVerts = segments + (wrap ? 0 : 1);
700: numPoints[0] = numVerts ; numPoints[1] = numCells;
701: PetscMalloc4(numCells+numVerts,&coneSize,numCells*2,&cones,numCells+numVerts,&coneOrientations,numVerts,&vertexCoords);
702: PetscArrayzero(coneOrientations,numCells+numVerts);
703: for (i = 0; i < numCells; ++i) { coneSize[i] = 2; }
704: for (i = 0; i < numVerts; ++i) { coneSize[numCells+i] = 0; }
705: for (i = 0; i < numCells; ++i) { cones[2*i] = numCells + i%numVerts; cones[2*i+1] = numCells + (i+1)%numVerts; }
706: for (i = 0; i < numVerts; ++i) { vertexCoords[i] = lower + (upper-lower)*((PetscReal)i/(PetscReal)numCells); }
707: DMPlexCreateFromDAG(dm,1,numPoints,coneSize,cones,coneOrientations,vertexCoords);
708: PetscFree4(coneSize,cones,coneOrientations,vertexCoords);
710: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dm_plex_separate_marker",&markerSeparate,NULL);
711: if (markerSeparate) { markerLeft = faceMarkerLeft; markerRight = faceMarkerRight;}
712: if (!wrap && rank == 0) {
713: DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);
714: DMSetLabelValue(dm,"marker",fStart,markerLeft);
715: DMSetLabelValue(dm,"marker",fEnd-1,markerRight);
716: DMSetLabelValue(dm,"Face Sets",fStart,faceMarkerLeft);
717: DMSetLabelValue(dm,"Face Sets",fEnd-1,faceMarkerRight);
718: }
719: if (wrap) {
720: L = upper - lower;
721: maxCell = (PetscReal)1.1*(L/(PetscReal)PetscMax(1,segments));
722: DMSetPeriodicity(dm,PETSC_TRUE,&maxCell,&L,&bd);
723: }
724: DMPlexSetRefinementUniform(dm, PETSC_TRUE);
725: return(0);
726: }
728: static PetscErrorCode DMPlexCreateBoxMesh_Simplex_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate)
729: {
730: DM boundary, vol;
731: PetscInt i;
736: for (i = 0; i < dim; ++i) if (periodicity[i] != DM_BOUNDARY_NONE) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity is not supported for simplex meshes");
737: DMCreate(PetscObjectComm((PetscObject) dm), &boundary);
738: DMSetType(boundary, DMPLEX);
739: DMPlexCreateBoxSurfaceMesh_Internal(boundary, dim, faces, lower, upper, PETSC_FALSE);
740: DMPlexGenerate(boundary, NULL, interpolate, &vol);
741: DMPlexReplace_Static(dm, &vol);
742: DMDestroy(&boundary);
743: return(0);
744: }
746: static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ)
747: {
748: DMLabel cutLabel = NULL;
749: PetscInt markerTop = 1, faceMarkerTop = 1;
750: PetscInt markerBottom = 1, faceMarkerBottom = 1;
751: PetscInt markerFront = 1, faceMarkerFront = 1;
752: PetscInt markerBack = 1, faceMarkerBack = 1;
753: PetscInt markerRight = 1, faceMarkerRight = 1;
754: PetscInt markerLeft = 1, faceMarkerLeft = 1;
755: PetscInt dim;
756: PetscBool markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE;
757: PetscMPIInt rank;
761: DMGetDimension(dm,&dim);
762: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
763: DMCreateLabel(dm,"marker");
764: DMCreateLabel(dm,"Face Sets");
765: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);
766: if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ||
767: bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ||
768: bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) {
770: if (cutMarker) {DMCreateLabel(dm, "periodic_cut"); DMGetLabel(dm, "periodic_cut", &cutLabel);}
771: }
772: switch (dim) {
773: case 2:
774: faceMarkerTop = 3;
775: faceMarkerBottom = 1;
776: faceMarkerRight = 2;
777: faceMarkerLeft = 4;
778: break;
779: case 3:
780: faceMarkerBottom = 1;
781: faceMarkerTop = 2;
782: faceMarkerFront = 3;
783: faceMarkerBack = 4;
784: faceMarkerRight = 5;
785: faceMarkerLeft = 6;
786: break;
787: default:
788: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D not supported",dim);
789: }
790: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);
791: if (markerSeparate) {
792: markerBottom = faceMarkerBottom;
793: markerTop = faceMarkerTop;
794: markerFront = faceMarkerFront;
795: markerBack = faceMarkerBack;
796: markerRight = faceMarkerRight;
797: markerLeft = faceMarkerLeft;
798: }
799: {
800: const PetscInt numXEdges = rank == 0 ? edges[0] : 0;
801: const PetscInt numYEdges = rank == 0 ? edges[1] : 0;
802: const PetscInt numZEdges = rank == 0 ? edges[2] : 0;
803: const PetscInt numXVertices = rank == 0 ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0;
804: const PetscInt numYVertices = rank == 0 ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0;
805: const PetscInt numZVertices = rank == 0 ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0;
806: const PetscInt numCells = numXEdges*numYEdges*numZEdges;
807: const PetscInt numXFaces = numYEdges*numZEdges;
808: const PetscInt numYFaces = numXEdges*numZEdges;
809: const PetscInt numZFaces = numXEdges*numYEdges;
810: const PetscInt numTotXFaces = numXVertices*numXFaces;
811: const PetscInt numTotYFaces = numYVertices*numYFaces;
812: const PetscInt numTotZFaces = numZVertices*numZFaces;
813: const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces;
814: const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices;
815: const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices;
816: const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices;
817: const PetscInt numVertices = numXVertices*numYVertices*numZVertices;
818: const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges;
819: const PetscInt firstVertex = (dim == 2) ? numFaces : numCells;
820: const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices;
821: const PetscInt firstYFace = firstXFace + numTotXFaces;
822: const PetscInt firstZFace = firstYFace + numTotYFaces;
823: const PetscInt firstXEdge = numCells + numFaces + numVertices;
824: const PetscInt firstYEdge = firstXEdge + numTotXEdges;
825: const PetscInt firstZEdge = firstYEdge + numTotYEdges;
826: Vec coordinates;
827: PetscSection coordSection;
828: PetscScalar *coords;
829: PetscInt coordSize;
830: PetscInt v, vx, vy, vz;
831: PetscInt c, f, fx, fy, fz, e, ex, ey, ez;
833: DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);
834: for (c = 0; c < numCells; c++) {
835: DMPlexSetConeSize(dm, c, 6);
836: }
837: for (f = firstXFace; f < firstXFace+numFaces; ++f) {
838: DMPlexSetConeSize(dm, f, 4);
839: }
840: for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
841: DMPlexSetConeSize(dm, e, 2);
842: }
843: DMSetUp(dm); /* Allocate space for cones */
844: /* Build cells */
845: for (fz = 0; fz < numZEdges; ++fz) {
846: for (fy = 0; fy < numYEdges; ++fy) {
847: for (fx = 0; fx < numXEdges; ++fx) {
848: PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx;
849: PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz;
850: PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices);
851: PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy;
852: PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices);
853: PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx;
854: PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices);
855: /* B, T, F, K, R, L */
856: PetscInt ornt[6] = {-2, 0, 0, -3, 0, -2}; /* ??? */
857: PetscInt cone[6];
859: /* no boundary twisting in 3D */
860: cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL;
861: DMPlexSetCone(dm, cell, cone);
862: DMPlexSetConeOrientation(dm, cell, ornt);
863: if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, cell, 2);}
864: if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, cell, 2);}
865: if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, cell, 2);}
866: }
867: }
868: }
869: /* Build x faces */
870: for (fz = 0; fz < numZEdges; ++fz) {
871: for (fy = 0; fy < numYEdges; ++fy) {
872: for (fx = 0; fx < numXVertices; ++fx) {
873: PetscInt face = firstXFace + (fz*numYEdges+fy) *numXVertices+fx;
874: PetscInt edgeL = firstZEdge + (fy *numXVertices+fx)*numZEdges + fz;
875: PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz;
876: PetscInt edgeB = firstYEdge + (fz *numXVertices+fx)*numYEdges + fy;
877: PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy;
878: PetscInt ornt[4] = {0, 0, -1, -1};
879: PetscInt cone[4];
881: if (dim == 3) {
882: /* markers */
883: if (bdX != DM_BOUNDARY_PERIODIC) {
884: if (fx == numXVertices-1) {
885: DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);
886: DMSetLabelValue(dm, "marker", face, markerRight);
887: }
888: else if (fx == 0) {
889: DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);
890: DMSetLabelValue(dm, "marker", face, markerLeft);
891: }
892: }
893: }
894: cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
895: DMPlexSetCone(dm, face, cone);
896: DMPlexSetConeOrientation(dm, face, ornt);
897: }
898: }
899: }
900: /* Build y faces */
901: for (fz = 0; fz < numZEdges; ++fz) {
902: for (fx = 0; fx < numXEdges; ++fx) {
903: for (fy = 0; fy < numYVertices; ++fy) {
904: PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy;
905: PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx)*numZEdges + fz;
906: PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz;
907: PetscInt edgeB = firstXEdge + (fz *numYVertices+fy)*numXEdges + fx;
908: PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx;
909: PetscInt ornt[4] = {0, 0, -1, -1};
910: PetscInt cone[4];
912: if (dim == 3) {
913: /* markers */
914: if (bdY != DM_BOUNDARY_PERIODIC) {
915: if (fy == numYVertices-1) {
916: DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);
917: DMSetLabelValue(dm, "marker", face, markerBack);
918: }
919: else if (fy == 0) {
920: DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);
921: DMSetLabelValue(dm, "marker", face, markerFront);
922: }
923: }
924: }
925: cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
926: DMPlexSetCone(dm, face, cone);
927: DMPlexSetConeOrientation(dm, face, ornt);
928: }
929: }
930: }
931: /* Build z faces */
932: for (fy = 0; fy < numYEdges; ++fy) {
933: for (fx = 0; fx < numXEdges; ++fx) {
934: for (fz = 0; fz < numZVertices; fz++) {
935: PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz;
936: PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx)*numYEdges + fy;
937: PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy;
938: PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy)*numXEdges + fx;
939: PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx;
940: PetscInt ornt[4] = {0, 0, -1, -1};
941: PetscInt cone[4];
943: if (dim == 2) {
944: if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -1;}
945: if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;}
946: if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, face, 2);}
947: if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {DMLabelSetValue(cutLabel, face, 2);}
948: } else {
949: /* markers */
950: if (bdZ != DM_BOUNDARY_PERIODIC) {
951: if (fz == numZVertices-1) {
952: DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);
953: DMSetLabelValue(dm, "marker", face, markerTop);
954: }
955: else if (fz == 0) {
956: DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);
957: DMSetLabelValue(dm, "marker", face, markerBottom);
958: }
959: }
960: }
961: cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
962: DMPlexSetCone(dm, face, cone);
963: DMPlexSetConeOrientation(dm, face, ornt);
964: }
965: }
966: }
967: /* Build Z edges*/
968: for (vy = 0; vy < numYVertices; vy++) {
969: for (vx = 0; vx < numXVertices; vx++) {
970: for (ez = 0; ez < numZEdges; ez++) {
971: const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez;
972: const PetscInt vertexB = firstVertex + (ez *numYVertices+vy)*numXVertices + vx;
973: const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx;
974: PetscInt cone[2];
976: if (dim == 3) {
977: if (bdX != DM_BOUNDARY_PERIODIC) {
978: if (vx == numXVertices-1) {
979: DMSetLabelValue(dm, "marker", edge, markerRight);
980: }
981: else if (vx == 0) {
982: DMSetLabelValue(dm, "marker", edge, markerLeft);
983: }
984: }
985: if (bdY != DM_BOUNDARY_PERIODIC) {
986: if (vy == numYVertices-1) {
987: DMSetLabelValue(dm, "marker", edge, markerBack);
988: }
989: else if (vy == 0) {
990: DMSetLabelValue(dm, "marker", edge, markerFront);
991: }
992: }
993: }
994: cone[0] = vertexB; cone[1] = vertexT;
995: DMPlexSetCone(dm, edge, cone);
996: }
997: }
998: }
999: /* Build Y edges*/
1000: for (vz = 0; vz < numZVertices; vz++) {
1001: for (vx = 0; vx < numXVertices; vx++) {
1002: for (ey = 0; ey < numYEdges; ey++) {
1003: const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx;
1004: const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey;
1005: const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx;
1006: const PetscInt vertexK = firstVertex + nextv;
1007: PetscInt cone[2];
1009: cone[0] = vertexF; cone[1] = vertexK;
1010: DMPlexSetCone(dm, edge, cone);
1011: if (dim == 2) {
1012: if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) {
1013: if (vx == numXVertices-1) {
1014: DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);
1015: DMSetLabelValue(dm, "marker", edge, markerRight);
1016: DMSetLabelValue(dm, "marker", cone[0], markerRight);
1017: if (ey == numYEdges-1) {
1018: DMSetLabelValue(dm, "marker", cone[1], markerRight);
1019: }
1020: } else if (vx == 0) {
1021: DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);
1022: DMSetLabelValue(dm, "marker", edge, markerLeft);
1023: DMSetLabelValue(dm, "marker", cone[0], markerLeft);
1024: if (ey == numYEdges-1) {
1025: DMSetLabelValue(dm, "marker", cone[1], markerLeft);
1026: }
1027: }
1028: } else {
1029: if (vx == 0 && cutLabel) {
1030: DMLabelSetValue(cutLabel, edge, 1);
1031: DMLabelSetValue(cutLabel, cone[0], 1);
1032: if (ey == numYEdges-1) {
1033: DMLabelSetValue(cutLabel, cone[1], 1);
1034: }
1035: }
1036: }
1037: } else {
1038: if (bdX != DM_BOUNDARY_PERIODIC) {
1039: if (vx == numXVertices-1) {
1040: DMSetLabelValue(dm, "marker", edge, markerRight);
1041: } else if (vx == 0) {
1042: DMSetLabelValue(dm, "marker", edge, markerLeft);
1043: }
1044: }
1045: if (bdZ != DM_BOUNDARY_PERIODIC) {
1046: if (vz == numZVertices-1) {
1047: DMSetLabelValue(dm, "marker", edge, markerTop);
1048: } else if (vz == 0) {
1049: DMSetLabelValue(dm, "marker", edge, markerBottom);
1050: }
1051: }
1052: }
1053: }
1054: }
1055: }
1056: /* Build X edges*/
1057: for (vz = 0; vz < numZVertices; vz++) {
1058: for (vy = 0; vy < numYVertices; vy++) {
1059: for (ex = 0; ex < numXEdges; ex++) {
1060: const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices;
1061: const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex;
1062: const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex;
1063: const PetscInt vertexR = firstVertex + nextv;
1064: PetscInt cone[2];
1066: cone[0] = vertexL; cone[1] = vertexR;
1067: DMPlexSetCone(dm, edge, cone);
1068: if (dim == 2) {
1069: if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) {
1070: if (vy == numYVertices-1) {
1071: DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);
1072: DMSetLabelValue(dm, "marker", edge, markerTop);
1073: DMSetLabelValue(dm, "marker", cone[0], markerTop);
1074: if (ex == numXEdges-1) {
1075: DMSetLabelValue(dm, "marker", cone[1], markerTop);
1076: }
1077: } else if (vy == 0) {
1078: DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);
1079: DMSetLabelValue(dm, "marker", edge, markerBottom);
1080: DMSetLabelValue(dm, "marker", cone[0], markerBottom);
1081: if (ex == numXEdges-1) {
1082: DMSetLabelValue(dm, "marker", cone[1], markerBottom);
1083: }
1084: }
1085: } else {
1086: if (vy == 0 && cutLabel) {
1087: DMLabelSetValue(cutLabel, edge, 1);
1088: DMLabelSetValue(cutLabel, cone[0], 1);
1089: if (ex == numXEdges-1) {
1090: DMLabelSetValue(cutLabel, cone[1], 1);
1091: }
1092: }
1093: }
1094: } else {
1095: if (bdY != DM_BOUNDARY_PERIODIC) {
1096: if (vy == numYVertices-1) {
1097: DMSetLabelValue(dm, "marker", edge, markerBack);
1098: }
1099: else if (vy == 0) {
1100: DMSetLabelValue(dm, "marker", edge, markerFront);
1101: }
1102: }
1103: if (bdZ != DM_BOUNDARY_PERIODIC) {
1104: if (vz == numZVertices-1) {
1105: DMSetLabelValue(dm, "marker", edge, markerTop);
1106: }
1107: else if (vz == 0) {
1108: DMSetLabelValue(dm, "marker", edge, markerBottom);
1109: }
1110: }
1111: }
1112: }
1113: }
1114: }
1115: DMPlexSymmetrize(dm);
1116: DMPlexStratify(dm);
1117: /* Build coordinates */
1118: DMGetCoordinateSection(dm, &coordSection);
1119: PetscSectionSetNumFields(coordSection, 1);
1120: PetscSectionSetFieldComponents(coordSection, 0, dim);
1121: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);
1122: for (v = firstVertex; v < firstVertex+numVertices; ++v) {
1123: PetscSectionSetDof(coordSection, v, dim);
1124: PetscSectionSetFieldDof(coordSection, v, 0, dim);
1125: }
1126: PetscSectionSetUp(coordSection);
1127: PetscSectionGetStorageSize(coordSection, &coordSize);
1128: VecCreate(PETSC_COMM_SELF, &coordinates);
1129: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1130: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1131: VecSetBlockSize(coordinates, dim);
1132: VecSetType(coordinates,VECSTANDARD);
1133: VecGetArray(coordinates, &coords);
1134: for (vz = 0; vz < numZVertices; ++vz) {
1135: for (vy = 0; vy < numYVertices; ++vy) {
1136: for (vx = 0; vx < numXVertices; ++vx) {
1137: coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
1138: coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
1139: if (dim == 3) {
1140: coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz;
1141: }
1142: }
1143: }
1144: }
1145: VecRestoreArray(coordinates, &coords);
1146: DMSetCoordinatesLocal(dm, coordinates);
1147: VecDestroy(&coordinates);
1148: }
1149: return(0);
1150: }
1152: static PetscErrorCode DMPlexCreateBoxMesh_Tensor_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[])
1153: {
1154: DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
1155: PetscInt fac[3] = {0, 0, 0}, d;
1161: DMSetDimension(dm, dim);
1162: for (d = 0; d < dim; ++d) {fac[d] = faces[d]; bdt[d] = periodicity[d];}
1163: DMPlexCreateCubeMesh_Internal(dm, lower, upper, fac, bdt[0], bdt[1], bdt[2]);
1164: if (periodicity[0] == DM_BOUNDARY_PERIODIC || periodicity[0] == DM_BOUNDARY_TWIST ||
1165: periodicity[1] == DM_BOUNDARY_PERIODIC || periodicity[1] == DM_BOUNDARY_TWIST ||
1166: (dim > 2 && (periodicity[2] == DM_BOUNDARY_PERIODIC || periodicity[2] == DM_BOUNDARY_TWIST))) {
1167: PetscReal L[3];
1168: PetscReal maxCell[3];
1170: for (d = 0; d < dim; ++d) {
1171: L[d] = upper[d] - lower[d];
1172: maxCell[d] = 1.1 * (L[d] / PetscMax(1, faces[d]));
1173: }
1174: DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, periodicity);
1175: }
1176: DMPlexSetRefinementUniform(dm, PETSC_TRUE);
1177: return(0);
1178: }
1180: static PetscErrorCode DMPlexCreateBoxMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate)
1181: {
1185: if (dim == 1) {DMPlexCreateLineMesh_Internal(dm, faces[0], lower[0], upper[0], periodicity[0]);}
1186: else if (simplex) {DMPlexCreateBoxMesh_Simplex_Internal(dm, dim, faces, lower, upper, periodicity, interpolate);}
1187: else {DMPlexCreateBoxMesh_Tensor_Internal(dm, dim, faces, lower, upper, periodicity);}
1188: if (!interpolate && dim > 1 && !simplex) {
1189: DM udm;
1191: DMPlexUninterpolate(dm, &udm);
1192: DMPlexCopyCoordinates(dm, udm);
1193: DMPlexReplace_Static(dm, &udm);
1194: }
1195: return(0);
1196: }
1198: /*@C
1199: DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices or tensor cells (hexahedra).
1201: Collective
1203: Input Parameters:
1204: + comm - The communicator for the DM object
1205: . dim - The spatial dimension
1206: . simplex - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells
1207: . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D
1208: . lower - The lower left corner, or NULL for (0, 0, 0)
1209: . upper - The upper right corner, or NULL for (1, 1, 1)
1210: . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE
1211: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
1213: Output Parameter:
1214: . dm - The DM object
1216: Note: If you want to customize this mesh using options, you just need to
1217: $ DMCreate(comm, &dm);
1218: $ DMSetType(dm, DMPLEX);
1219: $ DMSetFromOptions(dm);
1220: and use the options on the DMSetFromOptions() page.
1222: Here is the numbering returned for 2 faces in each direction for tensor cells:
1223: $ 10---17---11---18----12
1224: $ | | |
1225: $ | | |
1226: $ 20 2 22 3 24
1227: $ | | |
1228: $ | | |
1229: $ 7---15----8---16----9
1230: $ | | |
1231: $ | | |
1232: $ 19 0 21 1 23
1233: $ | | |
1234: $ | | |
1235: $ 4---13----5---14----6
1237: and for simplicial cells
1239: $ 14----8---15----9----16
1240: $ |\ 5 |\ 7 |
1241: $ | \ | \ |
1242: $ 13 2 14 3 15
1243: $ | 4 \ | 6 \ |
1244: $ | \ | \ |
1245: $ 11----6---12----7----13
1246: $ |\ |\ |
1247: $ | \ 1 | \ 3 |
1248: $ 10 0 11 1 12
1249: $ | 0 \ | 2 \ |
1250: $ | \ | \ |
1251: $ 8----4----9----5----10
1253: Level: beginner
1255: .seealso: DMSetFromOptions(), DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate()
1256: @*/
1257: PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm)
1258: {
1259: PetscInt fac[3] = {1, 1, 1};
1260: PetscReal low[3] = {0, 0, 0};
1261: PetscReal upp[3] = {1, 1, 1};
1262: DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
1266: DMCreate(comm,dm);
1267: DMSetType(*dm,DMPLEX);
1268: DMPlexCreateBoxMesh_Internal(*dm, dim, simplex, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt, interpolate);
1269: return(0);
1270: }
1272: static PetscErrorCode DMPlexCreateWedgeBoxMesh_Internal(DM dm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate)
1273: {
1274: DM bdm, vol;
1275: PetscReal normal[3] = {0., 0., 1.};
1276: PetscInt i;
1280: for (i = 0; i < 3; ++i) if (periodicity[i] != DM_BOUNDARY_NONE) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity not yet supported");
1281: DMCreate(PetscObjectComm((PetscObject) dm), &bdm);
1282: DMSetType(bdm, DMPLEX);
1283: DMSetDimension(bdm, 2);
1284: DMPlexCreateBoxMesh_Simplex_Internal(bdm, 2, faces, lower, upper, periodicity, interpolate);
1285: DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], orderHeight, normal, interpolate, &vol);
1286: DMDestroy(&bdm);
1287: DMPlexReplace_Static(dm, &vol);
1288: if (lower[2] != 0.0) {
1289: Vec v;
1290: PetscScalar *x;
1291: PetscInt cDim, n;
1293: DMGetCoordinatesLocal(dm, &v);
1294: VecGetBlockSize(v, &cDim);
1295: VecGetLocalSize(v, &n);
1296: VecGetArray(v, &x);
1297: x += cDim;
1298: for (i = 0; i < n; i += cDim) x[i] += lower[2];
1299: VecRestoreArray(v,&x);
1300: DMSetCoordinatesLocal(dm, v);
1301: }
1302: return(0);
1303: }
1305: /*@
1306: DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells.
1308: Collective
1310: Input Parameters:
1311: + comm - The communicator for the DM object
1312: . faces - Number of faces per dimension, or NULL for (1, 1, 1)
1313: . lower - The lower left corner, or NULL for (0, 0, 0)
1314: . upper - The upper right corner, or NULL for (1, 1, 1)
1315: . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE
1316: . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first
1317: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
1319: Output Parameter:
1320: . dm - The DM object
1322: Level: beginner
1324: .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1325: @*/
1326: PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm)
1327: {
1328: PetscInt fac[3] = {1, 1, 1};
1329: PetscReal low[3] = {0, 0, 0};
1330: PetscReal upp[3] = {1, 1, 1};
1331: DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
1335: DMCreate(comm,dm);
1336: DMSetType(*dm,DMPLEX);
1337: DMPlexCreateWedgeBoxMesh_Internal(*dm, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt, orderHeight, interpolate);
1338: return(0);
1339: }
1341: /*@C
1342: DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells.
1344: Collective on idm
1346: Input Parameters:
1347: + idm - The mesh to be extruded
1348: . layers - The number of layers, or PETSC_DETERMINE to use the default
1349: . height - The total height of the extrusion, or PETSC_DETERMINE to use the default
1350: . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first
1351: . extNormal - The normal direction in which the mesh should be extruded, or NULL to extrude using the surface normal
1352: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
1354: Output Parameter:
1355: . dm - The DM object
1357: Notes:
1358: The mesh created has prismatic cells, and the vertex ordering in the cone of the cell is that of the tensor prismatic cells. Not currently supported in Fortran.
1360: Options Database Keys:
1361: + -dm_plex_extrude_layers <k> - Sets the number of layers k
1362: . -dm_plex_extrude_height <h> - Sets the total height of the extrusion
1363: . -dm_plex_extrude_heights <h0,h1,...> - Sets the height of each layer
1364: . -dm_plex_extrude_order_height - If true, order cells by height first
1365: - -dm_plex_extrude_normal <n0,...,nd> - Sets the normal vector along which to extrude
1367: Level: advanced
1369: .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMSetType(), DMCreate()
1370: @*/
1371: PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool orderHeight, const PetscReal extNormal[], PetscBool interpolate, DM* dm)
1372: {
1373: PetscScalar *coordsB;
1374: const PetscScalar *coordsA;
1375: PetscReal *normals = NULL, *heights = NULL;
1376: PetscReal clNormal[3];
1377: Vec coordinatesA, coordinatesB;
1378: PetscSection coordSectionA, coordSectionB;
1379: PetscInt dim, cDim, cDimB, c, l, v, coordSize, *newCone, nl;
1380: PetscInt cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices;
1381: const char *prefix;
1382: PetscBool haveCLNormal, flg;
1383: PetscErrorCode ierr;
1390: DMGetDimension(idm, &dim);
1391: DMGetCoordinateDim(idm, &cDim);
1392: cDimB = cDim == dim ? cDim+1 : cDim;
1393: if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim);
1395: PetscObjectGetOptionsPrefix((PetscObject) idm, &prefix);
1396: if (layers < 0) layers = 1;
1397: PetscOptionsGetInt(NULL, prefix, "-dm_plex_extrude_layers", &layers, NULL);
1398: if (layers <= 0) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Number of layers %D must be positive", layers);
1399: if (height < 0.) height = 1.;
1400: PetscOptionsGetReal(NULL, prefix, "-dm_plex_extrude_height", &height, NULL);
1401: if (height <= 0.) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double) height);
1402: PetscMalloc1(layers, &heights);
1403: nl = layers;
1404: PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_heights", heights, &nl, &flg);
1405: if (flg) {
1406: if (!nl) SETERRQ(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one height for -dm_plex_extrude_heights");
1407: for (l = nl; l < layers; ++l) heights[l] = heights[l-1];
1408: for (l = 0; l < layers; ++l) if (heights[l] <= 0.) SETERRQ2(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height %g of layers %D must be positive", (double) heights[l], l);
1409: } else {
1410: for (l = 0; l < layers; ++l) heights[l] = height/layers;
1411: }
1412: PetscOptionsGetBool(NULL, prefix, "-dm_plex_extrude_order_height", &orderHeight, NULL);
1413: c = 3;
1414: PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_normal", clNormal, &c, &haveCLNormal);
1415: if (haveCLNormal && c != cDimB) SETERRQ2(PetscObjectComm((PetscObject)idm), PETSC_ERR_ARG_SIZ, "Input normal has size %D != %D extruded coordinate dimension", c, cDimB);
1417: DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);
1418: DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);
1419: numCells = (cEnd - cStart)*layers;
1420: numVertices = (vEnd - vStart)*(layers+1);
1421: DMCreate(PetscObjectComm((PetscObject)idm), dm);
1422: DMSetType(*dm, DMPLEX);
1423: DMSetDimension(*dm, dim+1);
1424: DMPlexSetChart(*dm, 0, numCells+numVertices);
1425: /* Must create the celltype label here so that we do not automatically try to compute the types */
1426: DMCreateLabel(*dm, "celltype");
1427: for (c = cStart, cellV = 0; c < cEnd; ++c) {
1428: DMPolytopeType ct, nct;
1429: PetscInt *closure = NULL;
1430: PetscInt closureSize, numCorners = 0;
1432: DMPlexGetCellType(idm, c, &ct);
1433: switch (ct) {
1434: case DM_POLYTOPE_SEGMENT: nct = DM_POLYTOPE_SEG_PRISM_TENSOR;break;
1435: case DM_POLYTOPE_TRIANGLE: nct = DM_POLYTOPE_TRI_PRISM_TENSOR;break;
1436: case DM_POLYTOPE_QUADRILATERAL: nct = DM_POLYTOPE_QUAD_PRISM_TENSOR;break;
1437: default: nct = DM_POLYTOPE_UNKNOWN;
1438: }
1439: DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);
1440: for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++;
1441: DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);
1442: for (l = 0; l < layers; ++l) {
1443: const PetscInt cell = orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart;
1445: DMPlexSetConeSize(*dm, cell, 2*numCorners);
1446: DMPlexSetCellType(*dm, cell, nct);
1447: }
1448: cellV = PetscMax(numCorners,cellV);
1449: }
1450: DMSetUp(*dm);
1452: if (dim != cDim && !(extNormal || haveCLNormal)) {PetscCalloc1(cDim*(vEnd - vStart), &normals);}
1453: PetscMalloc1(3*cellV,&newCone);
1454: for (c = cStart; c < cEnd; ++c) {
1455: PetscInt *closure = NULL;
1456: PetscInt closureSize, numCorners = 0, l;
1457: PetscReal normal[3] = {0, 0, 0};
1459: if (normals) {DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);}
1460: DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);
1461: for (v = 0; v < closureSize*2; v += 2) {
1462: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
1463: PetscInt d;
1465: newCone[numCorners++] = closure[v] - vStart;
1466: if (normals) {for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d];}
1467: }
1468: }
1469: DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);
1470: for (l = 0; l < layers; ++l) {
1471: PetscInt i;
1473: for (i = 0; i < numCorners; ++i) {
1474: newCone[ numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + numCells : l*(vEnd - vStart) + newCone[i] + numCells;
1475: newCone[2*numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells;
1476: }
1477: DMPlexSetCone(*dm, orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);
1478: }
1479: }
1480: DMPlexSymmetrize(*dm);
1481: DMPlexStratify(*dm);
1482: PetscFree(newCone);
1484: DMGetCoordinateSection(*dm, &coordSectionB);
1485: PetscSectionSetNumFields(coordSectionB, 1);
1486: PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);
1487: PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);
1488: for (v = numCells; v < numCells+numVertices; ++v) {
1489: PetscSectionSetDof(coordSectionB, v, cDimB);
1490: PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);
1491: DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1492: }
1493: PetscSectionSetUp(coordSectionB);
1494: PetscSectionGetStorageSize(coordSectionB, &coordSize);
1495: VecCreate(PETSC_COMM_SELF, &coordinatesB);
1496: PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
1497: VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);
1498: VecSetBlockSize(coordinatesB, cDimB);
1499: VecSetType(coordinatesB,VECSTANDARD);
1501: DMGetCoordinateSection(idm, &coordSectionA);
1502: DMGetCoordinatesLocal(idm, &coordinatesA);
1503: VecGetArray(coordinatesB, &coordsB);
1504: VecGetArrayRead(coordinatesA, &coordsA);
1505: for (v = vStart; v < vEnd; ++v) {
1506: const PetscScalar *cptr;
1507: PetscReal ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.};
1508: PetscReal normal[3];
1509: PetscReal norm;
1510: PetscInt offA, d, cDimA = cDim;
1512: if (normals) {for (d = 0; d < cDimB; ++d) normal[d] = normals[cDimB*(v - vStart)+d];}
1513: else if (haveCLNormal) {for (d = 0; d < cDimB; ++d) normal[d] = clNormal[d];}
1514: else if (extNormal) {for (d = 0; d < cDimB; ++d) normal[d] = extNormal[d];}
1515: else if (cDimB == 2) {for (d = 0; d < cDimB; ++d) normal[d] = ones2[d];}
1516: else if (cDimB == 3) {for (d = 0; d < cDimB; ++d) normal[d] = ones3[d];}
1517: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
1518: for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d];
1519: for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm);
1521: PetscSectionGetOffset(coordSectionA, v, &offA);
1522: cptr = coordsA + offA;
1523: for (l = 0; l <= layers; ++l) {
1524: PetscInt offB, d, newV;
1526: newV = orderHeight ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells;
1527: PetscSectionGetOffset(coordSectionB, newV, &offB);
1528: for (d = 0; d < cDimA; ++d) { coordsB[offB+d] = cptr[d]; }
1529: for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*heights[l-1] : 0.0; }
1530: cptr = coordsB + offB;
1531: cDimA = cDimB;
1532: }
1533: }
1534: VecRestoreArrayRead(coordinatesA, &coordsA);
1535: VecRestoreArray(coordinatesB, &coordsB);
1536: DMSetCoordinatesLocal(*dm, coordinatesB);
1537: VecDestroy(&coordinatesB);
1538: PetscFree(normals);
1539: PetscFree(heights);
1540: if (interpolate) {
1541: DM idm;
1543: DMPlexInterpolate(*dm, &idm);
1544: DMPlexCopyCoordinates(*dm, idm);
1545: DMDestroy(dm);
1546: *dm = idm;
1547: }
1548: return(0);
1549: }
1551: /*@C
1552: DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database.
1554: Logically Collective on dm
1556: Input Parameters:
1557: + dm - the DM context
1558: - prefix - the prefix to prepend to all option names
1560: Notes:
1561: A hyphen (-) must NOT be given at the beginning of the prefix name.
1562: The first character of all runtime options is AUTOMATICALLY the hyphen.
1564: Level: advanced
1566: .seealso: SNESSetFromOptions()
1567: @*/
1568: PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[])
1569: {
1570: DM_Plex *mesh = (DM_Plex *) dm->data;
1575: PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);
1576: PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);
1577: return(0);
1578: }
1580: /* Remap geometry to cylinder
1581: TODO: This only works for a single refinement, then it is broken
1583: Interior square: Linear interpolation is correct
1584: The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing
1585: such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance
1587: phi = arctan(y/x)
1588: d_close = sqrt(1/8 + 1/4 sin^2(phi))
1589: d_far = sqrt(1/2 + sin^2(phi))
1591: so we remap them using
1593: x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close)
1594: y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close)
1596: If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y.
1597: */
1598: static void snapToCylinder(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1599: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1600: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1601: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1602: {
1603: const PetscReal dis = 1.0/PetscSqrtReal(2.0);
1604: const PetscReal ds2 = 0.5*dis;
1606: if ((PetscAbsScalar(u[0]) <= ds2) && (PetscAbsScalar(u[1]) <= ds2)) {
1607: f0[0] = u[0];
1608: f0[1] = u[1];
1609: } else {
1610: PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc;
1612: x = PetscRealPart(u[0]);
1613: y = PetscRealPart(u[1]);
1614: phi = PetscAtan2Real(y, x);
1615: sinp = PetscSinReal(phi);
1616: cosp = PetscCosReal(phi);
1617: if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) {
1618: dc = PetscAbsReal(ds2/sinp);
1619: df = PetscAbsReal(dis/sinp);
1620: xc = ds2*x/PetscAbsReal(y);
1621: yc = ds2*PetscSignReal(y);
1622: } else {
1623: dc = PetscAbsReal(ds2/cosp);
1624: df = PetscAbsReal(dis/cosp);
1625: xc = ds2*PetscSignReal(x);
1626: yc = ds2*y/PetscAbsReal(x);
1627: }
1628: f0[0] = xc + (u[0] - xc)*(1.0 - dc)/(df - dc);
1629: f0[1] = yc + (u[1] - yc)*(1.0 - dc)/(df - dc);
1630: }
1631: f0[2] = u[2];
1632: }
1634: static PetscErrorCode DMPlexCreateHexCylinderMesh_Internal(DM dm, DMBoundaryType periodicZ)
1635: {
1636: const PetscInt dim = 3;
1637: PetscInt numCells, numVertices;
1638: PetscMPIInt rank;
1642: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1643: DMSetDimension(dm, dim);
1644: /* Create topology */
1645: {
1646: PetscInt cone[8], c;
1648: numCells = rank == 0 ? 5 : 0;
1649: numVertices = rank == 0 ? 16 : 0;
1650: if (periodicZ == DM_BOUNDARY_PERIODIC) {
1651: numCells *= 3;
1652: numVertices = rank == 0 ? 24 : 0;
1653: }
1654: DMPlexSetChart(dm, 0, numCells+numVertices);
1655: for (c = 0; c < numCells; c++) {DMPlexSetConeSize(dm, c, 8);}
1656: DMSetUp(dm);
1657: if (rank == 0) {
1658: if (periodicZ == DM_BOUNDARY_PERIODIC) {
1659: cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16;
1660: cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34;
1661: DMPlexSetCone(dm, 0, cone);
1662: cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23;
1663: cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */
1664: DMPlexSetCone(dm, 1, cone);
1665: cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17;
1666: cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38;
1667: DMPlexSetCone(dm, 2, cone);
1668: cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15;
1669: cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38;
1670: DMPlexSetCone(dm, 3, cone);
1671: cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23;
1672: cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31;
1673: DMPlexSetCone(dm, 4, cone);
1675: cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32;
1676: cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20;
1677: DMPlexSetCone(dm, 5, cone);
1678: cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36;
1679: cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21;
1680: DMPlexSetCone(dm, 6, cone);
1681: cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33;
1682: cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28;
1683: DMPlexSetCone(dm, 7, cone);
1684: cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31;
1685: cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28;
1686: DMPlexSetCone(dm, 8, cone);
1687: cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36;
1688: cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19;
1689: DMPlexSetCone(dm, 9, cone);
1691: cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22;
1692: cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18;
1693: DMPlexSetCone(dm, 10, cone);
1694: cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25;
1695: cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17;
1696: DMPlexSetCone(dm, 11, cone);
1697: cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21;
1698: cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27;
1699: DMPlexSetCone(dm, 12, cone);
1700: cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19;
1701: cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27;
1702: DMPlexSetCone(dm, 13, cone);
1703: cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25;
1704: cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15;
1705: DMPlexSetCone(dm, 14, cone);
1706: } else {
1707: cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6;
1708: cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10;
1709: DMPlexSetCone(dm, 0, cone);
1710: cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13;
1711: cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11;
1712: DMPlexSetCone(dm, 1, cone);
1713: cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7;
1714: cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18;
1715: DMPlexSetCone(dm, 2, cone);
1716: cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5;
1717: cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18;
1718: DMPlexSetCone(dm, 3, cone);
1719: cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13;
1720: cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9;
1721: DMPlexSetCone(dm, 4, cone);
1722: }
1723: }
1724: DMPlexSymmetrize(dm);
1725: DMPlexStratify(dm);
1726: }
1727: /* Create cube geometry */
1728: {
1729: Vec coordinates;
1730: PetscSection coordSection;
1731: PetscScalar *coords;
1732: PetscInt coordSize, v;
1733: const PetscReal dis = 1.0/PetscSqrtReal(2.0);
1734: const PetscReal ds2 = dis/2.0;
1736: /* Build coordinates */
1737: DMGetCoordinateSection(dm, &coordSection);
1738: PetscSectionSetNumFields(coordSection, 1);
1739: PetscSectionSetFieldComponents(coordSection, 0, dim);
1740: PetscSectionSetChart(coordSection, numCells, numCells+numVertices);
1741: for (v = numCells; v < numCells+numVertices; ++v) {
1742: PetscSectionSetDof(coordSection, v, dim);
1743: PetscSectionSetFieldDof(coordSection, v, 0, dim);
1744: }
1745: PetscSectionSetUp(coordSection);
1746: PetscSectionGetStorageSize(coordSection, &coordSize);
1747: VecCreate(PETSC_COMM_SELF, &coordinates);
1748: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1749: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1750: VecSetBlockSize(coordinates, dim);
1751: VecSetType(coordinates,VECSTANDARD);
1752: VecGetArray(coordinates, &coords);
1753: if (rank == 0) {
1754: coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0;
1755: coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0;
1756: coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0;
1757: coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0;
1758: coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0;
1759: coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0;
1760: coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0;
1761: coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0;
1762: coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0;
1763: coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0;
1764: coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0;
1765: coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0;
1766: coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0;
1767: coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0;
1768: coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0;
1769: coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0;
1770: if (periodicZ == DM_BOUNDARY_PERIODIC) {
1771: /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5;
1772: /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5;
1773: /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5;
1774: /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5;
1775: /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5;
1776: /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5;
1777: /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5;
1778: /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5;
1779: }
1780: }
1781: VecRestoreArray(coordinates, &coords);
1782: DMSetCoordinatesLocal(dm, coordinates);
1783: VecDestroy(&coordinates);
1784: }
1785: /* Create periodicity */
1786: if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) {
1787: PetscReal L[3];
1788: PetscReal maxCell[3];
1789: DMBoundaryType bdType[3];
1790: PetscReal lower[3] = {0.0, 0.0, 0.0};
1791: PetscReal upper[3] = {1.0, 1.0, 1.5};
1792: PetscInt i, numZCells = 3;
1794: bdType[0] = DM_BOUNDARY_NONE;
1795: bdType[1] = DM_BOUNDARY_NONE;
1796: bdType[2] = periodicZ;
1797: for (i = 0; i < dim; i++) {
1798: L[i] = upper[i] - lower[i];
1799: maxCell[i] = 1.1 * (L[i] / numZCells);
1800: }
1801: DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, bdType);
1802: }
1803: {
1804: DM cdm;
1805: PetscDS cds;
1806: PetscScalar c[2] = {1.0, 1.0};
1808: DMPlexCreateCoordinateSpace(dm, 1, snapToCylinder);
1809: DMGetCoordinateDM(dm, &cdm);
1810: DMGetDS(cdm, &cds);
1811: PetscDSSetConstants(cds, 2, c);
1812: }
1813: /* Wait for coordinate creation before doing in-place modification */
1814: DMPlexInterpolateInPlace_Internal(dm);
1815: return(0);
1816: }
1818: /*@
1819: DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra.
1821: Collective
1823: Input Parameters:
1824: + comm - The communicator for the DM object
1825: - periodicZ - The boundary type for the Z direction
1827: Output Parameter:
1828: . dm - The DM object
1830: Note:
1831: Here is the output numbering looking from the bottom of the cylinder:
1832: $ 17-----14
1833: $ | |
1834: $ | 2 |
1835: $ | |
1836: $ 17-----8-----7-----14
1837: $ | | | |
1838: $ | 3 | 0 | 1 |
1839: $ | | | |
1840: $ 19-----5-----6-----13
1841: $ | |
1842: $ | 4 |
1843: $ | |
1844: $ 19-----13
1845: $
1846: $ and up through the top
1847: $
1848: $ 18-----16
1849: $ | |
1850: $ | 2 |
1851: $ | |
1852: $ 18----10----11-----16
1853: $ | | | |
1854: $ | 3 | 0 | 1 |
1855: $ | | | |
1856: $ 20-----9----12-----15
1857: $ | |
1858: $ | 4 |
1859: $ | |
1860: $ 20-----15
1862: Level: beginner
1864: .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1865: @*/
1866: PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, DMBoundaryType periodicZ, DM *dm)
1867: {
1872: DMCreate(comm, dm);
1873: DMSetType(*dm, DMPLEX);
1874: DMPlexCreateHexCylinderMesh_Internal(*dm, periodicZ);
1875: return(0);
1876: }
1878: static PetscErrorCode DMPlexCreateWedgeCylinderMesh_Internal(DM dm, PetscInt n, PetscBool interpolate)
1879: {
1880: const PetscInt dim = 3;
1881: PetscInt numCells, numVertices, v;
1882: PetscMPIInt rank;
1886: if (n < 0) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n);
1887: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1888: DMSetDimension(dm, dim);
1889: /* Must create the celltype label here so that we do not automatically try to compute the types */
1890: DMCreateLabel(dm, "celltype");
1891: /* Create topology */
1892: {
1893: PetscInt cone[6], c;
1895: numCells = rank == 0 ? n : 0;
1896: numVertices = rank == 0 ? 2*(n+1) : 0;
1897: DMPlexSetChart(dm, 0, numCells+numVertices);
1898: for (c = 0; c < numCells; c++) {DMPlexSetConeSize(dm, c, 6);}
1899: DMSetUp(dm);
1900: for (c = 0; c < numCells; c++) {
1901: cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n;
1902: cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n;
1903: DMPlexSetCone(dm, c, cone);
1904: DMPlexSetCellType(dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR);
1905: }
1906: DMPlexSymmetrize(dm);
1907: DMPlexStratify(dm);
1908: }
1909: for (v = numCells; v < numCells+numVertices; ++v) {
1910: DMPlexSetCellType(dm, v, DM_POLYTOPE_POINT);
1911: }
1912: /* Create cylinder geometry */
1913: {
1914: Vec coordinates;
1915: PetscSection coordSection;
1916: PetscScalar *coords;
1917: PetscInt coordSize, c;
1919: /* Build coordinates */
1920: DMGetCoordinateSection(dm, &coordSection);
1921: PetscSectionSetNumFields(coordSection, 1);
1922: PetscSectionSetFieldComponents(coordSection, 0, dim);
1923: PetscSectionSetChart(coordSection, numCells, numCells+numVertices);
1924: for (v = numCells; v < numCells+numVertices; ++v) {
1925: PetscSectionSetDof(coordSection, v, dim);
1926: PetscSectionSetFieldDof(coordSection, v, 0, dim);
1927: }
1928: PetscSectionSetUp(coordSection);
1929: PetscSectionGetStorageSize(coordSection, &coordSize);
1930: VecCreate(PETSC_COMM_SELF, &coordinates);
1931: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1932: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1933: VecSetBlockSize(coordinates, dim);
1934: VecSetType(coordinates,VECSTANDARD);
1935: VecGetArray(coordinates, &coords);
1936: for (c = 0; c < numCells; c++) {
1937: coords[(c+0*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+2] = 1.0;
1938: coords[(c+1*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+2] = 0.0;
1939: }
1940: if (rank == 0) {
1941: coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0;
1942: coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0;
1943: }
1944: VecRestoreArray(coordinates, &coords);
1945: DMSetCoordinatesLocal(dm, coordinates);
1946: VecDestroy(&coordinates);
1947: }
1948: /* Interpolate */
1949: if (interpolate) {DMPlexInterpolateInPlace_Internal(dm);}
1950: return(0);
1951: }
1953: /*@
1954: DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges.
1956: Collective
1958: Input Parameters:
1959: + comm - The communicator for the DM object
1960: . n - The number of wedges around the origin
1961: - interpolate - Create edges and faces
1963: Output Parameter:
1964: . dm - The DM object
1966: Level: beginner
1968: .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1969: @*/
1970: PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm)
1971: {
1976: DMCreate(comm, dm);
1977: DMSetType(*dm, DMPLEX);
1978: DMPlexCreateWedgeCylinderMesh_Internal(*dm, n, interpolate);
1979: return(0);
1980: }
1982: PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[])
1983: {
1984: PetscReal prod = 0.0;
1985: PetscInt i;
1986: for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]);
1987: return PetscSqrtReal(prod);
1988: }
1989: PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[])
1990: {
1991: PetscReal prod = 0.0;
1992: PetscInt i;
1993: for (i = 0; i < dim; ++i) prod += x[i]*y[i];
1994: return prod;
1995: }
1997: /* The first constant is the sphere radius */
1998: static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1999: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2000: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2001: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
2002: {
2003: PetscReal r = PetscRealPart(constants[0]);
2004: PetscReal norm2 = 0.0, fac;
2005: PetscInt n = uOff[1] - uOff[0], d;
2007: for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d]));
2008: fac = r/PetscSqrtReal(norm2);
2009: for (d = 0; d < n; ++d) f0[d] = u[d]*fac;
2010: }
2012: static PetscErrorCode DMPlexCreateSphereMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, PetscReal R)
2013: {
2014: const PetscInt embedDim = dim+1;
2015: PetscSection coordSection;
2016: Vec coordinates;
2017: PetscScalar *coords;
2018: PetscReal *coordsIn;
2019: PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e;
2020: PetscMPIInt rank;
2021: PetscErrorCode ierr;
2025: DMSetDimension(dm, dim);
2026: DMSetCoordinateDim(dm, dim+1);
2027: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2028: switch (dim) {
2029: case 2:
2030: if (simplex) {
2031: const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI);
2032: const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius);
2033: const PetscInt degree = 5;
2034: PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)};
2035: PetscInt s[3] = {1, 1, 1};
2036: PetscInt cone[3];
2037: PetscInt *graph, p, i, j, k;
2039: vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius;
2040: numCells = rank == 0 ? 20 : 0;
2041: numVerts = rank == 0 ? 12 : 0;
2042: firstVertex = numCells;
2043: /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of
2045: (0, \pm 1/\phi+1, \pm \phi/\phi+1)
2047: where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge
2048: length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393.
2049: */
2050: /* Construct vertices */
2051: PetscCalloc1(numVerts * embedDim, &coordsIn);
2052: if (rank == 0) {
2053: for (p = 0, i = 0; p < embedDim; ++p) {
2054: for (s[1] = -1; s[1] < 2; s[1] += 2) {
2055: for (s[2] = -1; s[2] < 2; s[2] += 2) {
2056: for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim];
2057: ++i;
2058: }
2059: }
2060: }
2061: }
2062: /* Construct graph */
2063: PetscCalloc1(numVerts * numVerts, &graph);
2064: for (i = 0; i < numVerts; ++i) {
2065: for (j = 0, k = 0; j < numVerts; ++j) {
2066: if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;}
2067: }
2068: if (k != degree) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree);
2069: }
2070: /* Build Topology */
2071: DMPlexSetChart(dm, 0, numCells+numVerts);
2072: for (c = 0; c < numCells; c++) {
2073: DMPlexSetConeSize(dm, c, embedDim);
2074: }
2075: DMSetUp(dm); /* Allocate space for cones */
2076: /* Cells */
2077: for (i = 0, c = 0; i < numVerts; ++i) {
2078: for (j = 0; j < i; ++j) {
2079: for (k = 0; k < j; ++k) {
2080: if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) {
2081: cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k;
2082: /* Check orientation */
2083: {
2084: const PetscInt epsilon[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
2085: PetscReal normal[3];
2086: PetscInt e, f;
2088: for (d = 0; d < embedDim; ++d) {
2089: normal[d] = 0.0;
2090: for (e = 0; e < embedDim; ++e) {
2091: for (f = 0; f < embedDim; ++f) {
2092: normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]);
2093: }
2094: }
2095: }
2096: if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;}
2097: }
2098: DMPlexSetCone(dm, c++, cone);
2099: }
2100: }
2101: }
2102: }
2103: DMPlexSymmetrize(dm);
2104: DMPlexStratify(dm);
2105: PetscFree(graph);
2106: } else {
2107: /*
2108: 12-21--13
2109: | |
2110: 25 4 24
2111: | |
2112: 12-25--9-16--8-24--13
2113: | | | |
2114: 23 5 17 0 15 3 22
2115: | | | |
2116: 10-20--6-14--7-19--11
2117: | |
2118: 20 1 19
2119: | |
2120: 10-18--11
2121: | |
2122: 23 2 22
2123: | |
2124: 12-21--13
2125: */
2126: PetscInt cone[4], ornt[4];
2128: numCells = rank == 0 ? 6 : 0;
2129: numEdges = rank == 0 ? 12 : 0;
2130: numVerts = rank == 0 ? 8 : 0;
2131: firstVertex = numCells;
2132: firstEdge = numCells + numVerts;
2133: /* Build Topology */
2134: DMPlexSetChart(dm, 0, numCells+numEdges+numVerts);
2135: for (c = 0; c < numCells; c++) {
2136: DMPlexSetConeSize(dm, c, 4);
2137: }
2138: for (e = firstEdge; e < firstEdge+numEdges; ++e) {
2139: DMPlexSetConeSize(dm, e, 2);
2140: }
2141: DMSetUp(dm); /* Allocate space for cones */
2142: if (rank == 0) {
2143: /* Cell 0 */
2144: cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17;
2145: DMPlexSetCone(dm, 0, cone);
2146: ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0;
2147: DMPlexSetConeOrientation(dm, 0, ornt);
2148: /* Cell 1 */
2149: cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20;
2150: DMPlexSetCone(dm, 1, cone);
2151: ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0;
2152: DMPlexSetConeOrientation(dm, 1, ornt);
2153: /* Cell 2 */
2154: cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23;
2155: DMPlexSetCone(dm, 2, cone);
2156: ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0;
2157: DMPlexSetConeOrientation(dm, 2, ornt);
2158: /* Cell 3 */
2159: cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15;
2160: DMPlexSetCone(dm, 3, cone);
2161: ornt[0] = -1; ornt[1] = -1; ornt[2] = 0; ornt[3] = -1;
2162: DMPlexSetConeOrientation(dm, 3, ornt);
2163: /* Cell 4 */
2164: cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25;
2165: DMPlexSetCone(dm, 4, cone);
2166: ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = 0;
2167: DMPlexSetConeOrientation(dm, 4, ornt);
2168: /* Cell 5 */
2169: cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23;
2170: DMPlexSetCone(dm, 5, cone);
2171: ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = -1;
2172: DMPlexSetConeOrientation(dm, 5, ornt);
2173: /* Edges */
2174: cone[0] = 6; cone[1] = 7;
2175: DMPlexSetCone(dm, 14, cone);
2176: cone[0] = 7; cone[1] = 8;
2177: DMPlexSetCone(dm, 15, cone);
2178: cone[0] = 8; cone[1] = 9;
2179: DMPlexSetCone(dm, 16, cone);
2180: cone[0] = 9; cone[1] = 6;
2181: DMPlexSetCone(dm, 17, cone);
2182: cone[0] = 10; cone[1] = 11;
2183: DMPlexSetCone(dm, 18, cone);
2184: cone[0] = 11; cone[1] = 7;
2185: DMPlexSetCone(dm, 19, cone);
2186: cone[0] = 6; cone[1] = 10;
2187: DMPlexSetCone(dm, 20, cone);
2188: cone[0] = 12; cone[1] = 13;
2189: DMPlexSetCone(dm, 21, cone);
2190: cone[0] = 13; cone[1] = 11;
2191: DMPlexSetCone(dm, 22, cone);
2192: cone[0] = 10; cone[1] = 12;
2193: DMPlexSetCone(dm, 23, cone);
2194: cone[0] = 13; cone[1] = 8;
2195: DMPlexSetCone(dm, 24, cone);
2196: cone[0] = 12; cone[1] = 9;
2197: DMPlexSetCone(dm, 25, cone);
2198: }
2199: DMPlexSymmetrize(dm);
2200: DMPlexStratify(dm);
2201: /* Build coordinates */
2202: PetscCalloc1(numVerts * embedDim, &coordsIn);
2203: if (rank == 0) {
2204: coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R;
2205: coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R;
2206: coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R;
2207: coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R;
2208: coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R;
2209: coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R;
2210: coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R;
2211: coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R;
2212: }
2213: }
2214: break;
2215: case 3:
2216: if (simplex) {
2217: const PetscReal edgeLen = 1.0/PETSC_PHI;
2218: PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5};
2219: PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0};
2220: PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0};
2221: const PetscInt degree = 12;
2222: PetscInt s[4] = {1, 1, 1};
2223: PetscInt evenPerm[12][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 0, 3, 2}, {1, 2, 0, 3}, {1, 3, 2, 0},
2224: {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}};
2225: PetscInt cone[4];
2226: PetscInt *graph, p, i, j, k, l;
2228: vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R;
2229: vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R;
2230: vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R;
2231: numCells = rank == 0 ? 600 : 0;
2232: numVerts = rank == 0 ? 120 : 0;
2233: firstVertex = numCells;
2234: /* Use the 600-cell, which for a unit sphere has coordinates which are
2236: 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16
2237: (\pm 1, 0, 0, 0) all cyclic permutations 8
2238: 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96
2240: where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge
2241: length is then given by 1/\phi = 0.61803.
2243: http://buzzard.pugetsound.edu/sage-practice/ch03s03.html
2244: http://mathworld.wolfram.com/600-Cell.html
2245: */
2246: /* Construct vertices */
2247: PetscCalloc1(numVerts * embedDim, &coordsIn);
2248: i = 0;
2249: if (rank == 0) {
2250: for (s[0] = -1; s[0] < 2; s[0] += 2) {
2251: for (s[1] = -1; s[1] < 2; s[1] += 2) {
2252: for (s[2] = -1; s[2] < 2; s[2] += 2) {
2253: for (s[3] = -1; s[3] < 2; s[3] += 2) {
2254: for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d];
2255: ++i;
2256: }
2257: }
2258: }
2259: }
2260: for (p = 0; p < embedDim; ++p) {
2261: s[1] = s[2] = s[3] = 1;
2262: for (s[0] = -1; s[0] < 2; s[0] += 2) {
2263: for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim];
2264: ++i;
2265: }
2266: }
2267: for (p = 0; p < 12; ++p) {
2268: s[3] = 1;
2269: for (s[0] = -1; s[0] < 2; s[0] += 2) {
2270: for (s[1] = -1; s[1] < 2; s[1] += 2) {
2271: for (s[2] = -1; s[2] < 2; s[2] += 2) {
2272: for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]];
2273: ++i;
2274: }
2275: }
2276: }
2277: }
2278: }
2279: if (i != numVerts) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts);
2280: /* Construct graph */
2281: PetscCalloc1(numVerts * numVerts, &graph);
2282: for (i = 0; i < numVerts; ++i) {
2283: for (j = 0, k = 0; j < numVerts; ++j) {
2284: if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;}
2285: }
2286: if (k != degree) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree);
2287: }
2288: /* Build Topology */
2289: DMPlexSetChart(dm, 0, numCells+numVerts);
2290: for (c = 0; c < numCells; c++) {
2291: DMPlexSetConeSize(dm, c, embedDim);
2292: }
2293: DMSetUp(dm); /* Allocate space for cones */
2294: /* Cells */
2295: if (rank == 0) {
2296: for (i = 0, c = 0; i < numVerts; ++i) {
2297: for (j = 0; j < i; ++j) {
2298: for (k = 0; k < j; ++k) {
2299: for (l = 0; l < k; ++l) {
2300: if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] &&
2301: graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) {
2302: cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l;
2303: /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */
2304: {
2305: const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}},
2306: {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}},
2307: {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}},
2308: {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}},
2310: {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}},
2311: {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}},
2312: {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}},
2313: {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}},
2315: {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}},
2316: {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}},
2317: {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}},
2318: {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}},
2320: {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}},
2321: {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}},
2322: {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}},
2323: {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}};
2324: PetscReal normal[4];
2325: PetscInt e, f, g;
2327: for (d = 0; d < embedDim; ++d) {
2328: normal[d] = 0.0;
2329: for (e = 0; e < embedDim; ++e) {
2330: for (f = 0; f < embedDim; ++f) {
2331: for (g = 0; g < embedDim; ++g) {
2332: normal[d] += epsilon[d][e][f][g]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f])*(coordsIn[l*embedDim+f] - coordsIn[i*embedDim+f]);
2333: }
2334: }
2335: }
2336: }
2337: if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;}
2338: }
2339: DMPlexSetCone(dm, c++, cone);
2340: }
2341: }
2342: }
2343: }
2344: }
2345: }
2346: DMPlexSymmetrize(dm);
2347: DMPlexStratify(dm);
2348: PetscFree(graph);
2349: break;
2350: }
2351: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim);
2352: }
2353: /* Create coordinates */
2354: DMGetCoordinateSection(dm, &coordSection);
2355: PetscSectionSetNumFields(coordSection, 1);
2356: PetscSectionSetFieldComponents(coordSection, 0, embedDim);
2357: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);
2358: for (v = firstVertex; v < firstVertex+numVerts; ++v) {
2359: PetscSectionSetDof(coordSection, v, embedDim);
2360: PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
2361: }
2362: PetscSectionSetUp(coordSection);
2363: PetscSectionGetStorageSize(coordSection, &coordSize);
2364: VecCreate(PETSC_COMM_SELF, &coordinates);
2365: VecSetBlockSize(coordinates, embedDim);
2366: PetscObjectSetName((PetscObject) coordinates, "coordinates");
2367: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
2368: VecSetType(coordinates,VECSTANDARD);
2369: VecGetArray(coordinates, &coords);
2370: for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];}
2371: VecRestoreArray(coordinates, &coords);
2372: DMSetCoordinatesLocal(dm, coordinates);
2373: VecDestroy(&coordinates);
2374: PetscFree(coordsIn);
2375: {
2376: DM cdm;
2377: PetscDS cds;
2378: PetscScalar c = R;
2380: DMPlexCreateCoordinateSpace(dm, 1, snapToSphere);
2381: DMGetCoordinateDM(dm, &cdm);
2382: DMGetDS(cdm, &cds);
2383: PetscDSSetConstants(cds, 1, &c);
2384: }
2385: /* Wait for coordinate creation before doing in-place modification */
2386: if (simplex) {DMPlexInterpolateInPlace_Internal(dm);}
2387: return(0);
2388: }
2390: /*@
2391: DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d.
2393: Collective
2395: Input Parameters:
2396: + comm - The communicator for the DM object
2397: . dim - The dimension
2398: . simplex - Use simplices, or tensor product cells
2399: - R - The radius
2401: Output Parameter:
2402: . dm - The DM object
2404: Level: beginner
2406: .seealso: DMPlexCreateBallMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
2407: @*/
2408: PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm)
2409: {
2414: DMCreate(comm, dm);
2415: DMSetType(*dm, DMPLEX);
2416: DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R);
2417: return(0);
2418: }
2420: static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R)
2421: {
2422: DM sdm, vol;
2423: DMLabel bdlabel;
2427: DMCreate(PetscObjectComm((PetscObject) dm), &sdm);
2428: DMSetType(sdm, DMPLEX);
2429: PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_");
2430: DMPlexCreateSphereMesh_Internal(sdm, dim-1, PETSC_TRUE, R);
2431: DMSetFromOptions(sdm);
2432: DMViewFromOptions(sdm, NULL, "-dm_view");
2433: DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol);
2434: DMDestroy(&sdm);
2435: DMPlexReplace_Static(dm, &vol);
2436: DMCreateLabel(dm, "marker");
2437: DMGetLabel(dm, "marker", &bdlabel);
2438: DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel);
2439: DMPlexLabelComplete(dm, bdlabel);
2440: return(0);
2441: }
2443: /*@
2444: DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d.
2446: Collective
2448: Input Parameters:
2449: + comm - The communicator for the DM object
2450: . dim - The dimension
2451: - R - The radius
2453: Output Parameter:
2454: . dm - The DM object
2456: Options Database Keys:
2457: - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry
2459: Level: beginner
2461: .seealso: DMPlexCreateSphereMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
2462: @*/
2463: PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm)
2464: {
2468: DMCreate(comm, dm);
2469: DMSetType(*dm, DMPLEX);
2470: DMPlexCreateBallMesh_Internal(*dm, dim, R);
2471: return(0);
2472: }
2474: static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct)
2475: {
2479: switch (ct) {
2480: case DM_POLYTOPE_POINT:
2481: {
2482: PetscInt numPoints[1] = {1};
2483: PetscInt coneSize[1] = {0};
2484: PetscInt cones[1] = {0};
2485: PetscInt coneOrientations[1] = {0};
2486: PetscScalar vertexCoords[1] = {0.0};
2488: DMSetDimension(rdm, 0);
2489: DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2490: }
2491: break;
2492: case DM_POLYTOPE_SEGMENT:
2493: {
2494: PetscInt numPoints[2] = {2, 1};
2495: PetscInt coneSize[3] = {2, 0, 0};
2496: PetscInt cones[2] = {1, 2};
2497: PetscInt coneOrientations[2] = {0, 0};
2498: PetscScalar vertexCoords[2] = {-1.0, 1.0};
2500: DMSetDimension(rdm, 1);
2501: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2502: }
2503: break;
2504: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2505: {
2506: PetscInt numPoints[2] = {2, 1};
2507: PetscInt coneSize[3] = {2, 0, 0};
2508: PetscInt cones[2] = {1, 2};
2509: PetscInt coneOrientations[2] = {0, 0};
2510: PetscScalar vertexCoords[2] = {-1.0, 1.0};
2512: DMSetDimension(rdm, 1);
2513: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2514: }
2515: break;
2516: case DM_POLYTOPE_TRIANGLE:
2517: {
2518: PetscInt numPoints[2] = {3, 1};
2519: PetscInt coneSize[4] = {3, 0, 0, 0};
2520: PetscInt cones[3] = {1, 2, 3};
2521: PetscInt coneOrientations[3] = {0, 0, 0};
2522: PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0};
2524: DMSetDimension(rdm, 2);
2525: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2526: }
2527: break;
2528: case DM_POLYTOPE_QUADRILATERAL:
2529: {
2530: PetscInt numPoints[2] = {4, 1};
2531: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
2532: PetscInt cones[4] = {1, 2, 3, 4};
2533: PetscInt coneOrientations[4] = {0, 0, 0, 0};
2534: PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
2536: DMSetDimension(rdm, 2);
2537: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2538: }
2539: break;
2540: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2541: {
2542: PetscInt numPoints[2] = {4, 1};
2543: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
2544: PetscInt cones[4] = {1, 2, 3, 4};
2545: PetscInt coneOrientations[4] = {0, 0, 0, 0};
2546: PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0};
2548: DMSetDimension(rdm, 2);
2549: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2550: }
2551: break;
2552: case DM_POLYTOPE_TETRAHEDRON:
2553: {
2554: PetscInt numPoints[2] = {4, 1};
2555: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
2556: PetscInt cones[4] = {1, 2, 3, 4};
2557: PetscInt coneOrientations[4] = {0, 0, 0, 0};
2558: 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};
2560: DMSetDimension(rdm, 3);
2561: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2562: }
2563: break;
2564: case DM_POLYTOPE_HEXAHEDRON:
2565: {
2566: PetscInt numPoints[2] = {8, 1};
2567: PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0};
2568: PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8};
2569: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
2570: 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,
2571: -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
2573: DMSetDimension(rdm, 3);
2574: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2575: }
2576: break;
2577: case DM_POLYTOPE_TRI_PRISM:
2578: {
2579: PetscInt numPoints[2] = {6, 1};
2580: PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0};
2581: PetscInt cones[6] = {1, 2, 3, 4, 5, 6};
2582: PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
2583: PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0,
2584: -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0};
2586: DMSetDimension(rdm, 3);
2587: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2588: }
2589: break;
2590: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2591: {
2592: PetscInt numPoints[2] = {6, 1};
2593: PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0};
2594: PetscInt cones[6] = {1, 2, 3, 4, 5, 6};
2595: PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
2596: PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0,
2597: -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0};
2599: DMSetDimension(rdm, 3);
2600: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2601: }
2602: break;
2603: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2604: {
2605: PetscInt numPoints[2] = {8, 1};
2606: PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0};
2607: PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8};
2608: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
2609: 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,
2610: -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
2612: DMSetDimension(rdm, 3);
2613: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2614: }
2615: break;
2616: case DM_POLYTOPE_PYRAMID:
2617: {
2618: PetscInt numPoints[2] = {5, 1};
2619: PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0};
2620: PetscInt cones[5] = {1, 2, 3, 4, 5};
2621: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
2622: 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,
2623: 0.0, 0.0, 1.0};
2625: DMSetDimension(rdm, 3);
2626: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
2627: }
2628: break;
2629: default: SETERRQ1(PetscObjectComm((PetscObject) rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]);
2630: }
2631: {
2632: PetscInt Nv, v;
2634: /* Must create the celltype label here so that we do not automatically try to compute the types */
2635: DMCreateLabel(rdm, "celltype");
2636: DMPlexSetCellType(rdm, 0, ct);
2637: DMPlexGetChart(rdm, NULL, &Nv);
2638: for (v = 1; v < Nv; ++v) {DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);}
2639: }
2640: DMPlexInterpolateInPlace_Internal(rdm);
2641: PetscObjectSetName((PetscObject) rdm, DMPolytopeTypes[ct]);
2642: return(0);
2643: }
2645: /*@
2646: DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
2648: Collective
2650: Input Parameters:
2651: + comm - The communicator
2652: - ct - The cell type of the reference cell
2654: Output Parameter:
2655: . refdm - The reference cell
2657: Level: intermediate
2659: .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh()
2660: @*/
2661: PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm)
2662: {
2666: DMCreate(comm, refdm);
2667: DMSetType(*refdm, DMPLEX);
2668: DMPlexCreateReferenceCell_Internal(*refdm, ct);
2669: return(0);
2670: }
2672: static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[])
2673: {
2674: DM plex;
2675: DMLabel label;
2676: PetscBool hasLabel;
2680: DMHasLabel(dm, name, &hasLabel);
2681: if (hasLabel) return(0);
2682: DMCreateLabel(dm, name);
2683: DMGetLabel(dm, name, &label);
2684: DMConvert(dm, DMPLEX, &plex);
2685: DMPlexMarkBoundaryFaces(plex, 1, label);
2686: DMDestroy(&plex);
2687: return(0);
2688: }
2690: const char * const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "unknown", "DMPlexShape", "DM_SHAPE_", NULL};
2692: static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm)
2693: {
2694: DMPlexShape shape = DM_SHAPE_BOX;
2695: DMPolytopeType cell = DM_POLYTOPE_TRIANGLE;
2696: PetscInt dim = 2;
2697: PetscBool simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE;
2698: PetscBool flg, flg2, fflg, bdfflg;
2699: MPI_Comm comm;
2700: char filename[PETSC_MAX_PATH_LEN], bdFilename[PETSC_MAX_PATH_LEN];
2704: PetscObjectGetComm((PetscObject) dm, &comm);
2705: /* TODO Turn this into a registration interface */
2706: PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg);
2707: PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg);
2708: PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum) cell, (PetscEnum *) &cell, NULL);
2709: PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL);
2710: PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum) shape, (PetscEnum *) &shape, &flg);
2711: PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0);
2712: if ((dim < 0) || (dim > 3)) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D should be in [1, 3]", dim);
2713: PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg);
2714: PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg);
2715: PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg);
2716: PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2);
2717: if (flg || flg2) {DMSetBasicAdjacency(dm, adjCone, adjClosure);}
2719: switch (cell) {
2720: case DM_POLYTOPE_POINT:
2721: case DM_POLYTOPE_SEGMENT:
2722: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2723: case DM_POLYTOPE_TRIANGLE:
2724: case DM_POLYTOPE_QUADRILATERAL:
2725: case DM_POLYTOPE_TETRAHEDRON:
2726: case DM_POLYTOPE_HEXAHEDRON:
2727: *useCoordSpace = PETSC_TRUE;break;
2728: default: *useCoordSpace = PETSC_FALSE;break;
2729: }
2731: if (fflg) {
2732: DM dmnew;
2734: DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), filename, interpolate, &dmnew);
2735: DMPlexReplace_Static(dm, &dmnew);
2736: } else if (refDomain) {
2737: DMPlexCreateReferenceCell_Internal(dm, cell);
2738: } else if (bdfflg) {
2739: DM bdm, dmnew;
2741: DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), bdFilename, interpolate, &bdm);
2742: PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_");
2743: DMSetFromOptions(bdm);
2744: DMPlexGenerate(bdm, NULL, interpolate, &dmnew);
2745: DMDestroy(&bdm);
2746: DMPlexReplace_Static(dm, &dmnew);
2747: } else {
2748: PetscObjectSetName((PetscObject) dm, DMPlexShapes[shape]);
2749: switch (shape) {
2750: case DM_SHAPE_BOX:
2751: {
2752: PetscInt faces[3] = {0, 0, 0};
2753: PetscReal lower[3] = {0, 0, 0};
2754: PetscReal upper[3] = {1, 1, 1};
2755: DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
2756: PetscInt i, n;
2758: n = dim;
2759: for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4-dim);
2760: PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);
2761: n = 3;
2762: PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);
2763: if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim);
2764: n = 3;
2765: PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);
2766: if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim);
2767: n = 3;
2768: PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg);
2769: if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %D values, should have been %D", n, dim);
2770: switch (cell) {
2771: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2772: DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt, PETSC_FALSE, interpolate);
2773: break;
2774: default:
2775: DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate);
2776: break;
2777: }
2778: }
2779: break;
2780: case DM_SHAPE_BOX_SURFACE:
2781: {
2782: PetscInt faces[3] = {0, 0, 0};
2783: PetscReal lower[3] = {0, 0, 0};
2784: PetscReal upper[3] = {1, 1, 1};
2785: PetscInt i, n;
2787: n = dim+1;
2788: for (i = 0; i < dim+1; ++i) faces[i] = (dim+1 == 1 ? 1 : 4-(dim+1));
2789: PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);
2790: n = 3;
2791: PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);
2792: if (flg && (n != dim+1)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim+1);
2793: n = 3;
2794: PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);
2795: if (flg && (n != dim+1)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim+1);
2796: DMPlexCreateBoxSurfaceMesh_Internal(dm, dim+1, faces, lower, upper, interpolate);
2797: }
2798: break;
2799: case DM_SHAPE_SPHERE:
2800: {
2801: PetscReal R = 1.0;
2803: PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg);
2804: DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R);
2805: }
2806: break;
2807: case DM_SHAPE_BALL:
2808: {
2809: PetscReal R = 1.0;
2811: PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg);
2812: DMPlexCreateBallMesh_Internal(dm, dim, R);
2813: }
2814: break;
2815: case DM_SHAPE_CYLINDER:
2816: {
2817: DMBoundaryType bdt = DM_BOUNDARY_NONE;
2818: PetscInt Nw = 6;
2820: PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum) bdt, (PetscEnum *) &bdt, NULL);
2821: PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL);
2822: switch (cell) {
2823: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2824: DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate);
2825: break;
2826: default:
2827: DMPlexCreateHexCylinderMesh_Internal(dm, bdt);
2828: break;
2829: }
2830: }
2831: break;
2832: default: SETERRQ1(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]);
2833: }
2834: }
2835: DMPlexSetRefinementUniform(dm, PETSC_TRUE);
2836: return(0);
2837: }
2839: PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject, DM dm)
2840: {
2841: DM_Plex *mesh = (DM_Plex*) dm->data;
2842: PetscBool flg;
2843: char bdLabel[PETSC_MAX_PATH_LEN];
2847: /* Handle viewing */
2848: PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);
2849: PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0);
2850: PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);
2851: PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0);
2852: DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);
2853: if (flg) {PetscLogDefaultBegin();}
2854: /* Labeling */
2855: PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg);
2856: if (flg) {DMPlexCreateBoundaryLabel_Private(dm, bdLabel);}
2857: /* Point Location */
2858: PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);
2859: /* Partitioning and distribution */
2860: PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);
2861: /* Generation and remeshing */
2862: PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);
2863: /* Projection behavior */
2864: PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0);
2865: PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);
2866: /* Checking structure */
2867: {
2868: PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE;
2870: PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2);
2871: PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);
2872: if (all || (flg && flg2)) {DMPlexCheckSymmetry(dm);}
2873: PetscOptionsBool("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices (only for homogeneous simplex or tensor meshes)", "DMPlexCheckSkeleton", PETSC_FALSE, &flg, &flg2);
2874: if (all || (flg && flg2)) {DMPlexCheckSkeleton(dm, 0);}
2875: PetscOptionsBool("-dm_plex_check_faces", "Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type", "DMPlexCheckFaces", PETSC_FALSE, &flg, &flg2);
2876: if (all || (flg && flg2)) {DMPlexCheckFaces(dm, 0);}
2877: PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);
2878: if (all || (flg && flg2)) {DMPlexCheckGeometry(dm);}
2879: PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);
2880: if (all || (flg && flg2)) {DMPlexCheckPointSF(dm);}
2881: PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2);
2882: if (all || (flg && flg2)) {DMPlexCheckInterfaceCones(dm);}
2883: PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);
2884: if (flg && flg2) {DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);}
2885: }
2886: {
2887: PetscReal scale = 1.0;
2889: PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg);
2890: if (flg) {
2891: Vec coordinates, coordinatesLocal;
2893: DMGetCoordinates(dm, &coordinates);
2894: DMGetCoordinatesLocal(dm, &coordinatesLocal);
2895: VecScale(coordinates, scale);
2896: VecScale(coordinatesLocal, scale);
2897: }
2898: }
2899: PetscPartitionerSetFromOptions(mesh->partitioner);
2900: return(0);
2901: }
2903: static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm)
2904: {
2905: PetscReal volume = -1.0, extThickness = 1.0;
2906: PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim;
2907: PetscBool uniformOrig, created = PETSC_FALSE, uniform = PETSC_TRUE, distribute = PETSC_FALSE, interpolate = PETSC_TRUE, extColOrder = PETSC_TRUE, coordSpace = PETSC_TRUE, remap = PETSC_TRUE, ghostCells = PETSC_FALSE, isHierarchy, ignoreModel = PETSC_FALSE, flg;
2912: PetscOptionsHead(PetscOptionsObject,"DMPlex Options");
2913: /* Handle automatic creation */
2914: DMGetDimension(dm, &dim);
2915: if (dim < 0) {DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm);created = PETSC_TRUE;}
2916: /* Handle interpolation before distribution */
2917: PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg);
2918: if (flg) {
2919: DMPlexInterpolatedFlag interpolated;
2921: DMPlexIsInterpolated(dm, &interpolated);
2922: if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) {
2923: DM udm;
2925: DMPlexUninterpolate(dm, &udm);
2926: DMPlexReplace_Static(dm, &udm);
2927: } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) {
2928: DM idm;
2930: DMPlexInterpolate(dm, &idm);
2931: DMPlexReplace_Static(dm, &idm);
2932: }
2933: }
2934: /* Handle DMPlex refinement before distribution */
2935: PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg);
2936: if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;}
2937: DMPlexGetRefinementUniform(dm, &uniformOrig);
2938: PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0);
2939: PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);
2940: PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);
2941: if (flg) {DMPlexSetRefinementUniform(dm, uniform);}
2942: PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);
2943: if (flg) {
2944: DMPlexSetRefinementUniform(dm, PETSC_FALSE);
2945: DMPlexSetRefinementLimit(dm, volume);
2946: prerefine = PetscMax(prerefine, 1);
2947: }
2948: for (r = 0; r < prerefine; ++r) {
2949: DM rdm;
2950: PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc;
2952: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
2953: DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);
2954: DMPlexReplace_Static(dm, &rdm);
2955: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
2956: if (coordFunc && remap) {
2957: DMPlexRemapGeometry(dm, 0.0, coordFunc);
2958: ((DM_Plex*) dm->data)->coordFunc = coordFunc;
2959: }
2960: }
2961: DMPlexSetRefinementUniform(dm, uniformOrig);
2962: /* Handle DMPlex extrusion before distribution */
2963: PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, NULL);
2964: PetscOptionsBoundedInt("-dm_extrude_layers", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0);
2965: PetscOptionsReal("-dm_extrude_thickness", "The thickness of the layer to be extruded", "", extThickness, &extThickness, NULL);
2966: PetscOptionsBool("-dm_extrude_column_first", "Order the cells in a vertical column first", "", extColOrder, &extColOrder, NULL);
2967: if (extLayers) {
2968: DM edm;
2970: DMPlexExtrude(dm, extLayers, extThickness, extColOrder, NULL, interpolate, &edm);
2971: DMPlexReplace_Static(dm, &edm);
2972: }
2973: /* Handle DMPlex distribution */
2974: PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL);
2975: PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0);
2976: if (distribute) {
2977: DM pdm = NULL;
2978: PetscPartitioner part;
2980: DMPlexGetPartitioner(dm, &part);
2981: PetscPartitionerSetFromOptions(part);
2982: DMPlexDistribute(dm, overlap, NULL, &pdm);
2983: if (pdm) {
2984: DMPlexReplace_Static(dm, &pdm);
2985: }
2986: }
2987: /* Create coordinate space */
2988: if (created) {
2989: DM_Plex *mesh = (DM_Plex *) dm->data;
2990: PetscInt degree = 1;
2991: PetscBool periodic, flg;
2993: PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg);
2994: PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL);
2995: if (coordSpace) {DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc);}
2996: if (flg && !coordSpace) {
2997: DM cdm;
2998: PetscDS cds;
2999: PetscObject obj;
3000: PetscClassId id;
3002: DMGetCoordinateDM(dm, &cdm);
3003: DMGetDS(cdm, &cds);
3004: PetscDSGetDiscretization(cds, 0, &obj);
3005: PetscObjectGetClassId(obj, &id);
3006: if (id == PETSCFE_CLASSID) {
3007: PetscContainer dummy;
3009: PetscContainerCreate(PETSC_COMM_SELF, &dummy);
3010: PetscObjectSetName((PetscObject) dummy, "coordinates");
3011: DMSetField(cdm, 0, NULL, (PetscObject) dummy);
3012: PetscContainerDestroy(&dummy);
3013: DMClearDS(cdm);
3014: }
3015: mesh->coordFunc = NULL;
3016: }
3017: DMLocalizeCoordinates(dm);
3018: DMGetPeriodicity(dm, &periodic, NULL, NULL, NULL);
3019: if (periodic) {DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, NULL);}
3020: }
3021: /* Handle DMPlex refinement */
3022: remap = PETSC_TRUE;
3023: PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0);
3024: PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);
3025: PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0);
3026: if (refine) {DMPlexSetRefinementUniform(dm, PETSC_TRUE);}
3027: if (refine && isHierarchy) {
3028: DM *dms, coarseDM;
3030: DMGetCoarseDM(dm, &coarseDM);
3031: PetscObjectReference((PetscObject)coarseDM);
3032: PetscMalloc1(refine,&dms);
3033: DMRefineHierarchy(dm, refine, dms);
3034: /* Total hack since we do not pass in a pointer */
3035: DMPlexSwap_Static(dm, dms[refine-1]);
3036: if (refine == 1) {
3037: DMSetCoarseDM(dm, dms[0]);
3038: DMPlexSetRegularRefinement(dm, PETSC_TRUE);
3039: } else {
3040: DMSetCoarseDM(dm, dms[refine-2]);
3041: DMPlexSetRegularRefinement(dm, PETSC_TRUE);
3042: DMSetCoarseDM(dms[0], dms[refine-1]);
3043: DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);
3044: }
3045: DMSetCoarseDM(dms[refine-1], coarseDM);
3046: PetscObjectDereference((PetscObject)coarseDM);
3047: /* Free DMs */
3048: for (r = 0; r < refine; ++r) {
3049: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);
3050: DMDestroy(&dms[r]);
3051: }
3052: PetscFree(dms);
3053: } else {
3054: for (r = 0; r < refine; ++r) {
3055: DM rdm;
3056: PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc;
3058: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
3059: DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);
3060: /* Total hack since we do not pass in a pointer */
3061: DMPlexReplace_Static(dm, &rdm);
3062: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
3063: if (coordFunc && remap) {
3064: DMPlexRemapGeometry(dm, 0.0, coordFunc);
3065: ((DM_Plex*) dm->data)->coordFunc = coordFunc;
3066: }
3067: }
3068: }
3069: /* Handle DMPlex coarsening */
3070: PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0);
3071: PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0);
3072: if (coarsen && isHierarchy) {
3073: DM *dms;
3075: PetscMalloc1(coarsen, &dms);
3076: DMCoarsenHierarchy(dm, coarsen, dms);
3077: /* Free DMs */
3078: for (r = 0; r < coarsen; ++r) {
3079: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);
3080: DMDestroy(&dms[r]);
3081: }
3082: PetscFree(dms);
3083: } else {
3084: for (r = 0; r < coarsen; ++r) {
3085: DM cdm;
3086: PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc;
3088: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
3089: DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm);
3090: /* Total hack since we do not pass in a pointer */
3091: DMPlexReplace_Static(dm, &cdm);
3092: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
3093: if (coordFunc) {
3094: DMPlexRemapGeometry(dm, 0.0, coordFunc);
3095: ((DM_Plex*) dm->data)->coordFunc = coordFunc;
3096: }
3097: }
3098: }
3099: /* Handle ghost cells */
3100: PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL);
3101: if (ghostCells) {
3102: DM gdm;
3103: char lname[PETSC_MAX_PATH_LEN];
3105: lname[0] = '\0';
3106: PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg);
3107: DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm);
3108: DMPlexReplace_Static(dm, &gdm);
3109: }
3110: /* Handle */
3111: DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);
3112: PetscOptionsTail();
3113: return(0);
3114: }
3116: static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
3117: {
3121: DMCreateGlobalVector_Section_Private(dm,vec);
3122: /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
3123: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);
3124: VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);
3125: VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);
3126: VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);
3127: return(0);
3128: }
3130: static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
3131: {
3135: DMCreateLocalVector_Section_Private(dm,vec);
3136: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);
3137: VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);
3138: return(0);
3139: }
3141: static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
3142: {
3143: PetscInt depth, d;
3147: DMPlexGetDepth(dm, &depth);
3148: if (depth == 1) {
3149: DMGetDimension(dm, &d);
3150: if (dim == 0) {DMPlexGetDepthStratum(dm, dim, pStart, pEnd);}
3151: else if (dim == d) {DMPlexGetDepthStratum(dm, 1, pStart, pEnd);}
3152: else {*pStart = 0; *pEnd = 0;}
3153: } else {
3154: DMPlexGetDepthStratum(dm, dim, pStart, pEnd);
3155: }
3156: return(0);
3157: }
3159: static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[])
3160: {
3161: PetscSF sf;
3162: PetscInt niranks, njranks, n;
3163: const PetscMPIInt *iranks, *jranks;
3164: DM_Plex *data = (DM_Plex*) dm->data;
3165: PetscErrorCode ierr;
3168: DMGetPointSF(dm, &sf);
3169: if (!data->neighbors) {
3170: PetscSFSetUp(sf);
3171: PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);
3172: PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);
3173: PetscMalloc1(njranks + niranks + 1, &data->neighbors);
3174: PetscArraycpy(data->neighbors + 1, jranks, njranks);
3175: PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);
3176: n = njranks + niranks;
3177: PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);
3178: /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */
3179: PetscMPIIntCast(n, data->neighbors);
3180: }
3181: if (nranks) *nranks = data->neighbors[0];
3182: if (ranks) {
3183: if (data->neighbors[0]) *ranks = data->neighbors + 1;
3184: else *ranks = NULL;
3185: }
3186: return(0);
3187: }
3189: PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec);
3191: static PetscErrorCode DMInitialize_Plex(DM dm)
3192: {
3196: dm->ops->view = DMView_Plex;
3197: dm->ops->load = DMLoad_Plex;
3198: dm->ops->setfromoptions = DMSetFromOptions_Plex;
3199: dm->ops->clone = DMClone_Plex;
3200: dm->ops->setup = DMSetUp_Plex;
3201: dm->ops->createlocalsection = DMCreateLocalSection_Plex;
3202: dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex;
3203: dm->ops->createglobalvector = DMCreateGlobalVector_Plex;
3204: dm->ops->createlocalvector = DMCreateLocalVector_Plex;
3205: dm->ops->getlocaltoglobalmapping = NULL;
3206: dm->ops->createfieldis = NULL;
3207: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex;
3208: dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex;
3209: dm->ops->getcoloring = NULL;
3210: dm->ops->creatematrix = DMCreateMatrix_Plex;
3211: dm->ops->createinterpolation = DMCreateInterpolation_Plex;
3212: dm->ops->createmassmatrix = DMCreateMassMatrix_Plex;
3213: dm->ops->createinjection = DMCreateInjection_Plex;
3214: dm->ops->refine = DMRefine_Plex;
3215: dm->ops->coarsen = DMCoarsen_Plex;
3216: dm->ops->refinehierarchy = DMRefineHierarchy_Plex;
3217: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex;
3218: dm->ops->adaptlabel = DMAdaptLabel_Plex;
3219: dm->ops->adaptmetric = DMAdaptMetric_Plex;
3220: dm->ops->globaltolocalbegin = NULL;
3221: dm->ops->globaltolocalend = NULL;
3222: dm->ops->localtoglobalbegin = NULL;
3223: dm->ops->localtoglobalend = NULL;
3224: dm->ops->destroy = DMDestroy_Plex;
3225: dm->ops->createsubdm = DMCreateSubDM_Plex;
3226: dm->ops->createsuperdm = DMCreateSuperDM_Plex;
3227: dm->ops->getdimpoints = DMGetDimPoints_Plex;
3228: dm->ops->locatepoints = DMLocatePoints_Plex;
3229: dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex;
3230: dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex;
3231: dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex;
3232: dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex;
3233: dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex;
3234: dm->ops->computel2diff = DMComputeL2Diff_Plex;
3235: dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex;
3236: dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex;
3237: dm->ops->getneighbors = DMGetNeighbors_Plex;
3238: PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);
3239: PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);
3240: PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);
3241: PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);
3242: PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);
3243: PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex);
3244: return(0);
3245: }
3247: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
3248: {
3249: DM_Plex *mesh = (DM_Plex *) dm->data;
3253: mesh->refct++;
3254: (*newdm)->data = mesh;
3255: PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);
3256: DMInitialize_Plex(*newdm);
3257: return(0);
3258: }
3260: /*MC
3261: DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram.
3262: In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is
3263: specified by a PetscSection object. Ownership in the global representation is determined by
3264: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
3266: Options Database Keys:
3267: + -dm_refine_pre - Refine mesh before distribution
3268: + -dm_refine_uniform_pre - Choose uniform or generator-based refinement
3269: + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator
3270: . -dm_distribute - Distribute mesh across processes
3271: . -dm_distribute_overlap - Number of cells to overlap for distribution
3272: . -dm_refine - Refine mesh after distribution
3273: . -dm_plex_hash_location - Use grid hashing for point location
3274: . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash
3275: . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes
3276: . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing
3277: . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally
3278: . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement
3279: . -dm_plex_check_all - Perform all shecks below
3280: . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric
3281: . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices
3282: . -dm_plex_check_faces <celltype> - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type
3283: . -dm_plex_check_geometry - Check that cells have positive volume
3284: . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ
3285: . -dm_plex_view_scale <num> - Scale the TikZ
3286: - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices
3288: Level: intermediate
3290: .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
3291: M*/
3293: PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
3294: {
3295: DM_Plex *mesh;
3296: PetscInt unit;
3301: PetscNewLog(dm,&mesh);
3302: dm->data = mesh;
3304: mesh->refct = 1;
3305: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);
3306: mesh->maxConeSize = 0;
3307: mesh->cones = NULL;
3308: mesh->coneOrientations = NULL;
3309: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);
3310: mesh->maxSupportSize = 0;
3311: mesh->supports = NULL;
3312: mesh->refinementUniform = PETSC_TRUE;
3313: mesh->refinementLimit = -1.0;
3314: mesh->interpolated = DMPLEX_INTERPOLATED_INVALID;
3315: mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID;
3317: mesh->facesTmp = NULL;
3319: mesh->tetgenOpts = NULL;
3320: mesh->triangleOpts = NULL;
3321: PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);
3322: mesh->remeshBd = PETSC_FALSE;
3324: mesh->subpointMap = NULL;
3326: for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0;
3328: mesh->regularRefinement = PETSC_FALSE;
3329: mesh->depthState = -1;
3330: mesh->celltypeState = -1;
3331: mesh->globalVertexNumbers = NULL;
3332: mesh->globalCellNumbers = NULL;
3333: mesh->anchorSection = NULL;
3334: mesh->anchorIS = NULL;
3335: mesh->createanchors = NULL;
3336: mesh->computeanchormatrix = NULL;
3337: mesh->parentSection = NULL;
3338: mesh->parents = NULL;
3339: mesh->childIDs = NULL;
3340: mesh->childSection = NULL;
3341: mesh->children = NULL;
3342: mesh->referenceTree = NULL;
3343: mesh->getchildsymmetry = NULL;
3344: mesh->vtkCellHeight = 0;
3345: mesh->useAnchors = PETSC_FALSE;
3347: mesh->maxProjectionHeight = 0;
3349: mesh->neighbors = NULL;
3351: mesh->printSetValues = PETSC_FALSE;
3352: mesh->printFEM = 0;
3353: mesh->printTol = 1.0e-10;
3355: DMInitialize_Plex(dm);
3356: return(0);
3357: }
3359: /*@
3360: DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.
3362: Collective
3364: Input Parameter:
3365: . comm - The communicator for the DMPlex object
3367: Output Parameter:
3368: . mesh - The DMPlex object
3370: Level: beginner
3372: @*/
3373: PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
3374: {
3379: DMCreate(comm, mesh);
3380: DMSetType(*mesh, DMPLEX);
3381: return(0);
3382: }
3384: /*@C
3385: DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output)
3387: Input Parameters:
3388: + dm - The DM
3389: . numCells - The number of cells owned by this process
3390: . numVertices - The number of vertices to be owned by this process, or PETSC_DECIDE
3391: . NVertices - The global number of vertices, or PETSC_DETERMINE
3392: . numCorners - The number of vertices for each cell
3393: - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell
3395: Output Parameter:
3396: . vertexSF - (Optional) SF describing complete vertex ownership
3398: Notes:
3399: Two triangles sharing a face
3400: $
3401: $ 2
3402: $ / | \
3403: $ / | \
3404: $ / | \
3405: $ 0 0 | 1 3
3406: $ \ | /
3407: $ \ | /
3408: $ \ | /
3409: $ 1
3410: would have input
3411: $ numCells = 2, numVertices = 4
3412: $ cells = [0 1 2 1 3 2]
3413: $
3414: which would result in the DMPlex
3415: $
3416: $ 4
3417: $ / | \
3418: $ / | \
3419: $ / | \
3420: $ 2 0 | 1 5
3421: $ \ | /
3422: $ \ | /
3423: $ \ | /
3424: $ 3
3426: Vertices are implicitly numbered consecutively 0,...,NVertices.
3427: Each rank owns a chunk of numVertices consecutive vertices.
3428: If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout.
3429: If NVertices is PETSC_DETERMINE and numVertices is PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1.
3430: If only NVertices is PETSC_DETERMINE, it is computed as the sum of numVertices over all ranks.
3432: The cell distribution is arbitrary non-overlapping, independent of the vertex distribution.
3434: Not currently supported in Fortran.
3436: Level: advanced
3438: .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel()
3439: @*/
3440: PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF)
3441: {
3442: PetscSF sfPoint;
3443: PetscLayout layout;
3444: PetscInt numVerticesAdj, *verticesAdj, *cones, c, p, dim;
3445: PetscMPIInt rank, size;
3446: PetscErrorCode ierr;
3450: PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);
3451: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3452: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
3453: DMGetDimension(dm, &dim);
3454: /* Get/check global number of vertices */
3455: {
3456: PetscInt NVerticesInCells, i;
3457: const PetscInt len = numCells * numCorners;
3459: /* NVerticesInCells = max(cells) + 1 */
3460: NVerticesInCells = PETSC_MIN_INT;
3461: for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i];
3462: ++NVerticesInCells;
3463: MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
3465: if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells;
3466: else if (NVertices != PETSC_DECIDE && NVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified global number of vertices %D must be greater than or equal to the number of vertices in cells %D",NVertices,NVerticesInCells);
3467: }
3468: /* Count locally unique vertices */
3469: {
3470: PetscHSetI vhash;
3471: PetscInt off = 0;
3473: PetscHSetICreate(&vhash);
3474: for (c = 0; c < numCells; ++c) {
3475: for (p = 0; p < numCorners; ++p) {
3476: PetscHSetIAdd(vhash, cells[c*numCorners+p]);
3477: }
3478: }
3479: PetscHSetIGetSize(vhash, &numVerticesAdj);
3480: PetscMalloc1(numVerticesAdj, &verticesAdj);
3481: PetscHSetIGetElems(vhash, &off, verticesAdj);
3482: PetscHSetIDestroy(&vhash);
3483: if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj);
3484: }
3485: PetscSortInt(numVerticesAdj, verticesAdj);
3486: /* Create cones */
3487: DMPlexSetChart(dm, 0, numCells+numVerticesAdj);
3488: for (c = 0; c < numCells; ++c) {DMPlexSetConeSize(dm, c, numCorners);}
3489: DMSetUp(dm);
3490: DMPlexGetCones(dm,&cones);
3491: for (c = 0; c < numCells; ++c) {
3492: for (p = 0; p < numCorners; ++p) {
3493: const PetscInt gv = cells[c*numCorners+p];
3494: PetscInt lv;
3496: /* Positions within verticesAdj form 0-based local vertex numbering;
3497: we need to shift it by numCells to get correct DAG points (cells go first) */
3498: PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);
3499: if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv);
3500: cones[c*numCorners+p] = lv+numCells;
3501: }
3502: }
3503: /* Build point sf */
3504: PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout);
3505: PetscLayoutSetSize(layout, NVertices);
3506: PetscLayoutSetLocalSize(layout, numVertices);
3507: PetscLayoutSetBlockSize(layout, 1);
3508: PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint);
3509: PetscLayoutDestroy(&layout);
3510: PetscFree(verticesAdj);
3511: PetscObjectSetName((PetscObject) sfPoint, "point SF");
3512: if (dm->sf) {
3513: const char *prefix;
3515: PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix);
3516: PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix);
3517: }
3518: DMSetPointSF(dm, sfPoint);
3519: PetscSFDestroy(&sfPoint);
3520: if (vertexSF) {PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF");}
3521: /* Fill in the rest of the topology structure */
3522: DMPlexSymmetrize(dm);
3523: DMPlexStratify(dm);
3524: PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);
3525: return(0);
3526: }
3528: /*@C
3529: DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output)
3531: Input Parameters:
3532: + dm - The DM
3533: . spaceDim - The spatial dimension used for coordinates
3534: . sfVert - SF describing complete vertex ownership
3535: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
3537: Level: advanced
3539: Notes:
3540: Not currently supported in Fortran.
3542: .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel()
3543: @*/
3544: PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[])
3545: {
3546: PetscSection coordSection;
3547: Vec coordinates;
3548: PetscScalar *coords;
3549: PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd;
3553: PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);
3554: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3555: if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first.");
3556: DMSetCoordinateDim(dm, spaceDim);
3557: PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);
3558: if (vEnd - vStart != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Supplied sfVert has wrong number of leaves = %D != %D = vEnd - vStart",numVerticesAdj,vEnd - vStart);
3559: DMGetCoordinateSection(dm, &coordSection);
3560: PetscSectionSetNumFields(coordSection, 1);
3561: PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
3562: PetscSectionSetChart(coordSection, vStart, vEnd);
3563: for (v = vStart; v < vEnd; ++v) {
3564: PetscSectionSetDof(coordSection, v, spaceDim);
3565: PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
3566: }
3567: PetscSectionSetUp(coordSection);
3568: PetscSectionGetStorageSize(coordSection, &coordSize);
3569: VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
3570: VecSetBlockSize(coordinates, spaceDim);
3571: PetscObjectSetName((PetscObject) coordinates, "coordinates");
3572: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
3573: VecSetType(coordinates,VECSTANDARD);
3574: VecGetArray(coordinates, &coords);
3575: {
3576: MPI_Datatype coordtype;
3578: /* Need a temp buffer for coords if we have complex/single */
3579: MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);
3580: MPI_Type_commit(&coordtype);
3581: #if defined(PETSC_USE_COMPLEX)
3582: {
3583: PetscScalar *svertexCoords;
3584: PetscInt i;
3585: PetscMalloc1(numVertices*spaceDim,&svertexCoords);
3586: for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i];
3587: PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);
3588: PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);
3589: PetscFree(svertexCoords);
3590: }
3591: #else
3592: PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);
3593: PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);
3594: #endif
3595: MPI_Type_free(&coordtype);
3596: }
3597: VecRestoreArray(coordinates, &coords);
3598: DMSetCoordinatesLocal(dm, coordinates);
3599: VecDestroy(&coordinates);
3600: PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);
3601: return(0);
3602: }
3604: /*@
3605: DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output)
3607: Input Parameters:
3608: + comm - The communicator
3609: . dim - The topological dimension of the mesh
3610: . numCells - The number of cells owned by this process
3611: . numVertices - The number of vertices owned by this process, or PETSC_DECIDE
3612: . NVertices - The global number of vertices, or PETSC_DECIDE
3613: . numCorners - The number of vertices for each cell
3614: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
3615: . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell
3616: . spaceDim - The spatial dimension used for coordinates
3617: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
3619: Output Parameters:
3620: + dm - The DM
3621: - vertexSF - (Optional) SF describing complete vertex ownership
3623: Notes:
3624: This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(),
3625: DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel()
3627: See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters.
3628: See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters.
3630: Level: intermediate
3632: .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate()
3633: @*/
3634: PetscErrorCode DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm)
3635: {
3636: PetscSF sfVert;
3640: DMCreate(comm, dm);
3641: DMSetType(*dm, DMPLEX);
3644: DMSetDimension(*dm, dim);
3645: DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert);
3646: if (interpolate) {
3647: DM idm;
3649: DMPlexInterpolate(*dm, &idm);
3650: DMDestroy(dm);
3651: *dm = idm;
3652: }
3653: DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);
3654: if (vertexSF) *vertexSF = sfVert;
3655: else {PetscSFDestroy(&sfVert);}
3656: return(0);
3657: }
3659: /*@
3660: DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc()
3662: Level: deprecated
3664: .seealso: DMPlexCreateFromCellListParallelPetsc()
3665: @*/
3666: PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm)
3667: {
3669: PetscInt i;
3670: PetscInt *pintCells;
3673: if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt));
3674: if (sizeof(int) == sizeof(PetscInt)) {
3675: pintCells = (PetscInt *) cells;
3676: } else {
3677: PetscMalloc1(numCells*numCorners, &pintCells);
3678: for (i = 0; i < numCells*numCorners; i++) {
3679: pintCells[i] = (PetscInt) cells[i];
3680: }
3681: }
3682: DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, dm);
3683: if (sizeof(int) != sizeof(PetscInt)) {
3684: PetscFree(pintCells);
3685: }
3686: return(0);
3687: }
3689: /*@C
3690: DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output)
3692: Input Parameters:
3693: + dm - The DM
3694: . numCells - The number of cells owned by this process
3695: . numVertices - The number of vertices owned by this process, or PETSC_DETERMINE
3696: . numCorners - The number of vertices for each cell
3697: - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell
3699: Level: advanced
3701: Notes:
3702: Two triangles sharing a face
3703: $
3704: $ 2
3705: $ / | \
3706: $ / | \
3707: $ / | \
3708: $ 0 0 | 1 3
3709: $ \ | /
3710: $ \ | /
3711: $ \ | /
3712: $ 1
3713: would have input
3714: $ numCells = 2, numVertices = 4
3715: $ cells = [0 1 2 1 3 2]
3716: $
3717: which would result in the DMPlex
3718: $
3719: $ 4
3720: $ / | \
3721: $ / | \
3722: $ / | \
3723: $ 2 0 | 1 5
3724: $ \ | /
3725: $ \ | /
3726: $ \ | /
3727: $ 3
3729: If numVertices is PETSC_DETERMINE, it is computed by PETSc as the maximum vertex index in cells + 1.
3731: Not currently supported in Fortran.
3733: .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc()
3734: @*/
3735: PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[])
3736: {
3737: PetscInt *cones, c, p, dim;
3741: PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);
3742: DMGetDimension(dm, &dim);
3743: /* Get/check global number of vertices */
3744: {
3745: PetscInt NVerticesInCells, i;
3746: const PetscInt len = numCells * numCorners;
3748: /* NVerticesInCells = max(cells) + 1 */
3749: NVerticesInCells = PETSC_MIN_INT;
3750: for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i];
3751: ++NVerticesInCells;
3753: if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells;
3754: else if (numVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified number of vertices %D must be greater than or equal to the number of vertices in cells %D",numVertices,NVerticesInCells);
3755: }
3756: DMPlexSetChart(dm, 0, numCells+numVertices);
3757: for (c = 0; c < numCells; ++c) {
3758: DMPlexSetConeSize(dm, c, numCorners);
3759: }
3760: DMSetUp(dm);
3761: DMPlexGetCones(dm,&cones);
3762: for (c = 0; c < numCells; ++c) {
3763: for (p = 0; p < numCorners; ++p) {
3764: cones[c*numCorners+p] = cells[c*numCorners+p]+numCells;
3765: }
3766: }
3767: DMPlexSymmetrize(dm);
3768: DMPlexStratify(dm);
3769: PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);
3770: return(0);
3771: }
3773: /*@C
3774: DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output)
3776: Input Parameters:
3777: + dm - The DM
3778: . spaceDim - The spatial dimension used for coordinates
3779: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
3781: Level: advanced
3783: Notes:
3784: Not currently supported in Fortran.
3786: .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList()
3787: @*/
3788: PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[])
3789: {
3790: PetscSection coordSection;
3791: Vec coordinates;
3792: DM cdm;
3793: PetscScalar *coords;
3794: PetscInt v, vStart, vEnd, d;
3798: PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);
3799: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3800: if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first.");
3801: DMSetCoordinateDim(dm, spaceDim);
3802: DMGetCoordinateSection(dm, &coordSection);
3803: PetscSectionSetNumFields(coordSection, 1);
3804: PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
3805: PetscSectionSetChart(coordSection, vStart, vEnd);
3806: for (v = vStart; v < vEnd; ++v) {
3807: PetscSectionSetDof(coordSection, v, spaceDim);
3808: PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
3809: }
3810: PetscSectionSetUp(coordSection);
3812: DMGetCoordinateDM(dm, &cdm);
3813: DMCreateLocalVector(cdm, &coordinates);
3814: VecSetBlockSize(coordinates, spaceDim);
3815: PetscObjectSetName((PetscObject) coordinates, "coordinates");
3816: VecGetArrayWrite(coordinates, &coords);
3817: for (v = 0; v < vEnd-vStart; ++v) {
3818: for (d = 0; d < spaceDim; ++d) {
3819: coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
3820: }
3821: }
3822: VecRestoreArrayWrite(coordinates, &coords);
3823: DMSetCoordinatesLocal(dm, coordinates);
3824: VecDestroy(&coordinates);
3825: PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);
3826: return(0);
3827: }
3829: /*@
3830: DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input
3832: Collective on comm
3834: Input Parameters:
3835: + comm - The communicator
3836: . dim - The topological dimension of the mesh
3837: . numCells - The number of cells, only on process 0
3838: . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0
3839: . numCorners - The number of vertices for each cell, only on process 0
3840: . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
3841: . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0
3842: . spaceDim - The spatial dimension used for coordinates
3843: - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0
3845: Output Parameter:
3846: . dm - The DM, which only has points on process 0
3848: Notes:
3849: This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(),
3850: DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList()
3852: See DMPlexBuildFromCellList() for an example and details about the topology-related parameters.
3853: See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters.
3854: See DMPlexCreateFromCellListParallelPetsc() for parallel input
3856: Level: intermediate
3858: .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate()
3859: @*/
3860: PetscErrorCode DMPlexCreateFromCellListPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], DM *dm)
3861: {
3862: PetscMPIInt rank;
3866: if (!dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "This is not appropriate for 0-dimensional meshes. Consider either creating the DM using DMPlexCreateFromDAG(), by hand, or using DMSwarm.");
3867: MPI_Comm_rank(comm, &rank);
3868: DMCreate(comm, dm);
3869: DMSetType(*dm, DMPLEX);
3870: DMSetDimension(*dm, dim);
3871: if (!rank) {DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);}
3872: else {DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL);}
3873: if (interpolate) {
3874: DM idm;
3876: DMPlexInterpolate(*dm, &idm);
3877: DMDestroy(dm);
3878: *dm = idm;
3879: }
3880: if (!rank) {DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);}
3881: else {DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL);}
3882: return(0);
3883: }
3885: /*@
3886: DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc()
3888: Level: deprecated
3890: .seealso: DMPlexCreateFromCellListPetsc()
3891: @*/
3892: 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)
3893: {
3895: PetscInt i;
3896: PetscInt *pintCells;
3897: PetscReal *prealVC;
3900: if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt));
3901: if (sizeof(int) == sizeof(PetscInt)) {
3902: pintCells = (PetscInt *) cells;
3903: } else {
3904: PetscMalloc1(numCells*numCorners, &pintCells);
3905: for (i = 0; i < numCells*numCorners; i++) {
3906: pintCells[i] = (PetscInt) cells[i];
3907: }
3908: }
3909: if (sizeof(double) > sizeof(PetscReal)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of double %zd greater than size of PetscReal %zd. Reconfigure PETSc --with-precision=<higher precision>.", sizeof(double), sizeof(PetscReal));
3910: if (sizeof(double) == sizeof(PetscReal)) {
3911: prealVC = (PetscReal *) vertexCoords;
3912: } else {
3913: PetscMalloc1(numVertices*spaceDim, &prealVC);
3914: for (i = 0; i < numVertices*spaceDim; i++) {
3915: prealVC[i] = (PetscReal) vertexCoords[i];
3916: }
3917: }
3918: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);
3919: if (sizeof(int) != sizeof(PetscInt)) {
3920: PetscFree(pintCells);
3921: }
3922: if (sizeof(double) != sizeof(PetscReal)) {
3923: PetscFree(prealVC);
3924: }
3925: return(0);
3926: }
3928: /*@
3929: DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM
3931: Input Parameters:
3932: + dm - The empty DM object, usually from DMCreate() and DMSetDimension()
3933: . depth - The depth of the DAG
3934: . numPoints - Array of size depth + 1 containing the number of points at each depth
3935: . coneSize - The cone size of each point
3936: . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
3937: . coneOrientations - The orientation of each cone point
3938: - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim()
3940: Output Parameter:
3941: . dm - The DM
3943: Note: Two triangles sharing a face would have input
3944: $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
3945: $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0]
3946: $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0]
3947: $
3948: which would result in the DMPlex
3949: $
3950: $ 4
3951: $ / | \
3952: $ / | \
3953: $ / | \
3954: $ 2 0 | 1 5
3955: $ \ | /
3956: $ \ | /
3957: $ \ | /
3958: $ 3
3959: $
3960: $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc()
3962: Level: advanced
3964: .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate()
3965: @*/
3966: PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
3967: {
3968: Vec coordinates;
3969: PetscSection coordSection;
3970: PetscScalar *coords;
3971: PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off;
3975: DMGetDimension(dm, &dim);
3976: DMGetCoordinateDim(dm, &dimEmbed);
3977: if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim);
3978: for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
3979: DMPlexSetChart(dm, pStart, pEnd);
3980: for (p = pStart; p < pEnd; ++p) {
3981: DMPlexSetConeSize(dm, p, coneSize[p-pStart]);
3982: if (firstVertex < 0 && !coneSize[p - pStart]) {
3983: firstVertex = p - pStart;
3984: }
3985: }
3986: if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]);
3987: DMSetUp(dm); /* Allocate space for cones */
3988: for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
3989: DMPlexSetCone(dm, p, &cones[off]);
3990: DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);
3991: }
3992: DMPlexSymmetrize(dm);
3993: DMPlexStratify(dm);
3994: /* Build coordinates */
3995: DMGetCoordinateSection(dm, &coordSection);
3996: PetscSectionSetNumFields(coordSection, 1);
3997: PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);
3998: PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);
3999: for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
4000: PetscSectionSetDof(coordSection, v, dimEmbed);
4001: PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);
4002: }
4003: PetscSectionSetUp(coordSection);
4004: PetscSectionGetStorageSize(coordSection, &coordSize);
4005: VecCreate(PETSC_COMM_SELF, &coordinates);
4006: PetscObjectSetName((PetscObject) coordinates, "coordinates");
4007: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
4008: VecSetBlockSize(coordinates, dimEmbed);
4009: VecSetType(coordinates,VECSTANDARD);
4010: if (vertexCoords) {
4011: VecGetArray(coordinates, &coords);
4012: for (v = 0; v < numPoints[0]; ++v) {
4013: PetscInt off;
4015: PetscSectionGetOffset(coordSection, v+firstVertex, &off);
4016: for (d = 0; d < dimEmbed; ++d) {
4017: coords[off+d] = vertexCoords[v*dimEmbed+d];
4018: }
4019: }
4020: }
4021: VecRestoreArray(coordinates, &coords);
4022: DMSetCoordinatesLocal(dm, coordinates);
4023: VecDestroy(&coordinates);
4024: return(0);
4025: }
4027: /*@C
4028: DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file.
4030: + comm - The MPI communicator
4031: . filename - Name of the .dat file
4032: - interpolate - Create faces and edges in the mesh
4034: Output Parameter:
4035: . dm - The DM object representing the mesh
4037: Note: The format is the simplest possible:
4038: $ Ne
4039: $ v0 v1 ... vk
4040: $ Nv
4041: $ x y z marker
4043: Level: beginner
4045: .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
4046: @*/
4047: PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
4048: {
4049: DMLabel marker;
4050: PetscViewer viewer;
4051: Vec coordinates;
4052: PetscSection coordSection;
4053: PetscScalar *coords;
4054: char line[PETSC_MAX_PATH_LEN];
4055: PetscInt dim = 3, cdim = 3, coordSize, v, c, d;
4056: PetscMPIInt rank;
4057: int snum, Nv, Nc;
4058: PetscErrorCode ierr;
4061: MPI_Comm_rank(comm, &rank);
4062: PetscViewerCreate(comm, &viewer);
4063: PetscViewerSetType(viewer, PETSCVIEWERASCII);
4064: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
4065: PetscViewerFileSetName(viewer, filename);
4066: if (rank == 0) {
4067: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
4068: snum = sscanf(line, "%d %d", &Nc, &Nv);
4069: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
4070: } else {
4071: Nc = Nv = 0;
4072: }
4073: DMCreate(comm, dm);
4074: DMSetType(*dm, DMPLEX);
4075: DMPlexSetChart(*dm, 0, Nc+Nv);
4076: DMSetDimension(*dm, dim);
4077: DMSetCoordinateDim(*dm, cdim);
4078: /* Read topology */
4079: if (rank == 0) {
4080: PetscInt cone[8], corners = 8;
4081: int vbuf[8], v;
4083: for (c = 0; c < Nc; ++c) {DMPlexSetConeSize(*dm, c, corners);}
4084: DMSetUp(*dm);
4085: for (c = 0; c < Nc; ++c) {
4086: PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);
4087: snum = sscanf(line, "%d %d %d %d %d %d %d %d", &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);
4088: if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
4089: for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc;
4090: /* Hexahedra are inverted */
4091: {
4092: PetscInt tmp = cone[1];
4093: cone[1] = cone[3];
4094: cone[3] = tmp;
4095: }
4096: DMPlexSetCone(*dm, c, cone);
4097: }
4098: }
4099: DMPlexSymmetrize(*dm);
4100: DMPlexStratify(*dm);
4101: /* Read coordinates */
4102: DMGetCoordinateSection(*dm, &coordSection);
4103: PetscSectionSetNumFields(coordSection, 1);
4104: PetscSectionSetFieldComponents(coordSection, 0, cdim);
4105: PetscSectionSetChart(coordSection, Nc, Nc + Nv);
4106: for (v = Nc; v < Nc+Nv; ++v) {
4107: PetscSectionSetDof(coordSection, v, cdim);
4108: PetscSectionSetFieldDof(coordSection, v, 0, cdim);
4109: }
4110: PetscSectionSetUp(coordSection);
4111: PetscSectionGetStorageSize(coordSection, &coordSize);
4112: VecCreate(PETSC_COMM_SELF, &coordinates);
4113: PetscObjectSetName((PetscObject) coordinates, "coordinates");
4114: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
4115: VecSetBlockSize(coordinates, cdim);
4116: VecSetType(coordinates, VECSTANDARD);
4117: VecGetArray(coordinates, &coords);
4118: if (rank == 0) {
4119: double x[3];
4120: int val;
4122: DMCreateLabel(*dm, "marker");
4123: DMGetLabel(*dm, "marker", &marker);
4124: for (v = 0; v < Nv; ++v) {
4125: PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);
4126: snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val);
4127: if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
4128: for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d];
4129: if (val) {DMLabelSetValue(marker, v+Nc, val);}
4130: }
4131: }
4132: VecRestoreArray(coordinates, &coords);
4133: DMSetCoordinatesLocal(*dm, coordinates);
4134: VecDestroy(&coordinates);
4135: PetscViewerDestroy(&viewer);
4136: if (interpolate) {
4137: DM idm;
4138: DMLabel bdlabel;
4140: DMPlexInterpolate(*dm, &idm);
4141: DMDestroy(dm);
4142: *dm = idm;
4144: DMGetLabel(*dm, "marker", &bdlabel);
4145: DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);
4146: DMPlexLabelComplete(*dm, bdlabel);
4147: }
4148: return(0);
4149: }
4151: /*@C
4152: DMPlexCreateFromFile - This takes a filename and produces a DM
4154: Input Parameters:
4155: + comm - The communicator
4156: . filename - A file name
4157: - interpolate - Flag to create intermediate mesh pieces (edges, faces)
4159: Output Parameter:
4160: . dm - The DM
4162: Options Database Keys:
4163: . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5
4165: Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g.
4166: $ -dm_plex_create_viewer_hdf5_collective
4168: Level: beginner
4170: .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate()
4171: @*/
4172: PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
4173: {
4174: const char *extGmsh = ".msh";
4175: const char *extGmsh2 = ".msh2";
4176: const char *extGmsh4 = ".msh4";
4177: const char *extCGNS = ".cgns";
4178: const char *extExodus = ".exo";
4179: const char *extGenesis = ".gen";
4180: const char *extFluent = ".cas";
4181: const char *extHDF5 = ".h5";
4182: const char *extMed = ".med";
4183: const char *extPLY = ".ply";
4184: const char *extEGADSLite = ".egadslite";
4185: const char *extEGADS = ".egads";
4186: const char *extIGES = ".igs";
4187: const char *extSTEP = ".stp";
4188: const char *extCV = ".dat";
4189: size_t len;
4190: PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV;
4191: PetscMPIInt rank;
4197: DMInitializePackage();
4198: PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);
4199: MPI_Comm_rank(comm, &rank);
4200: PetscStrlen(filename, &len);
4201: if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path");
4202: PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);
4203: PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);
4204: PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);
4205: PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);
4206: PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);
4207: PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);
4208: PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);
4209: PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);
4210: PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);
4211: PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);
4212: PetscStrncmp(&filename[PetscMax(0,len-10)], extEGADSLite, 10, &isEGADSLite);
4213: PetscStrncmp(&filename[PetscMax(0,len-6)], extEGADS, 6, &isEGADS);
4214: PetscStrncmp(&filename[PetscMax(0,len-4)], extIGES, 4, &isIGES);
4215: PetscStrncmp(&filename[PetscMax(0,len-4)], extSTEP, 4, &isSTEP);
4216: PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);
4217: if (isGmsh || isGmsh2 || isGmsh4) {
4218: DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);
4219: } else if (isCGNS) {
4220: DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);
4221: } else if (isExodus || isGenesis) {
4222: DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);
4223: } else if (isFluent) {
4224: DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);
4225: } else if (isHDF5) {
4226: PetscBool load_hdf5_xdmf = PETSC_FALSE;
4227: PetscViewer viewer;
4229: /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */
4230: PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);
4231: PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);
4232: PetscViewerCreate(comm, &viewer);
4233: PetscViewerSetType(viewer, PETSCVIEWERHDF5);
4234: PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");
4235: PetscViewerSetFromOptions(viewer);
4236: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
4237: PetscViewerFileSetName(viewer, filename);
4238: DMCreate(comm, dm);
4239: DMSetType(*dm, DMPLEX);
4240: if (load_hdf5_xdmf) {PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);}
4241: DMLoad(*dm, viewer);
4242: if (load_hdf5_xdmf) {PetscViewerPopFormat(viewer);}
4243: PetscViewerDestroy(&viewer);
4245: if (interpolate) {
4246: DM idm;
4248: DMPlexInterpolate(*dm, &idm);
4249: DMDestroy(dm);
4250: *dm = idm;
4251: }
4252: } else if (isMed) {
4253: DMPlexCreateMedFromFile(comm, filename, interpolate, dm);
4254: } else if (isPLY) {
4255: DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);
4256: } else if (isEGADSLite || isEGADS || isIGES || isSTEP) {
4257: if (isEGADSLite) {DMPlexCreateEGADSLiteFromFile(comm, filename, dm);}
4258: else {DMPlexCreateEGADSFromFile(comm, filename, dm);}
4259: if (!interpolate) {
4260: DM udm;
4262: DMPlexUninterpolate(*dm, &udm);
4263: DMDestroy(dm);
4264: *dm = udm;
4265: }
4266: } else if (isCV) {
4267: DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);
4268: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename);
4269: PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);
4270: return(0);
4271: }