Actual source code: plexsubmesh.c
petsc-3.13.6 2020-09-29
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);
622: PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
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: DMPlexCreateSubpointIS(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) {
1742: ISDestroy(&subpointIS);
1743: return(0);
1744: }
1745: ISGetLocalSize(dimIS, &numPoints);
1746: ISGetIndices(dimIS, &points);
1747: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1748: const PetscInt *support;
1749: PetscInt supportSize, s;
1751: DMPlexGetSupportSize(dm, points[p], &supportSize);
1752: #if 0
1753: if (supportSize != 2) {
1754: const PetscInt *lp;
1755: PetscInt Nlp, pind;
1757: /* Check that for a cell with a single support face, that face is in the SF */
1758: /* THis check only works for the remote side. We would need root side information */
1759: PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1760: PetscFindInt(points[p], Nlp, lp, &pind);
1761: 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);
1762: }
1763: #endif
1764: DMPlexGetSupport(dm, points[p], &support);
1765: for (s = 0; s < supportSize; ++s) {
1766: const PetscInt *cone;
1767: PetscInt coneSize, c;
1768: PetscBool pos;
1770: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1771: if (pos) {DMLabelSetValue(label, support[s], rev*(shift+dim));}
1772: else {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1773: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1774: /* Put faces touching the fault in the label */
1775: DMPlexGetConeSize(dm, support[s], &coneSize);
1776: DMPlexGetCone(dm, support[s], &cone);
1777: for (c = 0; c < coneSize; ++c) {
1778: const PetscInt point = cone[c];
1780: DMLabelGetValue(label, point, &val);
1781: if (val == -1) {
1782: PetscInt *closure = NULL;
1783: PetscInt closureSize, cl;
1785: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1786: for (cl = 0; cl < closureSize*2; cl += 2) {
1787: const PetscInt clp = closure[cl];
1788: PetscInt bval = -1;
1790: DMLabelGetValue(label, clp, &val);
1791: if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1792: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1793: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1794: break;
1795: }
1796: }
1797: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1798: }
1799: }
1800: }
1801: }
1802: ISRestoreIndices(dimIS, &points);
1803: ISDestroy(&dimIS);
1804: if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1805: ISDestroy(&subpointIS);
1806: /* Mark boundary points as unsplit */
1807: if (blabel) {
1808: DMLabelGetStratumIS(blabel, 1, &dimIS);
1809: ISGetLocalSize(dimIS, &numPoints);
1810: ISGetIndices(dimIS, &points);
1811: for (p = 0; p < numPoints; ++p) {
1812: const PetscInt point = points[p];
1813: PetscInt val, bval;
1815: DMLabelGetValue(blabel, point, &bval);
1816: if (bval >= 0) {
1817: DMLabelGetValue(label, point, &val);
1818: if ((val < 0) || (val > dim)) {
1819: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1820: DMLabelClearValue(blabel, point, bval);
1821: }
1822: }
1823: }
1824: for (p = 0; p < numPoints; ++p) {
1825: const PetscInt point = points[p];
1826: PetscInt val, bval;
1828: DMLabelGetValue(blabel, point, &bval);
1829: if (bval >= 0) {
1830: const PetscInt *cone, *support;
1831: PetscInt coneSize, supportSize, s, valA, valB, valE;
1833: /* Mark as unsplit */
1834: DMLabelGetValue(label, point, &val);
1835: 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);
1836: DMLabelClearValue(label, point, val);
1837: DMLabelSetValue(label, point, shift2+val);
1838: /* Check for cross-edge
1839: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1840: if (val != 0) continue;
1841: DMPlexGetSupport(dm, point, &support);
1842: DMPlexGetSupportSize(dm, point, &supportSize);
1843: for (s = 0; s < supportSize; ++s) {
1844: DMPlexGetCone(dm, support[s], &cone);
1845: DMPlexGetConeSize(dm, support[s], &coneSize);
1846: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1847: DMLabelGetValue(blabel, cone[0], &valA);
1848: DMLabelGetValue(blabel, cone[1], &valB);
1849: DMLabelGetValue(blabel, support[s], &valE);
1850: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1851: }
1852: }
1853: }
1854: ISRestoreIndices(dimIS, &points);
1855: ISDestroy(&dimIS);
1856: }
1857: /* Search for other cells/faces/edges connected to the fault by a vertex */
1858: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1859: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1860: DMLabelGetStratumIS(label, 0, &dimIS);
1861: if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1862: if (dimIS && crossEdgeIS) {
1863: IS vertIS = dimIS;
1865: ISExpand(vertIS, crossEdgeIS, &dimIS);
1866: ISDestroy(&crossEdgeIS);
1867: ISDestroy(&vertIS);
1868: }
1869: if (!dimIS) {
1870: return(0);
1871: }
1872: ISGetLocalSize(dimIS, &numPoints);
1873: ISGetIndices(dimIS, &points);
1874: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1875: PetscInt *star = NULL;
1876: PetscInt starSize, s;
1877: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1879: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1880: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1881: while (again) {
1882: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1883: again = 0;
1884: for (s = 0; s < starSize*2; s += 2) {
1885: const PetscInt point = star[s];
1886: const PetscInt *cone;
1887: PetscInt coneSize, c;
1889: if ((point < cStart) || (point >= cEnd)) continue;
1890: DMLabelGetValue(label, point, &val);
1891: if (val != -1) continue;
1892: again = again == 1 ? 1 : 2;
1893: DMPlexGetConeSize(dm, point, &coneSize);
1894: DMPlexGetCone(dm, point, &cone);
1895: for (c = 0; c < coneSize; ++c) {
1896: DMLabelGetValue(label, cone[c], &val);
1897: if (val != -1) {
1898: const PetscInt *ccone;
1899: PetscInt cconeSize, cc, side;
1901: 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);
1902: if (val > 0) side = 1;
1903: else side = -1;
1904: DMLabelSetValue(label, point, side*(shift+dim));
1905: /* Mark cell faces which touch the fault */
1906: DMPlexGetConeSize(dm, point, &cconeSize);
1907: DMPlexGetCone(dm, point, &ccone);
1908: for (cc = 0; cc < cconeSize; ++cc) {
1909: PetscInt *closure = NULL;
1910: PetscInt closureSize, cl;
1912: DMLabelGetValue(label, ccone[cc], &val);
1913: if (val != -1) continue;
1914: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1915: for (cl = 0; cl < closureSize*2; cl += 2) {
1916: const PetscInt clp = closure[cl];
1918: DMLabelGetValue(label, clp, &val);
1919: if (val == -1) continue;
1920: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1921: break;
1922: }
1923: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1924: }
1925: again = 1;
1926: break;
1927: }
1928: }
1929: }
1930: }
1931: /* Classify the rest by cell membership */
1932: for (s = 0; s < starSize*2; s += 2) {
1933: const PetscInt point = star[s];
1935: DMLabelGetValue(label, point, &val);
1936: if (val == -1) {
1937: PetscInt *sstar = NULL;
1938: PetscInt sstarSize, ss;
1939: PetscBool marked = PETSC_FALSE, isHybrid;
1941: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1942: for (ss = 0; ss < sstarSize*2; ss += 2) {
1943: const PetscInt spoint = sstar[ss];
1945: if ((spoint < cStart) || (spoint >= cEnd)) continue;
1946: DMLabelGetValue(label, spoint, &val);
1947: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1948: DMLabelGetValue(depthLabel, point, &dep);
1949: if (val > 0) {
1950: DMLabelSetValue(label, point, shift+dep);
1951: } else {
1952: DMLabelSetValue(label, point, -(shift+dep));
1953: }
1954: marked = PETSC_TRUE;
1955: break;
1956: }
1957: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1958: DMPlexCellIsHybrid_Internal(dm, point, &isHybrid);
1959: if (!isHybrid && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1960: }
1961: }
1962: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1963: }
1964: ISRestoreIndices(dimIS, &points);
1965: ISDestroy(&dimIS);
1966: /* If any faces touching the fault divide cells on either side, split them */
1967: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1968: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1969: ISExpand(facePosIS, faceNegIS, &dimIS);
1970: ISDestroy(&facePosIS);
1971: ISDestroy(&faceNegIS);
1972: ISGetLocalSize(dimIS, &numPoints);
1973: ISGetIndices(dimIS, &points);
1974: for (p = 0; p < numPoints; ++p) {
1975: const PetscInt point = points[p];
1976: const PetscInt *support;
1977: PetscInt supportSize, valA, valB;
1979: DMPlexGetSupportSize(dm, point, &supportSize);
1980: if (supportSize != 2) continue;
1981: DMPlexGetSupport(dm, point, &support);
1982: DMLabelGetValue(label, support[0], &valA);
1983: DMLabelGetValue(label, support[1], &valB);
1984: if ((valA == -1) || (valB == -1)) continue;
1985: if (valA*valB > 0) continue;
1986: /* Split the face */
1987: DMLabelGetValue(label, point, &valA);
1988: DMLabelClearValue(label, point, valA);
1989: DMLabelSetValue(label, point, dim-1);
1990: /* Label its closure:
1991: unmarked: label as unsplit
1992: incident: relabel as split
1993: split: do nothing
1994: */
1995: {
1996: PetscInt *closure = NULL;
1997: PetscInt closureSize, cl;
1999: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2000: for (cl = 0; cl < closureSize*2; cl += 2) {
2001: DMLabelGetValue(label, closure[cl], &valA);
2002: if (valA == -1) { /* Mark as unsplit */
2003: DMLabelGetValue(depthLabel, closure[cl], &dep);
2004: DMLabelSetValue(label, closure[cl], shift2+dep);
2005: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
2006: DMLabelGetValue(depthLabel, closure[cl], &dep);
2007: DMLabelClearValue(label, closure[cl], valA);
2008: DMLabelSetValue(label, closure[cl], dep);
2009: }
2010: }
2011: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2012: }
2013: }
2014: ISRestoreIndices(dimIS, &points);
2015: ISDestroy(&dimIS);
2016: return(0);
2017: }
2019: /* Check that no cell have all vertices on the fault */
2020: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
2021: {
2022: IS subpointIS;
2023: const PetscInt *dmpoints;
2024: PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd;
2025: PetscErrorCode ierr;
2028: if (!label) return(0);
2029: DMLabelGetDefaultValue(label, &defaultValue);
2030: DMPlexCreateSubpointIS(subdm, &subpointIS);
2031: if (!subpointIS) return(0);
2032: DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
2033: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2034: ISGetIndices(subpointIS, &dmpoints);
2035: for (c = cStart; c < cEnd; ++c) {
2036: PetscBool invalidCell = PETSC_TRUE;
2037: PetscInt *closure = NULL;
2038: PetscInt closureSize, cl;
2040: DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2041: for (cl = 0; cl < closureSize*2; cl += 2) {
2042: PetscInt value = 0;
2044: if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
2045: DMLabelGetValue(label, closure[cl], &value);
2046: if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
2047: }
2048: DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2049: if (invalidCell) {
2050: ISRestoreIndices(subpointIS, &dmpoints);
2051: ISDestroy(&subpointIS);
2052: DMDestroy(&subdm);
2053: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
2054: }
2055: }
2056: ISRestoreIndices(subpointIS, &dmpoints);
2057: ISDestroy(&subpointIS);
2058: return(0);
2059: }
2062: /*@
2063: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
2065: Collective on dm
2067: Input Parameters:
2068: + dm - The original DM
2069: . label - The label specifying the interface vertices
2070: - bdlabel - The optional label specifying the interface boundary vertices
2072: Output Parameters:
2073: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
2074: . splitLabel - The label containing the split points, or NULL if no output is desired
2075: . dmInterface - The new interface DM, or NULL
2076: - dmHybrid - The new DM with cohesive cells
2078: Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points
2079: directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
2080: 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
2081: one vertex 3 on the surface would be 6 (101) and 3 (0) in hybridLabel. If an edge 9 from the negative side of the
2082: surface also hits vertex 3, it would be 9 (-101) in hybridLabel.
2084: The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
2085: mesh. The label value is +=100+dim for each point. For example, if two edges 10 and 14 in the hybrid resulting from
2086: splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.
2088: The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using
2089: DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh.
2091: Level: developer
2093: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
2094: @*/
2095: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2096: {
2097: DM idm;
2098: DMLabel subpointMap, hlabel, slabel = NULL;
2099: PetscInt dim;
2109: DMGetDimension(dm, &dim);
2110: DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2111: DMPlexCheckValidSubmesh_Private(dm, label, idm);
2112: DMPlexOrient(idm);
2113: DMPlexGetSubpointMap(idm, &subpointMap);
2114: DMLabelDuplicate(subpointMap, &hlabel);
2115: DMLabelClearStratum(hlabel, dim);
2116: if (splitLabel) {
2117: const char *name;
2118: char sname[PETSC_MAX_PATH_LEN];
2120: PetscObjectGetName((PetscObject) hlabel, &name);
2121: PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2122: PetscStrcat(sname, " split");
2123: DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2124: }
2125: DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2126: if (dmInterface) {*dmInterface = idm;}
2127: else {DMDestroy(&idm);}
2128: DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2129: if (hybridLabel) *hybridLabel = hlabel;
2130: else {DMLabelDestroy(&hlabel);}
2131: if (splitLabel) *splitLabel = slabel;
2132: return(0);
2133: }
2135: /* Here we need the explicit assumption that:
2137: For any marked cell, the marked vertices constitute a single face
2138: */
2139: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2140: {
2141: IS subvertexIS = NULL;
2142: const PetscInt *subvertices;
2143: PetscInt *pStart, *pEnd, pSize;
2144: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
2145: PetscErrorCode ierr;
2148: *numFaces = 0;
2149: *nFV = 0;
2150: DMPlexGetDepth(dm, &depth);
2151: DMGetDimension(dm, &dim);
2152: pSize = PetscMax(depth, dim) + 1;
2153: PetscMalloc2(pSize, &pStart, pSize, &pEnd);
2154: for (d = 0; d <= depth; ++d) {
2155: DMPlexGetSimplexOrBoxCells(dm, depth-d, &pStart[d], &pEnd[d]);
2156: }
2157: /* Loop over initial vertices and mark all faces in the collective star() */
2158: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
2159: if (subvertexIS) {
2160: ISGetSize(subvertexIS, &numSubVerticesInitial);
2161: ISGetIndices(subvertexIS, &subvertices);
2162: }
2163: for (v = 0; v < numSubVerticesInitial; ++v) {
2164: const PetscInt vertex = subvertices[v];
2165: PetscInt *star = NULL;
2166: PetscInt starSize, s, numCells = 0, c;
2168: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2169: for (s = 0; s < starSize*2; s += 2) {
2170: const PetscInt point = star[s];
2171: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2172: }
2173: for (c = 0; c < numCells; ++c) {
2174: const PetscInt cell = star[c];
2175: PetscInt *closure = NULL;
2176: PetscInt closureSize, cl;
2177: PetscInt cellLoc, numCorners = 0, faceSize = 0;
2179: DMLabelGetValue(subpointMap, cell, &cellLoc);
2180: if (cellLoc == 2) continue;
2181: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2182: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2183: for (cl = 0; cl < closureSize*2; cl += 2) {
2184: const PetscInt point = closure[cl];
2185: PetscInt vertexLoc;
2187: if ((point >= pStart[0]) && (point < pEnd[0])) {
2188: ++numCorners;
2189: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2190: if (vertexLoc == value) closure[faceSize++] = point;
2191: }
2192: }
2193: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2194: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2195: if (faceSize == *nFV) {
2196: const PetscInt *cells = NULL;
2197: PetscInt numCells, nc;
2199: ++(*numFaces);
2200: for (cl = 0; cl < faceSize; ++cl) {
2201: DMLabelSetValue(subpointMap, closure[cl], 0);
2202: }
2203: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2204: for (nc = 0; nc < numCells; ++nc) {
2205: DMLabelSetValue(subpointMap, cells[nc], 2);
2206: }
2207: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2208: }
2209: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2210: }
2211: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2212: }
2213: if (subvertexIS) {
2214: ISRestoreIndices(subvertexIS, &subvertices);
2215: }
2216: ISDestroy(&subvertexIS);
2217: PetscFree2(pStart, pEnd);
2218: return(0);
2219: }
2221: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2222: {
2223: IS subvertexIS = NULL;
2224: const PetscInt *subvertices;
2225: PetscInt *pStart, *pEnd;
2226: PetscInt dim, d, numSubVerticesInitial = 0, v;
2227: PetscErrorCode ierr;
2230: DMGetDimension(dm, &dim);
2231: PetscMalloc2(dim+1, &pStart, dim+1, &pEnd);
2232: for (d = 0; d <= dim; ++d) {
2233: DMPlexGetSimplexOrBoxCells(dm, dim-d, &pStart[d], &pEnd[d]);
2234: }
2235: /* Loop over initial vertices and mark all faces in the collective star() */
2236: if (vertexLabel) {
2237: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2238: if (subvertexIS) {
2239: ISGetSize(subvertexIS, &numSubVerticesInitial);
2240: ISGetIndices(subvertexIS, &subvertices);
2241: }
2242: }
2243: for (v = 0; v < numSubVerticesInitial; ++v) {
2244: const PetscInt vertex = subvertices[v];
2245: PetscInt *star = NULL;
2246: PetscInt starSize, s, numFaces = 0, f;
2248: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2249: for (s = 0; s < starSize*2; s += 2) {
2250: const PetscInt point = star[s];
2251: PetscInt faceLoc;
2253: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2254: if (markedFaces) {
2255: DMLabelGetValue(vertexLabel, point, &faceLoc);
2256: if (faceLoc < 0) continue;
2257: }
2258: star[numFaces++] = point;
2259: }
2260: }
2261: for (f = 0; f < numFaces; ++f) {
2262: const PetscInt face = star[f];
2263: PetscInt *closure = NULL;
2264: PetscInt closureSize, c;
2265: PetscInt faceLoc;
2267: DMLabelGetValue(subpointMap, face, &faceLoc);
2268: if (faceLoc == dim-1) continue;
2269: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2270: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2271: for (c = 0; c < closureSize*2; c += 2) {
2272: const PetscInt point = closure[c];
2273: PetscInt vertexLoc;
2275: if ((point >= pStart[0]) && (point < pEnd[0])) {
2276: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2277: if (vertexLoc != value) break;
2278: }
2279: }
2280: if (c == closureSize*2) {
2281: const PetscInt *support;
2282: PetscInt supportSize, s;
2284: for (c = 0; c < closureSize*2; c += 2) {
2285: const PetscInt point = closure[c];
2287: for (d = 0; d < dim; ++d) {
2288: if ((point >= pStart[d]) && (point < pEnd[d])) {
2289: DMLabelSetValue(subpointMap, point, d);
2290: break;
2291: }
2292: }
2293: }
2294: DMPlexGetSupportSize(dm, face, &supportSize);
2295: DMPlexGetSupport(dm, face, &support);
2296: for (s = 0; s < supportSize; ++s) {
2297: DMLabelSetValue(subpointMap, support[s], dim);
2298: }
2299: }
2300: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2301: }
2302: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2303: }
2304: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2305: ISDestroy(&subvertexIS);
2306: PetscFree2(pStart, pEnd);
2307: return(0);
2308: }
2310: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2311: {
2312: DMLabel label = NULL;
2313: const PetscInt *cone;
2314: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2315: PetscErrorCode ierr;
2318: *numFaces = 0;
2319: *nFV = 0;
2320: if (labelname) {DMGetLabel(dm, labelname, &label);}
2321: *subCells = NULL;
2322: DMGetDimension(dm, &dim);
2323: DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2324: if (cMax < 0) return(0);
2325: if (label) {
2326: for (c = cMax; c < cEnd; ++c) {
2327: PetscInt val;
2329: DMLabelGetValue(label, c, &val);
2330: if (val == value) {
2331: ++(*numFaces);
2332: DMPlexGetConeSize(dm, c, &coneSize);
2333: }
2334: }
2335: } else {
2336: *numFaces = cEnd - cMax;
2337: DMPlexGetConeSize(dm, cMax, &coneSize);
2338: }
2339: PetscMalloc1(*numFaces *2, subCells);
2340: if (!(*numFaces)) return(0);
2341: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2342: for (c = cMax; c < cEnd; ++c) {
2343: const PetscInt *cells;
2344: PetscInt numCells;
2346: if (label) {
2347: PetscInt val;
2349: DMLabelGetValue(label, c, &val);
2350: if (val != value) continue;
2351: }
2352: DMPlexGetCone(dm, c, &cone);
2353: for (p = 0; p < *nFV; ++p) {
2354: DMLabelSetValue(subpointMap, cone[p], 0);
2355: }
2356: /* Negative face */
2357: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2358: /* Not true in parallel
2359: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2360: for (p = 0; p < numCells; ++p) {
2361: DMLabelSetValue(subpointMap, cells[p], 2);
2362: (*subCells)[subc++] = cells[p];
2363: }
2364: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2365: /* Positive face is not included */
2366: }
2367: return(0);
2368: }
2370: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2371: {
2372: PetscInt *pStart, *pEnd;
2373: PetscInt dim, cMax, cEnd, c, d;
2377: DMGetDimension(dm, &dim);
2378: DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2379: if (cMax < 0) return(0);
2380: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2381: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2382: for (c = cMax; c < cEnd; ++c) {
2383: const PetscInt *cone;
2384: PetscInt *closure = NULL;
2385: PetscInt fconeSize, coneSize, closureSize, cl, val;
2387: if (label) {
2388: DMLabelGetValue(label, c, &val);
2389: if (val != value) continue;
2390: }
2391: DMPlexGetConeSize(dm, c, &coneSize);
2392: DMPlexGetCone(dm, c, &cone);
2393: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2394: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2395: /* Negative face */
2396: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2397: for (cl = 0; cl < closureSize*2; cl += 2) {
2398: const PetscInt point = closure[cl];
2400: for (d = 0; d <= dim; ++d) {
2401: if ((point >= pStart[d]) && (point < pEnd[d])) {
2402: DMLabelSetValue(subpointMap, point, d);
2403: break;
2404: }
2405: }
2406: }
2407: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2408: /* Cells -- positive face is not included */
2409: for (cl = 0; cl < 1; ++cl) {
2410: const PetscInt *support;
2411: PetscInt supportSize, s;
2413: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2414: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2415: DMPlexGetSupport(dm, cone[cl], &support);
2416: for (s = 0; s < supportSize; ++s) {
2417: DMLabelSetValue(subpointMap, support[s], dim);
2418: }
2419: }
2420: }
2421: PetscFree2(pStart, pEnd);
2422: return(0);
2423: }
2425: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2426: {
2427: MPI_Comm comm;
2428: PetscBool posOrient = PETSC_FALSE;
2429: const PetscInt debug = 0;
2430: PetscInt cellDim, faceSize, f;
2434: PetscObjectGetComm((PetscObject)dm,&comm);
2435: DMGetDimension(dm, &cellDim);
2436: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2438: if (cellDim == 1 && numCorners == 2) {
2439: /* Triangle */
2440: faceSize = numCorners-1;
2441: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2442: } else if (cellDim == 2 && numCorners == 3) {
2443: /* Triangle */
2444: faceSize = numCorners-1;
2445: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2446: } else if (cellDim == 3 && numCorners == 4) {
2447: /* Tetrahedron */
2448: faceSize = numCorners-1;
2449: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2450: } else if (cellDim == 1 && numCorners == 3) {
2451: /* Quadratic line */
2452: faceSize = 1;
2453: posOrient = PETSC_TRUE;
2454: } else if (cellDim == 2 && numCorners == 4) {
2455: /* Quads */
2456: faceSize = 2;
2457: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2458: posOrient = PETSC_TRUE;
2459: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2460: posOrient = PETSC_TRUE;
2461: } else {
2462: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2463: posOrient = PETSC_FALSE;
2464: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2465: }
2466: } else if (cellDim == 2 && numCorners == 6) {
2467: /* Quadratic triangle (I hate this) */
2468: /* Edges are determined by the first 2 vertices (corners of edges) */
2469: const PetscInt faceSizeTri = 3;
2470: PetscInt sortedIndices[3], i, iFace;
2471: PetscBool found = PETSC_FALSE;
2472: PetscInt faceVerticesTriSorted[9] = {
2473: 0, 3, 4, /* bottom */
2474: 1, 4, 5, /* right */
2475: 2, 3, 5, /* left */
2476: };
2477: PetscInt faceVerticesTri[9] = {
2478: 0, 3, 4, /* bottom */
2479: 1, 4, 5, /* right */
2480: 2, 5, 3, /* left */
2481: };
2483: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2484: PetscSortInt(faceSizeTri, sortedIndices);
2485: for (iFace = 0; iFace < 3; ++iFace) {
2486: const PetscInt ii = iFace*faceSizeTri;
2487: PetscInt fVertex, cVertex;
2489: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2490: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2491: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2492: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2493: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2494: faceVertices[fVertex] = origVertices[cVertex];
2495: break;
2496: }
2497: }
2498: }
2499: found = PETSC_TRUE;
2500: break;
2501: }
2502: }
2503: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2504: if (posOriented) *posOriented = PETSC_TRUE;
2505: return(0);
2506: } else if (cellDim == 2 && numCorners == 9) {
2507: /* Quadratic quad (I hate this) */
2508: /* Edges are determined by the first 2 vertices (corners of edges) */
2509: const PetscInt faceSizeQuad = 3;
2510: PetscInt sortedIndices[3], i, iFace;
2511: PetscBool found = PETSC_FALSE;
2512: PetscInt faceVerticesQuadSorted[12] = {
2513: 0, 1, 4, /* bottom */
2514: 1, 2, 5, /* right */
2515: 2, 3, 6, /* top */
2516: 0, 3, 7, /* left */
2517: };
2518: PetscInt faceVerticesQuad[12] = {
2519: 0, 1, 4, /* bottom */
2520: 1, 2, 5, /* right */
2521: 2, 3, 6, /* top */
2522: 3, 0, 7, /* left */
2523: };
2525: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2526: PetscSortInt(faceSizeQuad, sortedIndices);
2527: for (iFace = 0; iFace < 4; ++iFace) {
2528: const PetscInt ii = iFace*faceSizeQuad;
2529: PetscInt fVertex, cVertex;
2531: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2532: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2533: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2534: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2535: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2536: faceVertices[fVertex] = origVertices[cVertex];
2537: break;
2538: }
2539: }
2540: }
2541: found = PETSC_TRUE;
2542: break;
2543: }
2544: }
2545: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2546: if (posOriented) *posOriented = PETSC_TRUE;
2547: return(0);
2548: } else if (cellDim == 3 && numCorners == 8) {
2549: /* Hexes
2550: A hex is two oriented quads with the normal of the first
2551: pointing up at the second.
2553: 7---6
2554: /| /|
2555: 4---5 |
2556: | 1-|-2
2557: |/ |/
2558: 0---3
2560: Faces are determined by the first 4 vertices (corners of faces) */
2561: const PetscInt faceSizeHex = 4;
2562: PetscInt sortedIndices[4], i, iFace;
2563: PetscBool found = PETSC_FALSE;
2564: PetscInt faceVerticesHexSorted[24] = {
2565: 0, 1, 2, 3, /* bottom */
2566: 4, 5, 6, 7, /* top */
2567: 0, 3, 4, 5, /* front */
2568: 2, 3, 5, 6, /* right */
2569: 1, 2, 6, 7, /* back */
2570: 0, 1, 4, 7, /* left */
2571: };
2572: PetscInt faceVerticesHex[24] = {
2573: 1, 2, 3, 0, /* bottom */
2574: 4, 5, 6, 7, /* top */
2575: 0, 3, 5, 4, /* front */
2576: 3, 2, 6, 5, /* right */
2577: 2, 1, 7, 6, /* back */
2578: 1, 0, 4, 7, /* left */
2579: };
2581: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2582: PetscSortInt(faceSizeHex, sortedIndices);
2583: for (iFace = 0; iFace < 6; ++iFace) {
2584: const PetscInt ii = iFace*faceSizeHex;
2585: PetscInt fVertex, cVertex;
2587: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2588: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2589: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2590: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2591: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2592: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2593: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2594: faceVertices[fVertex] = origVertices[cVertex];
2595: break;
2596: }
2597: }
2598: }
2599: found = PETSC_TRUE;
2600: break;
2601: }
2602: }
2603: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2604: if (posOriented) *posOriented = PETSC_TRUE;
2605: return(0);
2606: } else if (cellDim == 3 && numCorners == 10) {
2607: /* Quadratic tet */
2608: /* Faces are determined by the first 3 vertices (corners of faces) */
2609: const PetscInt faceSizeTet = 6;
2610: PetscInt sortedIndices[6], i, iFace;
2611: PetscBool found = PETSC_FALSE;
2612: PetscInt faceVerticesTetSorted[24] = {
2613: 0, 1, 2, 6, 7, 8, /* bottom */
2614: 0, 3, 4, 6, 7, 9, /* front */
2615: 1, 4, 5, 7, 8, 9, /* right */
2616: 2, 3, 5, 6, 8, 9, /* left */
2617: };
2618: PetscInt faceVerticesTet[24] = {
2619: 0, 1, 2, 6, 7, 8, /* bottom */
2620: 0, 4, 3, 6, 7, 9, /* front */
2621: 1, 5, 4, 7, 8, 9, /* right */
2622: 2, 3, 5, 8, 6, 9, /* left */
2623: };
2625: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2626: PetscSortInt(faceSizeTet, sortedIndices);
2627: for (iFace=0; iFace < 4; ++iFace) {
2628: const PetscInt ii = iFace*faceSizeTet;
2629: PetscInt fVertex, cVertex;
2631: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2632: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2633: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2634: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2635: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2636: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2637: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2638: faceVertices[fVertex] = origVertices[cVertex];
2639: break;
2640: }
2641: }
2642: }
2643: found = PETSC_TRUE;
2644: break;
2645: }
2646: }
2647: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2648: if (posOriented) *posOriented = PETSC_TRUE;
2649: return(0);
2650: } else if (cellDim == 3 && numCorners == 27) {
2651: /* Quadratic hexes (I hate this)
2652: A hex is two oriented quads with the normal of the first
2653: pointing up at the second.
2655: 7---6
2656: /| /|
2657: 4---5 |
2658: | 3-|-2
2659: |/ |/
2660: 0---1
2662: Faces are determined by the first 4 vertices (corners of faces) */
2663: const PetscInt faceSizeQuadHex = 9;
2664: PetscInt sortedIndices[9], i, iFace;
2665: PetscBool found = PETSC_FALSE;
2666: PetscInt faceVerticesQuadHexSorted[54] = {
2667: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2668: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2669: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2670: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2671: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2672: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2673: };
2674: PetscInt faceVerticesQuadHex[54] = {
2675: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2676: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2677: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2678: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2679: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2680: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2681: };
2683: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2684: PetscSortInt(faceSizeQuadHex, sortedIndices);
2685: for (iFace = 0; iFace < 6; ++iFace) {
2686: const PetscInt ii = iFace*faceSizeQuadHex;
2687: PetscInt fVertex, cVertex;
2689: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2690: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2691: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2692: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2693: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2694: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2695: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2696: faceVertices[fVertex] = origVertices[cVertex];
2697: break;
2698: }
2699: }
2700: }
2701: found = PETSC_TRUE;
2702: break;
2703: }
2704: }
2705: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2706: if (posOriented) *posOriented = PETSC_TRUE;
2707: return(0);
2708: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2709: if (!posOrient) {
2710: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2711: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2712: } else {
2713: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2714: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2715: }
2716: if (posOriented) *posOriented = posOrient;
2717: return(0);
2718: }
2720: /*@
2721: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2722: in faceVertices. The orientation is such that the face normal points out of the cell
2724: Not collective
2726: Input Parameters:
2727: + dm - The original mesh
2728: . cell - The cell mesh point
2729: . faceSize - The number of vertices on the face
2730: . face - The face vertices
2731: . numCorners - The number of vertices on the cell
2732: . indices - Local numbering of face vertices in cell cone
2733: - origVertices - Original face vertices
2735: Output Parameter:
2736: + faceVertices - The face vertices properly oriented
2737: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2739: Level: developer
2741: .seealso: DMPlexGetCone()
2742: @*/
2743: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2744: {
2745: const PetscInt *cone = NULL;
2746: PetscInt coneSize, v, f, v2;
2747: PetscInt oppositeVertex = -1;
2748: PetscErrorCode ierr;
2751: DMPlexGetConeSize(dm, cell, &coneSize);
2752: DMPlexGetCone(dm, cell, &cone);
2753: for (v = 0, v2 = 0; v < coneSize; ++v) {
2754: PetscBool found = PETSC_FALSE;
2756: for (f = 0; f < faceSize; ++f) {
2757: if (face[f] == cone[v]) {
2758: found = PETSC_TRUE; break;
2759: }
2760: }
2761: if (found) {
2762: indices[v2] = v;
2763: origVertices[v2] = cone[v];
2764: ++v2;
2765: } else {
2766: oppositeVertex = v;
2767: }
2768: }
2769: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2770: return(0);
2771: }
2773: /*
2774: DMPlexInsertFace_Internal - Puts a face into the mesh
2776: Not collective
2778: Input Parameters:
2779: + dm - The DMPlex
2780: . numFaceVertex - The number of vertices in the face
2781: . faceVertices - The vertices in the face for dm
2782: . subfaceVertices - The vertices in the face for subdm
2783: . numCorners - The number of vertices in the cell
2784: . cell - A cell in dm containing the face
2785: . subcell - A cell in subdm containing the face
2786: . firstFace - First face in the mesh
2787: - newFacePoint - Next face in the mesh
2789: Output Parameters:
2790: . newFacePoint - Contains next face point number on input, updated on output
2792: Level: developer
2793: */
2794: 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)
2795: {
2796: MPI_Comm comm;
2797: DM_Plex *submesh = (DM_Plex*) subdm->data;
2798: const PetscInt *faces;
2799: PetscInt numFaces, coneSize;
2800: PetscErrorCode ierr;
2803: PetscObjectGetComm((PetscObject)dm,&comm);
2804: DMPlexGetConeSize(subdm, subcell, &coneSize);
2805: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2806: #if 0
2807: /* Cannot use this because support() has not been constructed yet */
2808: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2809: #else
2810: {
2811: PetscInt f;
2813: numFaces = 0;
2814: DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2815: for (f = firstFace; f < *newFacePoint; ++f) {
2816: PetscInt dof, off, d;
2818: PetscSectionGetDof(submesh->coneSection, f, &dof);
2819: PetscSectionGetOffset(submesh->coneSection, f, &off);
2820: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2821: for (d = 0; d < dof; ++d) {
2822: const PetscInt p = submesh->cones[off+d];
2823: PetscInt v;
2825: for (v = 0; v < numFaceVertices; ++v) {
2826: if (subfaceVertices[v] == p) break;
2827: }
2828: if (v == numFaceVertices) break;
2829: }
2830: if (d == dof) {
2831: numFaces = 1;
2832: ((PetscInt*) faces)[0] = f;
2833: }
2834: }
2835: }
2836: #endif
2837: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2838: else if (numFaces == 1) {
2839: /* Add the other cell neighbor for this face */
2840: DMPlexSetCone(subdm, subcell, faces);
2841: } else {
2842: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2843: PetscBool posOriented;
2845: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2846: origVertices = &orientedVertices[numFaceVertices];
2847: indices = &orientedVertices[numFaceVertices*2];
2848: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2849: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2850: /* TODO: I know that routine should return a permutation, not the indices */
2851: for (v = 0; v < numFaceVertices; ++v) {
2852: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2853: for (ov = 0; ov < numFaceVertices; ++ov) {
2854: if (orientedVertices[ov] == vertex) {
2855: orientedSubVertices[ov] = subvertex;
2856: break;
2857: }
2858: }
2859: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2860: }
2861: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2862: DMPlexSetCone(subdm, subcell, newFacePoint);
2863: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2864: ++(*newFacePoint);
2865: }
2866: #if 0
2867: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2868: #else
2869: DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2870: #endif
2871: return(0);
2872: }
2874: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2875: {
2876: MPI_Comm comm;
2877: DMLabel subpointMap;
2878: IS subvertexIS, subcellIS;
2879: const PetscInt *subVertices, *subCells;
2880: PetscInt numSubVertices, firstSubVertex, numSubCells;
2881: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2882: PetscInt vStart, vEnd, c, f;
2883: PetscErrorCode ierr;
2886: PetscObjectGetComm((PetscObject)dm,&comm);
2887: /* Create subpointMap which marks the submesh */
2888: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2889: DMPlexSetSubpointMap(subdm, subpointMap);
2890: DMLabelDestroy(&subpointMap);
2891: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2892: /* Setup chart */
2893: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2894: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2895: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2896: DMPlexSetVTKCellHeight(subdm, 1);
2897: /* Set cone sizes */
2898: firstSubVertex = numSubCells;
2899: firstSubFace = numSubCells+numSubVertices;
2900: newFacePoint = firstSubFace;
2901: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2902: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2903: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2904: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2905: for (c = 0; c < numSubCells; ++c) {
2906: DMPlexSetConeSize(subdm, c, 1);
2907: }
2908: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2909: DMPlexSetConeSize(subdm, f, nFV);
2910: }
2911: DMSetUp(subdm);
2912: /* Create face cones */
2913: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2914: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2915: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2916: for (c = 0; c < numSubCells; ++c) {
2917: const PetscInt cell = subCells[c];
2918: const PetscInt subcell = c;
2919: PetscInt *closure = NULL;
2920: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2922: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2923: for (cl = 0; cl < closureSize*2; cl += 2) {
2924: const PetscInt point = closure[cl];
2925: PetscInt subVertex;
2927: if ((point >= vStart) && (point < vEnd)) {
2928: ++numCorners;
2929: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2930: if (subVertex >= 0) {
2931: closure[faceSize] = point;
2932: subface[faceSize] = firstSubVertex+subVertex;
2933: ++faceSize;
2934: }
2935: }
2936: }
2937: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2938: if (faceSize == nFV) {
2939: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2940: }
2941: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2942: }
2943: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2944: DMPlexSymmetrize(subdm);
2945: DMPlexStratify(subdm);
2946: /* Build coordinates */
2947: {
2948: PetscSection coordSection, subCoordSection;
2949: Vec coordinates, subCoordinates;
2950: PetscScalar *coords, *subCoords;
2951: PetscInt numComp, coordSize, v;
2952: const char *name;
2954: DMGetCoordinateSection(dm, &coordSection);
2955: DMGetCoordinatesLocal(dm, &coordinates);
2956: DMGetCoordinateSection(subdm, &subCoordSection);
2957: PetscSectionSetNumFields(subCoordSection, 1);
2958: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2959: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2960: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2961: for (v = 0; v < numSubVertices; ++v) {
2962: const PetscInt vertex = subVertices[v];
2963: const PetscInt subvertex = firstSubVertex+v;
2964: PetscInt dof;
2966: PetscSectionGetDof(coordSection, vertex, &dof);
2967: PetscSectionSetDof(subCoordSection, subvertex, dof);
2968: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2969: }
2970: PetscSectionSetUp(subCoordSection);
2971: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2972: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2973: PetscObjectGetName((PetscObject)coordinates,&name);
2974: PetscObjectSetName((PetscObject)subCoordinates,name);
2975: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2976: VecSetType(subCoordinates,VECSTANDARD);
2977: if (coordSize) {
2978: VecGetArray(coordinates, &coords);
2979: VecGetArray(subCoordinates, &subCoords);
2980: for (v = 0; v < numSubVertices; ++v) {
2981: const PetscInt vertex = subVertices[v];
2982: const PetscInt subvertex = firstSubVertex+v;
2983: PetscInt dof, off, sdof, soff, d;
2985: PetscSectionGetDof(coordSection, vertex, &dof);
2986: PetscSectionGetOffset(coordSection, vertex, &off);
2987: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2988: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2989: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2990: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2991: }
2992: VecRestoreArray(coordinates, &coords);
2993: VecRestoreArray(subCoordinates, &subCoords);
2994: }
2995: DMSetCoordinatesLocal(subdm, subCoordinates);
2996: VecDestroy(&subCoordinates);
2997: }
2998: /* Cleanup */
2999: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3000: ISDestroy(&subvertexIS);
3001: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
3002: ISDestroy(&subcellIS);
3003: return(0);
3004: }
3006: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
3007: {
3008: PetscInt subPoint;
3011: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
3012: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
3013: }
3015: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
3016: {
3017: MPI_Comm comm;
3018: DMLabel subpointMap;
3019: IS *subpointIS;
3020: const PetscInt **subpoints;
3021: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3022: PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v;
3023: PetscMPIInt rank;
3024: PetscErrorCode ierr;
3027: PetscObjectGetComm((PetscObject)dm,&comm);
3028: MPI_Comm_rank(comm, &rank);
3029: /* Create subpointMap which marks the submesh */
3030: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3031: DMPlexSetSubpointMap(subdm, subpointMap);
3032: if (cellHeight) {
3033: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
3034: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
3035: } else {
3036: DMLabel depth;
3037: IS pointIS;
3038: const PetscInt *points;
3039: PetscInt numPoints=0;
3041: DMPlexGetDepthLabel(dm, &depth);
3042: DMLabelGetStratumIS(label, value, &pointIS);
3043: if (pointIS) {
3044: ISGetIndices(pointIS, &points);
3045: ISGetLocalSize(pointIS, &numPoints);
3046: }
3047: for (p = 0; p < numPoints; ++p) {
3048: PetscInt *closure = NULL;
3049: PetscInt closureSize, c, pdim;
3051: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3052: for (c = 0; c < closureSize*2; c += 2) {
3053: DMLabelGetValue(depth, closure[c], &pdim);
3054: DMLabelSetValue(subpointMap, closure[c], pdim);
3055: }
3056: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3057: }
3058: if (pointIS) {ISRestoreIndices(pointIS, &points);}
3059: ISDestroy(&pointIS);
3060: }
3061: /* Setup chart */
3062: DMGetDimension(dm, &dim);
3063: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
3064: for (d = 0; d <= dim; ++d) {
3065: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3066: totSubPoints += numSubPoints[d];
3067: }
3068: DMPlexSetChart(subdm, 0, totSubPoints);
3069: DMPlexSetVTKCellHeight(subdm, cellHeight);
3070: /* Set cone sizes */
3071: firstSubPoint[dim] = 0;
3072: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
3073: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
3074: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
3075: for (d = 0; d <= dim; ++d) {
3076: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3077: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
3078: }
3079: /* We do not want this label automatically computed, instead we compute it here */
3080: DMCreateLabel(subdm, "celltype");
3081: for (d = 0; d <= dim; ++d) {
3082: for (p = 0; p < numSubPoints[d]; ++p) {
3083: const PetscInt point = subpoints[d][p];
3084: const PetscInt subpoint = firstSubPoint[d] + p;
3085: const PetscInt *cone;
3086: PetscInt coneSize, coneSizeNew, c, val;
3087: DMPolytopeType ct;
3089: DMPlexGetConeSize(dm, point, &coneSize);
3090: DMPlexSetConeSize(subdm, subpoint, coneSize);
3091: DMPlexGetCellType(dm, point, &ct);
3092: DMPlexSetCellType(subdm, subpoint, ct);
3093: if (cellHeight && (d == dim)) {
3094: DMPlexGetCone(dm, point, &cone);
3095: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3096: DMLabelGetValue(subpointMap, cone[c], &val);
3097: if (val >= 0) coneSizeNew++;
3098: }
3099: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3100: DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST);
3101: }
3102: }
3103: }
3104: DMLabelDestroy(&subpointMap);
3105: DMSetUp(subdm);
3106: /* Set cones */
3107: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3108: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3109: for (d = 0; d <= dim; ++d) {
3110: for (p = 0; p < numSubPoints[d]; ++p) {
3111: const PetscInt point = subpoints[d][p];
3112: const PetscInt subpoint = firstSubPoint[d] + p;
3113: const PetscInt *cone, *ornt;
3114: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
3116: if (d == dim-1) {
3117: const PetscInt *support, *cone, *ornt;
3118: PetscInt supportSize, coneSize, s, subc;
3120: DMPlexGetSupport(dm, point, &support);
3121: DMPlexGetSupportSize(dm, point, &supportSize);
3122: for (s = 0; s < supportSize; ++s) {
3123: PetscBool isHybrid;
3125: DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);
3126: if (!isHybrid) continue;
3127: PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3128: if (subc >= 0) {
3129: const PetscInt ccell = subpoints[d+1][subc];
3131: DMPlexGetCone(dm, ccell, &cone);
3132: DMPlexGetConeSize(dm, ccell, &coneSize);
3133: DMPlexGetConeOrientation(dm, ccell, &ornt);
3134: for (c = 0; c < coneSize; ++c) {
3135: if (cone[c] == point) {
3136: fornt = ornt[c];
3137: break;
3138: }
3139: }
3140: break;
3141: }
3142: }
3143: }
3144: DMPlexGetConeSize(dm, point, &coneSize);
3145: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3146: DMPlexGetCone(dm, point, &cone);
3147: DMPlexGetConeOrientation(dm, point, &ornt);
3148: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3149: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3150: if (subc >= 0) {
3151: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3152: orntNew[coneSizeNew] = ornt[c];
3153: ++coneSizeNew;
3154: }
3155: }
3156: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3157: if (fornt < 0) {
3158: /* This should be replaced by a call to DMPlexReverseCell() */
3159: #if 0
3160: DMPlexReverseCell(subdm, subpoint);
3161: #else
3162: for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
3163: PetscInt faceSize, tmp;
3165: tmp = coneNew[c];
3166: coneNew[c] = coneNew[coneSizeNew-1-c];
3167: coneNew[coneSizeNew-1-c] = tmp;
3168: DMPlexGetConeSize(dm, cone[c], &faceSize);
3169: tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
3170: orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
3171: orntNew[coneSizeNew-1-c] = tmp;
3172: }
3173: }
3174: DMPlexSetCone(subdm, subpoint, coneNew);
3175: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3176: #endif
3177: }
3178: }
3179: PetscFree2(coneNew,orntNew);
3180: DMPlexSymmetrize(subdm);
3181: DMPlexStratify(subdm);
3182: /* Build coordinates */
3183: {
3184: PetscSection coordSection, subCoordSection;
3185: Vec coordinates, subCoordinates;
3186: PetscScalar *coords, *subCoords;
3187: PetscInt cdim, numComp, coordSize;
3188: const char *name;
3190: DMGetCoordinateDim(dm, &cdim);
3191: DMGetCoordinateSection(dm, &coordSection);
3192: DMGetCoordinatesLocal(dm, &coordinates);
3193: DMGetCoordinateSection(subdm, &subCoordSection);
3194: PetscSectionSetNumFields(subCoordSection, 1);
3195: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3196: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3197: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3198: for (v = 0; v < numSubPoints[0]; ++v) {
3199: const PetscInt vertex = subpoints[0][v];
3200: const PetscInt subvertex = firstSubPoint[0]+v;
3201: PetscInt dof;
3203: PetscSectionGetDof(coordSection, vertex, &dof);
3204: PetscSectionSetDof(subCoordSection, subvertex, dof);
3205: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3206: }
3207: PetscSectionSetUp(subCoordSection);
3208: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3209: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3210: PetscObjectGetName((PetscObject)coordinates,&name);
3211: PetscObjectSetName((PetscObject)subCoordinates,name);
3212: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3213: VecSetBlockSize(subCoordinates, cdim);
3214: VecSetType(subCoordinates,VECSTANDARD);
3215: VecGetArray(coordinates, &coords);
3216: VecGetArray(subCoordinates, &subCoords);
3217: for (v = 0; v < numSubPoints[0]; ++v) {
3218: const PetscInt vertex = subpoints[0][v];
3219: const PetscInt subvertex = firstSubPoint[0]+v;
3220: PetscInt dof, off, sdof, soff, d;
3222: PetscSectionGetDof(coordSection, vertex, &dof);
3223: PetscSectionGetOffset(coordSection, vertex, &off);
3224: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3225: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3226: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3227: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3228: }
3229: VecRestoreArray(coordinates, &coords);
3230: VecRestoreArray(subCoordinates, &subCoords);
3231: DMSetCoordinatesLocal(subdm, subCoordinates);
3232: VecDestroy(&subCoordinates);
3233: }
3234: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3235: {
3236: PetscSF sfPoint, sfPointSub;
3237: IS subpIS;
3238: const PetscSFNode *remotePoints;
3239: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3240: const PetscInt *localPoints, *subpoints;
3241: PetscInt *slocalPoints;
3242: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3243: PetscMPIInt rank;
3245: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3246: DMGetPointSF(dm, &sfPoint);
3247: DMGetPointSF(subdm, &sfPointSub);
3248: DMPlexGetChart(dm, &pStart, &pEnd);
3249: DMPlexGetChart(subdm, NULL, &numSubroots);
3250: DMPlexCreateSubpointIS(subdm, &subpIS);
3251: if (subpIS) {
3252: ISGetIndices(subpIS, &subpoints);
3253: ISGetLocalSize(subpIS, &numSubpoints);
3254: }
3255: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3256: if (numRoots >= 0) {
3257: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3258: for (p = 0; p < pEnd-pStart; ++p) {
3259: newLocalPoints[p].rank = -2;
3260: newLocalPoints[p].index = -2;
3261: }
3262: /* Set subleaves */
3263: for (l = 0; l < numLeaves; ++l) {
3264: const PetscInt point = localPoints[l];
3265: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3267: if (subpoint < 0) continue;
3268: newLocalPoints[point-pStart].rank = rank;
3269: newLocalPoints[point-pStart].index = subpoint;
3270: ++numSubleaves;
3271: }
3272: /* Must put in owned subpoints */
3273: for (p = pStart; p < pEnd; ++p) {
3274: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3276: if (subpoint < 0) {
3277: newOwners[p-pStart].rank = -3;
3278: newOwners[p-pStart].index = -3;
3279: } else {
3280: newOwners[p-pStart].rank = rank;
3281: newOwners[p-pStart].index = subpoint;
3282: }
3283: }
3284: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3285: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3286: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3287: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3288: PetscMalloc1(numSubleaves, &slocalPoints);
3289: PetscMalloc1(numSubleaves, &sremotePoints);
3290: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3291: const PetscInt point = localPoints[l];
3292: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3294: if (subpoint < 0) continue;
3295: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3296: slocalPoints[sl] = subpoint;
3297: sremotePoints[sl].rank = newLocalPoints[point].rank;
3298: sremotePoints[sl].index = newLocalPoints[point].index;
3299: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3300: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3301: ++sl;
3302: }
3303: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3304: PetscFree2(newLocalPoints,newOwners);
3305: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3306: }
3307: if (subpIS) {
3308: ISRestoreIndices(subpIS, &subpoints);
3309: ISDestroy(&subpIS);
3310: }
3311: }
3312: /* Cleanup */
3313: for (d = 0; d <= dim; ++d) {
3314: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3315: ISDestroy(&subpointIS[d]);
3316: }
3317: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3318: return(0);
3319: }
3321: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3322: {
3326: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3327: return(0);
3328: }
3330: /*@
3331: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3333: Input Parameters:
3334: + dm - The original mesh
3335: . vertexLabel - The DMLabel marking points contained in the surface
3336: . value - The label value to use
3337: - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked
3339: Output Parameter:
3340: . subdm - The surface mesh
3342: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3344: Level: developer
3346: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3347: @*/
3348: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3349: {
3350: DMPlexInterpolatedFlag interpolated;
3351: PetscInt dim, cdim;
3357: DMGetDimension(dm, &dim);
3358: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3359: DMSetType(*subdm, DMPLEX);
3360: DMSetDimension(*subdm, dim-1);
3361: DMGetCoordinateDim(dm, &cdim);
3362: DMSetCoordinateDim(*subdm, cdim);
3363: DMPlexIsInterpolated(dm, &interpolated);
3364: if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
3365: if (interpolated) {
3366: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3367: } else {
3368: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3369: }
3370: return(0);
3371: }
3373: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3374: {
3375: MPI_Comm comm;
3376: DMLabel subpointMap;
3377: IS subvertexIS;
3378: const PetscInt *subVertices;
3379: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3380: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3381: PetscInt c, f;
3382: PetscErrorCode ierr;
3385: PetscObjectGetComm((PetscObject)dm, &comm);
3386: /* Create subpointMap which marks the submesh */
3387: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3388: DMPlexSetSubpointMap(subdm, subpointMap);
3389: DMLabelDestroy(&subpointMap);
3390: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3391: /* Setup chart */
3392: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3393: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3394: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3395: DMPlexSetVTKCellHeight(subdm, 1);
3396: /* Set cone sizes */
3397: firstSubVertex = numSubCells;
3398: firstSubFace = numSubCells+numSubVertices;
3399: newFacePoint = firstSubFace;
3400: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3401: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3402: for (c = 0; c < numSubCells; ++c) {
3403: DMPlexSetConeSize(subdm, c, 1);
3404: }
3405: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3406: DMPlexSetConeSize(subdm, f, nFV);
3407: }
3408: DMSetUp(subdm);
3409: /* Create face cones */
3410: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3411: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3412: for (c = 0; c < numSubCells; ++c) {
3413: const PetscInt cell = subCells[c];
3414: const PetscInt subcell = c;
3415: const PetscInt *cone, *cells;
3416: PetscBool isHybrid;
3417: PetscInt numCells, subVertex, p, v;
3419: DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid);
3420: if (!isHybrid) continue;
3421: DMPlexGetCone(dm, cell, &cone);
3422: for (v = 0; v < nFV; ++v) {
3423: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3424: subface[v] = firstSubVertex+subVertex;
3425: }
3426: DMPlexSetCone(subdm, newFacePoint, subface);
3427: DMPlexSetCone(subdm, subcell, &newFacePoint);
3428: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3429: /* Not true in parallel
3430: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3431: for (p = 0; p < numCells; ++p) {
3432: PetscInt negsubcell;
3433: PetscBool isHybrid;
3435: DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid);
3436: if (isHybrid) continue;
3437: /* I know this is a crap search */
3438: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3439: if (subCells[negsubcell] == cells[p]) break;
3440: }
3441: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3442: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3443: }
3444: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3445: ++newFacePoint;
3446: }
3447: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3448: DMPlexSymmetrize(subdm);
3449: DMPlexStratify(subdm);
3450: /* Build coordinates */
3451: {
3452: PetscSection coordSection, subCoordSection;
3453: Vec coordinates, subCoordinates;
3454: PetscScalar *coords, *subCoords;
3455: PetscInt cdim, numComp, coordSize, v;
3456: const char *name;
3458: DMGetCoordinateDim(dm, &cdim);
3459: DMGetCoordinateSection(dm, &coordSection);
3460: DMGetCoordinatesLocal(dm, &coordinates);
3461: DMGetCoordinateSection(subdm, &subCoordSection);
3462: PetscSectionSetNumFields(subCoordSection, 1);
3463: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3464: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3465: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3466: for (v = 0; v < numSubVertices; ++v) {
3467: const PetscInt vertex = subVertices[v];
3468: const PetscInt subvertex = firstSubVertex+v;
3469: PetscInt dof;
3471: PetscSectionGetDof(coordSection, vertex, &dof);
3472: PetscSectionSetDof(subCoordSection, subvertex, dof);
3473: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3474: }
3475: PetscSectionSetUp(subCoordSection);
3476: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3477: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3478: PetscObjectGetName((PetscObject)coordinates,&name);
3479: PetscObjectSetName((PetscObject)subCoordinates,name);
3480: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3481: VecSetBlockSize(subCoordinates, cdim);
3482: VecSetType(subCoordinates,VECSTANDARD);
3483: VecGetArray(coordinates, &coords);
3484: VecGetArray(subCoordinates, &subCoords);
3485: for (v = 0; v < numSubVertices; ++v) {
3486: const PetscInt vertex = subVertices[v];
3487: const PetscInt subvertex = firstSubVertex+v;
3488: PetscInt dof, off, sdof, soff, d;
3490: PetscSectionGetDof(coordSection, vertex, &dof);
3491: PetscSectionGetOffset(coordSection, vertex, &off);
3492: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3493: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3494: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3495: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3496: }
3497: VecRestoreArray(coordinates, &coords);
3498: VecRestoreArray(subCoordinates, &subCoords);
3499: DMSetCoordinatesLocal(subdm, subCoordinates);
3500: VecDestroy(&subCoordinates);
3501: }
3502: /* Build SF */
3503: CHKMEMQ;
3504: {
3505: PetscSF sfPoint, sfPointSub;
3506: const PetscSFNode *remotePoints;
3507: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3508: const PetscInt *localPoints;
3509: PetscInt *slocalPoints;
3510: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3511: PetscMPIInt rank;
3513: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3514: DMGetPointSF(dm, &sfPoint);
3515: DMGetPointSF(subdm, &sfPointSub);
3516: DMPlexGetChart(dm, &pStart, &pEnd);
3517: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3518: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3519: if (numRoots >= 0) {
3520: /* Only vertices should be shared */
3521: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3522: for (p = 0; p < pEnd-pStart; ++p) {
3523: newLocalPoints[p].rank = -2;
3524: newLocalPoints[p].index = -2;
3525: }
3526: /* Set subleaves */
3527: for (l = 0; l < numLeaves; ++l) {
3528: const PetscInt point = localPoints[l];
3529: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3531: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3532: if (subPoint < 0) continue;
3533: newLocalPoints[point-pStart].rank = rank;
3534: newLocalPoints[point-pStart].index = subPoint;
3535: ++numSubLeaves;
3536: }
3537: /* Must put in owned subpoints */
3538: for (p = pStart; p < pEnd; ++p) {
3539: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3541: if (subPoint < 0) {
3542: newOwners[p-pStart].rank = -3;
3543: newOwners[p-pStart].index = -3;
3544: } else {
3545: newOwners[p-pStart].rank = rank;
3546: newOwners[p-pStart].index = subPoint;
3547: }
3548: }
3549: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3550: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3551: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3552: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3553: PetscMalloc1(numSubLeaves, &slocalPoints);
3554: PetscMalloc1(numSubLeaves, &sremotePoints);
3555: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3556: const PetscInt point = localPoints[l];
3557: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3559: if (subPoint < 0) continue;
3560: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3561: slocalPoints[sl] = subPoint;
3562: sremotePoints[sl].rank = newLocalPoints[point].rank;
3563: sremotePoints[sl].index = newLocalPoints[point].index;
3564: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3565: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3566: ++sl;
3567: }
3568: PetscFree2(newLocalPoints,newOwners);
3569: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3570: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3571: }
3572: }
3573: CHKMEMQ;
3574: /* Cleanup */
3575: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3576: ISDestroy(&subvertexIS);
3577: PetscFree(subCells);
3578: return(0);
3579: }
3581: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3582: {
3583: DMLabel label = NULL;
3587: if (labelname) {DMGetLabel(dm, labelname, &label);}
3588: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3589: return(0);
3590: }
3592: /*@C
3593: 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.
3595: Input Parameters:
3596: + dm - The original mesh
3597: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3598: . label - A label name, or NULL
3599: - value - A label value
3601: Output Parameter:
3602: . subdm - The surface mesh
3604: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3606: Level: developer
3608: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3609: @*/
3610: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3611: {
3612: PetscInt dim, cdim, depth;
3618: DMGetDimension(dm, &dim);
3619: DMPlexGetDepth(dm, &depth);
3620: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3621: DMSetType(*subdm, DMPLEX);
3622: DMSetDimension(*subdm, dim-1);
3623: DMGetCoordinateDim(dm, &cdim);
3624: DMSetCoordinateDim(*subdm, cdim);
3625: if (depth == dim) {
3626: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3627: } else {
3628: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3629: }
3630: return(0);
3631: }
3633: /*@
3634: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3636: Input Parameters:
3637: + dm - The original mesh
3638: . cellLabel - The DMLabel marking cells contained in the new mesh
3639: - value - The label value to use
3641: Output Parameter:
3642: . subdm - The new mesh
3644: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3646: Level: developer
3648: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3649: @*/
3650: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3651: {
3652: PetscInt dim;
3658: DMGetDimension(dm, &dim);
3659: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3660: DMSetType(*subdm, DMPLEX);
3661: DMSetDimension(*subdm, dim);
3662: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3663: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3664: return(0);
3665: }
3667: /*@
3668: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3670: Input Parameter:
3671: . dm - The submesh DM
3673: Output Parameter:
3674: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3676: Level: developer
3678: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3679: @*/
3680: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3681: {
3685: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3686: return(0);
3687: }
3689: /*@
3690: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3692: Input Parameters:
3693: + dm - The submesh DM
3694: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3696: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3698: Level: developer
3700: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3701: @*/
3702: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3703: {
3704: DM_Plex *mesh = (DM_Plex *) dm->data;
3705: DMLabel tmp;
3710: tmp = mesh->subpointMap;
3711: mesh->subpointMap = subpointMap;
3712: PetscObjectReference((PetscObject) mesh->subpointMap);
3713: DMLabelDestroy(&tmp);
3714: return(0);
3715: }
3717: /*@
3718: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
3720: Input Parameter:
3721: . dm - The submesh DM
3723: Output Parameter:
3724: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3726: Note: This IS is guaranteed to be sorted by the construction of the submesh
3728: Level: developer
3730: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3731: @*/
3732: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3733: {
3734: MPI_Comm comm;
3735: DMLabel subpointMap;
3736: IS is;
3737: const PetscInt *opoints;
3738: PetscInt *points, *depths;
3739: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3740: PetscErrorCode ierr;
3745: PetscObjectGetComm((PetscObject)dm,&comm);
3746: *subpointIS = NULL;
3747: DMPlexGetSubpointMap(dm, &subpointMap);
3748: DMPlexGetDepth(dm, &depth);
3749: if (subpointMap && depth >= 0) {
3750: DMPlexGetChart(dm, &pStart, &pEnd);
3751: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3752: DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3753: depths[0] = depth;
3754: depths[1] = 0;
3755: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3756: PetscMalloc1(pEnd, &points);
3757: for(d = 0, off = 0; d <= depth; ++d) {
3758: const PetscInt dep = depths[d];
3760: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3761: DMLabelGetStratumSize(subpointMap, dep, &n);
3762: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3763: if (n != depEnd-depStart) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d at depth %d should be %d", n, dep, depEnd-depStart);
3764: } else {
3765: if (!n) {
3766: if (d == 0) {
3767: /* Missing cells */
3768: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3769: } else {
3770: /* Missing faces */
3771: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3772: }
3773: }
3774: }
3775: if (n) {
3776: DMLabelGetStratumIS(subpointMap, dep, &is);
3777: ISGetIndices(is, &opoints);
3778: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3779: ISRestoreIndices(is, &opoints);
3780: ISDestroy(&is);
3781: }
3782: }
3783: DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3784: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3785: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3786: }
3787: return(0);
3788: }
3790: /*@
3791: DMGetEnclosureRelation - Get the relationship between dmA and dmB
3793: Input Parameters:
3794: + dmA - The first DM
3795: - dmB - The second DM
3797: Output Parameter:
3798: . rel - The relation of dmA to dmB
3800: Level: intermediate
3802: .seealso: DMPlexGetEnclosurePoint()
3803: @*/
3804: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
3805: {
3806: DM plexA, plexB, sdm;
3807: DMLabel spmap;
3808: PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB;
3813: *rel = DM_ENC_NONE;
3814: if (!dmA || !dmB) return(0);
3817: if (dmA == dmB) {*rel = DM_ENC_EQUALITY; return(0);}
3818: DMConvert(dmA, DMPLEX, &plexA);
3819: DMConvert(dmB, DMPLEX, &plexB);
3820: DMPlexGetChart(plexA, &pStartA, &pEndA);
3821: DMPlexGetChart(plexB, &pStartB, &pEndB);
3822: /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
3823: - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
3824: if ((pStartA == pStartB) && (pEndA == pEndB)) {
3825: *rel = DM_ENC_EQUALITY;
3826: goto end;
3827: }
3828: NpA = pEndA - pStartA;
3829: NpB = pEndB - pStartB;
3830: if (NpA == NpB) goto end;
3831: sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
3832: DMPlexGetSubpointMap(sdm, &spmap);
3833: if (!spmap) goto end;
3834: /* TODO Check the space mapped to by subpointMap is same size as dm */
3835: if (NpA > NpB) {
3836: *rel = DM_ENC_SUPERMESH;
3837: } else {
3838: *rel = DM_ENC_SUBMESH;
3839: }
3840: end:
3841: DMDestroy(&plexA);
3842: DMDestroy(&plexB);
3843: return(0);
3844: }
3846: /*@
3847: DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB
3849: Input Parameters:
3850: + dmA - The first DM
3851: . dmB - The second DM
3852: . etype - The type of enclosure relation that dmA has to dmB
3853: - pB - A point of dmB
3855: Output Parameter:
3856: . pA - The corresponding point of dmA
3858: Level: intermediate
3860: .seealso: DMGetEnclosureRelation()
3861: @*/
3862: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
3863: {
3864: DM sdm;
3865: DMLabel spmap;
3866: IS subpointIS;
3867: const PetscInt *subpoints;
3868: PetscInt numSubpoints;
3869: PetscErrorCode ierr;
3872: /* TODO Cache the IS, making it look like an index */
3873: switch (etype) {
3874: case DM_ENC_SUPERMESH:
3875: sdm = dmB;
3876: DMPlexGetSubpointMap(sdm, &spmap);
3877: DMPlexCreateSubpointIS(sdm, &subpointIS);
3878: ISGetIndices(subpointIS, &subpoints);
3879: *pA = subpoints[pB];
3880: ISRestoreIndices(subpointIS, &subpoints);
3881: ISDestroy(&subpointIS);
3882: break;
3883: case DM_ENC_SUBMESH:
3884: sdm = dmA;
3885: DMPlexGetSubpointMap(sdm, &spmap);
3886: DMPlexCreateSubpointIS(sdm, &subpointIS);
3887: ISGetLocalSize(subpointIS, &numSubpoints);
3888: ISGetIndices(subpointIS, &subpoints);
3889: PetscFindInt(pB, numSubpoints, subpoints, pA);
3890: if (*pA < 0) {
3891: DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
3892: DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
3893: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB);
3894: }
3895: ISRestoreIndices(subpointIS, &subpoints);
3896: ISDestroy(&subpointIS);
3897: break;
3898: case DM_ENC_EQUALITY:
3899: case DM_ENC_NONE:
3900: *pA = pB;break;
3901: case DM_ENC_UNKNOWN:
3902: {
3903: DMEnclosureType enc;
3905: DMGetEnclosureRelation(dmA, dmB, &enc);
3906: DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
3907: }
3908: break;
3909: default: SETERRQ1(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype);
3910: }
3911: return(0);
3912: }