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