Actual source code: plexsubmesh.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petscsf.h>
5: static PetscErrorCode DMPlexCellIsHybrid_Internal(DM dm, PetscInt p, PetscBool *isHybrid)
6: {
7: DMPolytopeType ct;
9: DMPlexGetCellType(dm, p, &ct);
10: switch (ct) {
11: case DM_POLYTOPE_POINT_PRISM_TENSOR:
12: case DM_POLYTOPE_SEG_PRISM_TENSOR:
13: case DM_POLYTOPE_TRI_PRISM_TENSOR:
14: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
15: *isHybrid = PETSC_TRUE;
16: default: *isHybrid = PETSC_FALSE;
17: }
18: return 0;
19: }
21: static PetscErrorCode DMPlexGetTensorPrismBounds_Internal(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
22: {
23: DMLabel ctLabel;
25: if (cStart) *cStart = -1;
26: if (cEnd) *cEnd = -1;
27: DMPlexGetCellTypeLabel(dm, &ctLabel);
28: switch (dim) {
29: case 1: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_POINT_PRISM_TENSOR, cStart, cEnd);break;
30: case 2: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_SEG_PRISM_TENSOR, cStart, cEnd);break;
31: case 3:
32: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_TRI_PRISM_TENSOR, cStart, cEnd);
33: if (*cStart < 0) DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_QUAD_PRISM_TENSOR, cStart, cEnd);
34: break;
35: default: return 0;
36: }
37: return 0;
38: }
40: static PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt val, PetscInt cellHeight, DMLabel label)
41: {
42: PetscInt fStart, fEnd, f;
44: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
45: for (f = fStart; f < fEnd; ++f) {
46: PetscInt supportSize;
48: DMPlexGetSupportSize(dm, f, &supportSize);
49: if (supportSize == 1) {
50: if (val < 0) {
51: PetscInt *closure = NULL;
52: PetscInt clSize, cl, cval;
54: DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
55: for (cl = 0; cl < clSize*2; cl += 2) {
56: DMLabelGetValue(label, closure[cl], &cval);
57: if (cval < 0) continue;
58: DMLabelSetValue(label, f, cval);
59: break;
60: }
61: if (cl == clSize*2) DMLabelSetValue(label, f, 1);
62: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
63: } else {
64: DMLabelSetValue(label, f, val);
65: }
66: }
67: }
68: return 0;
69: }
71: /*@
72: DMPlexMarkBoundaryFaces - Mark all faces on the boundary
74: Not Collective
76: Input Parameters:
77: + dm - The original DM
78: - val - The marker value, or PETSC_DETERMINE to use some value in the closure (or 1 if none are found)
80: Output Parameter:
81: . label - The DMLabel marking boundary faces with the given value
83: Level: developer
85: .seealso: DMLabelCreate(), DMCreateLabel()
86: @*/
87: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, PetscInt val, DMLabel label)
88: {
89: DMPlexInterpolatedFlag flg;
92: DMPlexIsInterpolated(dm, &flg);
94: DMPlexMarkBoundaryFaces_Internal(dm, val, 0, label);
95: return 0;
96: }
98: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
99: {
100: IS valueIS;
101: PetscSF sfPoint;
102: const PetscInt *values;
103: PetscInt numValues, v, cStart, cEnd, nroots;
105: DMLabelGetNumValues(label, &numValues);
106: DMLabelGetValueIS(label, &valueIS);
107: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
108: ISGetIndices(valueIS, &values);
109: for (v = 0; v < numValues; ++v) {
110: IS pointIS;
111: const PetscInt *points;
112: PetscInt numPoints, p;
114: DMLabelGetStratumSize(label, values[v], &numPoints);
115: DMLabelGetStratumIS(label, values[v], &pointIS);
116: ISGetIndices(pointIS, &points);
117: for (p = 0; p < numPoints; ++p) {
118: PetscInt q = points[p];
119: PetscInt *closure = NULL;
120: PetscInt closureSize, c;
122: if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
123: continue;
124: }
125: DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
126: for (c = 0; c < closureSize*2; c += 2) {
127: DMLabelSetValue(label, closure[c], values[v]);
128: }
129: DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
130: }
131: ISRestoreIndices(pointIS, &points);
132: ISDestroy(&pointIS);
133: }
134: ISRestoreIndices(valueIS, &values);
135: ISDestroy(&valueIS);
136: DMGetPointSF(dm, &sfPoint);
137: PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
138: if (nroots >= 0) {
139: DMLabel lblRoots, lblLeaves;
140: IS valueIS, pointIS;
141: const PetscInt *values;
142: PetscInt numValues, v;
144: /* Pull point contributions from remote leaves into local roots */
145: DMLabelGather(label, sfPoint, &lblLeaves);
146: DMLabelGetValueIS(lblLeaves, &valueIS);
147: ISGetLocalSize(valueIS, &numValues);
148: ISGetIndices(valueIS, &values);
149: for (v = 0; v < numValues; ++v) {
150: const PetscInt value = values[v];
152: DMLabelGetStratumIS(lblLeaves, value, &pointIS);
153: DMLabelInsertIS(label, pointIS, value);
154: ISDestroy(&pointIS);
155: }
156: ISRestoreIndices(valueIS, &values);
157: ISDestroy(&valueIS);
158: DMLabelDestroy(&lblLeaves);
159: /* Push point contributions from roots into remote leaves */
160: DMLabelDistribute(label, sfPoint, &lblRoots);
161: DMLabelGetValueIS(lblRoots, &valueIS);
162: ISGetLocalSize(valueIS, &numValues);
163: ISGetIndices(valueIS, &values);
164: for (v = 0; v < numValues; ++v) {
165: const PetscInt value = values[v];
167: DMLabelGetStratumIS(lblRoots, value, &pointIS);
168: DMLabelInsertIS(label, pointIS, value);
169: ISDestroy(&pointIS);
170: }
171: ISRestoreIndices(valueIS, &values);
172: ISDestroy(&valueIS);
173: DMLabelDestroy(&lblRoots);
174: }
175: return 0;
176: }
178: /*@
179: DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface
181: Input Parameters:
182: + dm - The DM
183: - label - A DMLabel marking the surface points
185: Output Parameter:
186: . label - A DMLabel marking all surface points in the transitive closure
188: Level: developer
190: .seealso: DMPlexLabelCohesiveComplete()
191: @*/
192: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
193: {
194: DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
195: return 0;
196: }
198: /*@
199: DMPlexLabelAddCells - Starting with a label marking points on a surface, we add a cell for each point
201: Input Parameters:
202: + dm - The DM
203: - label - A DMLabel marking the surface points
205: Output Parameter:
206: . label - A DMLabel incorporating cells
208: Level: developer
210: Note: The cells allow FEM boundary conditions to be applied using the cell geometry
212: .seealso: DMPlexLabelAddFaceCells(), DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
213: @*/
214: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
215: {
216: IS valueIS;
217: const PetscInt *values;
218: PetscInt numValues, v, cStart, cEnd;
220: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
221: DMLabelGetNumValues(label, &numValues);
222: DMLabelGetValueIS(label, &valueIS);
223: ISGetIndices(valueIS, &values);
224: for (v = 0; v < numValues; ++v) {
225: IS pointIS;
226: const PetscInt *points;
227: PetscInt numPoints, p;
229: DMLabelGetStratumSize(label, values[v], &numPoints);
230: DMLabelGetStratumIS(label, values[v], &pointIS);
231: ISGetIndices(pointIS, &points);
232: for (p = 0; p < numPoints; ++p) {
233: PetscInt *closure = NULL;
234: PetscInt closureSize, cl;
236: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
237: for (cl = closureSize-1; cl > 0; --cl) {
238: const PetscInt cell = closure[cl*2];
239: if ((cell >= cStart) && (cell < cEnd)) {DMLabelSetValue(label, cell, values[v]); break;}
240: }
241: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
242: }
243: ISRestoreIndices(pointIS, &points);
244: ISDestroy(&pointIS);
245: }
246: ISRestoreIndices(valueIS, &values);
247: ISDestroy(&valueIS);
248: return 0;
249: }
251: /*@
252: DMPlexLabelAddFaceCells - Starting with a label marking faces on a surface, we add a cell for each face
254: Input Parameters:
255: + dm - The DM
256: - label - A DMLabel marking the surface points
258: Output Parameter:
259: . label - A DMLabel incorporating cells
261: Level: developer
263: Note: The cells allow FEM boundary conditions to be applied using the cell geometry
265: .seealso: DMPlexLabelAddCells(), DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
266: @*/
267: PetscErrorCode DMPlexLabelAddFaceCells(DM dm, DMLabel label)
268: {
269: IS valueIS;
270: const PetscInt *values;
271: PetscInt numValues, v, cStart, cEnd, fStart, fEnd;
273: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
274: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
275: DMLabelGetNumValues(label, &numValues);
276: DMLabelGetValueIS(label, &valueIS);
277: ISGetIndices(valueIS, &values);
278: for (v = 0; v < numValues; ++v) {
279: IS pointIS;
280: const PetscInt *points;
281: PetscInt numPoints, p;
283: DMLabelGetStratumSize(label, values[v], &numPoints);
284: DMLabelGetStratumIS(label, values[v], &pointIS);
285: ISGetIndices(pointIS, &points);
286: for (p = 0; p < numPoints; ++p) {
287: const PetscInt face = points[p];
288: PetscInt *closure = NULL;
289: PetscInt closureSize, cl;
291: if ((face < fStart) || (face >= fEnd)) continue;
292: DMPlexGetTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
293: for (cl = closureSize-1; cl > 0; --cl) {
294: const PetscInt cell = closure[cl*2];
295: if ((cell >= cStart) && (cell < cEnd)) {DMLabelSetValue(label, cell, values[v]); break;}
296: }
297: DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
298: }
299: ISRestoreIndices(pointIS, &points);
300: ISDestroy(&pointIS);
301: }
302: ISRestoreIndices(valueIS, &values);
303: ISDestroy(&valueIS);
304: return 0;
305: }
307: /*@
308: DMPlexLabelClearCells - Remove cells from a label
310: Input Parameters:
311: + dm - The DM
312: - label - A DMLabel marking surface points and their adjacent cells
314: Output Parameter:
315: . label - A DMLabel without cells
317: Level: developer
319: Note: This undoes DMPlexLabelAddCells() or DMPlexLabelAddFaceCells()
321: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
322: @*/
323: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
324: {
325: IS valueIS;
326: const PetscInt *values;
327: PetscInt numValues, v, cStart, cEnd;
329: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
330: DMLabelGetNumValues(label, &numValues);
331: DMLabelGetValueIS(label, &valueIS);
332: ISGetIndices(valueIS, &values);
333: for (v = 0; v < numValues; ++v) {
334: IS pointIS;
335: const PetscInt *points;
336: PetscInt numPoints, p;
338: DMLabelGetStratumSize(label, values[v], &numPoints);
339: DMLabelGetStratumIS(label, values[v], &pointIS);
340: ISGetIndices(pointIS, &points);
341: for (p = 0; p < numPoints; ++p) {
342: PetscInt point = points[p];
344: if (point >= cStart && point < cEnd) {
345: DMLabelClearValue(label,point,values[v]);
346: }
347: }
348: ISRestoreIndices(pointIS, &points);
349: ISDestroy(&pointIS);
350: }
351: ISRestoreIndices(valueIS, &values);
352: ISDestroy(&valueIS);
353: return 0;
354: }
356: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
357: * index (skipping first, which is (0,0)) */
358: static inline PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
359: {
360: PetscInt d, off = 0;
362: /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
363: for (d = 0; d < depth; d++) {
364: PetscInt firstd = d;
365: PetscInt firstStart = depthShift[2*d];
366: PetscInt e;
368: for (e = d+1; e <= depth; e++) {
369: if (depthShift[2*e] < firstStart) {
370: firstd = e;
371: firstStart = depthShift[2*d];
372: }
373: }
374: if (firstd != d) {
375: PetscInt swap[2];
377: e = firstd;
378: swap[0] = depthShift[2*d];
379: swap[1] = depthShift[2*d+1];
380: depthShift[2*d] = depthShift[2*e];
381: depthShift[2*d+1] = depthShift[2*e+1];
382: depthShift[2*e] = swap[0];
383: depthShift[2*e+1] = swap[1];
384: }
385: }
386: /* convert (oldstart, added) to (oldstart, newstart) */
387: for (d = 0; d <= depth; d++) {
388: off += depthShift[2*d+1];
389: depthShift[2*d+1] = depthShift[2*d] + off;
390: }
391: return 0;
392: }
394: /* depthShift is a list of (old, new) pairs */
395: static inline PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
396: {
397: PetscInt d;
398: PetscInt newOff = 0;
400: for (d = 0; d <= depth; d++) {
401: if (p < depthShift[2*d]) return p + newOff;
402: else newOff = depthShift[2*d+1] - depthShift[2*d];
403: }
404: return p + newOff;
405: }
407: /* depthShift is a list of (old, new) pairs */
408: static inline PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
409: {
410: PetscInt d;
411: PetscInt newOff = 0;
413: for (d = 0; d <= depth; d++) {
414: if (p < depthShift[2*d+1]) return p + newOff;
415: else newOff = depthShift[2*d] - depthShift[2*d+1];
416: }
417: return p + newOff;
418: }
420: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
421: {
422: PetscInt depth = 0, d, pStart, pEnd, p;
423: DMLabel depthLabel;
425: DMPlexGetDepth(dm, &depth);
426: if (depth < 0) return 0;
427: /* Step 1: Expand chart */
428: DMPlexGetChart(dm, &pStart, &pEnd);
429: pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
430: DMPlexSetChart(dmNew, pStart, pEnd);
431: DMCreateLabel(dmNew,"depth");
432: DMPlexGetDepthLabel(dmNew,&depthLabel);
433: DMCreateLabel(dmNew, "celltype");
434: /* Step 2: Set cone and support sizes */
435: for (d = 0; d <= depth; ++d) {
436: PetscInt pStartNew, pEndNew;
437: IS pIS;
439: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
440: pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
441: pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
442: ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS);
443: DMLabelSetStratumIS(depthLabel, d, pIS);
444: ISDestroy(&pIS);
445: for (p = pStart; p < pEnd; ++p) {
446: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
447: PetscInt size;
448: DMPolytopeType ct;
450: DMPlexGetConeSize(dm, p, &size);
451: DMPlexSetConeSize(dmNew, newp, size);
452: DMPlexGetSupportSize(dm, p, &size);
453: DMPlexSetSupportSize(dmNew, newp, size);
454: DMPlexGetCellType(dm, p, &ct);
455: DMPlexSetCellType(dmNew, newp, ct);
456: }
457: }
458: return 0;
459: }
461: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
462: {
463: PetscInt *newpoints;
464: PetscInt depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;
466: DMPlexGetDepth(dm, &depth);
467: if (depth < 0) return 0;
468: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
469: DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
470: PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
471: /* Step 5: Set cones and supports */
472: DMPlexGetChart(dm, &pStart, &pEnd);
473: for (p = pStart; p < pEnd; ++p) {
474: const PetscInt *points = NULL, *orientations = NULL;
475: PetscInt size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
477: DMPlexGetConeSize(dm, p, &size);
478: DMPlexGetCone(dm, p, &points);
479: DMPlexGetConeOrientation(dm, p, &orientations);
480: for (i = 0; i < size; ++i) {
481: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
482: }
483: DMPlexSetCone(dmNew, newp, newpoints);
484: DMPlexSetConeOrientation(dmNew, newp, orientations);
485: DMPlexGetSupportSize(dm, p, &size);
486: DMPlexGetSupportSize(dmNew, newp, &sizeNew);
487: DMPlexGetSupport(dm, p, &points);
488: for (i = 0; i < size; ++i) {
489: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
490: }
491: for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
492: DMPlexSetSupport(dmNew, newp, newpoints);
493: }
494: PetscFree(newpoints);
495: return 0;
496: }
498: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
499: {
500: PetscSection coordSection, newCoordSection;
501: Vec coordinates, newCoordinates;
502: PetscScalar *coords, *newCoords;
503: PetscInt coordSize, sStart, sEnd;
504: PetscInt dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
505: PetscBool hasCells;
507: DMGetCoordinateDim(dm, &dim);
508: DMSetCoordinateDim(dmNew, dim);
509: DMPlexGetDepth(dm, &depth);
510: /* Step 8: Convert coordinates */
511: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
512: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
513: DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
514: DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
515: DMGetCoordinateSection(dm, &coordSection);
516: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
517: PetscSectionSetNumFields(newCoordSection, 1);
518: PetscSectionSetFieldComponents(newCoordSection, 0, dim);
519: PetscSectionGetChart(coordSection, &sStart, &sEnd);
520: hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
521: PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
522: if (hasCells) {
523: for (c = cStart; c < cEnd; ++c) {
524: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;
526: PetscSectionGetDof(coordSection, c, &dof);
527: PetscSectionSetDof(newCoordSection, cNew, dof);
528: PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
529: }
530: }
531: for (v = vStartNew; v < vEndNew; ++v) {
532: PetscSectionSetDof(newCoordSection, v, dim);
533: PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
534: }
535: PetscSectionSetUp(newCoordSection);
536: DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
537: PetscSectionGetStorageSize(newCoordSection, &coordSize);
538: VecCreate(PETSC_COMM_SELF, &newCoordinates);
539: PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
540: VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
541: VecSetBlockSize(newCoordinates, dim);
542: VecSetType(newCoordinates,VECSTANDARD);
543: DMSetCoordinatesLocal(dmNew, newCoordinates);
544: DMGetCoordinatesLocal(dm, &coordinates);
545: VecGetArray(coordinates, &coords);
546: VecGetArray(newCoordinates, &newCoords);
547: if (hasCells) {
548: for (c = cStart; c < cEnd; ++c) {
549: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;
551: PetscSectionGetDof(coordSection, c, &dof);
552: PetscSectionGetOffset(coordSection, c, &off);
553: PetscSectionGetOffset(newCoordSection, cNew, &noff);
554: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
555: }
556: }
557: for (v = vStart; v < vEnd; ++v) {
558: PetscInt dof, off, noff, d;
560: PetscSectionGetDof(coordSection, v, &dof);
561: PetscSectionGetOffset(coordSection, v, &off);
562: PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
563: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
564: }
565: VecRestoreArray(coordinates, &coords);
566: VecRestoreArray(newCoordinates, &newCoords);
567: VecDestroy(&newCoordinates);
568: PetscSectionDestroy(&newCoordSection);
569: return 0;
570: }
572: static PetscErrorCode DMPlexShiftSF_Single(DM dm, PetscInt depthShift[], PetscSF sf, PetscSF sfNew)
573: {
574: const PetscSFNode *remotePoints;
575: PetscSFNode *gremotePoints;
576: const PetscInt *localPoints;
577: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
578: PetscInt numRoots, numLeaves, l, pStart, pEnd, depth = 0, totShift = 0;
580: DMPlexGetDepth(dm, &depth);
581: DMPlexGetChart(dm, &pStart, &pEnd);
582: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
583: totShift = DMPlexShiftPoint_Internal(pEnd, depth, depthShift) - pEnd;
584: if (numRoots >= 0) {
585: PetscMalloc2(numRoots, &newLocation, pEnd-pStart, &newRemoteLocation);
586: for (l = 0; l < numRoots; ++l) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
587: PetscSFBcastBegin(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE);
588: PetscSFBcastEnd(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE);
589: PetscMalloc1(numLeaves, &glocalPoints);
590: PetscMalloc1(numLeaves, &gremotePoints);
591: for (l = 0; l < numLeaves; ++l) {
592: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
593: gremotePoints[l].rank = remotePoints[l].rank;
594: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
595: }
596: PetscFree2(newLocation, newRemoteLocation);
597: PetscSFSetGraph(sfNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
598: }
599: return 0;
600: }
602: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
603: {
604: PetscSF sfPoint, sfPointNew;
605: PetscBool useNatural;
607: /* Step 9: Convert pointSF */
608: DMGetPointSF(dm, &sfPoint);
609: DMGetPointSF(dmNew, &sfPointNew);
610: DMPlexShiftSF_Single(dm, depthShift, sfPoint, sfPointNew);
611: /* Step 9b: Convert naturalSF */
612: DMGetUseNatural(dm, &useNatural);
613: if (useNatural) {
614: PetscSF sfNat, sfNatNew;
616: DMSetUseNatural(dmNew, useNatural);
617: DMGetNaturalSF(dm, &sfNat);
618: DMGetNaturalSF(dmNew, &sfNatNew);
619: DMPlexShiftSF_Single(dm, depthShift, sfNat, sfNatNew);
620: }
621: return 0;
622: }
624: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
625: {
626: PetscInt depth = 0, numLabels, l;
628: DMPlexGetDepth(dm, &depth);
629: /* Step 10: Convert labels */
630: DMGetNumLabels(dm, &numLabels);
631: for (l = 0; l < numLabels; ++l) {
632: DMLabel label, newlabel;
633: const char *lname;
634: PetscBool isDepth, isDim;
635: IS valueIS;
636: const PetscInt *values;
637: PetscInt numValues, val;
639: DMGetLabelName(dm, l, &lname);
640: PetscStrcmp(lname, "depth", &isDepth);
641: if (isDepth) continue;
642: PetscStrcmp(lname, "dim", &isDim);
643: if (isDim) continue;
644: DMCreateLabel(dmNew, lname);
645: DMGetLabel(dm, lname, &label);
646: DMGetLabel(dmNew, lname, &newlabel);
647: DMLabelGetDefaultValue(label,&val);
648: DMLabelSetDefaultValue(newlabel,val);
649: DMLabelGetValueIS(label, &valueIS);
650: ISGetLocalSize(valueIS, &numValues);
651: ISGetIndices(valueIS, &values);
652: for (val = 0; val < numValues; ++val) {
653: IS pointIS;
654: const PetscInt *points;
655: PetscInt numPoints, p;
657: DMLabelGetStratumIS(label, values[val], &pointIS);
658: ISGetLocalSize(pointIS, &numPoints);
659: ISGetIndices(pointIS, &points);
660: for (p = 0; p < numPoints; ++p) {
661: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);
663: DMLabelSetValue(newlabel, newpoint, values[val]);
664: }
665: ISRestoreIndices(pointIS, &points);
666: ISDestroy(&pointIS);
667: }
668: ISRestoreIndices(valueIS, &values);
669: ISDestroy(&valueIS);
670: }
671: return 0;
672: }
674: static PetscErrorCode DMPlexCreateVTKLabel_Internal(DM dm, PetscBool createGhostLabel, DM dmNew)
675: {
676: PetscSF sfPoint;
677: DMLabel vtkLabel, ghostLabel = NULL;
678: const PetscSFNode *leafRemote;
679: const PetscInt *leafLocal;
680: PetscInt cellHeight, cStart, cEnd, c, fStart, fEnd, f, numLeaves, l;
681: PetscMPIInt rank;
683: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
684: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
685: DMGetPointSF(dm, &sfPoint);
686: DMPlexGetVTKCellHeight(dmNew, &cellHeight);
687: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
688: PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
689: DMCreateLabel(dmNew, "vtk");
690: DMGetLabel(dmNew, "vtk", &vtkLabel);
691: if (createGhostLabel) {
692: DMCreateLabel(dmNew, "ghost");
693: DMGetLabel(dmNew, "ghost", &ghostLabel);
694: }
695: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
696: for (; c < leafLocal[l] && c < cEnd; ++c) {
697: DMLabelSetValue(vtkLabel, c, 1);
698: }
699: if (leafLocal[l] >= cEnd) break;
700: if (leafRemote[l].rank == rank) {
701: DMLabelSetValue(vtkLabel, c, 1);
702: } else if (ghostLabel) {
703: DMLabelSetValue(ghostLabel, c, 2);
704: }
705: }
706: for (; c < cEnd; ++c) {
707: DMLabelSetValue(vtkLabel, c, 1);
708: }
709: if (ghostLabel) {
710: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
711: for (f = fStart; f < fEnd; ++f) {
712: PetscInt numCells;
714: DMPlexGetSupportSize(dmNew, f, &numCells);
715: if (numCells < 2) {
716: DMLabelSetValue(ghostLabel, f, 1);
717: } else {
718: const PetscInt *cells = NULL;
719: PetscInt vA, vB;
721: DMPlexGetSupport(dmNew, f, &cells);
722: DMLabelGetValue(vtkLabel, cells[0], &vA);
723: DMLabelGetValue(vtkLabel, cells[1], &vB);
724: if (vA != 1 && vB != 1) DMLabelSetValue(ghostLabel, f, 1);
725: }
726: }
727: }
728: return 0;
729: }
731: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
732: {
733: DM refTree;
734: PetscSection pSec;
735: PetscInt *parents, *childIDs;
737: DMPlexGetReferenceTree(dm,&refTree);
738: DMPlexSetReferenceTree(dmNew,refTree);
739: DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
740: if (pSec) {
741: PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
742: PetscInt *childIDsShifted;
743: PetscSection pSecShifted;
745: PetscSectionGetChart(pSec,&pStart,&pEnd);
746: DMPlexGetDepth(dm,&depth);
747: pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
748: pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
749: PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
750: PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
751: PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
752: for (p = pStartShifted; p < pEndShifted; p++) {
753: /* start off assuming no children */
754: PetscSectionSetDof(pSecShifted,p,0);
755: }
756: for (p = pStart; p < pEnd; p++) {
757: PetscInt dof;
758: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
760: PetscSectionGetDof(pSec,p,&dof);
761: PetscSectionSetDof(pSecShifted,pNew,dof);
762: }
763: PetscSectionSetUp(pSecShifted);
764: for (p = pStart; p < pEnd; p++) {
765: PetscInt dof;
766: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
768: PetscSectionGetDof(pSec,p,&dof);
769: if (dof) {
770: PetscInt off, offNew;
772: PetscSectionGetOffset(pSec,p,&off);
773: PetscSectionGetOffset(pSecShifted,pNew,&offNew);
774: parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
775: childIDsShifted[offNew] = childIDs[off];
776: }
777: }
778: DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
779: PetscFree2(parentsShifted,childIDsShifted);
780: PetscSectionDestroy(&pSecShifted);
781: }
782: return 0;
783: }
785: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
786: {
787: PetscSF sf;
788: IS valueIS;
789: const PetscInt *values, *leaves;
790: PetscInt *depthShift;
791: PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
792: PetscBool isper;
793: const PetscReal *maxCell, *L;
794: const DMBoundaryType *bd;
796: DMGetPointSF(dm, &sf);
797: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
798: nleaves = PetscMax(0, nleaves);
799: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
800: /* Count ghost cells */
801: DMLabelGetValueIS(label, &valueIS);
802: ISGetLocalSize(valueIS, &numFS);
803: ISGetIndices(valueIS, &values);
804: Ng = 0;
805: for (fs = 0; fs < numFS; ++fs) {
806: IS faceIS;
807: const PetscInt *faces;
808: PetscInt numFaces, f, numBdFaces = 0;
810: DMLabelGetStratumIS(label, values[fs], &faceIS);
811: ISGetLocalSize(faceIS, &numFaces);
812: ISGetIndices(faceIS, &faces);
813: for (f = 0; f < numFaces; ++f) {
814: PetscInt numChildren;
816: PetscFindInt(faces[f], nleaves, leaves, &loc);
817: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
818: /* non-local and ancestors points don't get to register ghosts */
819: if (loc >= 0 || numChildren) continue;
820: if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
821: }
822: Ng += numBdFaces;
823: ISRestoreIndices(faceIS, &faces);
824: ISDestroy(&faceIS);
825: }
826: DMPlexGetDepth(dm, &depth);
827: PetscMalloc1(2*(depth+1), &depthShift);
828: for (d = 0; d <= depth; d++) {
829: PetscInt dEnd;
831: DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
832: depthShift[2*d] = dEnd;
833: depthShift[2*d+1] = 0;
834: }
835: if (depth >= 0) depthShift[2*depth+1] = Ng;
836: DMPlexShiftPointSetUp_Internal(depth,depthShift);
837: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
838: /* Step 3: Set cone/support sizes for new points */
839: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
840: for (c = cEnd; c < cEnd + Ng; ++c) {
841: DMPlexSetConeSize(gdm, c, 1);
842: }
843: for (fs = 0; fs < numFS; ++fs) {
844: IS faceIS;
845: const PetscInt *faces;
846: PetscInt numFaces, f;
848: DMLabelGetStratumIS(label, values[fs], &faceIS);
849: ISGetLocalSize(faceIS, &numFaces);
850: ISGetIndices(faceIS, &faces);
851: for (f = 0; f < numFaces; ++f) {
852: PetscInt size, numChildren;
854: PetscFindInt(faces[f], nleaves, leaves, &loc);
855: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
856: if (loc >= 0 || numChildren) continue;
857: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
858: DMPlexGetSupportSize(dm, faces[f], &size);
860: DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
861: }
862: ISRestoreIndices(faceIS, &faces);
863: ISDestroy(&faceIS);
864: }
865: /* Step 4: Setup ghosted DM */
866: DMSetUp(gdm);
867: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
868: /* Step 6: Set cones and supports for new points */
869: ghostCell = cEnd;
870: for (fs = 0; fs < numFS; ++fs) {
871: IS faceIS;
872: const PetscInt *faces;
873: PetscInt numFaces, f;
875: DMLabelGetStratumIS(label, values[fs], &faceIS);
876: ISGetLocalSize(faceIS, &numFaces);
877: ISGetIndices(faceIS, &faces);
878: for (f = 0; f < numFaces; ++f) {
879: PetscInt newFace = faces[f] + Ng, numChildren;
881: PetscFindInt(faces[f], nleaves, leaves, &loc);
882: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
883: if (loc >= 0 || numChildren) continue;
884: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
885: DMPlexSetCone(gdm, ghostCell, &newFace);
886: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
887: ++ghostCell;
888: }
889: ISRestoreIndices(faceIS, &faces);
890: ISDestroy(&faceIS);
891: }
892: ISRestoreIndices(valueIS, &values);
893: ISDestroy(&valueIS);
894: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
895: DMPlexShiftSF_Internal(dm, depthShift, gdm);
896: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
897: DMPlexCreateVTKLabel_Internal(dm, PETSC_TRUE, gdm);
898: DMPlexShiftTree_Internal(dm, depthShift, gdm);
899: PetscFree(depthShift);
900: for (c = cEnd; c < cEnd + Ng; ++c) {
901: DMPlexSetCellType(gdm, c, DM_POLYTOPE_FV_GHOST);
902: }
903: /* Step 7: Periodicity */
904: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
905: DMSetPeriodicity(gdm, isper, maxCell, L, bd);
906: if (numGhostCells) *numGhostCells = Ng;
907: return 0;
908: }
910: /*@C
911: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
913: Collective on dm
915: Input Parameters:
916: + dm - The original DM
917: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
919: Output Parameters:
920: + numGhostCells - The number of ghost cells added to the DM
921: - dmGhosted - The new DM
923: Note: If no label exists of that name, one will be created marking all boundary faces
925: Level: developer
927: .seealso: DMCreate()
928: @*/
929: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
930: {
931: DM gdm;
932: DMLabel label;
933: const char *name = labelName ? labelName : "Face Sets";
934: PetscInt dim, Ng = 0;
935: PetscBool useCone, useClosure;
940: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
941: DMSetType(gdm, DMPLEX);
942: DMGetDimension(dm, &dim);
943: DMSetDimension(gdm, dim);
944: DMGetBasicAdjacency(dm, &useCone, &useClosure);
945: DMSetBasicAdjacency(gdm, useCone, useClosure);
946: DMGetLabel(dm, name, &label);
947: if (!label) {
948: /* Get label for boundary faces */
949: DMCreateLabel(dm, name);
950: DMGetLabel(dm, name, &label);
951: DMPlexMarkBoundaryFaces(dm, 1, label);
952: }
953: DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm);
954: DMCopyDisc(dm, gdm);
955: DMPlexCopy_Internal(dm, PETSC_TRUE, gdm);
956: gdm->setfromoptionscalled = dm->setfromoptionscalled;
957: if (numGhostCells) *numGhostCells = Ng;
958: *dmGhosted = gdm;
959: return 0;
960: }
962: /*
963: We are adding three kinds of points here:
964: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
965: Non-replicated: Points which exist on the fault, but are not replicated
966: Hybrid: Entirely new points, such as cohesive cells
968: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
969: each depth so that the new split/hybrid points can be inserted as a block.
970: */
971: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
972: {
973: MPI_Comm comm;
974: IS valueIS;
975: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
976: const PetscInt *values; /* List of depths for which we have replicated points */
977: IS *splitIS;
978: IS *unsplitIS;
979: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
980: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
981: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
982: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
983: const PetscInt **splitPoints; /* Replicated points for each depth */
984: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
985: PetscSection coordSection;
986: Vec coordinates;
987: PetscScalar *coords;
988: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
989: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
990: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
991: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
992: PetscInt *coneNew, *coneONew, *supportNew;
993: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
995: PetscObjectGetComm((PetscObject)dm,&comm);
996: DMGetDimension(dm, &dim);
997: DMPlexGetDepth(dm, &depth);
998: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
999: /* We do not want this label automatically computed, instead we compute it here */
1000: DMCreateLabel(sdm, "celltype");
1001: /* Count split points and add cohesive cells */
1002: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
1003: PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
1004: PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
1005: for (d = 0; d <= depth; ++d) {
1006: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
1007: DMPlexGetTensorPrismBounds_Internal(dm, d, &depthMax[d], NULL);
1008: depthEnd[d] = pMaxNew[d];
1009: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
1010: numSplitPoints[d] = 0;
1011: numUnsplitPoints[d] = 0;
1012: numHybridPoints[d] = 0;
1013: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
1014: splitPoints[d] = NULL;
1015: unsplitPoints[d] = NULL;
1016: splitIS[d] = NULL;
1017: unsplitIS[d] = NULL;
1018: /* we are shifting the existing hybrid points with the stratum behind them, so
1019: * the split comes at the end of the normal points, i.e., at depthMax[d] */
1020: depthShift[2*d] = depthMax[d];
1021: depthShift[2*d+1] = 0;
1022: }
1023: if (label) {
1024: DMLabelGetValueIS(label, &valueIS);
1025: ISGetLocalSize(valueIS, &numSP);
1026: ISGetIndices(valueIS, &values);
1027: }
1028: for (sp = 0; sp < numSP; ++sp) {
1029: const PetscInt dep = values[sp];
1031: if ((dep < 0) || (dep > depth)) continue;
1032: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
1033: if (splitIS[dep]) {
1034: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
1035: ISGetIndices(splitIS[dep], &splitPoints[dep]);
1036: }
1037: DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
1038: if (unsplitIS[dep]) {
1039: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
1040: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
1041: }
1042: }
1043: /* Calculate number of hybrid points */
1044: for (d = 1; d <= depth; ++d) numHybridPoints[d] = numSplitPoints[d-1] + numUnsplitPoints[d-1]; /* There is a hybrid cell/face/edge for every split face/edge/vertex */
1045: for (d = 0; d <= depth; ++d) depthShift[2*d+1] = numSplitPoints[d] + numHybridPoints[d];
1046: DMPlexShiftPointSetUp_Internal(depth,depthShift);
1047: /* the end of the points in this stratum that come before the new points:
1048: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
1049: * added points */
1050: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
1051: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
1052: /* Step 3: Set cone/support sizes for new points */
1053: for (dep = 0; dep <= depth; ++dep) {
1054: for (p = 0; p < numSplitPoints[dep]; ++p) {
1055: const PetscInt oldp = splitPoints[dep][p];
1056: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1057: const PetscInt splitp = p + pMaxNew[dep];
1058: const PetscInt *support;
1059: DMPolytopeType ct;
1060: PetscInt coneSize, supportSize, qf, qn, qp, e;
1062: DMPlexGetConeSize(dm, oldp, &coneSize);
1063: DMPlexSetConeSize(sdm, splitp, coneSize);
1064: DMPlexGetSupportSize(dm, oldp, &supportSize);
1065: DMPlexSetSupportSize(sdm, splitp, supportSize);
1066: DMPlexGetCellType(dm, oldp, &ct);
1067: DMPlexSetCellType(sdm, splitp, ct);
1068: if (dep == depth-1) {
1069: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1071: /* Add cohesive cells, they are prisms */
1072: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
1073: switch (coneSize) {
1074: case 2: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR);break;
1075: case 3: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR);break;
1076: case 4: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_QUAD_PRISM_TENSOR);break;
1077: }
1078: } else if (dep == 0) {
1079: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1081: DMPlexGetSupport(dm, oldp, &support);
1082: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1083: PetscInt val;
1085: DMLabelGetValue(label, support[e], &val);
1086: if (val == 1) ++qf;
1087: if ((val == 1) || (val == (shift + 1))) ++qn;
1088: if ((val == 1) || (val == -(shift + 1))) ++qp;
1089: }
1090: /* Split old vertex: Edges into original vertex and new cohesive edge */
1091: DMPlexSetSupportSize(sdm, newp, qn+1);
1092: /* Split new vertex: Edges into split vertex and new cohesive edge */
1093: DMPlexSetSupportSize(sdm, splitp, qp+1);
1094: /* Add hybrid edge */
1095: DMPlexSetConeSize(sdm, hybedge, 2);
1096: DMPlexSetSupportSize(sdm, hybedge, qf);
1097: DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1098: } else if (dep == dim-2) {
1099: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1101: DMPlexGetSupport(dm, oldp, &support);
1102: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1103: PetscInt val;
1105: DMLabelGetValue(label, support[e], &val);
1106: if (val == dim-1) ++qf;
1107: if ((val == dim-1) || (val == (shift + dim-1))) ++qn;
1108: if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
1109: }
1110: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1111: DMPlexSetSupportSize(sdm, newp, qn+1);
1112: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1113: DMPlexSetSupportSize(sdm, splitp, qp+1);
1114: /* Add hybrid face */
1115: DMPlexSetConeSize(sdm, hybface, 4);
1116: DMPlexSetSupportSize(sdm, hybface, qf);
1117: DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1118: }
1119: }
1120: }
1121: for (dep = 0; dep <= depth; ++dep) {
1122: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1123: const PetscInt oldp = unsplitPoints[dep][p];
1124: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1125: const PetscInt *support;
1126: PetscInt coneSize, supportSize, qf, e, s;
1128: DMPlexGetConeSize(dm, oldp, &coneSize);
1129: DMPlexGetSupportSize(dm, oldp, &supportSize);
1130: DMPlexGetSupport(dm, oldp, &support);
1131: if (dep == 0) {
1132: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1134: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1135: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1136: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1137: if (e >= 0) ++qf;
1138: }
1139: DMPlexSetSupportSize(sdm, newp, qf+2);
1140: /* Add hybrid edge */
1141: DMPlexSetConeSize(sdm, hybedge, 2);
1142: for (e = 0, qf = 0; e < supportSize; ++e) {
1143: PetscInt val;
1145: DMLabelGetValue(label, support[e], &val);
1146: /* Split and unsplit edges produce hybrid faces */
1147: if (val == 1) ++qf;
1148: if (val == (shift2 + 1)) ++qf;
1149: }
1150: DMPlexSetSupportSize(sdm, hybedge, qf);
1151: DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1152: } else if (dep == dim-2) {
1153: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1154: PetscInt val;
1156: for (e = 0, qf = 0; e < supportSize; ++e) {
1157: DMLabelGetValue(label, support[e], &val);
1158: if (val == dim-1) qf += 2;
1159: else ++qf;
1160: }
1161: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1162: DMPlexSetSupportSize(sdm, newp, qf+2);
1163: /* Add hybrid face */
1164: for (e = 0, qf = 0; e < supportSize; ++e) {
1165: DMLabelGetValue(label, support[e], &val);
1166: if (val == dim-1) ++qf;
1167: }
1168: DMPlexSetConeSize(sdm, hybface, 4);
1169: DMPlexSetSupportSize(sdm, hybface, qf);
1170: DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1171: }
1172: }
1173: }
1174: /* Step 4: Setup split DM */
1175: DMSetUp(sdm);
1176: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1177: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1178: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1179: /* Step 6: Set cones and supports for new points */
1180: for (dep = 0; dep <= depth; ++dep) {
1181: for (p = 0; p < numSplitPoints[dep]; ++p) {
1182: const PetscInt oldp = splitPoints[dep][p];
1183: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1184: const PetscInt splitp = p + pMaxNew[dep];
1185: const PetscInt *cone, *support, *ornt;
1186: DMPolytopeType ct;
1187: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1189: DMPlexGetCellType(dm, oldp, &ct);
1190: DMPlexGetConeSize(dm, oldp, &coneSize);
1191: DMPlexGetCone(dm, oldp, &cone);
1192: DMPlexGetConeOrientation(dm, oldp, &ornt);
1193: DMPlexGetSupportSize(dm, oldp, &supportSize);
1194: DMPlexGetSupport(dm, oldp, &support);
1195: if (dep == depth-1) {
1196: PetscBool hasUnsplit = PETSC_FALSE;
1197: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1198: const PetscInt *supportF;
1200: /* Split face: copy in old face to new face to start */
1201: DMPlexGetSupport(sdm, newp, &supportF);
1202: DMPlexSetSupport(sdm, splitp, supportF);
1203: /* Split old face: old vertices/edges in cone so no change */
1204: /* Split new face: new vertices/edges in cone */
1205: for (q = 0; q < coneSize; ++q) {
1206: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1207: if (v < 0) {
1208: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1210: coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1211: hasUnsplit = PETSC_TRUE;
1212: } else {
1213: coneNew[2+q] = v + pMaxNew[dep-1];
1214: if (dep > 1) {
1215: const PetscInt *econe;
1216: PetscInt econeSize, r, vs, vu;
1218: DMPlexGetConeSize(dm, cone[q], &econeSize);
1219: DMPlexGetCone(dm, cone[q], &econe);
1220: for (r = 0; r < econeSize; ++r) {
1221: PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);
1222: PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1223: if (vs >= 0) continue;
1225: hasUnsplit = PETSC_TRUE;
1226: }
1227: }
1228: }
1229: }
1230: DMPlexSetCone(sdm, splitp, &coneNew[2]);
1231: DMPlexSetConeOrientation(sdm, splitp, ornt);
1232: /* Face support */
1233: for (s = 0; s < supportSize; ++s) {
1234: PetscInt val;
1236: DMLabelGetValue(label, support[s], &val);
1237: if (val < 0) {
1238: /* Split old face: Replace negative side cell with cohesive cell */
1239: DMPlexInsertSupport(sdm, newp, s, hybcell);
1240: } else {
1241: /* Split new face: Replace positive side cell with cohesive cell */
1242: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1243: /* Get orientation for cohesive face */
1244: {
1245: const PetscInt *ncone, *nconeO;
1246: PetscInt nconeSize, nc;
1248: DMPlexGetConeSize(dm, support[s], &nconeSize);
1249: DMPlexGetCone(dm, support[s], &ncone);
1250: DMPlexGetConeOrientation(dm, support[s], &nconeO);
1251: for (nc = 0; nc < nconeSize; ++nc) {
1252: if (ncone[nc] == oldp) {
1253: coneONew[0] = nconeO[nc];
1254: break;
1255: }
1256: }
1258: }
1259: }
1260: }
1261: /* Cohesive cell: Old and new split face, then new cohesive faces */
1262: const PetscInt *arr = DMPolytopeTypeGetArrangment(ct, coneONew[0]);
1264: coneNew[0] = newp; /* Extracted negative side orientation above */
1265: coneNew[1] = splitp;
1266: coneONew[1] = coneONew[0];
1267: for (q = 0; q < coneSize; ++q) {
1268: /* Hybrid faces must follow order from oriented end face */
1269: const PetscInt qa = arr[q*2+0];
1270: const PetscInt qo = arr[q*2+1];
1271: DMPolytopeType ft = dep == 2 ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
1273: PetscFindInt(cone[qa], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1274: if (v < 0) {
1275: PetscFindInt(cone[qa], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1276: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1277: } else {
1278: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep];
1279: }
1280: coneONew[2+q] = DMPolytopeTypeComposeOrientation(ft, qo, ornt[qa]);
1281: }
1282: DMPlexSetCone(sdm, hybcell, coneNew);
1283: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1284: /* Label the hybrid cells on the boundary of the split */
1285: if (hasUnsplit) DMLabelSetValue(label, -hybcell, dim);
1286: } else if (dep == 0) {
1287: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1289: /* Split old vertex: Edges in old split faces and new cohesive edge */
1290: for (e = 0, qn = 0; e < supportSize; ++e) {
1291: PetscInt val;
1293: DMLabelGetValue(label, support[e], &val);
1294: if ((val == 1) || (val == (shift + 1))) {
1295: supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1296: }
1297: }
1298: supportNew[qn] = hybedge;
1299: DMPlexSetSupport(sdm, newp, supportNew);
1300: /* Split new vertex: Edges in new split faces and new cohesive edge */
1301: for (e = 0, qp = 0; e < supportSize; ++e) {
1302: PetscInt val, edge;
1304: DMLabelGetValue(label, support[e], &val);
1305: if (val == 1) {
1306: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1308: supportNew[qp++] = edge + pMaxNew[dep+1];
1309: } else if (val == -(shift + 1)) {
1310: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1311: }
1312: }
1313: supportNew[qp] = hybedge;
1314: DMPlexSetSupport(sdm, splitp, supportNew);
1315: /* Hybrid edge: Old and new split vertex */
1316: coneNew[0] = newp;
1317: coneNew[1] = splitp;
1318: DMPlexSetCone(sdm, hybedge, coneNew);
1319: for (e = 0, qf = 0; e < supportSize; ++e) {
1320: PetscInt val, edge;
1322: DMLabelGetValue(label, support[e], &val);
1323: if (val == 1) {
1324: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1326: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1327: }
1328: }
1329: DMPlexSetSupport(sdm, hybedge, supportNew);
1330: } else if (dep == dim-2) {
1331: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1333: /* Split old edge: old vertices in cone so no change */
1334: /* Split new edge: new vertices in cone */
1335: for (q = 0; q < coneSize; ++q) {
1336: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1337: if (v < 0) {
1338: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1340: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1341: } else {
1342: coneNew[q] = v + pMaxNew[dep-1];
1343: }
1344: }
1345: DMPlexSetCone(sdm, splitp, coneNew);
1346: /* Split old edge: Faces in positive side cells and old split faces */
1347: for (e = 0, q = 0; e < supportSize; ++e) {
1348: PetscInt val;
1350: DMLabelGetValue(label, support[e], &val);
1351: if (val == dim-1) {
1352: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1353: } else if (val == (shift + dim-1)) {
1354: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1355: }
1356: }
1357: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1358: DMPlexSetSupport(sdm, newp, supportNew);
1359: /* Split new edge: Faces in negative side cells and new split faces */
1360: for (e = 0, q = 0; e < supportSize; ++e) {
1361: PetscInt val, face;
1363: DMLabelGetValue(label, support[e], &val);
1364: if (val == dim-1) {
1365: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1367: supportNew[q++] = face + pMaxNew[dep+1];
1368: } else if (val == -(shift + dim-1)) {
1369: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1370: }
1371: }
1372: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1373: DMPlexSetSupport(sdm, splitp, supportNew);
1374: /* Hybrid face */
1375: coneNew[0] = newp;
1376: coneNew[1] = splitp;
1377: for (v = 0; v < coneSize; ++v) {
1378: PetscInt vertex;
1379: PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1380: if (vertex < 0) {
1381: PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1383: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1384: } else {
1385: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1386: }
1387: }
1388: DMPlexSetCone(sdm, hybface, coneNew);
1389: for (e = 0, qf = 0; e < supportSize; ++e) {
1390: PetscInt val, face;
1392: DMLabelGetValue(label, support[e], &val);
1393: if (val == dim-1) {
1394: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1396: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1397: }
1398: }
1399: DMPlexSetSupport(sdm, hybface, supportNew);
1400: }
1401: }
1402: }
1403: for (dep = 0; dep <= depth; ++dep) {
1404: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1405: const PetscInt oldp = unsplitPoints[dep][p];
1406: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1407: const PetscInt *cone, *support;
1408: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1410: DMPlexGetConeSize(dm, oldp, &coneSize);
1411: DMPlexGetCone(dm, oldp, &cone);
1412: DMPlexGetSupportSize(dm, oldp, &supportSize);
1413: DMPlexGetSupport(dm, oldp, &support);
1414: if (dep == 0) {
1415: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1417: /* Unsplit vertex */
1418: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1419: for (s = 0, q = 0; s < supportSize; ++s) {
1420: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1421: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1422: if (e >= 0) {
1423: supportNew[q++] = e + pMaxNew[dep+1];
1424: }
1425: }
1426: supportNew[q++] = hybedge;
1427: supportNew[q++] = hybedge;
1429: DMPlexSetSupport(sdm, newp, supportNew);
1430: /* Hybrid edge */
1431: coneNew[0] = newp;
1432: coneNew[1] = newp;
1433: DMPlexSetCone(sdm, hybedge, coneNew);
1434: for (e = 0, qf = 0; e < supportSize; ++e) {
1435: PetscInt val, edge;
1437: DMLabelGetValue(label, support[e], &val);
1438: if (val == 1) {
1439: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1441: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1442: } else if (val == (shift2 + 1)) {
1443: PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1445: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1446: }
1447: }
1448: DMPlexSetSupport(sdm, hybedge, supportNew);
1449: } else if (dep == dim-2) {
1450: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1452: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1453: for (f = 0, qf = 0; f < supportSize; ++f) {
1454: PetscInt val, face;
1456: DMLabelGetValue(label, support[f], &val);
1457: if (val == dim-1) {
1458: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1460: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1461: supportNew[qf++] = face + pMaxNew[dep+1];
1462: } else {
1463: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1464: }
1465: }
1466: supportNew[qf++] = hybface;
1467: supportNew[qf++] = hybface;
1468: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1470: DMPlexSetSupport(sdm, newp, supportNew);
1471: /* Add hybrid face */
1472: coneNew[0] = newp;
1473: coneNew[1] = newp;
1474: PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1476: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1477: PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1479: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1480: DMPlexSetCone(sdm, hybface, coneNew);
1481: for (f = 0, qf = 0; f < supportSize; ++f) {
1482: PetscInt val, face;
1484: DMLabelGetValue(label, support[f], &val);
1485: if (val == dim-1) {
1486: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1487: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1488: }
1489: }
1490: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1492: DMPlexSetSupport(sdm, hybface, supportNew);
1493: }
1494: }
1495: }
1496: /* Step 6b: Replace split points in negative side cones */
1497: for (sp = 0; sp < numSP; ++sp) {
1498: PetscInt dep = values[sp];
1499: IS pIS;
1500: PetscInt numPoints;
1501: const PetscInt *points;
1503: if (dep >= 0) continue;
1504: DMLabelGetStratumIS(label, dep, &pIS);
1505: if (!pIS) continue;
1506: dep = -dep - shift;
1507: ISGetLocalSize(pIS, &numPoints);
1508: ISGetIndices(pIS, &points);
1509: for (p = 0; p < numPoints; ++p) {
1510: const PetscInt oldp = points[p];
1511: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1512: const PetscInt *cone;
1513: PetscInt coneSize, c;
1514: /* PetscBool replaced = PETSC_FALSE; */
1516: /* Negative edge: replace split vertex */
1517: /* Negative cell: replace split face */
1518: DMPlexGetConeSize(sdm, newp, &coneSize);
1519: DMPlexGetCone(sdm, newp, &cone);
1520: for (c = 0; c < coneSize; ++c) {
1521: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1522: PetscInt csplitp, cp, val;
1524: DMLabelGetValue(label, coldp, &val);
1525: if (val == dep-1) {
1526: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1528: csplitp = pMaxNew[dep-1] + cp;
1529: DMPlexInsertCone(sdm, newp, c, csplitp);
1530: /* replaced = PETSC_TRUE; */
1531: }
1532: }
1533: /* Cells with only a vertex or edge on the submesh have no replacement */
1535: }
1536: ISRestoreIndices(pIS, &points);
1537: ISDestroy(&pIS);
1538: }
1539: /* Step 7: Coordinates */
1540: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1541: DMGetCoordinateSection(sdm, &coordSection);
1542: DMGetCoordinatesLocal(sdm, &coordinates);
1543: VecGetArray(coordinates, &coords);
1544: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1545: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1546: const PetscInt splitp = pMaxNew[0] + v;
1547: PetscInt dof, off, soff, d;
1549: PetscSectionGetDof(coordSection, newp, &dof);
1550: PetscSectionGetOffset(coordSection, newp, &off);
1551: PetscSectionGetOffset(coordSection, splitp, &soff);
1552: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1553: }
1554: VecRestoreArray(coordinates, &coords);
1555: /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1556: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1557: /* Step 9: Labels */
1558: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1559: DMPlexCreateVTKLabel_Internal(dm, PETSC_FALSE, sdm);
1560: DMGetNumLabels(sdm, &numLabels);
1561: for (dep = 0; dep <= depth; ++dep) {
1562: for (p = 0; p < numSplitPoints[dep]; ++p) {
1563: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1564: const PetscInt splitp = pMaxNew[dep] + p;
1565: PetscInt l;
1567: if (splitLabel) {
1568: const PetscInt val = 100 + dep;
1570: DMLabelSetValue(splitLabel, newp, val);
1571: DMLabelSetValue(splitLabel, splitp, -val);
1572: }
1573: for (l = 0; l < numLabels; ++l) {
1574: DMLabel mlabel;
1575: const char *lname;
1576: PetscInt val;
1577: PetscBool isDepth;
1579: DMGetLabelName(sdm, l, &lname);
1580: PetscStrcmp(lname, "depth", &isDepth);
1581: if (isDepth) continue;
1582: DMGetLabel(sdm, lname, &mlabel);
1583: DMLabelGetValue(mlabel, newp, &val);
1584: if (val >= 0) {
1585: DMLabelSetValue(mlabel, splitp, val);
1586: }
1587: }
1588: }
1589: }
1590: for (sp = 0; sp < numSP; ++sp) {
1591: const PetscInt dep = values[sp];
1593: if ((dep < 0) || (dep > depth)) continue;
1594: if (splitIS[dep]) ISRestoreIndices(splitIS[dep], &splitPoints[dep]);
1595: ISDestroy(&splitIS[dep]);
1596: if (unsplitIS[dep]) ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);
1597: ISDestroy(&unsplitIS[dep]);
1598: }
1599: if (label) {
1600: ISRestoreIndices(valueIS, &values);
1601: ISDestroy(&valueIS);
1602: }
1603: for (d = 0; d <= depth; ++d) {
1604: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1605: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1606: }
1607: PetscFree3(coneNew, coneONew, supportNew);
1608: PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1609: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1610: return 0;
1611: }
1613: /*@C
1614: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1616: Collective on dm
1618: Input Parameters:
1619: + dm - The original DM
1620: - label - The label specifying the boundary faces (this could be auto-generated)
1622: Output Parameters:
1623: + splitLabel - The label containing the split points, or NULL if no output is desired
1624: - dmSplit - The new DM
1626: Level: developer
1628: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1629: @*/
1630: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
1631: {
1632: DM sdm;
1633: PetscInt dim;
1637: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1638: DMSetType(sdm, DMPLEX);
1639: DMGetDimension(dm, &dim);
1640: DMSetDimension(sdm, dim);
1641: switch (dim) {
1642: case 2:
1643: case 3:
1644: DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
1645: break;
1646: default:
1647: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1648: }
1649: *dmSplit = sdm;
1650: return 0;
1651: }
1653: /* Returns the side of the surface for a given cell with a face on the surface */
1654: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1655: {
1656: const PetscInt *cone, *ornt;
1657: PetscInt dim, coneSize, c;
1659: *pos = PETSC_TRUE;
1660: DMGetDimension(dm, &dim);
1661: DMPlexGetConeSize(dm, cell, &coneSize);
1662: DMPlexGetCone(dm, cell, &cone);
1663: DMPlexGetConeOrientation(dm, cell, &ornt);
1664: for (c = 0; c < coneSize; ++c) {
1665: if (cone[c] == face) {
1666: PetscInt o = ornt[c];
1668: if (subdm) {
1669: const PetscInt *subcone, *subornt;
1670: PetscInt subpoint, subface, subconeSize, sc;
1672: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1673: PetscFindInt(face, numSubpoints, subpoints, &subface);
1674: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1675: DMPlexGetCone(subdm, subpoint, &subcone);
1676: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1677: for (sc = 0; sc < subconeSize; ++sc) {
1678: if (subcone[sc] == subface) {
1679: o = subornt[0];
1680: break;
1681: }
1682: }
1684: }
1685: if (o >= 0) *pos = PETSC_TRUE;
1686: else *pos = PETSC_FALSE;
1687: break;
1688: }
1689: }
1691: return 0;
1692: }
1694: /*@
1695: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1696: to complete the surface
1698: Input Parameters:
1699: + dm - The DM
1700: . label - A DMLabel marking the surface
1701: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1702: . flip - Flag to flip the submesh normal and replace points on the other side
1703: - subdm - The subDM associated with the label, or NULL
1705: Output Parameter:
1706: . label - A DMLabel marking all surface points
1708: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1710: Level: developer
1712: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1713: @*/
1714: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1715: {
1716: DMLabel depthLabel;
1717: IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1718: const PetscInt *points, *subpoints;
1719: const PetscInt rev = flip ? -1 : 1;
1720: PetscInt shift = 100, shift2 = 200, dim, depth, dep, cStart, cEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1722: DMPlexGetDepth(dm, &depth);
1723: DMGetDimension(dm, &dim);
1724: DMPlexGetDepthLabel(dm, &depthLabel);
1725: if (subdm) {
1726: DMPlexGetSubpointIS(subdm, &subpointIS);
1727: if (subpointIS) {
1728: ISGetLocalSize(subpointIS, &numSubpoints);
1729: ISGetIndices(subpointIS, &subpoints);
1730: }
1731: }
1732: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1733: DMLabelGetStratumIS(label, dim-1, &dimIS);
1734: if (!dimIS) return 0;
1735: ISGetLocalSize(dimIS, &numPoints);
1736: ISGetIndices(dimIS, &points);
1737: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1738: const PetscInt *support;
1739: PetscInt supportSize, s;
1741: DMPlexGetSupportSize(dm, points[p], &supportSize);
1742: #if 0
1743: if (supportSize != 2) {
1744: const PetscInt *lp;
1745: PetscInt Nlp, pind;
1747: /* Check that for a cell with a single support face, that face is in the SF */
1748: /* THis check only works for the remote side. We would need root side information */
1749: PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1750: PetscFindInt(points[p], Nlp, lp, &pind);
1752: }
1753: #endif
1754: DMPlexGetSupport(dm, points[p], &support);
1755: for (s = 0; s < supportSize; ++s) {
1756: const PetscInt *cone;
1757: PetscInt coneSize, c;
1758: PetscBool pos;
1760: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1761: if (pos) DMLabelSetValue(label, support[s], rev*(shift+dim));
1762: else DMLabelSetValue(label, support[s], -rev*(shift+dim));
1763: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1764: /* Put faces touching the fault in the label */
1765: DMPlexGetConeSize(dm, support[s], &coneSize);
1766: DMPlexGetCone(dm, support[s], &cone);
1767: for (c = 0; c < coneSize; ++c) {
1768: const PetscInt point = cone[c];
1770: DMLabelGetValue(label, point, &val);
1771: if (val == -1) {
1772: PetscInt *closure = NULL;
1773: PetscInt closureSize, cl;
1775: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1776: for (cl = 0; cl < closureSize*2; cl += 2) {
1777: const PetscInt clp = closure[cl];
1778: PetscInt bval = -1;
1780: DMLabelGetValue(label, clp, &val);
1781: if (blabel) DMLabelGetValue(blabel, clp, &bval);
1782: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1783: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1784: break;
1785: }
1786: }
1787: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1788: }
1789: }
1790: }
1791: }
1792: ISRestoreIndices(dimIS, &points);
1793: ISDestroy(&dimIS);
1794: if (subpointIS) ISRestoreIndices(subpointIS, &subpoints);
1795: /* Mark boundary points as unsplit */
1796: if (blabel) {
1797: DMLabelGetStratumIS(blabel, 1, &dimIS);
1798: ISGetLocalSize(dimIS, &numPoints);
1799: ISGetIndices(dimIS, &points);
1800: for (p = 0; p < numPoints; ++p) {
1801: const PetscInt point = points[p];
1802: PetscInt val, bval;
1804: DMLabelGetValue(blabel, point, &bval);
1805: if (bval >= 0) {
1806: DMLabelGetValue(label, point, &val);
1807: if ((val < 0) || (val > dim)) {
1808: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1809: DMLabelClearValue(blabel, point, bval);
1810: }
1811: }
1812: }
1813: for (p = 0; p < numPoints; ++p) {
1814: const PetscInt point = points[p];
1815: PetscInt val, bval;
1817: DMLabelGetValue(blabel, point, &bval);
1818: if (bval >= 0) {
1819: const PetscInt *cone, *support;
1820: PetscInt coneSize, supportSize, s, valA, valB, valE;
1822: /* Mark as unsplit */
1823: DMLabelGetValue(label, point, &val);
1825: DMLabelClearValue(label, point, val);
1826: DMLabelSetValue(label, point, shift2+val);
1827: /* Check for cross-edge
1828: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1829: if (val != 0) continue;
1830: DMPlexGetSupport(dm, point, &support);
1831: DMPlexGetSupportSize(dm, point, &supportSize);
1832: for (s = 0; s < supportSize; ++s) {
1833: DMPlexGetCone(dm, support[s], &cone);
1834: DMPlexGetConeSize(dm, support[s], &coneSize);
1836: DMLabelGetValue(blabel, cone[0], &valA);
1837: DMLabelGetValue(blabel, cone[1], &valB);
1838: DMLabelGetValue(blabel, support[s], &valE);
1839: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) DMLabelSetValue(blabel, support[s], 2);
1840: }
1841: }
1842: }
1843: ISRestoreIndices(dimIS, &points);
1844: ISDestroy(&dimIS);
1845: }
1846: /* Search for other cells/faces/edges connected to the fault by a vertex */
1847: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1848: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1849: DMLabelGetStratumIS(label, 0, &dimIS);
1850: /* TODO Why are we including cross edges here? Shouldn't they be in the star of boundary vertices? */
1851: if (blabel) DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);
1852: if (dimIS && crossEdgeIS) {
1853: IS vertIS = dimIS;
1855: ISExpand(vertIS, crossEdgeIS, &dimIS);
1856: ISDestroy(&crossEdgeIS);
1857: ISDestroy(&vertIS);
1858: }
1859: if (!dimIS) {
1860: return 0;
1861: }
1862: ISGetLocalSize(dimIS, &numPoints);
1863: ISGetIndices(dimIS, &points);
1864: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1865: PetscInt *star = NULL;
1866: PetscInt starSize, s;
1867: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1869: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1870: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1871: while (again) {
1873: again = 0;
1874: for (s = 0; s < starSize*2; s += 2) {
1875: const PetscInt point = star[s];
1876: const PetscInt *cone;
1877: PetscInt coneSize, c;
1879: if ((point < cStart) || (point >= cEnd)) continue;
1880: DMLabelGetValue(label, point, &val);
1881: if (val != -1) continue;
1882: again = again == 1 ? 1 : 2;
1883: DMPlexGetConeSize(dm, point, &coneSize);
1884: DMPlexGetCone(dm, point, &cone);
1885: for (c = 0; c < coneSize; ++c) {
1886: DMLabelGetValue(label, cone[c], &val);
1887: if (val != -1) {
1888: const PetscInt *ccone;
1889: PetscInt cconeSize, cc, side;
1892: if (val > 0) side = 1;
1893: else side = -1;
1894: DMLabelSetValue(label, point, side*(shift+dim));
1895: /* Mark cell faces which touch the fault */
1896: DMPlexGetConeSize(dm, point, &cconeSize);
1897: DMPlexGetCone(dm, point, &ccone);
1898: for (cc = 0; cc < cconeSize; ++cc) {
1899: PetscInt *closure = NULL;
1900: PetscInt closureSize, cl;
1902: DMLabelGetValue(label, ccone[cc], &val);
1903: if (val != -1) continue;
1904: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1905: for (cl = 0; cl < closureSize*2; cl += 2) {
1906: const PetscInt clp = closure[cl];
1908: DMLabelGetValue(label, clp, &val);
1909: if (val == -1) continue;
1910: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1911: break;
1912: }
1913: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1914: }
1915: again = 1;
1916: break;
1917: }
1918: }
1919: }
1920: }
1921: /* Classify the rest by cell membership */
1922: for (s = 0; s < starSize*2; s += 2) {
1923: const PetscInt point = star[s];
1925: DMLabelGetValue(label, point, &val);
1926: if (val == -1) {
1927: PetscInt *sstar = NULL;
1928: PetscInt sstarSize, ss;
1929: PetscBool marked = PETSC_FALSE, isHybrid;
1931: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1932: for (ss = 0; ss < sstarSize*2; ss += 2) {
1933: const PetscInt spoint = sstar[ss];
1935: if ((spoint < cStart) || (spoint >= cEnd)) continue;
1936: DMLabelGetValue(label, spoint, &val);
1938: DMLabelGetValue(depthLabel, point, &dep);
1939: if (val > 0) {
1940: DMLabelSetValue(label, point, shift+dep);
1941: } else {
1942: DMLabelSetValue(label, point, -(shift+dep));
1943: }
1944: marked = PETSC_TRUE;
1945: break;
1946: }
1947: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1948: DMPlexCellIsHybrid_Internal(dm, point, &isHybrid);
1950: }
1951: }
1952: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1953: }
1954: ISRestoreIndices(dimIS, &points);
1955: ISDestroy(&dimIS);
1956: /* If any faces touching the fault divide cells on either side, split them
1957: This only occurs without a surface boundary */
1958: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1959: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1960: ISExpand(facePosIS, faceNegIS, &dimIS);
1961: ISDestroy(&facePosIS);
1962: ISDestroy(&faceNegIS);
1963: ISGetLocalSize(dimIS, &numPoints);
1964: ISGetIndices(dimIS, &points);
1965: for (p = 0; p < numPoints; ++p) {
1966: const PetscInt point = points[p];
1967: const PetscInt *support;
1968: PetscInt supportSize, valA, valB;
1970: DMPlexGetSupportSize(dm, point, &supportSize);
1971: if (supportSize != 2) continue;
1972: DMPlexGetSupport(dm, point, &support);
1973: DMLabelGetValue(label, support[0], &valA);
1974: DMLabelGetValue(label, support[1], &valB);
1975: if ((valA == -1) || (valB == -1)) continue;
1976: if (valA*valB > 0) continue;
1977: /* Split the face */
1978: DMLabelGetValue(label, point, &valA);
1979: DMLabelClearValue(label, point, valA);
1980: DMLabelSetValue(label, point, dim-1);
1981: /* Label its closure:
1982: unmarked: label as unsplit
1983: incident: relabel as split
1984: split: do nothing
1985: */
1986: {
1987: PetscInt *closure = NULL;
1988: PetscInt closureSize, cl;
1990: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1991: for (cl = 0; cl < closureSize*2; cl += 2) {
1992: DMLabelGetValue(label, closure[cl], &valA);
1993: if (valA == -1) { /* Mark as unsplit */
1994: DMLabelGetValue(depthLabel, closure[cl], &dep);
1995: DMLabelSetValue(label, closure[cl], shift2+dep);
1996: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1997: DMLabelGetValue(depthLabel, closure[cl], &dep);
1998: DMLabelClearValue(label, closure[cl], valA);
1999: DMLabelSetValue(label, closure[cl], dep);
2000: }
2001: }
2002: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2003: }
2004: }
2005: ISRestoreIndices(dimIS, &points);
2006: ISDestroy(&dimIS);
2007: return 0;
2008: }
2010: /* Check that no cell have all vertices on the fault */
2011: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
2012: {
2013: IS subpointIS;
2014: const PetscInt *dmpoints;
2015: PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd;
2017: if (!label) return 0;
2018: DMLabelGetDefaultValue(label, &defaultValue);
2019: DMPlexGetSubpointIS(subdm, &subpointIS);
2020: if (!subpointIS) return 0;
2021: DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
2022: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2023: ISGetIndices(subpointIS, &dmpoints);
2024: for (c = cStart; c < cEnd; ++c) {
2025: PetscBool invalidCell = PETSC_TRUE;
2026: PetscInt *closure = NULL;
2027: PetscInt closureSize, cl;
2029: DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2030: for (cl = 0; cl < closureSize*2; cl += 2) {
2031: PetscInt value = 0;
2033: if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
2034: DMLabelGetValue(label, closure[cl], &value);
2035: if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
2036: }
2037: DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2038: if (invalidCell) {
2039: ISRestoreIndices(subpointIS, &dmpoints);
2040: ISDestroy(&subpointIS);
2041: DMDestroy(&subdm);
2042: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
2043: }
2044: }
2045: ISRestoreIndices(subpointIS, &dmpoints);
2046: return 0;
2047: }
2049: /*@
2050: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
2052: Collective on dm
2054: Input Parameters:
2055: + dm - The original DM
2056: . label - The label specifying the interface vertices
2057: - bdlabel - The optional label specifying the interface boundary vertices
2059: Output Parameters:
2060: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
2061: . splitLabel - The label containing the split points, or NULL if no output is desired
2062: . dmInterface - The new interface DM, or NULL
2063: - dmHybrid - The new DM with cohesive cells
2065: Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points
2066: directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
2067: 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
2068: one vertex 3 on the surface would be 6 (101) and 3 (0) in hybridLabel. If an edge 9 from the negative side of the
2069: surface also hits vertex 3, it would be 9 (-101) in hybridLabel.
2071: The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
2072: mesh. The label value is +=100+dim for each point. For example, if two edges 10 and 14 in the hybrid resulting from
2073: splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.
2075: The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using
2076: DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh.
2078: Level: developer
2080: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
2081: @*/
2082: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2083: {
2084: DM idm;
2085: DMLabel subpointMap, hlabel, slabel = NULL;
2086: PetscInt dim;
2094: DMGetDimension(dm, &dim);
2095: DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2096: DMPlexCheckValidSubmesh_Private(dm, label, idm);
2097: DMPlexOrient(idm);
2098: DMPlexGetSubpointMap(idm, &subpointMap);
2099: DMLabelDuplicate(subpointMap, &hlabel);
2100: DMLabelClearStratum(hlabel, dim);
2101: if (splitLabel) {
2102: const char *name;
2103: char sname[PETSC_MAX_PATH_LEN];
2105: PetscObjectGetName((PetscObject) hlabel, &name);
2106: PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2107: PetscStrcat(sname, " split");
2108: DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2109: }
2110: DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2111: if (dmInterface) {*dmInterface = idm;}
2112: else DMDestroy(&idm);
2113: DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2114: DMPlexCopy_Internal(dm, PETSC_TRUE, *dmHybrid);
2115: if (hybridLabel) *hybridLabel = hlabel;
2116: else DMLabelDestroy(&hlabel);
2117: if (splitLabel) *splitLabel = slabel;
2118: {
2119: DM cdm;
2120: DMLabel ctLabel;
2122: /* We need to somehow share the celltype label with the coordinate dm */
2123: DMGetCoordinateDM(*dmHybrid, &cdm);
2124: DMPlexGetCellTypeLabel(*dmHybrid, &ctLabel);
2125: DMSetLabel(cdm, ctLabel);
2126: }
2127: return 0;
2128: }
2130: /* Here we need the explicit assumption that:
2132: For any marked cell, the marked vertices constitute a single face
2133: */
2134: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2135: {
2136: IS subvertexIS = NULL;
2137: const PetscInt *subvertices;
2138: PetscInt *pStart, *pEnd, pSize;
2139: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
2141: *numFaces = 0;
2142: *nFV = 0;
2143: DMPlexGetDepth(dm, &depth);
2144: DMGetDimension(dm, &dim);
2145: pSize = PetscMax(depth, dim) + 1;
2146: PetscMalloc2(pSize, &pStart, pSize, &pEnd);
2147: for (d = 0; d <= depth; ++d) {
2148: DMPlexGetSimplexOrBoxCells(dm, depth-d, &pStart[d], &pEnd[d]);
2149: }
2150: /* Loop over initial vertices and mark all faces in the collective star() */
2151: if (vertexLabel) DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2152: if (subvertexIS) {
2153: ISGetSize(subvertexIS, &numSubVerticesInitial);
2154: ISGetIndices(subvertexIS, &subvertices);
2155: }
2156: for (v = 0; v < numSubVerticesInitial; ++v) {
2157: const PetscInt vertex = subvertices[v];
2158: PetscInt *star = NULL;
2159: PetscInt starSize, s, numCells = 0, c;
2161: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2162: for (s = 0; s < starSize*2; s += 2) {
2163: const PetscInt point = star[s];
2164: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2165: }
2166: for (c = 0; c < numCells; ++c) {
2167: const PetscInt cell = star[c];
2168: PetscInt *closure = NULL;
2169: PetscInt closureSize, cl;
2170: PetscInt cellLoc, numCorners = 0, faceSize = 0;
2172: DMLabelGetValue(subpointMap, cell, &cellLoc);
2173: if (cellLoc == 2) continue;
2175: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2176: for (cl = 0; cl < closureSize*2; cl += 2) {
2177: const PetscInt point = closure[cl];
2178: PetscInt vertexLoc;
2180: if ((point >= pStart[0]) && (point < pEnd[0])) {
2181: ++numCorners;
2182: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2183: if (vertexLoc == value) closure[faceSize++] = point;
2184: }
2185: }
2186: if (!(*nFV)) DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);
2188: if (faceSize == *nFV) {
2189: const PetscInt *cells = NULL;
2190: PetscInt numCells, nc;
2192: ++(*numFaces);
2193: for (cl = 0; cl < faceSize; ++cl) {
2194: DMLabelSetValue(subpointMap, closure[cl], 0);
2195: }
2196: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2197: for (nc = 0; nc < numCells; ++nc) {
2198: DMLabelSetValue(subpointMap, cells[nc], 2);
2199: }
2200: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2201: }
2202: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2203: }
2204: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2205: }
2206: if (subvertexIS) {
2207: ISRestoreIndices(subvertexIS, &subvertices);
2208: }
2209: ISDestroy(&subvertexIS);
2210: PetscFree2(pStart, pEnd);
2211: return 0;
2212: }
2214: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2215: {
2216: IS subvertexIS = NULL;
2217: const PetscInt *subvertices;
2218: PetscInt *pStart, *pEnd;
2219: PetscInt dim, d, numSubVerticesInitial = 0, v;
2221: DMGetDimension(dm, &dim);
2222: PetscMalloc2(dim+1, &pStart, dim+1, &pEnd);
2223: for (d = 0; d <= dim; ++d) {
2224: DMPlexGetSimplexOrBoxCells(dm, dim-d, &pStart[d], &pEnd[d]);
2225: }
2226: /* Loop over initial vertices and mark all faces in the collective star() */
2227: if (vertexLabel) {
2228: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2229: if (subvertexIS) {
2230: ISGetSize(subvertexIS, &numSubVerticesInitial);
2231: ISGetIndices(subvertexIS, &subvertices);
2232: }
2233: }
2234: for (v = 0; v < numSubVerticesInitial; ++v) {
2235: const PetscInt vertex = subvertices[v];
2236: PetscInt *star = NULL;
2237: PetscInt starSize, s, numFaces = 0, f;
2239: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2240: for (s = 0; s < starSize*2; s += 2) {
2241: const PetscInt point = star[s];
2242: PetscInt faceLoc;
2244: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2245: if (markedFaces) {
2246: DMLabelGetValue(vertexLabel, point, &faceLoc);
2247: if (faceLoc < 0) continue;
2248: }
2249: star[numFaces++] = point;
2250: }
2251: }
2252: for (f = 0; f < numFaces; ++f) {
2253: const PetscInt face = star[f];
2254: PetscInt *closure = NULL;
2255: PetscInt closureSize, c;
2256: PetscInt faceLoc;
2258: DMLabelGetValue(subpointMap, face, &faceLoc);
2259: if (faceLoc == dim-1) continue;
2261: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2262: for (c = 0; c < closureSize*2; c += 2) {
2263: const PetscInt point = closure[c];
2264: PetscInt vertexLoc;
2266: if ((point >= pStart[0]) && (point < pEnd[0])) {
2267: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2268: if (vertexLoc != value) break;
2269: }
2270: }
2271: if (c == closureSize*2) {
2272: const PetscInt *support;
2273: PetscInt supportSize, s;
2275: for (c = 0; c < closureSize*2; c += 2) {
2276: const PetscInt point = closure[c];
2278: for (d = 0; d < dim; ++d) {
2279: if ((point >= pStart[d]) && (point < pEnd[d])) {
2280: DMLabelSetValue(subpointMap, point, d);
2281: break;
2282: }
2283: }
2284: }
2285: DMPlexGetSupportSize(dm, face, &supportSize);
2286: DMPlexGetSupport(dm, face, &support);
2287: for (s = 0; s < supportSize; ++s) {
2288: DMLabelSetValue(subpointMap, support[s], dim);
2289: }
2290: }
2291: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2292: }
2293: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2294: }
2295: if (subvertexIS) ISRestoreIndices(subvertexIS, &subvertices);
2296: ISDestroy(&subvertexIS);
2297: PetscFree2(pStart, pEnd);
2298: return 0;
2299: }
2301: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2302: {
2303: DMLabel label = NULL;
2304: const PetscInt *cone;
2305: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2307: *numFaces = 0;
2308: *nFV = 0;
2309: if (labelname) DMGetLabel(dm, labelname, &label);
2310: *subCells = NULL;
2311: DMGetDimension(dm, &dim);
2312: DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2313: if (cMax < 0) return 0;
2314: if (label) {
2315: for (c = cMax; c < cEnd; ++c) {
2316: PetscInt val;
2318: DMLabelGetValue(label, c, &val);
2319: if (val == value) {
2320: ++(*numFaces);
2321: DMPlexGetConeSize(dm, c, &coneSize);
2322: }
2323: }
2324: } else {
2325: *numFaces = cEnd - cMax;
2326: DMPlexGetConeSize(dm, cMax, &coneSize);
2327: }
2328: PetscMalloc1(*numFaces *2, subCells);
2329: if (!(*numFaces)) return 0;
2330: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2331: for (c = cMax; c < cEnd; ++c) {
2332: const PetscInt *cells;
2333: PetscInt numCells;
2335: if (label) {
2336: PetscInt val;
2338: DMLabelGetValue(label, c, &val);
2339: if (val != value) continue;
2340: }
2341: DMPlexGetCone(dm, c, &cone);
2342: for (p = 0; p < *nFV; ++p) {
2343: DMLabelSetValue(subpointMap, cone[p], 0);
2344: }
2345: /* Negative face */
2346: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2347: /* Not true in parallel
2349: for (p = 0; p < numCells; ++p) {
2350: DMLabelSetValue(subpointMap, cells[p], 2);
2351: (*subCells)[subc++] = cells[p];
2352: }
2353: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2354: /* Positive face is not included */
2355: }
2356: return 0;
2357: }
2359: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2360: {
2361: PetscInt *pStart, *pEnd;
2362: PetscInt dim, cMax, cEnd, c, d;
2364: DMGetDimension(dm, &dim);
2365: DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2366: if (cMax < 0) return 0;
2367: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2368: for (d = 0; d <= dim; ++d) DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2369: for (c = cMax; c < cEnd; ++c) {
2370: const PetscInt *cone;
2371: PetscInt *closure = NULL;
2372: PetscInt fconeSize, coneSize, closureSize, cl, val;
2374: if (label) {
2375: DMLabelGetValue(label, c, &val);
2376: if (val != value) continue;
2377: }
2378: DMPlexGetConeSize(dm, c, &coneSize);
2379: DMPlexGetCone(dm, c, &cone);
2380: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2382: /* Negative face */
2383: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2384: for (cl = 0; cl < closureSize*2; cl += 2) {
2385: const PetscInt point = closure[cl];
2387: for (d = 0; d <= dim; ++d) {
2388: if ((point >= pStart[d]) && (point < pEnd[d])) {
2389: DMLabelSetValue(subpointMap, point, d);
2390: break;
2391: }
2392: }
2393: }
2394: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2395: /* Cells -- positive face is not included */
2396: for (cl = 0; cl < 1; ++cl) {
2397: const PetscInt *support;
2398: PetscInt supportSize, s;
2400: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2402: DMPlexGetSupport(dm, cone[cl], &support);
2403: for (s = 0; s < supportSize; ++s) {
2404: DMLabelSetValue(subpointMap, support[s], dim);
2405: }
2406: }
2407: }
2408: PetscFree2(pStart, pEnd);
2409: return 0;
2410: }
2412: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2413: {
2414: MPI_Comm comm;
2415: PetscBool posOrient = PETSC_FALSE;
2416: const PetscInt debug = 0;
2417: PetscInt cellDim, faceSize, f;
2419: PetscObjectGetComm((PetscObject)dm,&comm);
2420: DMGetDimension(dm, &cellDim);
2421: if (debug) PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);
2423: if (cellDim == 1 && numCorners == 2) {
2424: /* Triangle */
2425: faceSize = numCorners-1;
2426: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2427: } else if (cellDim == 2 && numCorners == 3) {
2428: /* Triangle */
2429: faceSize = numCorners-1;
2430: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2431: } else if (cellDim == 3 && numCorners == 4) {
2432: /* Tetrahedron */
2433: faceSize = numCorners-1;
2434: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2435: } else if (cellDim == 1 && numCorners == 3) {
2436: /* Quadratic line */
2437: faceSize = 1;
2438: posOrient = PETSC_TRUE;
2439: } else if (cellDim == 2 && numCorners == 4) {
2440: /* Quads */
2441: faceSize = 2;
2442: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2443: posOrient = PETSC_TRUE;
2444: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2445: posOrient = PETSC_TRUE;
2446: } else {
2447: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2448: posOrient = PETSC_FALSE;
2449: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2450: }
2451: } else if (cellDim == 2 && numCorners == 6) {
2452: /* Quadratic triangle (I hate this) */
2453: /* Edges are determined by the first 2 vertices (corners of edges) */
2454: const PetscInt faceSizeTri = 3;
2455: PetscInt sortedIndices[3], i, iFace;
2456: PetscBool found = PETSC_FALSE;
2457: PetscInt faceVerticesTriSorted[9] = {
2458: 0, 3, 4, /* bottom */
2459: 1, 4, 5, /* right */
2460: 2, 3, 5, /* left */
2461: };
2462: PetscInt faceVerticesTri[9] = {
2463: 0, 3, 4, /* bottom */
2464: 1, 4, 5, /* right */
2465: 2, 5, 3, /* left */
2466: };
2468: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2469: PetscSortInt(faceSizeTri, sortedIndices);
2470: for (iFace = 0; iFace < 3; ++iFace) {
2471: const PetscInt ii = iFace*faceSizeTri;
2472: PetscInt fVertex, cVertex;
2474: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2475: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2476: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2477: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2478: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2479: faceVertices[fVertex] = origVertices[cVertex];
2480: break;
2481: }
2482: }
2483: }
2484: found = PETSC_TRUE;
2485: break;
2486: }
2487: }
2489: if (posOriented) *posOriented = PETSC_TRUE;
2490: return 0;
2491: } else if (cellDim == 2 && numCorners == 9) {
2492: /* Quadratic quad (I hate this) */
2493: /* Edges are determined by the first 2 vertices (corners of edges) */
2494: const PetscInt faceSizeQuad = 3;
2495: PetscInt sortedIndices[3], i, iFace;
2496: PetscBool found = PETSC_FALSE;
2497: PetscInt faceVerticesQuadSorted[12] = {
2498: 0, 1, 4, /* bottom */
2499: 1, 2, 5, /* right */
2500: 2, 3, 6, /* top */
2501: 0, 3, 7, /* left */
2502: };
2503: PetscInt faceVerticesQuad[12] = {
2504: 0, 1, 4, /* bottom */
2505: 1, 2, 5, /* right */
2506: 2, 3, 6, /* top */
2507: 3, 0, 7, /* left */
2508: };
2510: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2511: PetscSortInt(faceSizeQuad, sortedIndices);
2512: for (iFace = 0; iFace < 4; ++iFace) {
2513: const PetscInt ii = iFace*faceSizeQuad;
2514: PetscInt fVertex, cVertex;
2516: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2517: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2518: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2519: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2520: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2521: faceVertices[fVertex] = origVertices[cVertex];
2522: break;
2523: }
2524: }
2525: }
2526: found = PETSC_TRUE;
2527: break;
2528: }
2529: }
2531: if (posOriented) *posOriented = PETSC_TRUE;
2532: return 0;
2533: } else if (cellDim == 3 && numCorners == 8) {
2534: /* Hexes
2535: A hex is two oriented quads with the normal of the first
2536: pointing up at the second.
2538: 7---6
2539: /| /|
2540: 4---5 |
2541: | 1-|-2
2542: |/ |/
2543: 0---3
2545: Faces are determined by the first 4 vertices (corners of faces) */
2546: const PetscInt faceSizeHex = 4;
2547: PetscInt sortedIndices[4], i, iFace;
2548: PetscBool found = PETSC_FALSE;
2549: PetscInt faceVerticesHexSorted[24] = {
2550: 0, 1, 2, 3, /* bottom */
2551: 4, 5, 6, 7, /* top */
2552: 0, 3, 4, 5, /* front */
2553: 2, 3, 5, 6, /* right */
2554: 1, 2, 6, 7, /* back */
2555: 0, 1, 4, 7, /* left */
2556: };
2557: PetscInt faceVerticesHex[24] = {
2558: 1, 2, 3, 0, /* bottom */
2559: 4, 5, 6, 7, /* top */
2560: 0, 3, 5, 4, /* front */
2561: 3, 2, 6, 5, /* right */
2562: 2, 1, 7, 6, /* back */
2563: 1, 0, 4, 7, /* left */
2564: };
2566: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2567: PetscSortInt(faceSizeHex, sortedIndices);
2568: for (iFace = 0; iFace < 6; ++iFace) {
2569: const PetscInt ii = iFace*faceSizeHex;
2570: PetscInt fVertex, cVertex;
2572: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2573: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2574: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2575: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2576: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2577: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2578: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2579: faceVertices[fVertex] = origVertices[cVertex];
2580: break;
2581: }
2582: }
2583: }
2584: found = PETSC_TRUE;
2585: break;
2586: }
2587: }
2589: if (posOriented) *posOriented = PETSC_TRUE;
2590: return 0;
2591: } else if (cellDim == 3 && numCorners == 10) {
2592: /* Quadratic tet */
2593: /* Faces are determined by the first 3 vertices (corners of faces) */
2594: const PetscInt faceSizeTet = 6;
2595: PetscInt sortedIndices[6], i, iFace;
2596: PetscBool found = PETSC_FALSE;
2597: PetscInt faceVerticesTetSorted[24] = {
2598: 0, 1, 2, 6, 7, 8, /* bottom */
2599: 0, 3, 4, 6, 7, 9, /* front */
2600: 1, 4, 5, 7, 8, 9, /* right */
2601: 2, 3, 5, 6, 8, 9, /* left */
2602: };
2603: PetscInt faceVerticesTet[24] = {
2604: 0, 1, 2, 6, 7, 8, /* bottom */
2605: 0, 4, 3, 6, 7, 9, /* front */
2606: 1, 5, 4, 7, 8, 9, /* right */
2607: 2, 3, 5, 8, 6, 9, /* left */
2608: };
2610: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2611: PetscSortInt(faceSizeTet, sortedIndices);
2612: for (iFace=0; iFace < 4; ++iFace) {
2613: const PetscInt ii = iFace*faceSizeTet;
2614: PetscInt fVertex, cVertex;
2616: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2617: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2618: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2619: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2620: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2621: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2622: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2623: faceVertices[fVertex] = origVertices[cVertex];
2624: break;
2625: }
2626: }
2627: }
2628: found = PETSC_TRUE;
2629: break;
2630: }
2631: }
2633: if (posOriented) *posOriented = PETSC_TRUE;
2634: return 0;
2635: } else if (cellDim == 3 && numCorners == 27) {
2636: /* Quadratic hexes (I hate this)
2637: A hex is two oriented quads with the normal of the first
2638: pointing up at the second.
2640: 7---6
2641: /| /|
2642: 4---5 |
2643: | 3-|-2
2644: |/ |/
2645: 0---1
2647: Faces are determined by the first 4 vertices (corners of faces) */
2648: const PetscInt faceSizeQuadHex = 9;
2649: PetscInt sortedIndices[9], i, iFace;
2650: PetscBool found = PETSC_FALSE;
2651: PetscInt faceVerticesQuadHexSorted[54] = {
2652: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2653: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2654: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2655: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2656: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2657: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2658: };
2659: PetscInt faceVerticesQuadHex[54] = {
2660: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2661: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2662: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2663: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2664: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2665: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2666: };
2668: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2669: PetscSortInt(faceSizeQuadHex, sortedIndices);
2670: for (iFace = 0; iFace < 6; ++iFace) {
2671: const PetscInt ii = iFace*faceSizeQuadHex;
2672: PetscInt fVertex, cVertex;
2674: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2675: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2676: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2677: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2678: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2679: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2680: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2681: faceVertices[fVertex] = origVertices[cVertex];
2682: break;
2683: }
2684: }
2685: }
2686: found = PETSC_TRUE;
2687: break;
2688: }
2689: }
2691: if (posOriented) *posOriented = PETSC_TRUE;
2692: return 0;
2693: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2694: if (!posOrient) {
2695: if (debug) PetscPrintf(comm, " Reversing initial face orientation\n");
2696: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2697: } else {
2698: if (debug) PetscPrintf(comm, " Keeping initial face orientation\n");
2699: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2700: }
2701: if (posOriented) *posOriented = posOrient;
2702: return 0;
2703: }
2705: /*@
2706: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2707: in faceVertices. The orientation is such that the face normal points out of the cell
2709: Not collective
2711: Input Parameters:
2712: + dm - The original mesh
2713: . cell - The cell mesh point
2714: . faceSize - The number of vertices on the face
2715: . face - The face vertices
2716: . numCorners - The number of vertices on the cell
2717: . indices - Local numbering of face vertices in cell cone
2718: - origVertices - Original face vertices
2720: Output Parameters:
2721: + faceVertices - The face vertices properly oriented
2722: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2724: Level: developer
2726: .seealso: DMPlexGetCone()
2727: @*/
2728: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2729: {
2730: const PetscInt *cone = NULL;
2731: PetscInt coneSize, v, f, v2;
2732: PetscInt oppositeVertex = -1;
2734: DMPlexGetConeSize(dm, cell, &coneSize);
2735: DMPlexGetCone(dm, cell, &cone);
2736: for (v = 0, v2 = 0; v < coneSize; ++v) {
2737: PetscBool found = PETSC_FALSE;
2739: for (f = 0; f < faceSize; ++f) {
2740: if (face[f] == cone[v]) {
2741: found = PETSC_TRUE; break;
2742: }
2743: }
2744: if (found) {
2745: indices[v2] = v;
2746: origVertices[v2] = cone[v];
2747: ++v2;
2748: } else {
2749: oppositeVertex = v;
2750: }
2751: }
2752: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2753: return 0;
2754: }
2756: /*
2757: DMPlexInsertFace_Internal - Puts a face into the mesh
2759: Not collective
2761: Input Parameters:
2762: + dm - The DMPlex
2763: . numFaceVertex - The number of vertices in the face
2764: . faceVertices - The vertices in the face for dm
2765: . subfaceVertices - The vertices in the face for subdm
2766: . numCorners - The number of vertices in the cell
2767: . cell - A cell in dm containing the face
2768: . subcell - A cell in subdm containing the face
2769: . firstFace - First face in the mesh
2770: - newFacePoint - Next face in the mesh
2772: Output Parameters:
2773: . newFacePoint - Contains next face point number on input, updated on output
2775: Level: developer
2776: */
2777: static PetscErrorCode DMPlexInsertFace_Internal(DM dm, DM subdm, PetscInt numFaceVertices, const PetscInt faceVertices[], const PetscInt subfaceVertices[], PetscInt numCorners, PetscInt cell, PetscInt subcell, PetscInt firstFace, PetscInt *newFacePoint)
2778: {
2779: MPI_Comm comm;
2780: DM_Plex *submesh = (DM_Plex*) subdm->data;
2781: const PetscInt *faces;
2782: PetscInt numFaces, coneSize;
2784: PetscObjectGetComm((PetscObject)dm,&comm);
2785: DMPlexGetConeSize(subdm, subcell, &coneSize);
2787: #if 0
2788: /* Cannot use this because support() has not been constructed yet */
2789: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2790: #else
2791: {
2792: PetscInt f;
2794: numFaces = 0;
2795: DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2796: for (f = firstFace; f < *newFacePoint; ++f) {
2797: PetscInt dof, off, d;
2799: PetscSectionGetDof(submesh->coneSection, f, &dof);
2800: PetscSectionGetOffset(submesh->coneSection, f, &off);
2801: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2802: for (d = 0; d < dof; ++d) {
2803: const PetscInt p = submesh->cones[off+d];
2804: PetscInt v;
2806: for (v = 0; v < numFaceVertices; ++v) {
2807: if (subfaceVertices[v] == p) break;
2808: }
2809: if (v == numFaceVertices) break;
2810: }
2811: if (d == dof) {
2812: numFaces = 1;
2813: ((PetscInt*) faces)[0] = f;
2814: }
2815: }
2816: }
2817: #endif
2819: else if (numFaces == 1) {
2820: /* Add the other cell neighbor for this face */
2821: DMPlexSetCone(subdm, subcell, faces);
2822: } else {
2823: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2824: PetscBool posOriented;
2826: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2827: origVertices = &orientedVertices[numFaceVertices];
2828: indices = &orientedVertices[numFaceVertices*2];
2829: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2830: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2831: /* TODO: I know that routine should return a permutation, not the indices */
2832: for (v = 0; v < numFaceVertices; ++v) {
2833: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2834: for (ov = 0; ov < numFaceVertices; ++ov) {
2835: if (orientedVertices[ov] == vertex) {
2836: orientedSubVertices[ov] = subvertex;
2837: break;
2838: }
2839: }
2841: }
2842: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2843: DMPlexSetCone(subdm, subcell, newFacePoint);
2844: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2845: ++(*newFacePoint);
2846: }
2847: #if 0
2848: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2849: #else
2850: DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2851: #endif
2852: return 0;
2853: }
2855: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2856: {
2857: MPI_Comm comm;
2858: DMLabel subpointMap;
2859: IS subvertexIS, subcellIS;
2860: const PetscInt *subVertices, *subCells;
2861: PetscInt numSubVertices, firstSubVertex, numSubCells;
2862: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2863: PetscInt vStart, vEnd, c, f;
2865: PetscObjectGetComm((PetscObject)dm,&comm);
2866: /* Create subpointMap which marks the submesh */
2867: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2868: DMPlexSetSubpointMap(subdm, subpointMap);
2869: DMLabelDestroy(&subpointMap);
2870: if (vertexLabel) DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);
2871: /* Setup chart */
2872: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2873: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2874: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2875: DMPlexSetVTKCellHeight(subdm, 1);
2876: /* Set cone sizes */
2877: firstSubVertex = numSubCells;
2878: firstSubFace = numSubCells+numSubVertices;
2879: newFacePoint = firstSubFace;
2880: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2881: if (subvertexIS) ISGetIndices(subvertexIS, &subVertices);
2882: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2883: if (subcellIS) ISGetIndices(subcellIS, &subCells);
2884: for (c = 0; c < numSubCells; ++c) {
2885: DMPlexSetConeSize(subdm, c, 1);
2886: }
2887: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2888: DMPlexSetConeSize(subdm, f, nFV);
2889: }
2890: DMSetUp(subdm);
2891: /* Create face cones */
2892: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2893: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2894: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2895: for (c = 0; c < numSubCells; ++c) {
2896: const PetscInt cell = subCells[c];
2897: const PetscInt subcell = c;
2898: PetscInt *closure = NULL;
2899: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2901: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2902: for (cl = 0; cl < closureSize*2; cl += 2) {
2903: const PetscInt point = closure[cl];
2904: PetscInt subVertex;
2906: if ((point >= vStart) && (point < vEnd)) {
2907: ++numCorners;
2908: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2909: if (subVertex >= 0) {
2910: closure[faceSize] = point;
2911: subface[faceSize] = firstSubVertex+subVertex;
2912: ++faceSize;
2913: }
2914: }
2915: }
2917: if (faceSize == nFV) {
2918: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2919: }
2920: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2921: }
2922: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2923: DMPlexSymmetrize(subdm);
2924: DMPlexStratify(subdm);
2925: /* Build coordinates */
2926: {
2927: PetscSection coordSection, subCoordSection;
2928: Vec coordinates, subCoordinates;
2929: PetscScalar *coords, *subCoords;
2930: PetscInt numComp, coordSize, v;
2931: const char *name;
2933: DMGetCoordinateSection(dm, &coordSection);
2934: DMGetCoordinatesLocal(dm, &coordinates);
2935: DMGetCoordinateSection(subdm, &subCoordSection);
2936: PetscSectionSetNumFields(subCoordSection, 1);
2937: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2938: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2939: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2940: for (v = 0; v < numSubVertices; ++v) {
2941: const PetscInt vertex = subVertices[v];
2942: const PetscInt subvertex = firstSubVertex+v;
2943: PetscInt dof;
2945: PetscSectionGetDof(coordSection, vertex, &dof);
2946: PetscSectionSetDof(subCoordSection, subvertex, dof);
2947: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2948: }
2949: PetscSectionSetUp(subCoordSection);
2950: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2951: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2952: PetscObjectGetName((PetscObject)coordinates,&name);
2953: PetscObjectSetName((PetscObject)subCoordinates,name);
2954: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2955: VecSetType(subCoordinates,VECSTANDARD);
2956: if (coordSize) {
2957: VecGetArray(coordinates, &coords);
2958: VecGetArray(subCoordinates, &subCoords);
2959: for (v = 0; v < numSubVertices; ++v) {
2960: const PetscInt vertex = subVertices[v];
2961: const PetscInt subvertex = firstSubVertex+v;
2962: PetscInt dof, off, sdof, soff, d;
2964: PetscSectionGetDof(coordSection, vertex, &dof);
2965: PetscSectionGetOffset(coordSection, vertex, &off);
2966: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2967: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2969: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2970: }
2971: VecRestoreArray(coordinates, &coords);
2972: VecRestoreArray(subCoordinates, &subCoords);
2973: }
2974: DMSetCoordinatesLocal(subdm, subCoordinates);
2975: VecDestroy(&subCoordinates);
2976: }
2977: /* Cleanup */
2978: if (subvertexIS) ISRestoreIndices(subvertexIS, &subVertices);
2979: ISDestroy(&subvertexIS);
2980: if (subcellIS) ISRestoreIndices(subcellIS, &subCells);
2981: ISDestroy(&subcellIS);
2982: return 0;
2983: }
2985: static inline PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2986: {
2987: PetscInt subPoint;
2990: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2991: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2992: }
2994: static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPoints[], const PetscInt *subpoints[], const PetscInt firstSubPoint[], DM subdm)
2995: {
2996: PetscInt Nl, l, d;
2998: DMGetNumLabels(dm, &Nl);
2999: for (l = 0; l < Nl; ++l) {
3000: DMLabel label, newlabel;
3001: const char *lname;
3002: PetscBool isDepth, isDim, isCelltype, isVTK;
3003: IS valueIS;
3004: const PetscInt *values;
3005: PetscInt Nv, v;
3007: DMGetLabelName(dm, l, &lname);
3008: PetscStrcmp(lname, "depth", &isDepth);
3009: PetscStrcmp(lname, "dim", &isDim);
3010: PetscStrcmp(lname, "celltype", &isCelltype);
3011: PetscStrcmp(lname, "vtk", &isVTK);
3012: if (isDepth || isDim || isCelltype || isVTK) continue;
3013: DMCreateLabel(subdm, lname);
3014: DMGetLabel(dm, lname, &label);
3015: DMGetLabel(subdm, lname, &newlabel);
3016: DMLabelGetDefaultValue(label, &v);
3017: DMLabelSetDefaultValue(newlabel, v);
3018: DMLabelGetValueIS(label, &valueIS);
3019: ISGetLocalSize(valueIS, &Nv);
3020: ISGetIndices(valueIS, &values);
3021: for (v = 0; v < Nv; ++v) {
3022: IS pointIS;
3023: const PetscInt *points;
3024: PetscInt Np, p;
3026: DMLabelGetStratumIS(label, values[v], &pointIS);
3027: ISGetLocalSize(pointIS, &Np);
3028: ISGetIndices(pointIS, &points);
3029: for (p = 0; p < Np; ++p) {
3030: const PetscInt point = points[p];
3031: PetscInt subp;
3033: DMPlexGetPointDepth(dm, point, &d);
3034: subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]);
3035: if (subp >= 0) DMLabelSetValue(newlabel, subp, values[v]);
3036: }
3037: ISRestoreIndices(pointIS, &points);
3038: ISDestroy(&pointIS);
3039: }
3040: ISRestoreIndices(valueIS, &values);
3041: ISDestroy(&valueIS);
3042: }
3043: return 0;
3044: }
3046: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
3047: {
3048: MPI_Comm comm;
3049: DMLabel subpointMap;
3050: IS *subpointIS;
3051: const PetscInt **subpoints;
3052: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3053: PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v;
3054: PetscMPIInt rank;
3056: PetscObjectGetComm((PetscObject)dm,&comm);
3057: MPI_Comm_rank(comm, &rank);
3058: /* Create subpointMap which marks the submesh */
3059: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3060: DMPlexSetSubpointMap(subdm, subpointMap);
3061: if (cellHeight) {
3062: if (isCohesive) DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);
3063: else DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);
3064: } else {
3065: DMLabel depth;
3066: IS pointIS;
3067: const PetscInt *points;
3068: PetscInt numPoints=0;
3070: DMPlexGetDepthLabel(dm, &depth);
3071: DMLabelGetStratumIS(label, value, &pointIS);
3072: if (pointIS) {
3073: ISGetIndices(pointIS, &points);
3074: ISGetLocalSize(pointIS, &numPoints);
3075: }
3076: for (p = 0; p < numPoints; ++p) {
3077: PetscInt *closure = NULL;
3078: PetscInt closureSize, c, pdim;
3080: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3081: for (c = 0; c < closureSize*2; c += 2) {
3082: DMLabelGetValue(depth, closure[c], &pdim);
3083: DMLabelSetValue(subpointMap, closure[c], pdim);
3084: }
3085: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3086: }
3087: if (pointIS) ISRestoreIndices(pointIS, &points);
3088: ISDestroy(&pointIS);
3089: }
3090: /* Setup chart */
3091: DMGetDimension(dm, &dim);
3092: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
3093: for (d = 0; d <= dim; ++d) {
3094: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3095: totSubPoints += numSubPoints[d];
3096: }
3097: DMPlexSetChart(subdm, 0, totSubPoints);
3098: DMPlexSetVTKCellHeight(subdm, cellHeight);
3099: /* Set cone sizes */
3100: firstSubPoint[dim] = 0;
3101: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
3102: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
3103: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
3104: for (d = 0; d <= dim; ++d) {
3105: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3106: if (subpointIS[d]) ISGetIndices(subpointIS[d], &subpoints[d]);
3107: }
3108: /* We do not want this label automatically computed, instead we compute it here */
3109: DMCreateLabel(subdm, "celltype");
3110: for (d = 0; d <= dim; ++d) {
3111: for (p = 0; p < numSubPoints[d]; ++p) {
3112: const PetscInt point = subpoints[d][p];
3113: const PetscInt subpoint = firstSubPoint[d] + p;
3114: const PetscInt *cone;
3115: PetscInt coneSize;
3117: DMPlexGetConeSize(dm, point, &coneSize);
3118: if (cellHeight && (d == dim)) {
3119: PetscInt coneSizeNew, c, val;
3121: DMPlexGetCone(dm, point, &cone);
3122: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3123: DMLabelGetValue(subpointMap, cone[c], &val);
3124: if (val >= 0) coneSizeNew++;
3125: }
3126: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3127: DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST);
3128: } else {
3129: DMPolytopeType ct;
3131: DMPlexSetConeSize(subdm, subpoint, coneSize);
3132: DMPlexGetCellType(dm, point, &ct);
3133: DMPlexSetCellType(subdm, subpoint, ct);
3134: }
3135: }
3136: }
3137: DMLabelDestroy(&subpointMap);
3138: DMSetUp(subdm);
3139: /* Set cones */
3140: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3141: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3142: for (d = 0; d <= dim; ++d) {
3143: for (p = 0; p < numSubPoints[d]; ++p) {
3144: const PetscInt point = subpoints[d][p];
3145: const PetscInt subpoint = firstSubPoint[d] + p;
3146: const PetscInt *cone, *ornt;
3147: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
3149: if (d == dim-1) {
3150: const PetscInt *support, *cone, *ornt;
3151: PetscInt supportSize, coneSize, s, subc;
3153: DMPlexGetSupport(dm, point, &support);
3154: DMPlexGetSupportSize(dm, point, &supportSize);
3155: for (s = 0; s < supportSize; ++s) {
3156: PetscBool isHybrid;
3158: DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);
3159: if (!isHybrid) continue;
3160: PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3161: if (subc >= 0) {
3162: const PetscInt ccell = subpoints[d+1][subc];
3164: DMPlexGetCone(dm, ccell, &cone);
3165: DMPlexGetConeSize(dm, ccell, &coneSize);
3166: DMPlexGetConeOrientation(dm, ccell, &ornt);
3167: for (c = 0; c < coneSize; ++c) {
3168: if (cone[c] == point) {
3169: fornt = ornt[c];
3170: break;
3171: }
3172: }
3173: break;
3174: }
3175: }
3176: }
3177: DMPlexGetConeSize(dm, point, &coneSize);
3178: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3179: DMPlexGetCone(dm, point, &cone);
3180: DMPlexGetConeOrientation(dm, point, &ornt);
3181: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3182: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3183: if (subc >= 0) {
3184: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3185: orntNew[coneSizeNew] = ornt[c];
3186: ++coneSizeNew;
3187: }
3188: }
3190: DMPlexSetCone(subdm, subpoint, coneNew);
3191: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3192: if (fornt < 0) DMPlexOrientPoint(subdm, subpoint, fornt);
3193: }
3194: }
3195: PetscFree2(coneNew,orntNew);
3196: DMPlexSymmetrize(subdm);
3197: DMPlexStratify(subdm);
3198: /* Build coordinates */
3199: {
3200: PetscSection coordSection, subCoordSection;
3201: Vec coordinates, subCoordinates;
3202: PetscScalar *coords, *subCoords;
3203: PetscInt cdim, numComp, coordSize;
3204: const char *name;
3206: DMGetCoordinateDim(dm, &cdim);
3207: DMGetCoordinateSection(dm, &coordSection);
3208: DMGetCoordinatesLocal(dm, &coordinates);
3209: DMGetCoordinateSection(subdm, &subCoordSection);
3210: PetscSectionSetNumFields(subCoordSection, 1);
3211: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3212: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3213: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3214: for (v = 0; v < numSubPoints[0]; ++v) {
3215: const PetscInt vertex = subpoints[0][v];
3216: const PetscInt subvertex = firstSubPoint[0]+v;
3217: PetscInt dof;
3219: PetscSectionGetDof(coordSection, vertex, &dof);
3220: PetscSectionSetDof(subCoordSection, subvertex, dof);
3221: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3222: }
3223: PetscSectionSetUp(subCoordSection);
3224: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3225: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3226: PetscObjectGetName((PetscObject)coordinates,&name);
3227: PetscObjectSetName((PetscObject)subCoordinates,name);
3228: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3229: VecSetBlockSize(subCoordinates, cdim);
3230: VecSetType(subCoordinates,VECSTANDARD);
3231: VecGetArray(coordinates, &coords);
3232: VecGetArray(subCoordinates, &subCoords);
3233: for (v = 0; v < numSubPoints[0]; ++v) {
3234: const PetscInt vertex = subpoints[0][v];
3235: const PetscInt subvertex = firstSubPoint[0]+v;
3236: PetscInt dof, off, sdof, soff, d;
3238: PetscSectionGetDof(coordSection, vertex, &dof);
3239: PetscSectionGetOffset(coordSection, vertex, &off);
3240: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3241: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3243: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3244: }
3245: VecRestoreArray(coordinates, &coords);
3246: VecRestoreArray(subCoordinates, &subCoords);
3247: DMSetCoordinatesLocal(subdm, subCoordinates);
3248: VecDestroy(&subCoordinates);
3249: }
3250: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3251: {
3252: PetscSF sfPoint, sfPointSub;
3253: IS subpIS;
3254: const PetscSFNode *remotePoints;
3255: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3256: const PetscInt *localPoints, *subpoints;
3257: PetscInt *slocalPoints;
3258: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3259: PetscMPIInt rank;
3261: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3262: DMGetPointSF(dm, &sfPoint);
3263: DMGetPointSF(subdm, &sfPointSub);
3264: DMPlexGetChart(dm, &pStart, &pEnd);
3265: DMPlexGetChart(subdm, NULL, &numSubroots);
3266: DMPlexGetSubpointIS(subdm, &subpIS);
3267: if (subpIS) {
3268: ISGetIndices(subpIS, &subpoints);
3269: ISGetLocalSize(subpIS, &numSubpoints);
3270: }
3271: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3272: if (numRoots >= 0) {
3273: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3274: for (p = 0; p < pEnd-pStart; ++p) {
3275: newLocalPoints[p].rank = -2;
3276: newLocalPoints[p].index = -2;
3277: }
3278: /* Set subleaves */
3279: for (l = 0; l < numLeaves; ++l) {
3280: const PetscInt point = localPoints[l];
3281: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3283: if (subpoint < 0) continue;
3284: newLocalPoints[point-pStart].rank = rank;
3285: newLocalPoints[point-pStart].index = subpoint;
3286: ++numSubleaves;
3287: }
3288: /* Must put in owned subpoints */
3289: for (p = pStart; p < pEnd; ++p) {
3290: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3292: if (subpoint < 0) {
3293: newOwners[p-pStart].rank = -3;
3294: newOwners[p-pStart].index = -3;
3295: } else {
3296: newOwners[p-pStart].rank = rank;
3297: newOwners[p-pStart].index = subpoint;
3298: }
3299: }
3300: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3301: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3302: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
3303: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
3304: PetscMalloc1(numSubleaves, &slocalPoints);
3305: PetscMalloc1(numSubleaves, &sremotePoints);
3306: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3307: const PetscInt point = localPoints[l];
3308: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3310: if (subpoint < 0) continue;
3311: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3312: slocalPoints[sl] = subpoint;
3313: sremotePoints[sl].rank = newLocalPoints[point].rank;
3314: sremotePoints[sl].index = newLocalPoints[point].index;
3317: ++sl;
3318: }
3320: PetscFree2(newLocalPoints,newOwners);
3321: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3322: }
3323: if (subpIS) {
3324: ISRestoreIndices(subpIS, &subpoints);
3325: }
3326: }
3327: /* Filter labels */
3328: DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm);
3329: /* Cleanup */
3330: for (d = 0; d <= dim; ++d) {
3331: if (subpointIS[d]) ISRestoreIndices(subpointIS[d], &subpoints[d]);
3332: ISDestroy(&subpointIS[d]);
3333: }
3334: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3335: return 0;
3336: }
3338: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3339: {
3340: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3341: return 0;
3342: }
3344: /*@
3345: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3347: Input Parameters:
3348: + dm - The original mesh
3349: . vertexLabel - The DMLabel marking points contained in the surface
3350: . value - The label value to use
3351: - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked
3353: Output Parameter:
3354: . subdm - The surface mesh
3356: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3358: Level: developer
3360: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3361: @*/
3362: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3363: {
3364: DMPlexInterpolatedFlag interpolated;
3365: PetscInt dim, cdim;
3369: DMGetDimension(dm, &dim);
3370: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3371: DMSetType(*subdm, DMPLEX);
3372: DMSetDimension(*subdm, dim-1);
3373: DMGetCoordinateDim(dm, &cdim);
3374: DMSetCoordinateDim(*subdm, cdim);
3375: DMPlexIsInterpolated(dm, &interpolated);
3377: if (interpolated) {
3378: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3379: } else {
3380: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3381: }
3382: DMPlexCopy_Internal(dm, PETSC_TRUE, *subdm);
3383: return 0;
3384: }
3386: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3387: {
3388: MPI_Comm comm;
3389: DMLabel subpointMap;
3390: IS subvertexIS;
3391: const PetscInt *subVertices;
3392: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3393: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3394: PetscInt c, f;
3396: PetscObjectGetComm((PetscObject)dm, &comm);
3397: /* Create subpointMap which marks the submesh */
3398: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3399: DMPlexSetSubpointMap(subdm, subpointMap);
3400: DMLabelDestroy(&subpointMap);
3401: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3402: /* Setup chart */
3403: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3404: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3405: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3406: DMPlexSetVTKCellHeight(subdm, 1);
3407: /* Set cone sizes */
3408: firstSubVertex = numSubCells;
3409: firstSubFace = numSubCells+numSubVertices;
3410: newFacePoint = firstSubFace;
3411: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3412: if (subvertexIS) ISGetIndices(subvertexIS, &subVertices);
3413: for (c = 0; c < numSubCells; ++c) {
3414: DMPlexSetConeSize(subdm, c, 1);
3415: }
3416: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3417: DMPlexSetConeSize(subdm, f, nFV);
3418: }
3419: DMSetUp(subdm);
3420: /* Create face cones */
3421: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3422: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3423: for (c = 0; c < numSubCells; ++c) {
3424: const PetscInt cell = subCells[c];
3425: const PetscInt subcell = c;
3426: const PetscInt *cone, *cells;
3427: PetscBool isHybrid;
3428: PetscInt numCells, subVertex, p, v;
3430: DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid);
3431: if (!isHybrid) continue;
3432: DMPlexGetCone(dm, cell, &cone);
3433: for (v = 0; v < nFV; ++v) {
3434: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3435: subface[v] = firstSubVertex+subVertex;
3436: }
3437: DMPlexSetCone(subdm, newFacePoint, subface);
3438: DMPlexSetCone(subdm, subcell, &newFacePoint);
3439: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3440: /* Not true in parallel
3442: for (p = 0; p < numCells; ++p) {
3443: PetscInt negsubcell;
3444: PetscBool isHybrid;
3446: DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid);
3447: if (isHybrid) continue;
3448: /* I know this is a crap search */
3449: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3450: if (subCells[negsubcell] == cells[p]) break;
3451: }
3453: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3454: }
3455: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3456: ++newFacePoint;
3457: }
3458: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3459: DMPlexSymmetrize(subdm);
3460: DMPlexStratify(subdm);
3461: /* Build coordinates */
3462: {
3463: PetscSection coordSection, subCoordSection;
3464: Vec coordinates, subCoordinates;
3465: PetscScalar *coords, *subCoords;
3466: PetscInt cdim, numComp, coordSize, v;
3467: const char *name;
3469: DMGetCoordinateDim(dm, &cdim);
3470: DMGetCoordinateSection(dm, &coordSection);
3471: DMGetCoordinatesLocal(dm, &coordinates);
3472: DMGetCoordinateSection(subdm, &subCoordSection);
3473: PetscSectionSetNumFields(subCoordSection, 1);
3474: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3475: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3476: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3477: for (v = 0; v < numSubVertices; ++v) {
3478: const PetscInt vertex = subVertices[v];
3479: const PetscInt subvertex = firstSubVertex+v;
3480: PetscInt dof;
3482: PetscSectionGetDof(coordSection, vertex, &dof);
3483: PetscSectionSetDof(subCoordSection, subvertex, dof);
3484: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3485: }
3486: PetscSectionSetUp(subCoordSection);
3487: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3488: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3489: PetscObjectGetName((PetscObject)coordinates,&name);
3490: PetscObjectSetName((PetscObject)subCoordinates,name);
3491: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3492: VecSetBlockSize(subCoordinates, cdim);
3493: VecSetType(subCoordinates,VECSTANDARD);
3494: VecGetArray(coordinates, &coords);
3495: VecGetArray(subCoordinates, &subCoords);
3496: for (v = 0; v < numSubVertices; ++v) {
3497: const PetscInt vertex = subVertices[v];
3498: const PetscInt subvertex = firstSubVertex+v;
3499: PetscInt dof, off, sdof, soff, d;
3501: PetscSectionGetDof(coordSection, vertex, &dof);
3502: PetscSectionGetOffset(coordSection, vertex, &off);
3503: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3504: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3506: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3507: }
3508: VecRestoreArray(coordinates, &coords);
3509: VecRestoreArray(subCoordinates, &subCoords);
3510: DMSetCoordinatesLocal(subdm, subCoordinates);
3511: VecDestroy(&subCoordinates);
3512: }
3513: /* Build SF */
3514: CHKMEMQ;
3515: {
3516: PetscSF sfPoint, sfPointSub;
3517: const PetscSFNode *remotePoints;
3518: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3519: const PetscInt *localPoints;
3520: PetscInt *slocalPoints;
3521: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3522: PetscMPIInt rank;
3524: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3525: DMGetPointSF(dm, &sfPoint);
3526: DMGetPointSF(subdm, &sfPointSub);
3527: DMPlexGetChart(dm, &pStart, &pEnd);
3528: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3529: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3530: if (numRoots >= 0) {
3531: /* Only vertices should be shared */
3532: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3533: for (p = 0; p < pEnd-pStart; ++p) {
3534: newLocalPoints[p].rank = -2;
3535: newLocalPoints[p].index = -2;
3536: }
3537: /* Set subleaves */
3538: for (l = 0; l < numLeaves; ++l) {
3539: const PetscInt point = localPoints[l];
3540: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3543: if (subPoint < 0) continue;
3544: newLocalPoints[point-pStart].rank = rank;
3545: newLocalPoints[point-pStart].index = subPoint;
3546: ++numSubLeaves;
3547: }
3548: /* Must put in owned subpoints */
3549: for (p = pStart; p < pEnd; ++p) {
3550: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3552: if (subPoint < 0) {
3553: newOwners[p-pStart].rank = -3;
3554: newOwners[p-pStart].index = -3;
3555: } else {
3556: newOwners[p-pStart].rank = rank;
3557: newOwners[p-pStart].index = subPoint;
3558: }
3559: }
3560: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3561: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3562: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
3563: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
3564: PetscMalloc1(numSubLeaves, &slocalPoints);
3565: PetscMalloc1(numSubLeaves, &sremotePoints);
3566: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3567: const PetscInt point = localPoints[l];
3568: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3570: if (subPoint < 0) continue;
3571: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3572: slocalPoints[sl] = subPoint;
3573: sremotePoints[sl].rank = newLocalPoints[point].rank;
3574: sremotePoints[sl].index = newLocalPoints[point].index;
3577: ++sl;
3578: }
3579: PetscFree2(newLocalPoints,newOwners);
3581: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3582: }
3583: }
3584: CHKMEMQ;
3585: /* Cleanup */
3586: if (subvertexIS) ISRestoreIndices(subvertexIS, &subVertices);
3587: ISDestroy(&subvertexIS);
3588: PetscFree(subCells);
3589: return 0;
3590: }
3592: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3593: {
3594: DMLabel label = NULL;
3596: if (labelname) DMGetLabel(dm, labelname, &label);
3597: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3598: return 0;
3599: }
3601: /*@C
3602: DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a Label can be given to restrict the cells.
3604: Input Parameters:
3605: + dm - The original mesh
3606: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3607: . label - A label name, or NULL
3608: - value - A label value
3610: Output Parameter:
3611: . subdm - The surface mesh
3613: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3615: Level: developer
3617: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3618: @*/
3619: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3620: {
3621: PetscInt dim, cdim, depth;
3625: DMGetDimension(dm, &dim);
3626: DMPlexGetDepth(dm, &depth);
3627: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3628: DMSetType(*subdm, DMPLEX);
3629: DMSetDimension(*subdm, dim-1);
3630: DMGetCoordinateDim(dm, &cdim);
3631: DMSetCoordinateDim(*subdm, cdim);
3632: if (depth == dim) {
3633: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3634: } else {
3635: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3636: }
3637: DMPlexCopy_Internal(dm, PETSC_TRUE, *subdm);
3638: return 0;
3639: }
3641: /*@
3642: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3644: Input Parameters:
3645: + dm - The original mesh
3646: . cellLabel - The DMLabel marking cells contained in the new mesh
3647: - value - The label value to use
3649: Output Parameter:
3650: . subdm - The new mesh
3652: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3654: Level: developer
3656: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3657: @*/
3658: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3659: {
3660: PetscInt dim;
3664: DMGetDimension(dm, &dim);
3665: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3666: DMSetType(*subdm, DMPLEX);
3667: DMSetDimension(*subdm, dim);
3668: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3669: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3670: DMPlexCopy_Internal(dm, PETSC_TRUE, *subdm);
3671: return 0;
3672: }
3674: /*@
3675: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3677: Input Parameter:
3678: . dm - The submesh DM
3680: Output Parameter:
3681: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3683: Level: developer
3685: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS()
3686: @*/
3687: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3688: {
3691: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3692: return 0;
3693: }
3695: /*@
3696: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3698: Input Parameters:
3699: + dm - The submesh DM
3700: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3702: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3704: Level: developer
3706: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS()
3707: @*/
3708: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3709: {
3710: DM_Plex *mesh = (DM_Plex *) dm->data;
3711: DMLabel tmp;
3714: tmp = mesh->subpointMap;
3715: mesh->subpointMap = subpointMap;
3716: PetscObjectReference((PetscObject) mesh->subpointMap);
3717: DMLabelDestroy(&tmp);
3718: return 0;
3719: }
3721: static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS)
3722: {
3723: DMLabel spmap;
3724: PetscInt depth, d;
3726: DMPlexGetSubpointMap(dm, &spmap);
3727: DMPlexGetDepth(dm, &depth);
3728: if (spmap && depth >= 0) {
3729: DM_Plex *mesh = (DM_Plex *) dm->data;
3730: PetscInt *points, *depths;
3731: PetscInt pStart, pEnd, p, off;
3733: DMPlexGetChart(dm, &pStart, &pEnd);
3735: PetscMalloc1(pEnd, &points);
3736: DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3737: depths[0] = depth;
3738: depths[1] = 0;
3739: for (d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3740: for (d = 0, off = 0; d <= depth; ++d) {
3741: const PetscInt dep = depths[d];
3742: PetscInt depStart, depEnd, n;
3744: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3745: DMLabelGetStratumSize(spmap, dep, &n);
3746: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3748: } else {
3749: if (!n) {
3750: if (d == 0) {
3751: /* Missing cells */
3752: for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3753: } else {
3754: /* Missing faces */
3755: for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3756: }
3757: }
3758: }
3759: if (n) {
3760: IS is;
3761: const PetscInt *opoints;
3763: DMLabelGetStratumIS(spmap, dep, &is);
3764: ISGetIndices(is, &opoints);
3765: for (p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3766: ISRestoreIndices(is, &opoints);
3767: ISDestroy(&is);
3768: }
3769: }
3770: DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3772: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3773: PetscObjectStateGet((PetscObject) spmap, &mesh->subpointState);
3774: }
3775: return 0;
3776: }
3778: /*@
3779: DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data
3781: Input Parameter:
3782: . dm - The submesh DM
3784: Output Parameter:
3785: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3787: Note: This IS is guaranteed to be sorted by the construction of the submesh
3789: Level: developer
3791: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3792: @*/
3793: PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS)
3794: {
3795: DM_Plex *mesh = (DM_Plex *) dm->data;
3796: DMLabel spmap;
3797: PetscObjectState state;
3801: DMPlexGetSubpointMap(dm, &spmap);
3802: PetscObjectStateGet((PetscObject) spmap, &state);
3803: if (state != mesh->subpointState || !mesh->subpointIS) DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS);
3804: *subpointIS = mesh->subpointIS;
3805: return 0;
3806: }
3808: /*@
3809: DMGetEnclosureRelation - Get the relationship between dmA and dmB
3811: Input Parameters:
3812: + dmA - The first DM
3813: - dmB - The second DM
3815: Output Parameter:
3816: . rel - The relation of dmA to dmB
3818: Level: intermediate
3820: .seealso: DMGetEnclosurePoint()
3821: @*/
3822: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
3823: {
3824: DM plexA, plexB, sdm;
3825: DMLabel spmap;
3826: PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB;
3829: *rel = DM_ENC_NONE;
3830: if (!dmA || !dmB) return 0;
3833: if (dmA == dmB) {*rel = DM_ENC_EQUALITY; return 0;}
3834: DMConvert(dmA, DMPLEX, &plexA);
3835: DMConvert(dmB, DMPLEX, &plexB);
3836: DMPlexGetChart(plexA, &pStartA, &pEndA);
3837: DMPlexGetChart(plexB, &pStartB, &pEndB);
3838: /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
3839: - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
3840: if ((pStartA == pStartB) && (pEndA == pEndB)) {
3841: *rel = DM_ENC_EQUALITY;
3842: goto end;
3843: }
3844: NpA = pEndA - pStartA;
3845: NpB = pEndB - pStartB;
3846: if (NpA == NpB) goto end;
3847: sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
3848: DMPlexGetSubpointMap(sdm, &spmap);
3849: if (!spmap) goto end;
3850: /* TODO Check the space mapped to by subpointMap is same size as dm */
3851: if (NpA > NpB) {
3852: *rel = DM_ENC_SUPERMESH;
3853: } else {
3854: *rel = DM_ENC_SUBMESH;
3855: }
3856: end:
3857: DMDestroy(&plexA);
3858: DMDestroy(&plexB);
3859: return 0;
3860: }
3862: /*@
3863: DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB
3865: Input Parameters:
3866: + dmA - The first DM
3867: . dmB - The second DM
3868: . etype - The type of enclosure relation that dmA has to dmB
3869: - pB - A point of dmB
3871: Output Parameter:
3872: . pA - The corresponding point of dmA
3874: Level: intermediate
3876: .seealso: DMGetEnclosureRelation()
3877: @*/
3878: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
3879: {
3880: DM sdm;
3881: IS subpointIS;
3882: const PetscInt *subpoints;
3883: PetscInt numSubpoints;
3885: /* TODO Cache the IS, making it look like an index */
3886: switch (etype) {
3887: case DM_ENC_SUPERMESH:
3888: sdm = dmB;
3889: DMPlexGetSubpointIS(sdm, &subpointIS);
3890: ISGetIndices(subpointIS, &subpoints);
3891: *pA = subpoints[pB];
3892: ISRestoreIndices(subpointIS, &subpoints);
3893: break;
3894: case DM_ENC_SUBMESH:
3895: sdm = dmA;
3896: DMPlexGetSubpointIS(sdm, &subpointIS);
3897: ISGetLocalSize(subpointIS, &numSubpoints);
3898: ISGetIndices(subpointIS, &subpoints);
3899: PetscFindInt(pB, numSubpoints, subpoints, pA);
3900: if (*pA < 0) {
3901: DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
3902: DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
3903: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB);
3904: }
3905: ISRestoreIndices(subpointIS, &subpoints);
3906: break;
3907: case DM_ENC_EQUALITY:
3908: case DM_ENC_NONE:
3909: *pA = pB;break;
3910: case DM_ENC_UNKNOWN:
3911: {
3912: DMEnclosureType enc;
3914: DMGetEnclosureRelation(dmA, dmB, &enc);
3915: DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
3916: }
3917: break;
3918: default: SETERRQ(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype);
3919: }
3920: return 0;
3921: }