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 Parameter:
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: DMCopyBoundary(dm, gdm);
962: DMCopyDisc(dm, gdm);
963: gdm->setfromoptionscalled = dm->setfromoptionscalled;
964: if (numGhostCells) *numGhostCells = Ng;
965: *dmGhosted = gdm;
966: return(0);
967: }
969: /*
970: We are adding three kinds of points here:
971: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
972: Non-replicated: Points which exist on the fault, but are not replicated
973: Hybrid: Entirely new points, such as cohesive cells
975: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
976: each depth so that the new split/hybrid points can be inserted as a block.
977: */
978: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
979: {
980: MPI_Comm comm;
981: IS valueIS;
982: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
983: const PetscInt *values; /* List of depths for which we have replicated points */
984: IS *splitIS;
985: IS *unsplitIS;
986: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
987: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
988: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
989: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
990: const PetscInt **splitPoints; /* Replicated points for each depth */
991: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
992: PetscSection coordSection;
993: Vec coordinates;
994: PetscScalar *coords;
995: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
996: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
997: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
998: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
999: PetscInt *coneNew, *coneONew, *supportNew;
1000: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
1001: PetscErrorCode ierr;
1004: PetscObjectGetComm((PetscObject)dm,&comm);
1005: DMGetDimension(dm, &dim);
1006: DMPlexGetDepth(dm, &depth);
1007: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1008: /* We do not want this label automatically computed, instead we compute it here */
1009: DMCreateLabel(sdm, "celltype");
1010: /* Count split points and add cohesive cells */
1011: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
1012: PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
1013: PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
1014: for (d = 0; d <= depth; ++d) {
1015: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
1016: DMPlexGetTensorPrismBounds_Internal(dm, d, &depthMax[d], NULL);
1017: depthEnd[d] = pMaxNew[d];
1018: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
1019: numSplitPoints[d] = 0;
1020: numUnsplitPoints[d] = 0;
1021: numHybridPoints[d] = 0;
1022: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
1023: splitPoints[d] = NULL;
1024: unsplitPoints[d] = NULL;
1025: splitIS[d] = NULL;
1026: unsplitIS[d] = NULL;
1027: /* we are shifting the existing hybrid points with the stratum behind them, so
1028: * the split comes at the end of the normal points, i.e., at depthMax[d] */
1029: depthShift[2*d] = depthMax[d];
1030: depthShift[2*d+1] = 0;
1031: }
1032: if (label) {
1033: DMLabelGetValueIS(label, &valueIS);
1034: ISGetLocalSize(valueIS, &numSP);
1035: ISGetIndices(valueIS, &values);
1036: }
1037: for (sp = 0; sp < numSP; ++sp) {
1038: const PetscInt dep = values[sp];
1040: if ((dep < 0) || (dep > depth)) continue;
1041: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
1042: if (splitIS[dep]) {
1043: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
1044: ISGetIndices(splitIS[dep], &splitPoints[dep]);
1045: }
1046: DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
1047: if (unsplitIS[dep]) {
1048: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
1049: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
1050: }
1051: }
1052: /* Calculate number of hybrid points */
1053: 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 */
1054: for (d = 0; d <= depth; ++d) depthShift[2*d+1] = numSplitPoints[d] + numHybridPoints[d];
1055: DMPlexShiftPointSetUp_Internal(depth,depthShift);
1056: /* the end of the points in this stratum that come before the new points:
1057: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
1058: * added points */
1059: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
1060: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
1061: /* Step 3: Set cone/support sizes for new points */
1062: for (dep = 0; dep <= depth; ++dep) {
1063: for (p = 0; p < numSplitPoints[dep]; ++p) {
1064: const PetscInt oldp = splitPoints[dep][p];
1065: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1066: const PetscInt splitp = p + pMaxNew[dep];
1067: const PetscInt *support;
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: if (dep == depth-1) {
1075: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1077: /* Add cohesive cells, they are prisms */
1078: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
1079: switch (coneSize) {
1080: case 2: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR);break;
1081: case 3: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR);break;
1082: case 4: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_QUAD_PRISM_TENSOR);break;
1083: }
1084: } else if (dep == 0) {
1085: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1087: DMPlexGetSupport(dm, oldp, &support);
1088: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1089: PetscInt val;
1091: DMLabelGetValue(label, support[e], &val);
1092: if (val == 1) ++qf;
1093: if ((val == 1) || (val == (shift + 1))) ++qn;
1094: if ((val == 1) || (val == -(shift + 1))) ++qp;
1095: }
1096: /* Split old vertex: Edges into original vertex and new cohesive edge */
1097: DMPlexSetSupportSize(sdm, newp, qn+1);
1098: /* Split new vertex: Edges into split vertex and new cohesive edge */
1099: DMPlexSetSupportSize(sdm, splitp, qp+1);
1100: /* Add hybrid edge */
1101: DMPlexSetConeSize(sdm, hybedge, 2);
1102: DMPlexSetSupportSize(sdm, hybedge, qf);
1103: DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1104: } else if (dep == dim-2) {
1105: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1107: DMPlexGetSupport(dm, oldp, &support);
1108: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1109: PetscInt val;
1111: DMLabelGetValue(label, support[e], &val);
1112: if (val == dim-1) ++qf;
1113: if ((val == dim-1) || (val == (shift + dim-1))) ++qn;
1114: if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
1115: }
1116: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1117: DMPlexSetSupportSize(sdm, newp, qn+1);
1118: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1119: DMPlexSetSupportSize(sdm, splitp, qp+1);
1120: /* Add hybrid face */
1121: DMPlexSetConeSize(sdm, hybface, 4);
1122: DMPlexSetSupportSize(sdm, hybface, qf);
1123: DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1124: }
1125: }
1126: }
1127: for (dep = 0; dep <= depth; ++dep) {
1128: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1129: const PetscInt oldp = unsplitPoints[dep][p];
1130: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1131: const PetscInt *support;
1132: PetscInt coneSize, supportSize, qf, e, s;
1134: DMPlexGetConeSize(dm, oldp, &coneSize);
1135: DMPlexGetSupportSize(dm, oldp, &supportSize);
1136: DMPlexGetSupport(dm, oldp, &support);
1137: if (dep == 0) {
1138: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1140: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1141: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1142: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1143: if (e >= 0) ++qf;
1144: }
1145: DMPlexSetSupportSize(sdm, newp, qf+2);
1146: /* Add hybrid edge */
1147: DMPlexSetConeSize(sdm, hybedge, 2);
1148: for (e = 0, qf = 0; e < supportSize; ++e) {
1149: PetscInt val;
1151: DMLabelGetValue(label, support[e], &val);
1152: /* Split and unsplit edges produce hybrid faces */
1153: if (val == 1) ++qf;
1154: if (val == (shift2 + 1)) ++qf;
1155: }
1156: DMPlexSetSupportSize(sdm, hybedge, qf);
1157: DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1158: } else if (dep == dim-2) {
1159: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1160: PetscInt val;
1162: for (e = 0, qf = 0; e < supportSize; ++e) {
1163: DMLabelGetValue(label, support[e], &val);
1164: if (val == dim-1) qf += 2;
1165: else ++qf;
1166: }
1167: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1168: DMPlexSetSupportSize(sdm, newp, qf+2);
1169: /* Add hybrid face */
1170: for (e = 0, qf = 0; e < supportSize; ++e) {
1171: DMLabelGetValue(label, support[e], &val);
1172: if (val == dim-1) ++qf;
1173: }
1174: DMPlexSetConeSize(sdm, hybface, 4);
1175: DMPlexSetSupportSize(sdm, hybface, qf);
1176: DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1177: }
1178: }
1179: }
1180: /* Step 4: Setup split DM */
1181: DMSetUp(sdm);
1182: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1183: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1184: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1185: /* Step 6: Set cones and supports for new points */
1186: for (dep = 0; dep <= depth; ++dep) {
1187: for (p = 0; p < numSplitPoints[dep]; ++p) {
1188: const PetscInt oldp = splitPoints[dep][p];
1189: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1190: const PetscInt splitp = p + pMaxNew[dep];
1191: const PetscInt *cone, *support, *ornt;
1192: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1194: DMPlexGetConeSize(dm, oldp, &coneSize);
1195: DMPlexGetCone(dm, oldp, &cone);
1196: DMPlexGetConeOrientation(dm, oldp, &ornt);
1197: DMPlexGetSupportSize(dm, oldp, &supportSize);
1198: DMPlexGetSupport(dm, oldp, &support);
1199: if (dep == depth-1) {
1200: PetscBool hasUnsplit = PETSC_FALSE;
1201: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1202: const PetscInt *supportF;
1204: /* Split face: copy in old face to new face to start */
1205: DMPlexGetSupport(sdm, newp, &supportF);
1206: DMPlexSetSupport(sdm, splitp, supportF);
1207: /* Split old face: old vertices/edges in cone so no change */
1208: /* Split new face: new vertices/edges in cone */
1209: for (q = 0; q < coneSize; ++q) {
1210: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1211: if (v < 0) {
1212: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1213: 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);
1214: coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1215: hasUnsplit = PETSC_TRUE;
1216: } else {
1217: coneNew[2+q] = v + pMaxNew[dep-1];
1218: if (dep > 1) {
1219: const PetscInt *econe;
1220: PetscInt econeSize, r, vs, vu;
1222: DMPlexGetConeSize(dm, cone[q], &econeSize);
1223: DMPlexGetCone(dm, cone[q], &econe);
1224: for (r = 0; r < econeSize; ++r) {
1225: PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);
1226: PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1227: if (vs >= 0) continue;
1228: 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);
1229: hasUnsplit = PETSC_TRUE;
1230: }
1231: }
1232: }
1233: }
1234: DMPlexSetCone(sdm, splitp, &coneNew[2]);
1235: DMPlexSetConeOrientation(sdm, splitp, ornt);
1236: /* Face support */
1237: for (s = 0; s < supportSize; ++s) {
1238: PetscInt val;
1240: DMLabelGetValue(label, support[s], &val);
1241: if (val < 0) {
1242: /* Split old face: Replace negative side cell with cohesive cell */
1243: DMPlexInsertSupport(sdm, newp, s, hybcell);
1244: } else {
1245: /* Split new face: Replace positive side cell with cohesive cell */
1246: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1247: /* Get orientation for cohesive face */
1248: {
1249: const PetscInt *ncone, *nconeO;
1250: PetscInt nconeSize, nc;
1252: DMPlexGetConeSize(dm, support[s], &nconeSize);
1253: DMPlexGetCone(dm, support[s], &ncone);
1254: DMPlexGetConeOrientation(dm, support[s], &nconeO);
1255: for (nc = 0; nc < nconeSize; ++nc) {
1256: if (ncone[nc] == oldp) {
1257: coneONew[0] = nconeO[nc];
1258: break;
1259: }
1260: }
1261: if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1262: }
1263: }
1264: }
1265: /* Cohesive cell: Old and new split face, then new cohesive faces */
1266: coneNew[0] = newp; /* Extracted negative side orientation above */
1267: coneNew[1] = splitp;
1268: coneONew[1] = coneONew[0];
1269: for (q = 0; q < coneSize; ++q) {
1270: /* Hybrid faces must follow order from oriented end face */
1271: const PetscInt o = coneONew[0];
1272: const PetscInt qo = o < 0 ? (-(o+1)+coneSize-q)%coneSize : (q+o)%coneSize;
1274: PetscFindInt(cone[qo], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1275: if (v < 0) {
1276: PetscFindInt(cone[qo], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1277: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1278: } else {
1279: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep];
1280: }
1281: coneONew[2+q] = ((o < 0) + (ornt[qo] < 0))%2 ? -1 : 0;
1282: }
1283: DMPlexSetCone(sdm, hybcell, coneNew);
1284: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1285: /* Label the hybrid cells on the boundary of the split */
1286: if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1287: } else if (dep == 0) {
1288: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1290: /* Split old vertex: Edges in old split faces and new cohesive edge */
1291: for (e = 0, qn = 0; e < supportSize; ++e) {
1292: PetscInt val;
1294: DMLabelGetValue(label, support[e], &val);
1295: if ((val == 1) || (val == (shift + 1))) {
1296: supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1297: }
1298: }
1299: supportNew[qn] = hybedge;
1300: DMPlexSetSupport(sdm, newp, supportNew);
1301: /* Split new vertex: Edges in new split faces and new cohesive edge */
1302: for (e = 0, qp = 0; e < supportSize; ++e) {
1303: PetscInt val, edge;
1305: DMLabelGetValue(label, support[e], &val);
1306: if (val == 1) {
1307: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1308: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1309: supportNew[qp++] = edge + pMaxNew[dep+1];
1310: } else if (val == -(shift + 1)) {
1311: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1312: }
1313: }
1314: supportNew[qp] = hybedge;
1315: DMPlexSetSupport(sdm, splitp, supportNew);
1316: /* Hybrid edge: Old and new split vertex */
1317: coneNew[0] = newp;
1318: coneNew[1] = splitp;
1319: DMPlexSetCone(sdm, hybedge, coneNew);
1320: for (e = 0, qf = 0; e < supportSize; ++e) {
1321: PetscInt val, edge;
1323: DMLabelGetValue(label, support[e], &val);
1324: if (val == 1) {
1325: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1326: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1327: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1328: }
1329: }
1330: DMPlexSetSupport(sdm, hybedge, supportNew);
1331: } else if (dep == dim-2) {
1332: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1334: /* Split old edge: old vertices in cone so no change */
1335: /* Split new edge: new vertices in cone */
1336: for (q = 0; q < coneSize; ++q) {
1337: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1338: if (v < 0) {
1339: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1340: 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);
1341: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1342: } else {
1343: coneNew[q] = v + pMaxNew[dep-1];
1344: }
1345: }
1346: DMPlexSetCone(sdm, splitp, coneNew);
1347: /* Split old edge: Faces in positive side cells and old split faces */
1348: for (e = 0, q = 0; e < supportSize; ++e) {
1349: PetscInt val;
1351: DMLabelGetValue(label, support[e], &val);
1352: if (val == dim-1) {
1353: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1354: } else if (val == (shift + dim-1)) {
1355: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1356: }
1357: }
1358: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1359: DMPlexSetSupport(sdm, newp, supportNew);
1360: /* Split new edge: Faces in negative side cells and new split faces */
1361: for (e = 0, q = 0; e < supportSize; ++e) {
1362: PetscInt val, face;
1364: DMLabelGetValue(label, support[e], &val);
1365: if (val == dim-1) {
1366: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1367: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1368: supportNew[q++] = face + pMaxNew[dep+1];
1369: } else if (val == -(shift + dim-1)) {
1370: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1371: }
1372: }
1373: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1374: DMPlexSetSupport(sdm, splitp, supportNew);
1375: /* Hybrid face */
1376: coneNew[0] = newp;
1377: coneNew[1] = splitp;
1378: for (v = 0; v < coneSize; ++v) {
1379: PetscInt vertex;
1380: PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1381: if (vertex < 0) {
1382: PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1383: 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);
1384: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1385: } else {
1386: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1387: }
1388: }
1389: DMPlexSetCone(sdm, hybface, coneNew);
1390: for (e = 0, qf = 0; e < supportSize; ++e) {
1391: PetscInt val, face;
1393: DMLabelGetValue(label, support[e], &val);
1394: if (val == dim-1) {
1395: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1396: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1397: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1398: }
1399: }
1400: DMPlexSetSupport(sdm, hybface, supportNew);
1401: }
1402: }
1403: }
1404: for (dep = 0; dep <= depth; ++dep) {
1405: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1406: const PetscInt oldp = unsplitPoints[dep][p];
1407: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1408: const PetscInt *cone, *support, *ornt;
1409: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1411: DMPlexGetConeSize(dm, oldp, &coneSize);
1412: DMPlexGetCone(dm, oldp, &cone);
1413: DMPlexGetConeOrientation(dm, oldp, &ornt);
1414: DMPlexGetSupportSize(dm, oldp, &supportSize);
1415: DMPlexGetSupport(dm, oldp, &support);
1416: if (dep == 0) {
1417: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1419: /* Unsplit vertex */
1420: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1421: for (s = 0, q = 0; s < supportSize; ++s) {
1422: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1423: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1424: if (e >= 0) {
1425: supportNew[q++] = e + pMaxNew[dep+1];
1426: }
1427: }
1428: supportNew[q++] = hybedge;
1429: supportNew[q++] = hybedge;
1430: if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1431: DMPlexSetSupport(sdm, newp, supportNew);
1432: /* Hybrid edge */
1433: coneNew[0] = newp;
1434: coneNew[1] = newp;
1435: DMPlexSetCone(sdm, hybedge, coneNew);
1436: for (e = 0, qf = 0; e < supportSize; ++e) {
1437: PetscInt val, edge;
1439: DMLabelGetValue(label, support[e], &val);
1440: if (val == 1) {
1441: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1442: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1443: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1444: } else if (val == (shift2 + 1)) {
1445: PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1446: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1447: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1448: }
1449: }
1450: DMPlexSetSupport(sdm, hybedge, supportNew);
1451: } else if (dep == dim-2) {
1452: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1454: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1455: for (f = 0, qf = 0; f < supportSize; ++f) {
1456: PetscInt val, face;
1458: DMLabelGetValue(label, support[f], &val);
1459: if (val == dim-1) {
1460: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1461: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1462: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1463: supportNew[qf++] = face + pMaxNew[dep+1];
1464: } else {
1465: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1466: }
1467: }
1468: supportNew[qf++] = hybface;
1469: supportNew[qf++] = hybface;
1470: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1471: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1472: DMPlexSetSupport(sdm, newp, supportNew);
1473: /* Add hybrid face */
1474: coneNew[0] = newp;
1475: coneNew[1] = newp;
1476: PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1477: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1478: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1479: PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1480: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1481: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1482: DMPlexSetCone(sdm, hybface, coneNew);
1483: for (f = 0, qf = 0; f < supportSize; ++f) {
1484: PetscInt val, face;
1486: DMLabelGetValue(label, support[f], &val);
1487: if (val == dim-1) {
1488: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1489: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1490: }
1491: }
1492: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1493: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1494: DMPlexSetSupport(sdm, hybface, supportNew);
1495: }
1496: }
1497: }
1498: /* Step 6b: Replace split points in negative side cones */
1499: for (sp = 0; sp < numSP; ++sp) {
1500: PetscInt dep = values[sp];
1501: IS pIS;
1502: PetscInt numPoints;
1503: const PetscInt *points;
1505: if (dep >= 0) continue;
1506: DMLabelGetStratumIS(label, dep, &pIS);
1507: if (!pIS) continue;
1508: dep = -dep - shift;
1509: ISGetLocalSize(pIS, &numPoints);
1510: ISGetIndices(pIS, &points);
1511: for (p = 0; p < numPoints; ++p) {
1512: const PetscInt oldp = points[p];
1513: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1514: const PetscInt *cone;
1515: PetscInt coneSize, c;
1516: /* PetscBool replaced = PETSC_FALSE; */
1518: /* Negative edge: replace split vertex */
1519: /* Negative cell: replace split face */
1520: DMPlexGetConeSize(sdm, newp, &coneSize);
1521: DMPlexGetCone(sdm, newp, &cone);
1522: for (c = 0; c < coneSize; ++c) {
1523: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1524: PetscInt csplitp, cp, val;
1526: DMLabelGetValue(label, coldp, &val);
1527: if (val == dep-1) {
1528: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1529: if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1530: csplitp = pMaxNew[dep-1] + cp;
1531: DMPlexInsertCone(sdm, newp, c, csplitp);
1532: /* replaced = PETSC_TRUE; */
1533: }
1534: }
1535: /* Cells with only a vertex or edge on the submesh have no replacement */
1536: /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1537: }
1538: ISRestoreIndices(pIS, &points);
1539: ISDestroy(&pIS);
1540: }
1541: /* Step 7: Coordinates */
1542: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1543: DMGetCoordinateSection(sdm, &coordSection);
1544: DMGetCoordinatesLocal(sdm, &coordinates);
1545: VecGetArray(coordinates, &coords);
1546: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1547: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1548: const PetscInt splitp = pMaxNew[0] + v;
1549: PetscInt dof, off, soff, d;
1551: PetscSectionGetDof(coordSection, newp, &dof);
1552: PetscSectionGetOffset(coordSection, newp, &off);
1553: PetscSectionGetOffset(coordSection, splitp, &soff);
1554: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1555: }
1556: VecRestoreArray(coordinates, &coords);
1557: /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1558: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1559: /* Step 9: Labels */
1560: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1561: DMGetNumLabels(sdm, &numLabels);
1562: for (dep = 0; dep <= depth; ++dep) {
1563: for (p = 0; p < numSplitPoints[dep]; ++p) {
1564: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1565: const PetscInt splitp = pMaxNew[dep] + p;
1566: PetscInt l;
1568: if (splitLabel) {
1569: const PetscInt val = 100 + dep;
1571: DMLabelSetValue(splitLabel, newp, val);
1572: DMLabelSetValue(splitLabel, splitp, -val);
1573: }
1574: for (l = 0; l < numLabels; ++l) {
1575: DMLabel mlabel;
1576: const char *lname;
1577: PetscInt val;
1578: PetscBool isDepth;
1580: DMGetLabelName(sdm, l, &lname);
1581: PetscStrcmp(lname, "depth", &isDepth);
1582: if (isDepth) continue;
1583: DMGetLabel(sdm, lname, &mlabel);
1584: DMLabelGetValue(mlabel, newp, &val);
1585: if (val >= 0) {
1586: DMLabelSetValue(mlabel, splitp, val);
1587: }
1588: }
1589: }
1590: }
1591: for (sp = 0; sp < numSP; ++sp) {
1592: const PetscInt dep = values[sp];
1594: if ((dep < 0) || (dep > depth)) continue;
1595: if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1596: ISDestroy(&splitIS[dep]);
1597: if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1598: ISDestroy(&unsplitIS[dep]);
1599: }
1600: if (label) {
1601: ISRestoreIndices(valueIS, &values);
1602: ISDestroy(&valueIS);
1603: }
1604: for (d = 0; d <= depth; ++d) {
1605: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1606: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1607: }
1608: PetscFree3(coneNew, coneONew, supportNew);
1609: PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1610: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1611: return(0);
1612: }
1614: /*@C
1615: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1617: Collective on dm
1619: Input Parameters:
1620: + dm - The original DM
1621: - label - The label specifying the boundary faces (this could be auto-generated)
1623: Output Parameters:
1624: + splitLabel - The label containing the split points, or NULL if no output is desired
1625: - dmSplit - The new DM
1627: Level: developer
1629: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1630: @*/
1631: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
1632: {
1633: DM sdm;
1634: PetscInt dim;
1640: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1641: DMSetType(sdm, DMPLEX);
1642: DMGetDimension(dm, &dim);
1643: DMSetDimension(sdm, dim);
1644: switch (dim) {
1645: case 2:
1646: case 3:
1647: DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
1648: break;
1649: default:
1650: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1651: }
1652: *dmSplit = sdm;
1653: return(0);
1654: }
1656: /* Returns the side of the surface for a given cell with a face on the surface */
1657: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1658: {
1659: const PetscInt *cone, *ornt;
1660: PetscInt dim, coneSize, c;
1661: PetscErrorCode ierr;
1664: *pos = PETSC_TRUE;
1665: DMGetDimension(dm, &dim);
1666: DMPlexGetConeSize(dm, cell, &coneSize);
1667: DMPlexGetCone(dm, cell, &cone);
1668: DMPlexGetConeOrientation(dm, cell, &ornt);
1669: for (c = 0; c < coneSize; ++c) {
1670: if (cone[c] == face) {
1671: PetscInt o = ornt[c];
1673: if (subdm) {
1674: const PetscInt *subcone, *subornt;
1675: PetscInt subpoint, subface, subconeSize, sc;
1677: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1678: PetscFindInt(face, numSubpoints, subpoints, &subface);
1679: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1680: DMPlexGetCone(subdm, subpoint, &subcone);
1681: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1682: for (sc = 0; sc < subconeSize; ++sc) {
1683: if (subcone[sc] == subface) {
1684: o = subornt[0];
1685: break;
1686: }
1687: }
1688: 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);
1689: }
1690: if (o >= 0) *pos = PETSC_TRUE;
1691: else *pos = PETSC_FALSE;
1692: break;
1693: }
1694: }
1695: 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);
1696: return(0);
1697: }
1699: /*@
1700: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1701: to complete the surface
1703: Input Parameters:
1704: + dm - The DM
1705: . label - A DMLabel marking the surface
1706: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1707: . flip - Flag to flip the submesh normal and replace points on the other side
1708: - subdm - The subDM associated with the label, or NULL
1710: Output Parameter:
1711: . label - A DMLabel marking all surface points
1713: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1715: Level: developer
1717: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1718: @*/
1719: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1720: {
1721: DMLabel depthLabel;
1722: IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1723: const PetscInt *points, *subpoints;
1724: const PetscInt rev = flip ? -1 : 1;
1725: PetscInt shift = 100, shift2 = 200, dim, depth, dep, cStart, cEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1726: PetscErrorCode ierr;
1729: DMPlexGetDepth(dm, &depth);
1730: DMGetDimension(dm, &dim);
1731: DMPlexGetDepthLabel(dm, &depthLabel);
1732: if (subdm) {
1733: DMPlexGetSubpointIS(subdm, &subpointIS);
1734: if (subpointIS) {
1735: ISGetLocalSize(subpointIS, &numSubpoints);
1736: ISGetIndices(subpointIS, &subpoints);
1737: }
1738: }
1739: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1740: DMLabelGetStratumIS(label, dim-1, &dimIS);
1741: if (!dimIS) return(0);
1742: ISGetLocalSize(dimIS, &numPoints);
1743: ISGetIndices(dimIS, &points);
1744: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1745: const PetscInt *support;
1746: PetscInt supportSize, s;
1748: DMPlexGetSupportSize(dm, points[p], &supportSize);
1749: #if 0
1750: if (supportSize != 2) {
1751: const PetscInt *lp;
1752: PetscInt Nlp, pind;
1754: /* Check that for a cell with a single support face, that face is in the SF */
1755: /* THis check only works for the remote side. We would need root side information */
1756: PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1757: PetscFindInt(points[p], Nlp, lp, &pind);
1758: 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);
1759: }
1760: #endif
1761: DMPlexGetSupport(dm, points[p], &support);
1762: for (s = 0; s < supportSize; ++s) {
1763: const PetscInt *cone;
1764: PetscInt coneSize, c;
1765: PetscBool pos;
1767: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1768: if (pos) {DMLabelSetValue(label, support[s], rev*(shift+dim));}
1769: else {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1770: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1771: /* Put faces touching the fault in the label */
1772: DMPlexGetConeSize(dm, support[s], &coneSize);
1773: DMPlexGetCone(dm, support[s], &cone);
1774: for (c = 0; c < coneSize; ++c) {
1775: const PetscInt point = cone[c];
1777: DMLabelGetValue(label, point, &val);
1778: if (val == -1) {
1779: PetscInt *closure = NULL;
1780: PetscInt closureSize, cl;
1782: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1783: for (cl = 0; cl < closureSize*2; cl += 2) {
1784: const PetscInt clp = closure[cl];
1785: PetscInt bval = -1;
1787: DMLabelGetValue(label, clp, &val);
1788: if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1789: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1790: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1791: break;
1792: }
1793: }
1794: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1795: }
1796: }
1797: }
1798: }
1799: ISRestoreIndices(dimIS, &points);
1800: ISDestroy(&dimIS);
1801: if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1802: /* Mark boundary points as unsplit */
1803: if (blabel) {
1804: DMLabelGetStratumIS(blabel, 1, &dimIS);
1805: ISGetLocalSize(dimIS, &numPoints);
1806: ISGetIndices(dimIS, &points);
1807: for (p = 0; p < numPoints; ++p) {
1808: const PetscInt point = points[p];
1809: PetscInt val, bval;
1811: DMLabelGetValue(blabel, point, &bval);
1812: if (bval >= 0) {
1813: DMLabelGetValue(label, point, &val);
1814: if ((val < 0) || (val > dim)) {
1815: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1816: DMLabelClearValue(blabel, point, bval);
1817: }
1818: }
1819: }
1820: for (p = 0; p < numPoints; ++p) {
1821: const PetscInt point = points[p];
1822: PetscInt val, bval;
1824: DMLabelGetValue(blabel, point, &bval);
1825: if (bval >= 0) {
1826: const PetscInt *cone, *support;
1827: PetscInt coneSize, supportSize, s, valA, valB, valE;
1829: /* Mark as unsplit */
1830: DMLabelGetValue(label, point, &val);
1831: 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);
1832: DMLabelClearValue(label, point, val);
1833: DMLabelSetValue(label, point, shift2+val);
1834: /* Check for cross-edge
1835: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1836: if (val != 0) continue;
1837: DMPlexGetSupport(dm, point, &support);
1838: DMPlexGetSupportSize(dm, point, &supportSize);
1839: for (s = 0; s < supportSize; ++s) {
1840: DMPlexGetCone(dm, support[s], &cone);
1841: DMPlexGetConeSize(dm, support[s], &coneSize);
1842: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1843: DMLabelGetValue(blabel, cone[0], &valA);
1844: DMLabelGetValue(blabel, cone[1], &valB);
1845: DMLabelGetValue(blabel, support[s], &valE);
1846: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1847: }
1848: }
1849: }
1850: ISRestoreIndices(dimIS, &points);
1851: ISDestroy(&dimIS);
1852: }
1853: /* Search for other cells/faces/edges connected to the fault by a vertex */
1854: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1855: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1856: DMLabelGetStratumIS(label, 0, &dimIS);
1857: if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1858: if (dimIS && crossEdgeIS) {
1859: IS vertIS = dimIS;
1861: ISExpand(vertIS, crossEdgeIS, &dimIS);
1862: ISDestroy(&crossEdgeIS);
1863: ISDestroy(&vertIS);
1864: }
1865: if (!dimIS) {
1866: return(0);
1867: }
1868: ISGetLocalSize(dimIS, &numPoints);
1869: ISGetIndices(dimIS, &points);
1870: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1871: PetscInt *star = NULL;
1872: PetscInt starSize, s;
1873: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1875: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1876: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1877: while (again) {
1878: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1879: again = 0;
1880: for (s = 0; s < starSize*2; s += 2) {
1881: const PetscInt point = star[s];
1882: const PetscInt *cone;
1883: PetscInt coneSize, c;
1885: if ((point < cStart) || (point >= cEnd)) continue;
1886: DMLabelGetValue(label, point, &val);
1887: if (val != -1) continue;
1888: again = again == 1 ? 1 : 2;
1889: DMPlexGetConeSize(dm, point, &coneSize);
1890: DMPlexGetCone(dm, point, &cone);
1891: for (c = 0; c < coneSize; ++c) {
1892: DMLabelGetValue(label, cone[c], &val);
1893: if (val != -1) {
1894: const PetscInt *ccone;
1895: PetscInt cconeSize, cc, side;
1897: 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);
1898: if (val > 0) side = 1;
1899: else side = -1;
1900: DMLabelSetValue(label, point, side*(shift+dim));
1901: /* Mark cell faces which touch the fault */
1902: DMPlexGetConeSize(dm, point, &cconeSize);
1903: DMPlexGetCone(dm, point, &ccone);
1904: for (cc = 0; cc < cconeSize; ++cc) {
1905: PetscInt *closure = NULL;
1906: PetscInt closureSize, cl;
1908: DMLabelGetValue(label, ccone[cc], &val);
1909: if (val != -1) continue;
1910: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1911: for (cl = 0; cl < closureSize*2; cl += 2) {
1912: const PetscInt clp = closure[cl];
1914: DMLabelGetValue(label, clp, &val);
1915: if (val == -1) continue;
1916: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1917: break;
1918: }
1919: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1920: }
1921: again = 1;
1922: break;
1923: }
1924: }
1925: }
1926: }
1927: /* Classify the rest by cell membership */
1928: for (s = 0; s < starSize*2; s += 2) {
1929: const PetscInt point = star[s];
1931: DMLabelGetValue(label, point, &val);
1932: if (val == -1) {
1933: PetscInt *sstar = NULL;
1934: PetscInt sstarSize, ss;
1935: PetscBool marked = PETSC_FALSE, isHybrid;
1937: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1938: for (ss = 0; ss < sstarSize*2; ss += 2) {
1939: const PetscInt spoint = sstar[ss];
1941: if ((spoint < cStart) || (spoint >= cEnd)) continue;
1942: DMLabelGetValue(label, spoint, &val);
1943: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1944: DMLabelGetValue(depthLabel, point, &dep);
1945: if (val > 0) {
1946: DMLabelSetValue(label, point, shift+dep);
1947: } else {
1948: DMLabelSetValue(label, point, -(shift+dep));
1949: }
1950: marked = PETSC_TRUE;
1951: break;
1952: }
1953: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1954: DMPlexCellIsHybrid_Internal(dm, point, &isHybrid);
1955: if (!isHybrid && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1956: }
1957: }
1958: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1959: }
1960: ISRestoreIndices(dimIS, &points);
1961: ISDestroy(&dimIS);
1962: /* If any faces touching the fault divide cells on either side, split them */
1963: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1964: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1965: ISExpand(facePosIS, faceNegIS, &dimIS);
1966: ISDestroy(&facePosIS);
1967: ISDestroy(&faceNegIS);
1968: ISGetLocalSize(dimIS, &numPoints);
1969: ISGetIndices(dimIS, &points);
1970: for (p = 0; p < numPoints; ++p) {
1971: const PetscInt point = points[p];
1972: const PetscInt *support;
1973: PetscInt supportSize, valA, valB;
1975: DMPlexGetSupportSize(dm, point, &supportSize);
1976: if (supportSize != 2) continue;
1977: DMPlexGetSupport(dm, point, &support);
1978: DMLabelGetValue(label, support[0], &valA);
1979: DMLabelGetValue(label, support[1], &valB);
1980: if ((valA == -1) || (valB == -1)) continue;
1981: if (valA*valB > 0) continue;
1982: /* Split the face */
1983: DMLabelGetValue(label, point, &valA);
1984: DMLabelClearValue(label, point, valA);
1985: DMLabelSetValue(label, point, dim-1);
1986: /* Label its closure:
1987: unmarked: label as unsplit
1988: incident: relabel as split
1989: split: do nothing
1990: */
1991: {
1992: PetscInt *closure = NULL;
1993: PetscInt closureSize, cl;
1995: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1996: for (cl = 0; cl < closureSize*2; cl += 2) {
1997: DMLabelGetValue(label, closure[cl], &valA);
1998: if (valA == -1) { /* Mark as unsplit */
1999: DMLabelGetValue(depthLabel, closure[cl], &dep);
2000: DMLabelSetValue(label, closure[cl], shift2+dep);
2001: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
2002: DMLabelGetValue(depthLabel, closure[cl], &dep);
2003: DMLabelClearValue(label, closure[cl], valA);
2004: DMLabelSetValue(label, closure[cl], dep);
2005: }
2006: }
2007: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2008: }
2009: }
2010: ISRestoreIndices(dimIS, &points);
2011: ISDestroy(&dimIS);
2012: return(0);
2013: }
2015: /* Check that no cell have all vertices on the fault */
2016: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
2017: {
2018: IS subpointIS;
2019: const PetscInt *dmpoints;
2020: PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd;
2021: PetscErrorCode ierr;
2024: if (!label) return(0);
2025: DMLabelGetDefaultValue(label, &defaultValue);
2026: DMPlexGetSubpointIS(subdm, &subpointIS);
2027: if (!subpointIS) return(0);
2028: DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
2029: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2030: ISGetIndices(subpointIS, &dmpoints);
2031: for (c = cStart; c < cEnd; ++c) {
2032: PetscBool invalidCell = PETSC_TRUE;
2033: PetscInt *closure = NULL;
2034: PetscInt closureSize, cl;
2036: DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2037: for (cl = 0; cl < closureSize*2; cl += 2) {
2038: PetscInt value = 0;
2040: if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
2041: DMLabelGetValue(label, closure[cl], &value);
2042: if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
2043: }
2044: DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2045: if (invalidCell) {
2046: ISRestoreIndices(subpointIS, &dmpoints);
2047: ISDestroy(&subpointIS);
2048: DMDestroy(&subdm);
2049: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
2050: }
2051: }
2052: ISRestoreIndices(subpointIS, &dmpoints);
2053: return(0);
2054: }
2057: /*@
2058: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
2060: Collective on dm
2062: Input Parameters:
2063: + dm - The original DM
2064: . label - The label specifying the interface vertices
2065: - bdlabel - The optional label specifying the interface boundary vertices
2067: Output Parameters:
2068: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
2069: . splitLabel - The label containing the split points, or NULL if no output is desired
2070: . dmInterface - The new interface DM, or NULL
2071: - dmHybrid - The new DM with cohesive cells
2073: Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points
2074: directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
2075: 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
2076: 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
2077: surface also hits vertex 3, it would be 9 (-101) in hybridLabel.
2079: The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
2080: mesh. The label value is +=100+dim for each point. For example, if two edges 10 and 14 in the hybrid resulting from
2081: splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.
2083: The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using
2084: DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh.
2086: Level: developer
2088: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
2089: @*/
2090: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2091: {
2092: DM idm;
2093: DMLabel subpointMap, hlabel, slabel = NULL;
2094: PetscInt dim;
2104: DMGetDimension(dm, &dim);
2105: DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2106: DMPlexCheckValidSubmesh_Private(dm, label, idm);
2107: DMPlexOrient(idm);
2108: DMPlexGetSubpointMap(idm, &subpointMap);
2109: DMLabelDuplicate(subpointMap, &hlabel);
2110: DMLabelClearStratum(hlabel, dim);
2111: if (splitLabel) {
2112: const char *name;
2113: char sname[PETSC_MAX_PATH_LEN];
2115: PetscObjectGetName((PetscObject) hlabel, &name);
2116: PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2117: PetscStrcat(sname, " split");
2118: DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2119: }
2120: DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2121: if (dmInterface) {*dmInterface = idm;}
2122: else {DMDestroy(&idm);}
2123: DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2124: if (hybridLabel) *hybridLabel = hlabel;
2125: else {DMLabelDestroy(&hlabel);}
2126: if (splitLabel) *splitLabel = slabel;
2127: return(0);
2128: }
2130: /* Here we need the explicit assumption that:
2132: For any marked cell, the marked vertices constitute a single face
2133: */
2134: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2135: {
2136: IS subvertexIS = NULL;
2137: const PetscInt *subvertices;
2138: PetscInt *pStart, *pEnd, pSize;
2139: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
2140: PetscErrorCode ierr;
2143: *numFaces = 0;
2144: *nFV = 0;
2145: DMPlexGetDepth(dm, &depth);
2146: DMGetDimension(dm, &dim);
2147: pSize = PetscMax(depth, dim) + 1;
2148: PetscMalloc2(pSize, &pStart, pSize, &pEnd);
2149: for (d = 0; d <= depth; ++d) {
2150: DMPlexGetSimplexOrBoxCells(dm, depth-d, &pStart[d], &pEnd[d]);
2151: }
2152: /* Loop over initial vertices and mark all faces in the collective star() */
2153: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
2154: if (subvertexIS) {
2155: ISGetSize(subvertexIS, &numSubVerticesInitial);
2156: ISGetIndices(subvertexIS, &subvertices);
2157: }
2158: for (v = 0; v < numSubVerticesInitial; ++v) {
2159: const PetscInt vertex = subvertices[v];
2160: PetscInt *star = NULL;
2161: PetscInt starSize, s, numCells = 0, c;
2163: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2164: for (s = 0; s < starSize*2; s += 2) {
2165: const PetscInt point = star[s];
2166: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2167: }
2168: for (c = 0; c < numCells; ++c) {
2169: const PetscInt cell = star[c];
2170: PetscInt *closure = NULL;
2171: PetscInt closureSize, cl;
2172: PetscInt cellLoc, numCorners = 0, faceSize = 0;
2174: DMLabelGetValue(subpointMap, cell, &cellLoc);
2175: if (cellLoc == 2) continue;
2176: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2177: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2178: for (cl = 0; cl < closureSize*2; cl += 2) {
2179: const PetscInt point = closure[cl];
2180: PetscInt vertexLoc;
2182: if ((point >= pStart[0]) && (point < pEnd[0])) {
2183: ++numCorners;
2184: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2185: if (vertexLoc == value) closure[faceSize++] = point;
2186: }
2187: }
2188: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2189: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2190: if (faceSize == *nFV) {
2191: const PetscInt *cells = NULL;
2192: PetscInt numCells, nc;
2194: ++(*numFaces);
2195: for (cl = 0; cl < faceSize; ++cl) {
2196: DMLabelSetValue(subpointMap, closure[cl], 0);
2197: }
2198: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2199: for (nc = 0; nc < numCells; ++nc) {
2200: DMLabelSetValue(subpointMap, cells[nc], 2);
2201: }
2202: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2203: }
2204: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2205: }
2206: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2207: }
2208: if (subvertexIS) {
2209: ISRestoreIndices(subvertexIS, &subvertices);
2210: }
2211: ISDestroy(&subvertexIS);
2212: PetscFree2(pStart, pEnd);
2213: return(0);
2214: }
2216: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2217: {
2218: IS subvertexIS = NULL;
2219: const PetscInt *subvertices;
2220: PetscInt *pStart, *pEnd;
2221: PetscInt dim, d, numSubVerticesInitial = 0, v;
2222: PetscErrorCode ierr;
2225: DMGetDimension(dm, &dim);
2226: PetscMalloc2(dim+1, &pStart, dim+1, &pEnd);
2227: for (d = 0; d <= dim; ++d) {
2228: DMPlexGetSimplexOrBoxCells(dm, dim-d, &pStart[d], &pEnd[d]);
2229: }
2230: /* Loop over initial vertices and mark all faces in the collective star() */
2231: if (vertexLabel) {
2232: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2233: if (subvertexIS) {
2234: ISGetSize(subvertexIS, &numSubVerticesInitial);
2235: ISGetIndices(subvertexIS, &subvertices);
2236: }
2237: }
2238: for (v = 0; v < numSubVerticesInitial; ++v) {
2239: const PetscInt vertex = subvertices[v];
2240: PetscInt *star = NULL;
2241: PetscInt starSize, s, numFaces = 0, f;
2243: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2244: for (s = 0; s < starSize*2; s += 2) {
2245: const PetscInt point = star[s];
2246: PetscInt faceLoc;
2248: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2249: if (markedFaces) {
2250: DMLabelGetValue(vertexLabel, point, &faceLoc);
2251: if (faceLoc < 0) continue;
2252: }
2253: star[numFaces++] = point;
2254: }
2255: }
2256: for (f = 0; f < numFaces; ++f) {
2257: const PetscInt face = star[f];
2258: PetscInt *closure = NULL;
2259: PetscInt closureSize, c;
2260: PetscInt faceLoc;
2262: DMLabelGetValue(subpointMap, face, &faceLoc);
2263: if (faceLoc == dim-1) continue;
2264: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2265: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2266: for (c = 0; c < closureSize*2; c += 2) {
2267: const PetscInt point = closure[c];
2268: PetscInt vertexLoc;
2270: if ((point >= pStart[0]) && (point < pEnd[0])) {
2271: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2272: if (vertexLoc != value) break;
2273: }
2274: }
2275: if (c == closureSize*2) {
2276: const PetscInt *support;
2277: PetscInt supportSize, s;
2279: for (c = 0; c < closureSize*2; c += 2) {
2280: const PetscInt point = closure[c];
2282: for (d = 0; d < dim; ++d) {
2283: if ((point >= pStart[d]) && (point < pEnd[d])) {
2284: DMLabelSetValue(subpointMap, point, d);
2285: break;
2286: }
2287: }
2288: }
2289: DMPlexGetSupportSize(dm, face, &supportSize);
2290: DMPlexGetSupport(dm, face, &support);
2291: for (s = 0; s < supportSize; ++s) {
2292: DMLabelSetValue(subpointMap, support[s], dim);
2293: }
2294: }
2295: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2296: }
2297: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2298: }
2299: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2300: ISDestroy(&subvertexIS);
2301: PetscFree2(pStart, pEnd);
2302: return(0);
2303: }
2305: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2306: {
2307: DMLabel label = NULL;
2308: const PetscInt *cone;
2309: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2310: PetscErrorCode ierr;
2313: *numFaces = 0;
2314: *nFV = 0;
2315: if (labelname) {DMGetLabel(dm, labelname, &label);}
2316: *subCells = NULL;
2317: DMGetDimension(dm, &dim);
2318: DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2319: if (cMax < 0) return(0);
2320: if (label) {
2321: for (c = cMax; c < cEnd; ++c) {
2322: PetscInt val;
2324: DMLabelGetValue(label, c, &val);
2325: if (val == value) {
2326: ++(*numFaces);
2327: DMPlexGetConeSize(dm, c, &coneSize);
2328: }
2329: }
2330: } else {
2331: *numFaces = cEnd - cMax;
2332: DMPlexGetConeSize(dm, cMax, &coneSize);
2333: }
2334: PetscMalloc1(*numFaces *2, subCells);
2335: if (!(*numFaces)) return(0);
2336: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2337: for (c = cMax; c < cEnd; ++c) {
2338: const PetscInt *cells;
2339: PetscInt numCells;
2341: if (label) {
2342: PetscInt val;
2344: DMLabelGetValue(label, c, &val);
2345: if (val != value) continue;
2346: }
2347: DMPlexGetCone(dm, c, &cone);
2348: for (p = 0; p < *nFV; ++p) {
2349: DMLabelSetValue(subpointMap, cone[p], 0);
2350: }
2351: /* Negative face */
2352: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2353: /* Not true in parallel
2354: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2355: for (p = 0; p < numCells; ++p) {
2356: DMLabelSetValue(subpointMap, cells[p], 2);
2357: (*subCells)[subc++] = cells[p];
2358: }
2359: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2360: /* Positive face is not included */
2361: }
2362: return(0);
2363: }
2365: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2366: {
2367: PetscInt *pStart, *pEnd;
2368: PetscInt dim, cMax, cEnd, c, d;
2372: DMGetDimension(dm, &dim);
2373: DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2374: if (cMax < 0) return(0);
2375: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2376: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2377: for (c = cMax; c < cEnd; ++c) {
2378: const PetscInt *cone;
2379: PetscInt *closure = NULL;
2380: PetscInt fconeSize, coneSize, closureSize, cl, val;
2382: if (label) {
2383: DMLabelGetValue(label, c, &val);
2384: if (val != value) continue;
2385: }
2386: DMPlexGetConeSize(dm, c, &coneSize);
2387: DMPlexGetCone(dm, c, &cone);
2388: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2389: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2390: /* Negative face */
2391: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2392: for (cl = 0; cl < closureSize*2; cl += 2) {
2393: const PetscInt point = closure[cl];
2395: for (d = 0; d <= dim; ++d) {
2396: if ((point >= pStart[d]) && (point < pEnd[d])) {
2397: DMLabelSetValue(subpointMap, point, d);
2398: break;
2399: }
2400: }
2401: }
2402: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2403: /* Cells -- positive face is not included */
2404: for (cl = 0; cl < 1; ++cl) {
2405: const PetscInt *support;
2406: PetscInt supportSize, s;
2408: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2409: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2410: DMPlexGetSupport(dm, cone[cl], &support);
2411: for (s = 0; s < supportSize; ++s) {
2412: DMLabelSetValue(subpointMap, support[s], dim);
2413: }
2414: }
2415: }
2416: PetscFree2(pStart, pEnd);
2417: return(0);
2418: }
2420: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2421: {
2422: MPI_Comm comm;
2423: PetscBool posOrient = PETSC_FALSE;
2424: const PetscInt debug = 0;
2425: PetscInt cellDim, faceSize, f;
2429: PetscObjectGetComm((PetscObject)dm,&comm);
2430: DMGetDimension(dm, &cellDim);
2431: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2433: if (cellDim == 1 && numCorners == 2) {
2434: /* Triangle */
2435: faceSize = numCorners-1;
2436: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2437: } else if (cellDim == 2 && numCorners == 3) {
2438: /* Triangle */
2439: faceSize = numCorners-1;
2440: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2441: } else if (cellDim == 3 && numCorners == 4) {
2442: /* Tetrahedron */
2443: faceSize = numCorners-1;
2444: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2445: } else if (cellDim == 1 && numCorners == 3) {
2446: /* Quadratic line */
2447: faceSize = 1;
2448: posOrient = PETSC_TRUE;
2449: } else if (cellDim == 2 && numCorners == 4) {
2450: /* Quads */
2451: faceSize = 2;
2452: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2453: posOrient = PETSC_TRUE;
2454: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2455: posOrient = PETSC_TRUE;
2456: } else {
2457: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2458: posOrient = PETSC_FALSE;
2459: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2460: }
2461: } else if (cellDim == 2 && numCorners == 6) {
2462: /* Quadratic triangle (I hate this) */
2463: /* Edges are determined by the first 2 vertices (corners of edges) */
2464: const PetscInt faceSizeTri = 3;
2465: PetscInt sortedIndices[3], i, iFace;
2466: PetscBool found = PETSC_FALSE;
2467: PetscInt faceVerticesTriSorted[9] = {
2468: 0, 3, 4, /* bottom */
2469: 1, 4, 5, /* right */
2470: 2, 3, 5, /* left */
2471: };
2472: PetscInt faceVerticesTri[9] = {
2473: 0, 3, 4, /* bottom */
2474: 1, 4, 5, /* right */
2475: 2, 5, 3, /* left */
2476: };
2478: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2479: PetscSortInt(faceSizeTri, sortedIndices);
2480: for (iFace = 0; iFace < 3; ++iFace) {
2481: const PetscInt ii = iFace*faceSizeTri;
2482: PetscInt fVertex, cVertex;
2484: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2485: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2486: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2487: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2488: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2489: faceVertices[fVertex] = origVertices[cVertex];
2490: break;
2491: }
2492: }
2493: }
2494: found = PETSC_TRUE;
2495: break;
2496: }
2497: }
2498: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2499: if (posOriented) *posOriented = PETSC_TRUE;
2500: return(0);
2501: } else if (cellDim == 2 && numCorners == 9) {
2502: /* Quadratic quad (I hate this) */
2503: /* Edges are determined by the first 2 vertices (corners of edges) */
2504: const PetscInt faceSizeQuad = 3;
2505: PetscInt sortedIndices[3], i, iFace;
2506: PetscBool found = PETSC_FALSE;
2507: PetscInt faceVerticesQuadSorted[12] = {
2508: 0, 1, 4, /* bottom */
2509: 1, 2, 5, /* right */
2510: 2, 3, 6, /* top */
2511: 0, 3, 7, /* left */
2512: };
2513: PetscInt faceVerticesQuad[12] = {
2514: 0, 1, 4, /* bottom */
2515: 1, 2, 5, /* right */
2516: 2, 3, 6, /* top */
2517: 3, 0, 7, /* left */
2518: };
2520: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2521: PetscSortInt(faceSizeQuad, sortedIndices);
2522: for (iFace = 0; iFace < 4; ++iFace) {
2523: const PetscInt ii = iFace*faceSizeQuad;
2524: PetscInt fVertex, cVertex;
2526: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2527: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2528: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2529: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2530: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2531: faceVertices[fVertex] = origVertices[cVertex];
2532: break;
2533: }
2534: }
2535: }
2536: found = PETSC_TRUE;
2537: break;
2538: }
2539: }
2540: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2541: if (posOriented) *posOriented = PETSC_TRUE;
2542: return(0);
2543: } else if (cellDim == 3 && numCorners == 8) {
2544: /* Hexes
2545: A hex is two oriented quads with the normal of the first
2546: pointing up at the second.
2548: 7---6
2549: /| /|
2550: 4---5 |
2551: | 1-|-2
2552: |/ |/
2553: 0---3
2555: Faces are determined by the first 4 vertices (corners of faces) */
2556: const PetscInt faceSizeHex = 4;
2557: PetscInt sortedIndices[4], i, iFace;
2558: PetscBool found = PETSC_FALSE;
2559: PetscInt faceVerticesHexSorted[24] = {
2560: 0, 1, 2, 3, /* bottom */
2561: 4, 5, 6, 7, /* top */
2562: 0, 3, 4, 5, /* front */
2563: 2, 3, 5, 6, /* right */
2564: 1, 2, 6, 7, /* back */
2565: 0, 1, 4, 7, /* left */
2566: };
2567: PetscInt faceVerticesHex[24] = {
2568: 1, 2, 3, 0, /* bottom */
2569: 4, 5, 6, 7, /* top */
2570: 0, 3, 5, 4, /* front */
2571: 3, 2, 6, 5, /* right */
2572: 2, 1, 7, 6, /* back */
2573: 1, 0, 4, 7, /* left */
2574: };
2576: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2577: PetscSortInt(faceSizeHex, sortedIndices);
2578: for (iFace = 0; iFace < 6; ++iFace) {
2579: const PetscInt ii = iFace*faceSizeHex;
2580: PetscInt fVertex, cVertex;
2582: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2583: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2584: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2585: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2586: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2587: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2588: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2589: faceVertices[fVertex] = origVertices[cVertex];
2590: break;
2591: }
2592: }
2593: }
2594: found = PETSC_TRUE;
2595: break;
2596: }
2597: }
2598: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2599: if (posOriented) *posOriented = PETSC_TRUE;
2600: return(0);
2601: } else if (cellDim == 3 && numCorners == 10) {
2602: /* Quadratic tet */
2603: /* Faces are determined by the first 3 vertices (corners of faces) */
2604: const PetscInt faceSizeTet = 6;
2605: PetscInt sortedIndices[6], i, iFace;
2606: PetscBool found = PETSC_FALSE;
2607: PetscInt faceVerticesTetSorted[24] = {
2608: 0, 1, 2, 6, 7, 8, /* bottom */
2609: 0, 3, 4, 6, 7, 9, /* front */
2610: 1, 4, 5, 7, 8, 9, /* right */
2611: 2, 3, 5, 6, 8, 9, /* left */
2612: };
2613: PetscInt faceVerticesTet[24] = {
2614: 0, 1, 2, 6, 7, 8, /* bottom */
2615: 0, 4, 3, 6, 7, 9, /* front */
2616: 1, 5, 4, 7, 8, 9, /* right */
2617: 2, 3, 5, 8, 6, 9, /* left */
2618: };
2620: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2621: PetscSortInt(faceSizeTet, sortedIndices);
2622: for (iFace=0; iFace < 4; ++iFace) {
2623: const PetscInt ii = iFace*faceSizeTet;
2624: PetscInt fVertex, cVertex;
2626: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2627: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2628: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2629: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2630: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2631: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2632: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2633: faceVertices[fVertex] = origVertices[cVertex];
2634: break;
2635: }
2636: }
2637: }
2638: found = PETSC_TRUE;
2639: break;
2640: }
2641: }
2642: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2643: if (posOriented) *posOriented = PETSC_TRUE;
2644: return(0);
2645: } else if (cellDim == 3 && numCorners == 27) {
2646: /* Quadratic hexes (I hate this)
2647: A hex is two oriented quads with the normal of the first
2648: pointing up at the second.
2650: 7---6
2651: /| /|
2652: 4---5 |
2653: | 3-|-2
2654: |/ |/
2655: 0---1
2657: Faces are determined by the first 4 vertices (corners of faces) */
2658: const PetscInt faceSizeQuadHex = 9;
2659: PetscInt sortedIndices[9], i, iFace;
2660: PetscBool found = PETSC_FALSE;
2661: PetscInt faceVerticesQuadHexSorted[54] = {
2662: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2663: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2664: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2665: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2666: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2667: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2668: };
2669: PetscInt faceVerticesQuadHex[54] = {
2670: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2671: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2672: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2673: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2674: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2675: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2676: };
2678: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2679: PetscSortInt(faceSizeQuadHex, sortedIndices);
2680: for (iFace = 0; iFace < 6; ++iFace) {
2681: const PetscInt ii = iFace*faceSizeQuadHex;
2682: PetscInt fVertex, cVertex;
2684: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2685: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2686: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2687: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2688: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2689: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2690: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2691: faceVertices[fVertex] = origVertices[cVertex];
2692: break;
2693: }
2694: }
2695: }
2696: found = PETSC_TRUE;
2697: break;
2698: }
2699: }
2700: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2701: if (posOriented) *posOriented = PETSC_TRUE;
2702: return(0);
2703: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2704: if (!posOrient) {
2705: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2706: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2707: } else {
2708: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2709: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2710: }
2711: if (posOriented) *posOriented = posOrient;
2712: return(0);
2713: }
2715: /*@
2716: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2717: in faceVertices. The orientation is such that the face normal points out of the cell
2719: Not collective
2721: Input Parameters:
2722: + dm - The original mesh
2723: . cell - The cell mesh point
2724: . faceSize - The number of vertices on the face
2725: . face - The face vertices
2726: . numCorners - The number of vertices on the cell
2727: . indices - Local numbering of face vertices in cell cone
2728: - origVertices - Original face vertices
2730: Output Parameter:
2731: + faceVertices - The face vertices properly oriented
2732: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2734: Level: developer
2736: .seealso: DMPlexGetCone()
2737: @*/
2738: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2739: {
2740: const PetscInt *cone = NULL;
2741: PetscInt coneSize, v, f, v2;
2742: PetscInt oppositeVertex = -1;
2743: PetscErrorCode ierr;
2746: DMPlexGetConeSize(dm, cell, &coneSize);
2747: DMPlexGetCone(dm, cell, &cone);
2748: for (v = 0, v2 = 0; v < coneSize; ++v) {
2749: PetscBool found = PETSC_FALSE;
2751: for (f = 0; f < faceSize; ++f) {
2752: if (face[f] == cone[v]) {
2753: found = PETSC_TRUE; break;
2754: }
2755: }
2756: if (found) {
2757: indices[v2] = v;
2758: origVertices[v2] = cone[v];
2759: ++v2;
2760: } else {
2761: oppositeVertex = v;
2762: }
2763: }
2764: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2765: return(0);
2766: }
2768: /*
2769: DMPlexInsertFace_Internal - Puts a face into the mesh
2771: Not collective
2773: Input Parameters:
2774: + dm - The DMPlex
2775: . numFaceVertex - The number of vertices in the face
2776: . faceVertices - The vertices in the face for dm
2777: . subfaceVertices - The vertices in the face for subdm
2778: . numCorners - The number of vertices in the cell
2779: . cell - A cell in dm containing the face
2780: . subcell - A cell in subdm containing the face
2781: . firstFace - First face in the mesh
2782: - newFacePoint - Next face in the mesh
2784: Output Parameters:
2785: . newFacePoint - Contains next face point number on input, updated on output
2787: Level: developer
2788: */
2789: 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)
2790: {
2791: MPI_Comm comm;
2792: DM_Plex *submesh = (DM_Plex*) subdm->data;
2793: const PetscInt *faces;
2794: PetscInt numFaces, coneSize;
2795: PetscErrorCode ierr;
2798: PetscObjectGetComm((PetscObject)dm,&comm);
2799: DMPlexGetConeSize(subdm, subcell, &coneSize);
2800: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2801: #if 0
2802: /* Cannot use this because support() has not been constructed yet */
2803: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2804: #else
2805: {
2806: PetscInt f;
2808: numFaces = 0;
2809: DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2810: for (f = firstFace; f < *newFacePoint; ++f) {
2811: PetscInt dof, off, d;
2813: PetscSectionGetDof(submesh->coneSection, f, &dof);
2814: PetscSectionGetOffset(submesh->coneSection, f, &off);
2815: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2816: for (d = 0; d < dof; ++d) {
2817: const PetscInt p = submesh->cones[off+d];
2818: PetscInt v;
2820: for (v = 0; v < numFaceVertices; ++v) {
2821: if (subfaceVertices[v] == p) break;
2822: }
2823: if (v == numFaceVertices) break;
2824: }
2825: if (d == dof) {
2826: numFaces = 1;
2827: ((PetscInt*) faces)[0] = f;
2828: }
2829: }
2830: }
2831: #endif
2832: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2833: else if (numFaces == 1) {
2834: /* Add the other cell neighbor for this face */
2835: DMPlexSetCone(subdm, subcell, faces);
2836: } else {
2837: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2838: PetscBool posOriented;
2840: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2841: origVertices = &orientedVertices[numFaceVertices];
2842: indices = &orientedVertices[numFaceVertices*2];
2843: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2844: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2845: /* TODO: I know that routine should return a permutation, not the indices */
2846: for (v = 0; v < numFaceVertices; ++v) {
2847: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2848: for (ov = 0; ov < numFaceVertices; ++ov) {
2849: if (orientedVertices[ov] == vertex) {
2850: orientedSubVertices[ov] = subvertex;
2851: break;
2852: }
2853: }
2854: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2855: }
2856: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2857: DMPlexSetCone(subdm, subcell, newFacePoint);
2858: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2859: ++(*newFacePoint);
2860: }
2861: #if 0
2862: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2863: #else
2864: DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2865: #endif
2866: return(0);
2867: }
2869: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2870: {
2871: MPI_Comm comm;
2872: DMLabel subpointMap;
2873: IS subvertexIS, subcellIS;
2874: const PetscInt *subVertices, *subCells;
2875: PetscInt numSubVertices, firstSubVertex, numSubCells;
2876: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2877: PetscInt vStart, vEnd, c, f;
2878: PetscErrorCode ierr;
2881: PetscObjectGetComm((PetscObject)dm,&comm);
2882: /* Create subpointMap which marks the submesh */
2883: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2884: DMPlexSetSubpointMap(subdm, subpointMap);
2885: DMLabelDestroy(&subpointMap);
2886: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2887: /* Setup chart */
2888: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2889: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2890: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2891: DMPlexSetVTKCellHeight(subdm, 1);
2892: /* Set cone sizes */
2893: firstSubVertex = numSubCells;
2894: firstSubFace = numSubCells+numSubVertices;
2895: newFacePoint = firstSubFace;
2896: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2897: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2898: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2899: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2900: for (c = 0; c < numSubCells; ++c) {
2901: DMPlexSetConeSize(subdm, c, 1);
2902: }
2903: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2904: DMPlexSetConeSize(subdm, f, nFV);
2905: }
2906: DMSetUp(subdm);
2907: /* Create face cones */
2908: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2909: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2910: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2911: for (c = 0; c < numSubCells; ++c) {
2912: const PetscInt cell = subCells[c];
2913: const PetscInt subcell = c;
2914: PetscInt *closure = NULL;
2915: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2917: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2918: for (cl = 0; cl < closureSize*2; cl += 2) {
2919: const PetscInt point = closure[cl];
2920: PetscInt subVertex;
2922: if ((point >= vStart) && (point < vEnd)) {
2923: ++numCorners;
2924: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2925: if (subVertex >= 0) {
2926: closure[faceSize] = point;
2927: subface[faceSize] = firstSubVertex+subVertex;
2928: ++faceSize;
2929: }
2930: }
2931: }
2932: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2933: if (faceSize == nFV) {
2934: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2935: }
2936: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2937: }
2938: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2939: DMPlexSymmetrize(subdm);
2940: DMPlexStratify(subdm);
2941: /* Build coordinates */
2942: {
2943: PetscSection coordSection, subCoordSection;
2944: Vec coordinates, subCoordinates;
2945: PetscScalar *coords, *subCoords;
2946: PetscInt numComp, coordSize, v;
2947: const char *name;
2949: DMGetCoordinateSection(dm, &coordSection);
2950: DMGetCoordinatesLocal(dm, &coordinates);
2951: DMGetCoordinateSection(subdm, &subCoordSection);
2952: PetscSectionSetNumFields(subCoordSection, 1);
2953: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2954: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2955: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2956: for (v = 0; v < numSubVertices; ++v) {
2957: const PetscInt vertex = subVertices[v];
2958: const PetscInt subvertex = firstSubVertex+v;
2959: PetscInt dof;
2961: PetscSectionGetDof(coordSection, vertex, &dof);
2962: PetscSectionSetDof(subCoordSection, subvertex, dof);
2963: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2964: }
2965: PetscSectionSetUp(subCoordSection);
2966: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2967: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2968: PetscObjectGetName((PetscObject)coordinates,&name);
2969: PetscObjectSetName((PetscObject)subCoordinates,name);
2970: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2971: VecSetType(subCoordinates,VECSTANDARD);
2972: if (coordSize) {
2973: VecGetArray(coordinates, &coords);
2974: VecGetArray(subCoordinates, &subCoords);
2975: for (v = 0; v < numSubVertices; ++v) {
2976: const PetscInt vertex = subVertices[v];
2977: const PetscInt subvertex = firstSubVertex+v;
2978: PetscInt dof, off, sdof, soff, d;
2980: PetscSectionGetDof(coordSection, vertex, &dof);
2981: PetscSectionGetOffset(coordSection, vertex, &off);
2982: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2983: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2984: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2985: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2986: }
2987: VecRestoreArray(coordinates, &coords);
2988: VecRestoreArray(subCoordinates, &subCoords);
2989: }
2990: DMSetCoordinatesLocal(subdm, subCoordinates);
2991: VecDestroy(&subCoordinates);
2992: }
2993: /* Cleanup */
2994: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2995: ISDestroy(&subvertexIS);
2996: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2997: ISDestroy(&subcellIS);
2998: return(0);
2999: }
3001: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
3002: {
3003: PetscInt subPoint;
3006: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
3007: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
3008: }
3010: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
3011: {
3012: MPI_Comm comm;
3013: DMLabel subpointMap;
3014: IS *subpointIS;
3015: const PetscInt **subpoints;
3016: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3017: PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v;
3018: PetscMPIInt rank;
3019: PetscErrorCode ierr;
3022: PetscObjectGetComm((PetscObject)dm,&comm);
3023: MPI_Comm_rank(comm, &rank);
3024: /* Create subpointMap which marks the submesh */
3025: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3026: DMPlexSetSubpointMap(subdm, subpointMap);
3027: if (cellHeight) {
3028: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
3029: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
3030: } else {
3031: DMLabel depth;
3032: IS pointIS;
3033: const PetscInt *points;
3034: PetscInt numPoints=0;
3036: DMPlexGetDepthLabel(dm, &depth);
3037: DMLabelGetStratumIS(label, value, &pointIS);
3038: if (pointIS) {
3039: ISGetIndices(pointIS, &points);
3040: ISGetLocalSize(pointIS, &numPoints);
3041: }
3042: for (p = 0; p < numPoints; ++p) {
3043: PetscInt *closure = NULL;
3044: PetscInt closureSize, c, pdim;
3046: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3047: for (c = 0; c < closureSize*2; c += 2) {
3048: DMLabelGetValue(depth, closure[c], &pdim);
3049: DMLabelSetValue(subpointMap, closure[c], pdim);
3050: }
3051: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3052: }
3053: if (pointIS) {ISRestoreIndices(pointIS, &points);}
3054: ISDestroy(&pointIS);
3055: }
3056: /* Setup chart */
3057: DMGetDimension(dm, &dim);
3058: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
3059: for (d = 0; d <= dim; ++d) {
3060: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3061: totSubPoints += numSubPoints[d];
3062: }
3063: DMPlexSetChart(subdm, 0, totSubPoints);
3064: DMPlexSetVTKCellHeight(subdm, cellHeight);
3065: /* Set cone sizes */
3066: firstSubPoint[dim] = 0;
3067: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
3068: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
3069: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
3070: for (d = 0; d <= dim; ++d) {
3071: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3072: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
3073: }
3074: /* We do not want this label automatically computed, instead we compute it here */
3075: DMCreateLabel(subdm, "celltype");
3076: for (d = 0; d <= dim; ++d) {
3077: for (p = 0; p < numSubPoints[d]; ++p) {
3078: const PetscInt point = subpoints[d][p];
3079: const PetscInt subpoint = firstSubPoint[d] + p;
3080: const PetscInt *cone;
3081: PetscInt coneSize, coneSizeNew, c, val;
3082: DMPolytopeType ct;
3084: DMPlexGetConeSize(dm, point, &coneSize);
3085: DMPlexSetConeSize(subdm, subpoint, coneSize);
3086: DMPlexGetCellType(dm, point, &ct);
3087: DMPlexSetCellType(subdm, subpoint, ct);
3088: if (cellHeight && (d == dim)) {
3089: DMPlexGetCone(dm, point, &cone);
3090: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3091: DMLabelGetValue(subpointMap, cone[c], &val);
3092: if (val >= 0) coneSizeNew++;
3093: }
3094: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3095: DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST);
3096: }
3097: }
3098: }
3099: DMLabelDestroy(&subpointMap);
3100: DMSetUp(subdm);
3101: /* Set cones */
3102: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3103: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3104: for (d = 0; d <= dim; ++d) {
3105: for (p = 0; p < numSubPoints[d]; ++p) {
3106: const PetscInt point = subpoints[d][p];
3107: const PetscInt subpoint = firstSubPoint[d] + p;
3108: const PetscInt *cone, *ornt;
3109: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
3111: if (d == dim-1) {
3112: const PetscInt *support, *cone, *ornt;
3113: PetscInt supportSize, coneSize, s, subc;
3115: DMPlexGetSupport(dm, point, &support);
3116: DMPlexGetSupportSize(dm, point, &supportSize);
3117: for (s = 0; s < supportSize; ++s) {
3118: PetscBool isHybrid;
3120: DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);
3121: if (!isHybrid) continue;
3122: PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3123: if (subc >= 0) {
3124: const PetscInt ccell = subpoints[d+1][subc];
3126: DMPlexGetCone(dm, ccell, &cone);
3127: DMPlexGetConeSize(dm, ccell, &coneSize);
3128: DMPlexGetConeOrientation(dm, ccell, &ornt);
3129: for (c = 0; c < coneSize; ++c) {
3130: if (cone[c] == point) {
3131: fornt = ornt[c];
3132: break;
3133: }
3134: }
3135: break;
3136: }
3137: }
3138: }
3139: DMPlexGetConeSize(dm, point, &coneSize);
3140: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3141: DMPlexGetCone(dm, point, &cone);
3142: DMPlexGetConeOrientation(dm, point, &ornt);
3143: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3144: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3145: if (subc >= 0) {
3146: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3147: orntNew[coneSizeNew] = ornt[c];
3148: ++coneSizeNew;
3149: }
3150: }
3151: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3152: if (fornt < 0) {
3153: /* This should be replaced by a call to DMPlexReverseCell() */
3154: #if 0
3155: DMPlexReverseCell(subdm, subpoint);
3156: #else
3157: for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
3158: PetscInt faceSize, tmp;
3160: tmp = coneNew[c];
3161: coneNew[c] = coneNew[coneSizeNew-1-c];
3162: coneNew[coneSizeNew-1-c] = tmp;
3163: DMPlexGetConeSize(dm, cone[c], &faceSize);
3164: tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
3165: orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
3166: orntNew[coneSizeNew-1-c] = tmp;
3167: }
3168: }
3169: DMPlexSetCone(subdm, subpoint, coneNew);
3170: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3171: #endif
3172: }
3173: }
3174: PetscFree2(coneNew,orntNew);
3175: DMPlexSymmetrize(subdm);
3176: DMPlexStratify(subdm);
3177: /* Build coordinates */
3178: {
3179: PetscSection coordSection, subCoordSection;
3180: Vec coordinates, subCoordinates;
3181: PetscScalar *coords, *subCoords;
3182: PetscInt cdim, numComp, coordSize;
3183: const char *name;
3185: DMGetCoordinateDim(dm, &cdim);
3186: DMGetCoordinateSection(dm, &coordSection);
3187: DMGetCoordinatesLocal(dm, &coordinates);
3188: DMGetCoordinateSection(subdm, &subCoordSection);
3189: PetscSectionSetNumFields(subCoordSection, 1);
3190: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3191: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3192: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3193: for (v = 0; v < numSubPoints[0]; ++v) {
3194: const PetscInt vertex = subpoints[0][v];
3195: const PetscInt subvertex = firstSubPoint[0]+v;
3196: PetscInt dof;
3198: PetscSectionGetDof(coordSection, vertex, &dof);
3199: PetscSectionSetDof(subCoordSection, subvertex, dof);
3200: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3201: }
3202: PetscSectionSetUp(subCoordSection);
3203: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3204: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3205: PetscObjectGetName((PetscObject)coordinates,&name);
3206: PetscObjectSetName((PetscObject)subCoordinates,name);
3207: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3208: VecSetBlockSize(subCoordinates, cdim);
3209: VecSetType(subCoordinates,VECSTANDARD);
3210: VecGetArray(coordinates, &coords);
3211: VecGetArray(subCoordinates, &subCoords);
3212: for (v = 0; v < numSubPoints[0]; ++v) {
3213: const PetscInt vertex = subpoints[0][v];
3214: const PetscInt subvertex = firstSubPoint[0]+v;
3215: PetscInt dof, off, sdof, soff, d;
3217: PetscSectionGetDof(coordSection, vertex, &dof);
3218: PetscSectionGetOffset(coordSection, vertex, &off);
3219: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3220: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3221: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3222: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3223: }
3224: VecRestoreArray(coordinates, &coords);
3225: VecRestoreArray(subCoordinates, &subCoords);
3226: DMSetCoordinatesLocal(subdm, subCoordinates);
3227: VecDestroy(&subCoordinates);
3228: }
3229: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3230: {
3231: PetscSF sfPoint, sfPointSub;
3232: IS subpIS;
3233: const PetscSFNode *remotePoints;
3234: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3235: const PetscInt *localPoints, *subpoints;
3236: PetscInt *slocalPoints;
3237: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3238: PetscMPIInt rank;
3240: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3241: DMGetPointSF(dm, &sfPoint);
3242: DMGetPointSF(subdm, &sfPointSub);
3243: DMPlexGetChart(dm, &pStart, &pEnd);
3244: DMPlexGetChart(subdm, NULL, &numSubroots);
3245: DMPlexGetSubpointIS(subdm, &subpIS);
3246: if (subpIS) {
3247: ISGetIndices(subpIS, &subpoints);
3248: ISGetLocalSize(subpIS, &numSubpoints);
3249: }
3250: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3251: if (numRoots >= 0) {
3252: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3253: for (p = 0; p < pEnd-pStart; ++p) {
3254: newLocalPoints[p].rank = -2;
3255: newLocalPoints[p].index = -2;
3256: }
3257: /* Set subleaves */
3258: for (l = 0; l < numLeaves; ++l) {
3259: const PetscInt point = localPoints[l];
3260: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3262: if (subpoint < 0) continue;
3263: newLocalPoints[point-pStart].rank = rank;
3264: newLocalPoints[point-pStart].index = subpoint;
3265: ++numSubleaves;
3266: }
3267: /* Must put in owned subpoints */
3268: for (p = pStart; p < pEnd; ++p) {
3269: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3271: if (subpoint < 0) {
3272: newOwners[p-pStart].rank = -3;
3273: newOwners[p-pStart].index = -3;
3274: } else {
3275: newOwners[p-pStart].rank = rank;
3276: newOwners[p-pStart].index = subpoint;
3277: }
3278: }
3279: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3280: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3281: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
3282: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
3283: PetscMalloc1(numSubleaves, &slocalPoints);
3284: PetscMalloc1(numSubleaves, &sremotePoints);
3285: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3286: const PetscInt point = localPoints[l];
3287: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3289: if (subpoint < 0) continue;
3290: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3291: slocalPoints[sl] = subpoint;
3292: sremotePoints[sl].rank = newLocalPoints[point].rank;
3293: sremotePoints[sl].index = newLocalPoints[point].index;
3294: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3295: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3296: ++sl;
3297: }
3298: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3299: PetscFree2(newLocalPoints,newOwners);
3300: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3301: }
3302: if (subpIS) {
3303: ISRestoreIndices(subpIS, &subpoints);
3304: }
3305: }
3306: /* Cleanup */
3307: for (d = 0; d <= dim; ++d) {
3308: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3309: ISDestroy(&subpointIS[d]);
3310: }
3311: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3312: return(0);
3313: }
3315: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3316: {
3320: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3321: return(0);
3322: }
3324: /*@
3325: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3327: Input Parameters:
3328: + dm - The original mesh
3329: . vertexLabel - The DMLabel marking points contained in the surface
3330: . value - The label value to use
3331: - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked
3333: Output Parameter:
3334: . subdm - The surface mesh
3336: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3338: Level: developer
3340: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3341: @*/
3342: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3343: {
3344: DMPlexInterpolatedFlag interpolated;
3345: PetscInt dim, cdim;
3351: DMGetDimension(dm, &dim);
3352: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3353: DMSetType(*subdm, DMPLEX);
3354: DMSetDimension(*subdm, dim-1);
3355: DMGetCoordinateDim(dm, &cdim);
3356: DMSetCoordinateDim(*subdm, cdim);
3357: DMPlexIsInterpolated(dm, &interpolated);
3358: if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
3359: if (interpolated) {
3360: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3361: } else {
3362: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3363: }
3364: return(0);
3365: }
3367: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3368: {
3369: MPI_Comm comm;
3370: DMLabel subpointMap;
3371: IS subvertexIS;
3372: const PetscInt *subVertices;
3373: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3374: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3375: PetscInt c, f;
3376: PetscErrorCode ierr;
3379: PetscObjectGetComm((PetscObject)dm, &comm);
3380: /* Create subpointMap which marks the submesh */
3381: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3382: DMPlexSetSubpointMap(subdm, subpointMap);
3383: DMLabelDestroy(&subpointMap);
3384: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3385: /* Setup chart */
3386: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3387: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3388: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3389: DMPlexSetVTKCellHeight(subdm, 1);
3390: /* Set cone sizes */
3391: firstSubVertex = numSubCells;
3392: firstSubFace = numSubCells+numSubVertices;
3393: newFacePoint = firstSubFace;
3394: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3395: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3396: for (c = 0; c < numSubCells; ++c) {
3397: DMPlexSetConeSize(subdm, c, 1);
3398: }
3399: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3400: DMPlexSetConeSize(subdm, f, nFV);
3401: }
3402: DMSetUp(subdm);
3403: /* Create face cones */
3404: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3405: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3406: for (c = 0; c < numSubCells; ++c) {
3407: const PetscInt cell = subCells[c];
3408: const PetscInt subcell = c;
3409: const PetscInt *cone, *cells;
3410: PetscBool isHybrid;
3411: PetscInt numCells, subVertex, p, v;
3413: DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid);
3414: if (!isHybrid) continue;
3415: DMPlexGetCone(dm, cell, &cone);
3416: for (v = 0; v < nFV; ++v) {
3417: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3418: subface[v] = firstSubVertex+subVertex;
3419: }
3420: DMPlexSetCone(subdm, newFacePoint, subface);
3421: DMPlexSetCone(subdm, subcell, &newFacePoint);
3422: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3423: /* Not true in parallel
3424: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3425: for (p = 0; p < numCells; ++p) {
3426: PetscInt negsubcell;
3427: PetscBool isHybrid;
3429: DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid);
3430: if (isHybrid) continue;
3431: /* I know this is a crap search */
3432: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3433: if (subCells[negsubcell] == cells[p]) break;
3434: }
3435: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3436: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3437: }
3438: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3439: ++newFacePoint;
3440: }
3441: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3442: DMPlexSymmetrize(subdm);
3443: DMPlexStratify(subdm);
3444: /* Build coordinates */
3445: {
3446: PetscSection coordSection, subCoordSection;
3447: Vec coordinates, subCoordinates;
3448: PetscScalar *coords, *subCoords;
3449: PetscInt cdim, numComp, coordSize, v;
3450: const char *name;
3452: DMGetCoordinateDim(dm, &cdim);
3453: DMGetCoordinateSection(dm, &coordSection);
3454: DMGetCoordinatesLocal(dm, &coordinates);
3455: DMGetCoordinateSection(subdm, &subCoordSection);
3456: PetscSectionSetNumFields(subCoordSection, 1);
3457: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3458: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3459: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3460: for (v = 0; v < numSubVertices; ++v) {
3461: const PetscInt vertex = subVertices[v];
3462: const PetscInt subvertex = firstSubVertex+v;
3463: PetscInt dof;
3465: PetscSectionGetDof(coordSection, vertex, &dof);
3466: PetscSectionSetDof(subCoordSection, subvertex, dof);
3467: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3468: }
3469: PetscSectionSetUp(subCoordSection);
3470: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3471: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3472: PetscObjectGetName((PetscObject)coordinates,&name);
3473: PetscObjectSetName((PetscObject)subCoordinates,name);
3474: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3475: VecSetBlockSize(subCoordinates, cdim);
3476: VecSetType(subCoordinates,VECSTANDARD);
3477: VecGetArray(coordinates, &coords);
3478: VecGetArray(subCoordinates, &subCoords);
3479: for (v = 0; v < numSubVertices; ++v) {
3480: const PetscInt vertex = subVertices[v];
3481: const PetscInt subvertex = firstSubVertex+v;
3482: PetscInt dof, off, sdof, soff, d;
3484: PetscSectionGetDof(coordSection, vertex, &dof);
3485: PetscSectionGetOffset(coordSection, vertex, &off);
3486: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3487: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3488: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3489: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3490: }
3491: VecRestoreArray(coordinates, &coords);
3492: VecRestoreArray(subCoordinates, &subCoords);
3493: DMSetCoordinatesLocal(subdm, subCoordinates);
3494: VecDestroy(&subCoordinates);
3495: }
3496: /* Build SF */
3497: CHKMEMQ;
3498: {
3499: PetscSF sfPoint, sfPointSub;
3500: const PetscSFNode *remotePoints;
3501: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3502: const PetscInt *localPoints;
3503: PetscInt *slocalPoints;
3504: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3505: PetscMPIInt rank;
3507: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3508: DMGetPointSF(dm, &sfPoint);
3509: DMGetPointSF(subdm, &sfPointSub);
3510: DMPlexGetChart(dm, &pStart, &pEnd);
3511: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3512: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3513: if (numRoots >= 0) {
3514: /* Only vertices should be shared */
3515: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3516: for (p = 0; p < pEnd-pStart; ++p) {
3517: newLocalPoints[p].rank = -2;
3518: newLocalPoints[p].index = -2;
3519: }
3520: /* Set subleaves */
3521: for (l = 0; l < numLeaves; ++l) {
3522: const PetscInt point = localPoints[l];
3523: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3525: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3526: if (subPoint < 0) continue;
3527: newLocalPoints[point-pStart].rank = rank;
3528: newLocalPoints[point-pStart].index = subPoint;
3529: ++numSubLeaves;
3530: }
3531: /* Must put in owned subpoints */
3532: for (p = pStart; p < pEnd; ++p) {
3533: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3535: if (subPoint < 0) {
3536: newOwners[p-pStart].rank = -3;
3537: newOwners[p-pStart].index = -3;
3538: } else {
3539: newOwners[p-pStart].rank = rank;
3540: newOwners[p-pStart].index = subPoint;
3541: }
3542: }
3543: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3544: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3545: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
3546: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);
3547: PetscMalloc1(numSubLeaves, &slocalPoints);
3548: PetscMalloc1(numSubLeaves, &sremotePoints);
3549: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3550: const PetscInt point = localPoints[l];
3551: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3553: if (subPoint < 0) continue;
3554: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3555: slocalPoints[sl] = subPoint;
3556: sremotePoints[sl].rank = newLocalPoints[point].rank;
3557: sremotePoints[sl].index = newLocalPoints[point].index;
3558: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3559: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3560: ++sl;
3561: }
3562: PetscFree2(newLocalPoints,newOwners);
3563: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3564: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3565: }
3566: }
3567: CHKMEMQ;
3568: /* Cleanup */
3569: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3570: ISDestroy(&subvertexIS);
3571: PetscFree(subCells);
3572: return(0);
3573: }
3575: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3576: {
3577: DMLabel label = NULL;
3581: if (labelname) {DMGetLabel(dm, labelname, &label);}
3582: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3583: return(0);
3584: }
3586: /*@C
3587: 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.
3589: Input Parameters:
3590: + dm - The original mesh
3591: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3592: . label - A label name, or NULL
3593: - value - A label value
3595: Output Parameter:
3596: . subdm - The surface mesh
3598: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3600: Level: developer
3602: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3603: @*/
3604: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3605: {
3606: PetscInt dim, cdim, depth;
3612: DMGetDimension(dm, &dim);
3613: DMPlexGetDepth(dm, &depth);
3614: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3615: DMSetType(*subdm, DMPLEX);
3616: DMSetDimension(*subdm, dim-1);
3617: DMGetCoordinateDim(dm, &cdim);
3618: DMSetCoordinateDim(*subdm, cdim);
3619: if (depth == dim) {
3620: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3621: } else {
3622: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3623: }
3624: return(0);
3625: }
3627: /*@
3628: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3630: Input Parameters:
3631: + dm - The original mesh
3632: . cellLabel - The DMLabel marking cells contained in the new mesh
3633: - value - The label value to use
3635: Output Parameter:
3636: . subdm - The new mesh
3638: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3640: Level: developer
3642: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3643: @*/
3644: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3645: {
3646: PetscInt dim;
3652: DMGetDimension(dm, &dim);
3653: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3654: DMSetType(*subdm, DMPLEX);
3655: DMSetDimension(*subdm, dim);
3656: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3657: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3658: return(0);
3659: }
3661: /*@
3662: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3664: Input Parameter:
3665: . dm - The submesh DM
3667: Output Parameter:
3668: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3670: Level: developer
3672: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS()
3673: @*/
3674: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3675: {
3679: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3680: return(0);
3681: }
3683: /*@
3684: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3686: Input Parameters:
3687: + dm - The submesh DM
3688: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3690: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3692: Level: developer
3694: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS()
3695: @*/
3696: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3697: {
3698: DM_Plex *mesh = (DM_Plex *) dm->data;
3699: DMLabel tmp;
3704: tmp = mesh->subpointMap;
3705: mesh->subpointMap = subpointMap;
3706: PetscObjectReference((PetscObject) mesh->subpointMap);
3707: DMLabelDestroy(&tmp);
3708: return(0);
3709: }
3711: static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS)
3712: {
3713: DMLabel spmap;
3714: PetscInt depth, d;
3718: DMPlexGetSubpointMap(dm, &spmap);
3719: DMPlexGetDepth(dm, &depth);
3720: if (spmap && depth >= 0) {
3721: DM_Plex *mesh = (DM_Plex *) dm->data;
3722: PetscInt *points, *depths;
3723: PetscInt pStart, pEnd, p, off;
3725: DMPlexGetChart(dm, &pStart, &pEnd);
3726: if (pStart) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3727: PetscMalloc1(pEnd, &points);
3728: DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3729: depths[0] = depth;
3730: depths[1] = 0;
3731: for (d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3732: for (d = 0, off = 0; d <= depth; ++d) {
3733: const PetscInt dep = depths[d];
3734: PetscInt depStart, depEnd, n;
3736: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3737: DMLabelGetStratumSize(spmap, dep, &n);
3738: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3739: 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);
3740: } else {
3741: if (!n) {
3742: if (d == 0) {
3743: /* Missing cells */
3744: for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3745: } else {
3746: /* Missing faces */
3747: for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3748: }
3749: }
3750: }
3751: if (n) {
3752: IS is;
3753: const PetscInt *opoints;
3755: DMLabelGetStratumIS(spmap, dep, &is);
3756: ISGetIndices(is, &opoints);
3757: for (p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3758: ISRestoreIndices(is, &opoints);
3759: ISDestroy(&is);
3760: }
3761: }
3762: DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3763: if (off != pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3764: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3765: PetscObjectStateGet((PetscObject) spmap, &mesh->subpointState);
3766: }
3767: return(0);
3768: }
3770: /*@
3771: DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data
3773: Input Parameter:
3774: . dm - The submesh DM
3776: Output Parameter:
3777: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3779: Note: This IS is guaranteed to be sorted by the construction of the submesh
3781: Level: developer
3783: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3784: @*/
3785: PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS)
3786: {
3787: DM_Plex *mesh = (DM_Plex *) dm->data;
3788: DMLabel spmap;
3789: PetscObjectState state;
3790: PetscErrorCode ierr;
3795: DMPlexGetSubpointMap(dm, &spmap);
3796: PetscObjectStateGet((PetscObject) spmap, &state);
3797: if (state != mesh->subpointState || !mesh->subpointIS) {DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS);}
3798: *subpointIS = mesh->subpointIS;
3799: return(0);
3800: }
3802: /*@
3803: DMGetEnclosureRelation - Get the relationship between dmA and dmB
3805: Input Parameters:
3806: + dmA - The first DM
3807: - dmB - The second DM
3809: Output Parameter:
3810: . rel - The relation of dmA to dmB
3812: Level: intermediate
3814: .seealso: DMGetEnclosurePoint()
3815: @*/
3816: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
3817: {
3818: DM plexA, plexB, sdm;
3819: DMLabel spmap;
3820: PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB;
3825: *rel = DM_ENC_NONE;
3826: if (!dmA || !dmB) return(0);
3829: if (dmA == dmB) {*rel = DM_ENC_EQUALITY; return(0);}
3830: DMConvert(dmA, DMPLEX, &plexA);
3831: DMConvert(dmB, DMPLEX, &plexB);
3832: DMPlexGetChart(plexA, &pStartA, &pEndA);
3833: DMPlexGetChart(plexB, &pStartB, &pEndB);
3834: /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
3835: - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
3836: if ((pStartA == pStartB) && (pEndA == pEndB)) {
3837: *rel = DM_ENC_EQUALITY;
3838: goto end;
3839: }
3840: NpA = pEndA - pStartA;
3841: NpB = pEndB - pStartB;
3842: if (NpA == NpB) goto end;
3843: sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
3844: DMPlexGetSubpointMap(sdm, &spmap);
3845: if (!spmap) goto end;
3846: /* TODO Check the space mapped to by subpointMap is same size as dm */
3847: if (NpA > NpB) {
3848: *rel = DM_ENC_SUPERMESH;
3849: } else {
3850: *rel = DM_ENC_SUBMESH;
3851: }
3852: end:
3853: DMDestroy(&plexA);
3854: DMDestroy(&plexB);
3855: return(0);
3856: }
3858: /*@
3859: DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB
3861: Input Parameters:
3862: + dmA - The first DM
3863: . dmB - The second DM
3864: . etype - The type of enclosure relation that dmA has to dmB
3865: - pB - A point of dmB
3867: Output Parameter:
3868: . pA - The corresponding point of dmA
3870: Level: intermediate
3872: .seealso: DMGetEnclosureRelation()
3873: @*/
3874: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
3875: {
3876: DM sdm;
3877: IS subpointIS;
3878: const PetscInt *subpoints;
3879: PetscInt numSubpoints;
3880: PetscErrorCode ierr;
3883: /* TODO Cache the IS, making it look like an index */
3884: switch (etype) {
3885: case DM_ENC_SUPERMESH:
3886: sdm = dmB;
3887: DMPlexGetSubpointIS(sdm, &subpointIS);
3888: ISGetIndices(subpointIS, &subpoints);
3889: *pA = subpoints[pB];
3890: ISRestoreIndices(subpointIS, &subpoints);
3891: break;
3892: case DM_ENC_SUBMESH:
3893: sdm = dmA;
3894: DMPlexGetSubpointIS(sdm, &subpointIS);
3895: ISGetLocalSize(subpointIS, &numSubpoints);
3896: ISGetIndices(subpointIS, &subpoints);
3897: PetscFindInt(pB, numSubpoints, subpoints, pA);
3898: if (*pA < 0) {
3899: DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
3900: DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
3901: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB);
3902: }
3903: ISRestoreIndices(subpointIS, &subpoints);
3904: break;
3905: case DM_ENC_EQUALITY:
3906: case DM_ENC_NONE:
3907: *pA = pB;break;
3908: case DM_ENC_UNKNOWN:
3909: {
3910: DMEnclosureType enc;
3912: DMGetEnclosureRelation(dmA, dmB, &enc);
3913: DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
3914: }
3915: break;
3916: default: SETERRQ1(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype);
3917: }
3918: return(0);
3919: }