Actual source code: plexsubmesh.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petscsf.h>
5: static PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt cellHeight, DMLabel label)
6: {
7: PetscInt fStart, fEnd, f;
12: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
13: for (f = fStart; f < fEnd; ++f) {
14: PetscInt supportSize;
16: DMPlexGetSupportSize(dm, f, &supportSize);
17: if (supportSize == 1) {DMLabelSetValue(label, f, 1);}
18: }
19: return(0);
20: }
22: /*@
23: DMPlexMarkBoundaryFaces - Mark all faces on the boundary
25: Not Collective
27: Input Parameter:
28: . dm - The original DM
30: Output Parameter:
31: . label - The DMLabel marking boundary faces with value 1
33: Level: developer
35: .seealso: DMLabelCreate(), DMCreateLabel()
36: @*/
37: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
38: {
42: DMPlexMarkBoundaryFaces_Internal(dm, 0, label);
43: return(0);
44: }
46: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
47: {
48: IS valueIS;
49: PetscSF sfPoint;
50: const PetscInt *values;
51: PetscInt numValues, v, cStart, cEnd, nroots;
52: PetscErrorCode ierr;
55: DMLabelGetNumValues(label, &numValues);
56: DMLabelGetValueIS(label, &valueIS);
57: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
58: ISGetIndices(valueIS, &values);
59: for (v = 0; v < numValues; ++v) {
60: IS pointIS;
61: const PetscInt *points;
62: PetscInt numPoints, p;
64: DMLabelGetStratumSize(label, values[v], &numPoints);
65: DMLabelGetStratumIS(label, values[v], &pointIS);
66: ISGetIndices(pointIS, &points);
67: for (p = 0; p < numPoints; ++p) {
68: PetscInt q = points[p];
69: PetscInt *closure = NULL;
70: PetscInt closureSize, c;
72: if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
73: continue;
74: }
75: DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
76: for (c = 0; c < closureSize*2; c += 2) {
77: DMLabelSetValue(label, closure[c], values[v]);
78: }
79: DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
80: }
81: ISRestoreIndices(pointIS, &points);
82: ISDestroy(&pointIS);
83: }
84: ISRestoreIndices(valueIS, &values);
85: ISDestroy(&valueIS);
86: DMGetPointSF(dm, &sfPoint);
87: PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
88: if (nroots >= 0) {
89: DMLabel lblRoots, lblLeaves;
90: IS valueIS, pointIS;
91: const PetscInt *values;
92: PetscInt numValues, v;
93: PetscErrorCode ierr;
95: DMGetPointSF(dm, &sfPoint);
96: /* Pull point contributions from remote leaves into local roots */
97: DMLabelGather(label, sfPoint, &lblLeaves);
98: DMLabelGetValueIS(lblLeaves, &valueIS);
99: ISGetLocalSize(valueIS, &numValues);
100: ISGetIndices(valueIS, &values);
101: for (v = 0; v < numValues; ++v) {
102: const PetscInt value = values[v];
103:
104: DMLabelGetStratumIS(lblLeaves, value, &pointIS);
105: DMLabelInsertIS(label, pointIS, value);
106: ISDestroy(&pointIS);
107: }
108: ISRestoreIndices(valueIS, &values);
109: ISDestroy(&valueIS);
110: DMLabelDestroy(&lblLeaves);
111: /* Push point contributions from roots into remote leaves */
112: DMLabelDistribute(label, sfPoint, &lblRoots);
113: DMLabelGetValueIS(lblRoots, &valueIS);
114: ISGetLocalSize(valueIS, &numValues);
115: ISGetIndices(valueIS, &values);
116: for (v = 0; v < numValues; ++v) {
117: const PetscInt value = values[v];
118:
119: DMLabelGetStratumIS(lblRoots, value, &pointIS);
120: DMLabelInsertIS(label, pointIS, value);
121: ISDestroy(&pointIS);
122: }
123: ISRestoreIndices(valueIS, &values);
124: ISDestroy(&valueIS);
125: DMLabelDestroy(&lblRoots);
126: }
127: return(0);
128: }
130: /*@
131: DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface
133: Input Parameters:
134: + dm - The DM
135: - label - A DMLabel marking the surface points
137: Output Parameter:
138: . label - A DMLabel marking all surface points in the transitive closure
140: Level: developer
142: .seealso: DMPlexLabelCohesiveComplete()
143: @*/
144: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
145: {
149: DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
150: return(0);
151: }
153: /*@
154: DMPlexLabelAddCells - Starting with a label marking faces on a surface, we add a cell for each face
156: Input Parameters:
157: + dm - The DM
158: - label - A DMLabel marking the surface points
160: Output Parameter:
161: . label - A DMLabel incorporating cells
163: Level: developer
165: Note: The cells allow FEM boundary conditions to be applied using the cell geometry
167: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
168: @*/
169: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
170: {
171: IS valueIS;
172: const PetscInt *values;
173: PetscInt numValues, v, cStart, cEnd, cEndInterior;
174: PetscErrorCode ierr;
177: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
178: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
179: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
180: DMLabelGetNumValues(label, &numValues);
181: DMLabelGetValueIS(label, &valueIS);
182: ISGetIndices(valueIS, &values);
183: for (v = 0; v < numValues; ++v) {
184: IS pointIS;
185: const PetscInt *points;
186: PetscInt numPoints, p;
188: DMLabelGetStratumSize(label, values[v], &numPoints);
189: DMLabelGetStratumIS(label, values[v], &pointIS);
190: ISGetIndices(pointIS, &points);
191: for (p = 0; p < numPoints; ++p) {
192: PetscInt *closure = NULL;
193: PetscInt closureSize, point, cl;
195: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
196: for (cl = closureSize-1; cl > 0; --cl) {
197: point = closure[cl*2];
198: if ((point >= cStart) && (point < cEnd)) {DMLabelSetValue(label, point, values[v]); break;}
199: }
200: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
201: }
202: ISRestoreIndices(pointIS, &points);
203: ISDestroy(&pointIS);
204: }
205: ISRestoreIndices(valueIS, &values);
206: ISDestroy(&valueIS);
207: return(0);
208: }
210: /*@
211: DMPlexLabelClearCells - Remove cells from a label
213: Input Parameters:
214: + dm - The DM
215: - label - A DMLabel marking surface points and their adjacent cells
217: Output Parameter:
218: . label - A DMLabel without cells
220: Level: developer
222: Note: This undoes DMPlexLabelAddCells()
224: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
225: @*/
226: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
227: {
228: IS valueIS;
229: const PetscInt *values;
230: PetscInt numValues, v, cStart, cEnd, cEndInterior;
231: PetscErrorCode ierr;
234: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
235: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
236: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
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 point = points[p];
251: if (point >= cStart && point < cEnd) {
252: DMLabelClearValue(label,point,values[v]);
253: }
254: }
255: ISRestoreIndices(pointIS, &points);
256: ISDestroy(&pointIS);
257: }
258: ISRestoreIndices(valueIS, &values);
259: ISDestroy(&valueIS);
260: return(0);
261: }
263: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
264: * index (skipping first, which is (0,0)) */
265: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
266: {
267: PetscInt d, off = 0;
270: /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
271: for (d = 0; d < depth; d++) {
272: PetscInt firstd = d;
273: PetscInt firstStart = depthShift[2*d];
274: PetscInt e;
276: for (e = d+1; e <= depth; e++) {
277: if (depthShift[2*e] < firstStart) {
278: firstd = e;
279: firstStart = depthShift[2*d];
280: }
281: }
282: if (firstd != d) {
283: PetscInt swap[2];
285: e = firstd;
286: swap[0] = depthShift[2*d];
287: swap[1] = depthShift[2*d+1];
288: depthShift[2*d] = depthShift[2*e];
289: depthShift[2*d+1] = depthShift[2*e+1];
290: depthShift[2*e] = swap[0];
291: depthShift[2*e+1] = swap[1];
292: }
293: }
294: /* convert (oldstart, added) to (oldstart, newstart) */
295: for (d = 0; d <= depth; d++) {
296: off += depthShift[2*d+1];
297: depthShift[2*d+1] = depthShift[2*d] + off;
298: }
299: return(0);
300: }
302: /* depthShift is a list of (old, new) pairs */
303: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
304: {
305: PetscInt d;
306: PetscInt newOff = 0;
308: for (d = 0; d <= depth; d++) {
309: if (p < depthShift[2*d]) return p + newOff;
310: else newOff = depthShift[2*d+1] - depthShift[2*d];
311: }
312: return p + newOff;
313: }
315: /* depthShift is a list of (old, new) pairs */
316: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
317: {
318: PetscInt d;
319: PetscInt newOff = 0;
321: for (d = 0; d <= depth; d++) {
322: if (p < depthShift[2*d+1]) return p + newOff;
323: else newOff = depthShift[2*d] - depthShift[2*d+1];
324: }
325: return p + newOff;
326: }
328: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
329: {
330: PetscInt depth = 0, d, pStart, pEnd, p;
334: DMPlexGetDepth(dm, &depth);
335: if (depth < 0) return(0);
336: /* Step 1: Expand chart */
337: DMPlexGetChart(dm, &pStart, &pEnd);
338: pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
339: DMPlexSetChart(dmNew, pStart, pEnd);
340: /* Step 2: Set cone and support sizes */
341: for (d = 0; d <= depth; ++d) {
342: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
343: for (p = pStart; p < pEnd; ++p) {
344: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
345: PetscInt size;
347: DMPlexGetConeSize(dm, p, &size);
348: DMPlexSetConeSize(dmNew, newp, size);
349: DMPlexGetSupportSize(dm, p, &size);
350: DMPlexSetSupportSize(dmNew, newp, size);
351: }
352: }
353: return(0);
354: }
356: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
357: {
358: PetscInt *newpoints;
359: PetscInt depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;
363: DMPlexGetDepth(dm, &depth);
364: if (depth < 0) return(0);
365: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
366: DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
367: PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
368: /* Step 5: Set cones and supports */
369: DMPlexGetChart(dm, &pStart, &pEnd);
370: for (p = pStart; p < pEnd; ++p) {
371: const PetscInt *points = NULL, *orientations = NULL;
372: PetscInt size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
374: DMPlexGetConeSize(dm, p, &size);
375: DMPlexGetCone(dm, p, &points);
376: DMPlexGetConeOrientation(dm, p, &orientations);
377: for (i = 0; i < size; ++i) {
378: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
379: }
380: DMPlexSetCone(dmNew, newp, newpoints);
381: DMPlexSetConeOrientation(dmNew, newp, orientations);
382: DMPlexGetSupportSize(dm, p, &size);
383: DMPlexGetSupportSize(dmNew, newp, &sizeNew);
384: DMPlexGetSupport(dm, p, &points);
385: for (i = 0; i < size; ++i) {
386: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
387: }
388: for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
389: DMPlexSetSupport(dmNew, newp, newpoints);
390: }
391: PetscFree(newpoints);
392: return(0);
393: }
395: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
396: {
397: PetscSection coordSection, newCoordSection;
398: Vec coordinates, newCoordinates;
399: PetscScalar *coords, *newCoords;
400: PetscInt coordSize, sStart, sEnd;
401: PetscInt dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
402: PetscBool hasCells;
406: DMGetCoordinateDim(dm, &dim);
407: DMSetCoordinateDim(dmNew, dim);
408: DMPlexGetDepth(dm, &depth);
409: /* Step 8: Convert coordinates */
410: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
411: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
412: DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
413: DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
414: DMGetCoordinateSection(dm, &coordSection);
415: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
416: PetscSectionSetNumFields(newCoordSection, 1);
417: PetscSectionSetFieldComponents(newCoordSection, 0, dim);
418: PetscSectionGetChart(coordSection, &sStart, &sEnd);
419: hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
420: PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
421: if (hasCells) {
422: for (c = cStart; c < cEnd; ++c) {
423: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;
425: PetscSectionGetDof(coordSection, c, &dof);
426: PetscSectionSetDof(newCoordSection, cNew, dof);
427: PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
428: }
429: }
430: for (v = vStartNew; v < vEndNew; ++v) {
431: PetscSectionSetDof(newCoordSection, v, dim);
432: PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
433: }
434: PetscSectionSetUp(newCoordSection);
435: DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
436: PetscSectionGetStorageSize(newCoordSection, &coordSize);
437: VecCreate(PETSC_COMM_SELF, &newCoordinates);
438: PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
439: VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
440: VecSetBlockSize(newCoordinates, dim);
441: VecSetType(newCoordinates,VECSTANDARD);
442: DMSetCoordinatesLocal(dmNew, newCoordinates);
443: DMGetCoordinatesLocal(dm, &coordinates);
444: VecGetArray(coordinates, &coords);
445: VecGetArray(newCoordinates, &newCoords);
446: if (hasCells) {
447: for (c = cStart; c < cEnd; ++c) {
448: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;
450: PetscSectionGetDof(coordSection, c, &dof);
451: PetscSectionGetOffset(coordSection, c, &off);
452: PetscSectionGetOffset(newCoordSection, cNew, &noff);
453: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
454: }
455: }
456: for (v = vStart; v < vEnd; ++v) {
457: PetscInt dof, off, noff, d;
459: PetscSectionGetDof(coordSection, v, &dof);
460: PetscSectionGetOffset(coordSection, v, &off);
461: PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
462: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
463: }
464: VecRestoreArray(coordinates, &coords);
465: VecRestoreArray(newCoordinates, &newCoords);
466: VecDestroy(&newCoordinates);
467: PetscSectionDestroy(&newCoordSection);
468: return(0);
469: }
471: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
472: {
473: PetscInt depth = 0;
474: PetscSF sfPoint, sfPointNew;
475: const PetscSFNode *remotePoints;
476: PetscSFNode *gremotePoints;
477: const PetscInt *localPoints;
478: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
479: PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
480: PetscErrorCode ierr;
483: DMPlexGetDepth(dm, &depth);
484: /* Step 9: Convert pointSF */
485: DMGetPointSF(dm, &sfPoint);
486: DMGetPointSF(dmNew, &sfPointNew);
487: DMPlexGetChart(dm, &pStart, &pEnd);
488: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
489: totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
490: if (numRoots >= 0) {
491: PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
492: for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
493: PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
494: PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
495: PetscMalloc1(numLeaves, &glocalPoints);
496: PetscMalloc1(numLeaves, &gremotePoints);
497: for (l = 0; l < numLeaves; ++l) {
498: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
499: gremotePoints[l].rank = remotePoints[l].rank;
500: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
501: }
502: PetscFree2(newLocation,newRemoteLocation);
503: PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
504: }
505: return(0);
506: }
508: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
509: {
510: PetscSF sfPoint;
511: DMLabel vtkLabel, ghostLabel;
512: const PetscSFNode *leafRemote;
513: const PetscInt *leafLocal;
514: PetscInt depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
515: PetscMPIInt rank;
516: PetscErrorCode ierr;
519: DMPlexGetDepth(dm, &depth);
520: /* Step 10: Convert labels */
521: DMGetNumLabels(dm, &numLabels);
522: for (l = 0; l < numLabels; ++l) {
523: DMLabel label, newlabel;
524: const char *lname;
525: PetscBool isDepth;
526: IS valueIS;
527: const PetscInt *values;
528: PetscInt numValues, val;
530: DMGetLabelName(dm, l, &lname);
531: PetscStrcmp(lname, "depth", &isDepth);
532: if (isDepth) continue;
533: DMCreateLabel(dmNew, lname);
534: DMGetLabel(dm, lname, &label);
535: DMGetLabel(dmNew, lname, &newlabel);
536: DMLabelGetDefaultValue(label,&val);
537: DMLabelSetDefaultValue(newlabel,val);
538: DMLabelGetValueIS(label, &valueIS);
539: ISGetLocalSize(valueIS, &numValues);
540: ISGetIndices(valueIS, &values);
541: for (val = 0; val < numValues; ++val) {
542: IS pointIS;
543: const PetscInt *points;
544: PetscInt numPoints, p;
546: DMLabelGetStratumIS(label, values[val], &pointIS);
547: ISGetLocalSize(pointIS, &numPoints);
548: ISGetIndices(pointIS, &points);
549: for (p = 0; p < numPoints; ++p) {
550: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);
552: DMLabelSetValue(newlabel, newpoint, values[val]);
553: }
554: ISRestoreIndices(pointIS, &points);
555: ISDestroy(&pointIS);
556: }
557: ISRestoreIndices(valueIS, &values);
558: ISDestroy(&valueIS);
559: }
560: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
561: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
562: DMGetPointSF(dm, &sfPoint);
563: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
564: PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
565: DMCreateLabel(dmNew, "vtk");
566: DMCreateLabel(dmNew, "ghost");
567: DMGetLabel(dmNew, "vtk", &vtkLabel);
568: DMGetLabel(dmNew, "ghost", &ghostLabel);
569: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
570: for (; c < leafLocal[l] && c < cEnd; ++c) {
571: DMLabelSetValue(vtkLabel, c, 1);
572: }
573: if (leafLocal[l] >= cEnd) break;
574: if (leafRemote[l].rank == rank) {
575: DMLabelSetValue(vtkLabel, c, 1);
576: } else {
577: DMLabelSetValue(ghostLabel, c, 2);
578: }
579: }
580: for (; c < cEnd; ++c) {
581: DMLabelSetValue(vtkLabel, c, 1);
582: }
583: if (0) {
584: DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
585: }
586: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
587: for (f = fStart; f < fEnd; ++f) {
588: PetscInt numCells;
590: DMPlexGetSupportSize(dmNew, f, &numCells);
591: if (numCells < 2) {
592: DMLabelSetValue(ghostLabel, f, 1);
593: } else {
594: const PetscInt *cells = NULL;
595: PetscInt vA, vB;
597: DMPlexGetSupport(dmNew, f, &cells);
598: DMLabelGetValue(vtkLabel, cells[0], &vA);
599: DMLabelGetValue(vtkLabel, cells[1], &vB);
600: if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
601: }
602: }
603: if (0) {
604: DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
605: }
606: return(0);
607: }
609: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
610: {
611: DM refTree;
612: PetscSection pSec;
613: PetscInt *parents, *childIDs;
617: DMPlexGetReferenceTree(dm,&refTree);
618: DMPlexSetReferenceTree(dmNew,refTree);
619: DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
620: if (pSec) {
621: PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
622: PetscInt *childIDsShifted;
623: PetscSection pSecShifted;
625: PetscSectionGetChart(pSec,&pStart,&pEnd);
626: DMPlexGetDepth(dm,&depth);
627: pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
628: pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
629: PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
630: PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
631: PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
632: for (p = pStartShifted; p < pEndShifted; p++) {
633: /* start off assuming no children */
634: PetscSectionSetDof(pSecShifted,p,0);
635: }
636: for (p = pStart; p < pEnd; p++) {
637: PetscInt dof;
638: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
640: PetscSectionGetDof(pSec,p,&dof);
641: PetscSectionSetDof(pSecShifted,pNew,dof);
642: }
643: PetscSectionSetUp(pSecShifted);
644: for (p = pStart; p < pEnd; p++) {
645: PetscInt dof;
646: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
648: PetscSectionGetDof(pSec,p,&dof);
649: if (dof) {
650: PetscInt off, offNew;
652: PetscSectionGetOffset(pSec,p,&off);
653: PetscSectionGetOffset(pSecShifted,pNew,&offNew);
654: parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
655: childIDsShifted[offNew] = childIDs[off];
656: }
657: }
658: DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
659: PetscFree2(parentsShifted,childIDsShifted);
660: PetscSectionDestroy(&pSecShifted);
661: }
662: return(0);
663: }
665: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
666: {
667: PetscSF sf;
668: IS valueIS;
669: const PetscInt *values, *leaves;
670: PetscInt *depthShift;
671: PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
672: PetscBool isper;
673: const PetscReal *maxCell, *L;
674: const DMBoundaryType *bd;
675: PetscErrorCode ierr;
678: DMGetPointSF(dm, &sf);
679: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
680: nleaves = PetscMax(0, nleaves);
681: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
682: /* Count ghost cells */
683: DMLabelGetValueIS(label, &valueIS);
684: ISGetLocalSize(valueIS, &numFS);
685: ISGetIndices(valueIS, &values);
686: Ng = 0;
687: for (fs = 0; fs < numFS; ++fs) {
688: IS faceIS;
689: const PetscInt *faces;
690: PetscInt numFaces, f, numBdFaces = 0;
692: DMLabelGetStratumIS(label, values[fs], &faceIS);
693: ISGetLocalSize(faceIS, &numFaces);
694: ISGetIndices(faceIS, &faces);
695: for (f = 0; f < numFaces; ++f) {
696: PetscInt numChildren;
698: PetscFindInt(faces[f], nleaves, leaves, &loc);
699: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
700: /* non-local and ancestors points don't get to register ghosts */
701: if (loc >= 0 || numChildren) continue;
702: if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
703: }
704: Ng += numBdFaces;
705: ISDestroy(&faceIS);
706: }
707: DMPlexGetDepth(dm, &depth);
708: PetscMalloc1(2*(depth+1), &depthShift);
709: for (d = 0; d <= depth; d++) {
710: PetscInt dEnd;
712: DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
713: depthShift[2*d] = dEnd;
714: depthShift[2*d+1] = 0;
715: }
716: if (depth >= 0) depthShift[2*depth+1] = Ng;
717: DMPlexShiftPointSetUp_Internal(depth,depthShift);
718: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
719: /* Step 3: Set cone/support sizes for new points */
720: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
721: DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
722: for (c = cEnd; c < cEnd + Ng; ++c) {
723: DMPlexSetConeSize(gdm, c, 1);
724: }
725: for (fs = 0; fs < numFS; ++fs) {
726: IS faceIS;
727: const PetscInt *faces;
728: PetscInt numFaces, f;
730: DMLabelGetStratumIS(label, values[fs], &faceIS);
731: ISGetLocalSize(faceIS, &numFaces);
732: ISGetIndices(faceIS, &faces);
733: for (f = 0; f < numFaces; ++f) {
734: PetscInt size, numChildren;
736: PetscFindInt(faces[f], nleaves, leaves, &loc);
737: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
738: if (loc >= 0 || numChildren) continue;
739: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
740: DMPlexGetSupportSize(dm, faces[f], &size);
741: if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
742: DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
743: }
744: ISRestoreIndices(faceIS, &faces);
745: ISDestroy(&faceIS);
746: }
747: /* Step 4: Setup ghosted DM */
748: DMSetUp(gdm);
749: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
750: /* Step 6: Set cones and supports for new points */
751: ghostCell = cEnd;
752: for (fs = 0; fs < numFS; ++fs) {
753: IS faceIS;
754: const PetscInt *faces;
755: PetscInt numFaces, f;
757: DMLabelGetStratumIS(label, values[fs], &faceIS);
758: ISGetLocalSize(faceIS, &numFaces);
759: ISGetIndices(faceIS, &faces);
760: for (f = 0; f < numFaces; ++f) {
761: PetscInt newFace = faces[f] + Ng, numChildren;
763: PetscFindInt(faces[f], nleaves, leaves, &loc);
764: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
765: if (loc >= 0 || numChildren) continue;
766: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
767: DMPlexSetCone(gdm, ghostCell, &newFace);
768: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
769: ++ghostCell;
770: }
771: ISRestoreIndices(faceIS, &faces);
772: ISDestroy(&faceIS);
773: }
774: ISRestoreIndices(valueIS, &values);
775: ISDestroy(&valueIS);
776: /* Step 7: Stratify */
777: DMPlexStratify(gdm);
778: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
779: DMPlexShiftSF_Internal(dm, depthShift, gdm);
780: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
781: DMPlexShiftTree_Internal(dm, depthShift, gdm);
782: PetscFree(depthShift);
783: /* Step 7: Periodicity */
784: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
785: DMSetPeriodicity(gdm, isper, maxCell, L, bd);
786: if (numGhostCells) *numGhostCells = Ng;
787: return(0);
788: }
790: /*@C
791: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
793: Collective on dm
795: Input Parameters:
796: + dm - The original DM
797: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
799: Output Parameters:
800: + numGhostCells - The number of ghost cells added to the DM
801: - dmGhosted - The new DM
803: Note: If no label exists of that name, one will be created marking all boundary faces
805: Level: developer
807: .seealso: DMCreate()
808: @*/
809: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
810: {
811: DM gdm;
812: DMLabel label;
813: const char *name = labelName ? labelName : "Face Sets";
814: PetscInt dim;
815: PetscBool flag;
822: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
823: DMSetType(gdm, DMPLEX);
824: DMGetDimension(dm, &dim);
825: DMSetDimension(gdm, dim);
826: DMPlexGetAdjacencyUseCone(dm, &flag);
827: DMPlexSetAdjacencyUseCone(gdm, flag);
828: DMPlexGetAdjacencyUseClosure(dm, &flag);
829: DMPlexSetAdjacencyUseClosure(gdm, flag);
830: DMGetLabel(dm, name, &label);
831: if (!label) {
832: /* Get label for boundary faces */
833: DMCreateLabel(dm, name);
834: DMGetLabel(dm, name, &label);
835: DMPlexMarkBoundaryFaces(dm, label);
836: }
837: DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
838: DMCopyBoundary(dm, gdm);
839: *dmGhosted = gdm;
840: return(0);
841: }
843: /*
844: We are adding three kinds of points here:
845: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
846: Non-replicated: Points which exist on the fault, but are not replicated
847: Hybrid: Entirely new points, such as cohesive cells
849: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
850: each depth so that the new split/hybrid points can be inserted as a block.
851: */
852: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
853: {
854: MPI_Comm comm;
855: IS valueIS;
856: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
857: const PetscInt *values; /* List of depths for which we have replicated points */
858: IS *splitIS;
859: IS *unsplitIS;
860: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
861: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
862: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
863: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
864: const PetscInt **splitPoints; /* Replicated points for each depth */
865: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
866: PetscSection coordSection;
867: Vec coordinates;
868: PetscScalar *coords;
869: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
870: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
871: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
872: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
873: PetscInt *coneNew, *coneONew, *supportNew;
874: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
875: PetscErrorCode ierr;
878: PetscObjectGetComm((PetscObject)dm,&comm);
879: DMGetDimension(dm, &dim);
880: DMPlexGetDepth(dm, &depth);
881: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
882: /* Count split points and add cohesive cells */
883: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
884: PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
885: PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
886: DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
887: for (d = 0; d <= depth; ++d) {
888: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
889: depthEnd[d] = pMaxNew[d];
890: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
891: numSplitPoints[d] = 0;
892: numUnsplitPoints[d] = 0;
893: numHybridPoints[d] = 0;
894: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
895: splitPoints[d] = NULL;
896: unsplitPoints[d] = NULL;
897: splitIS[d] = NULL;
898: unsplitIS[d] = NULL;
899: /* we are shifting the existing hybrid points with the stratum behind them, so
900: * the split comes at the end of the normal points, i.e., at depthMax[d] */
901: depthShift[2*d] = depthMax[d];
902: depthShift[2*d+1] = 0;
903: }
904: if (label) {
905: DMLabelGetValueIS(label, &valueIS);
906: ISGetLocalSize(valueIS, &numSP);
907: ISGetIndices(valueIS, &values);
908: }
909: for (sp = 0; sp < numSP; ++sp) {
910: const PetscInt dep = values[sp];
912: if ((dep < 0) || (dep > depth)) continue;
913: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
914: if (splitIS[dep]) {
915: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
916: ISGetIndices(splitIS[dep], &splitPoints[dep]);
917: }
918: DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
919: if (unsplitIS[dep]) {
920: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
921: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
922: }
923: }
924: /* Calculate number of hybrid points */
925: 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 */
926: for (d = 0; d <= depth; ++d) depthShift[2*d+1] = numSplitPoints[d] + numHybridPoints[d];
927: DMPlexShiftPointSetUp_Internal(depth,depthShift);
928: /* the end of the points in this stratum that come before the new points:
929: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
930: * added points */
931: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
932: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
933: /* Step 3: Set cone/support sizes for new points */
934: for (dep = 0; dep <= depth; ++dep) {
935: for (p = 0; p < numSplitPoints[dep]; ++p) {
936: const PetscInt oldp = splitPoints[dep][p];
937: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
938: const PetscInt splitp = p + pMaxNew[dep];
939: const PetscInt *support;
940: PetscInt coneSize, supportSize, qf, qn, qp, e;
942: DMPlexGetConeSize(dm, oldp, &coneSize);
943: DMPlexSetConeSize(sdm, splitp, coneSize);
944: DMPlexGetSupportSize(dm, oldp, &supportSize);
945: DMPlexSetSupportSize(sdm, splitp, supportSize);
946: if (dep == depth-1) {
947: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
949: /* Add cohesive cells, they are prisms */
950: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
951: } else if (dep == 0) {
952: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
954: DMPlexGetSupport(dm, oldp, &support);
955: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
956: PetscInt val;
958: DMLabelGetValue(label, support[e], &val);
959: if (val == 1) ++qf;
960: if ((val == 1) || (val == (shift + 1))) ++qn;
961: if ((val == 1) || (val == -(shift + 1))) ++qp;
962: }
963: /* Split old vertex: Edges into original vertex and new cohesive edge */
964: DMPlexSetSupportSize(sdm, newp, qn+1);
965: /* Split new vertex: Edges into split vertex and new cohesive edge */
966: DMPlexSetSupportSize(sdm, splitp, qp+1);
967: /* Add hybrid edge */
968: DMPlexSetConeSize(sdm, hybedge, 2);
969: DMPlexSetSupportSize(sdm, hybedge, qf);
970: } else if (dep == dim-2) {
971: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
973: DMPlexGetSupport(dm, oldp, &support);
974: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
975: PetscInt val;
977: DMLabelGetValue(label, support[e], &val);
978: if (val == dim-1) ++qf;
979: if ((val == dim-1) || (val == (shift + dim-1))) ++qn;
980: if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
981: }
982: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
983: DMPlexSetSupportSize(sdm, newp, qn+1);
984: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
985: DMPlexSetSupportSize(sdm, splitp, qp+1);
986: /* Add hybrid face */
987: DMPlexSetConeSize(sdm, hybface, 4);
988: DMPlexSetSupportSize(sdm, hybface, qf);
989: }
990: }
991: }
992: for (dep = 0; dep <= depth; ++dep) {
993: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
994: const PetscInt oldp = unsplitPoints[dep][p];
995: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
996: const PetscInt *support;
997: PetscInt coneSize, supportSize, qf, e, s;
999: DMPlexGetConeSize(dm, oldp, &coneSize);
1000: DMPlexGetSupportSize(dm, oldp, &supportSize);
1001: DMPlexGetSupport(dm, oldp, &support);
1002: if (dep == 0) {
1003: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1005: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1006: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1007: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1008: if (e >= 0) ++qf;
1009: }
1010: DMPlexSetSupportSize(sdm, newp, qf+2);
1011: /* Add hybrid edge */
1012: DMPlexSetConeSize(sdm, hybedge, 2);
1013: for (e = 0, qf = 0; e < supportSize; ++e) {
1014: PetscInt val;
1016: DMLabelGetValue(label, support[e], &val);
1017: /* Split and unsplit edges produce hybrid faces */
1018: if (val == 1) ++qf;
1019: if (val == (shift2 + 1)) ++qf;
1020: }
1021: DMPlexSetSupportSize(sdm, hybedge, qf);
1022: } else if (dep == dim-2) {
1023: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1024: PetscInt val;
1026: for (e = 0, qf = 0; e < supportSize; ++e) {
1027: DMLabelGetValue(label, support[e], &val);
1028: if (val == dim-1) qf += 2;
1029: else ++qf;
1030: }
1031: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1032: DMPlexSetSupportSize(sdm, newp, qf+2);
1033: /* Add hybrid face */
1034: for (e = 0, qf = 0; e < supportSize; ++e) {
1035: DMLabelGetValue(label, support[e], &val);
1036: if (val == dim-1) ++qf;
1037: }
1038: DMPlexSetConeSize(sdm, hybface, 4);
1039: DMPlexSetSupportSize(sdm, hybface, qf);
1040: }
1041: }
1042: }
1043: /* Step 4: Setup split DM */
1044: DMSetUp(sdm);
1045: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1046: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1047: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1048: /* Step 6: Set cones and supports for new points */
1049: for (dep = 0; dep <= depth; ++dep) {
1050: for (p = 0; p < numSplitPoints[dep]; ++p) {
1051: const PetscInt oldp = splitPoints[dep][p];
1052: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1053: const PetscInt splitp = p + pMaxNew[dep];
1054: const PetscInt *cone, *support, *ornt;
1055: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1057: DMPlexGetConeSize(dm, oldp, &coneSize);
1058: DMPlexGetCone(dm, oldp, &cone);
1059: DMPlexGetConeOrientation(dm, oldp, &ornt);
1060: DMPlexGetSupportSize(dm, oldp, &supportSize);
1061: DMPlexGetSupport(dm, oldp, &support);
1062: if (dep == depth-1) {
1063: PetscBool hasUnsplit = PETSC_FALSE;
1064: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1065: const PetscInt *supportF;
1067: /* Split face: copy in old face to new face to start */
1068: DMPlexGetSupport(sdm, newp, &supportF);
1069: DMPlexSetSupport(sdm, splitp, supportF);
1070: /* Split old face: old vertices/edges in cone so no change */
1071: /* Split new face: new vertices/edges in cone */
1072: for (q = 0; q < coneSize; ++q) {
1073: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1074: if (v < 0) {
1075: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1076: 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);
1077: coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1078: hasUnsplit = PETSC_TRUE;
1079: } else {
1080: coneNew[2+q] = v + pMaxNew[dep-1];
1081: if (dep > 1) {
1082: const PetscInt *econe;
1083: PetscInt econeSize, r, vs, vu;
1085: DMPlexGetConeSize(dm, cone[q], &econeSize);
1086: DMPlexGetCone(dm, cone[q], &econe);
1087: for (r = 0; r < econeSize; ++r) {
1088: PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);
1089: PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1090: if (vs >= 0) continue;
1091: 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);
1092: hasUnsplit = PETSC_TRUE;
1093: }
1094: }
1095: }
1096: }
1097: DMPlexSetCone(sdm, splitp, &coneNew[2]);
1098: DMPlexSetConeOrientation(sdm, splitp, ornt);
1099: /* Face support */
1100: for (s = 0; s < supportSize; ++s) {
1101: PetscInt val;
1103: DMLabelGetValue(label, support[s], &val);
1104: if (val < 0) {
1105: /* Split old face: Replace negative side cell with cohesive cell */
1106: DMPlexInsertSupport(sdm, newp, s, hybcell);
1107: } else {
1108: /* Split new face: Replace positive side cell with cohesive cell */
1109: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1110: /* Get orientation for cohesive face */
1111: {
1112: const PetscInt *ncone, *nconeO;
1113: PetscInt nconeSize, nc;
1115: DMPlexGetConeSize(dm, support[s], &nconeSize);
1116: DMPlexGetCone(dm, support[s], &ncone);
1117: DMPlexGetConeOrientation(dm, support[s], &nconeO);
1118: for (nc = 0; nc < nconeSize; ++nc) {
1119: if (ncone[nc] == oldp) {
1120: coneONew[0] = nconeO[nc];
1121: break;
1122: }
1123: }
1124: if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1125: }
1126: }
1127: }
1128: /* Cohesive cell: Old and new split face, then new cohesive faces */
1129: coneNew[0] = newp; /* Extracted negative side orientation above */
1130: coneNew[1] = splitp;
1131: coneONew[1] = coneONew[0];
1132: for (q = 0; q < coneSize; ++q) {
1133: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1134: if (v < 0) {
1135: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1136: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1137: coneONew[2+q] = 0;
1138: } else {
1139: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep];
1140: }
1141: coneONew[2+q] = 0;
1142: }
1143: DMPlexSetCone(sdm, hybcell, coneNew);
1144: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1145: /* Label the hybrid cells on the boundary of the split */
1146: if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1147: } else if (dep == 0) {
1148: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1150: /* Split old vertex: Edges in old split faces and new cohesive edge */
1151: for (e = 0, qn = 0; e < supportSize; ++e) {
1152: PetscInt val;
1154: DMLabelGetValue(label, support[e], &val);
1155: if ((val == 1) || (val == (shift + 1))) {
1156: supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1157: }
1158: }
1159: supportNew[qn] = hybedge;
1160: DMPlexSetSupport(sdm, newp, supportNew);
1161: /* Split new vertex: Edges in new split faces and new cohesive edge */
1162: for (e = 0, qp = 0; e < supportSize; ++e) {
1163: PetscInt val, edge;
1165: DMLabelGetValue(label, support[e], &val);
1166: if (val == 1) {
1167: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1168: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1169: supportNew[qp++] = edge + pMaxNew[dep+1];
1170: } else if (val == -(shift + 1)) {
1171: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1172: }
1173: }
1174: supportNew[qp] = hybedge;
1175: DMPlexSetSupport(sdm, splitp, supportNew);
1176: /* Hybrid edge: Old and new split vertex */
1177: coneNew[0] = newp;
1178: coneNew[1] = splitp;
1179: DMPlexSetCone(sdm, hybedge, coneNew);
1180: for (e = 0, qf = 0; e < supportSize; ++e) {
1181: PetscInt val, edge;
1183: DMLabelGetValue(label, support[e], &val);
1184: if (val == 1) {
1185: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1186: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1187: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1188: }
1189: }
1190: DMPlexSetSupport(sdm, hybedge, supportNew);
1191: } else if (dep == dim-2) {
1192: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1194: /* Split old edge: old vertices in cone so no change */
1195: /* Split new edge: new vertices in cone */
1196: for (q = 0; q < coneSize; ++q) {
1197: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1198: if (v < 0) {
1199: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1200: 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);
1201: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1202: } else {
1203: coneNew[q] = v + pMaxNew[dep-1];
1204: }
1205: }
1206: DMPlexSetCone(sdm, splitp, coneNew);
1207: /* Split old edge: Faces in positive side cells and old split faces */
1208: for (e = 0, q = 0; e < supportSize; ++e) {
1209: PetscInt val;
1211: DMLabelGetValue(label, support[e], &val);
1212: if (val == dim-1) {
1213: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1214: } else if (val == (shift + dim-1)) {
1215: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1216: }
1217: }
1218: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1219: DMPlexSetSupport(sdm, newp, supportNew);
1220: /* Split new edge: Faces in negative side cells and new split faces */
1221: for (e = 0, q = 0; e < supportSize; ++e) {
1222: PetscInt val, face;
1224: DMLabelGetValue(label, support[e], &val);
1225: if (val == dim-1) {
1226: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1227: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1228: supportNew[q++] = face + pMaxNew[dep+1];
1229: } else if (val == -(shift + dim-1)) {
1230: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1231: }
1232: }
1233: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1234: DMPlexSetSupport(sdm, splitp, supportNew);
1235: /* Hybrid face */
1236: coneNew[0] = newp;
1237: coneNew[1] = splitp;
1238: for (v = 0; v < coneSize; ++v) {
1239: PetscInt vertex;
1240: PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1241: if (vertex < 0) {
1242: PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1243: 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);
1244: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1245: } else {
1246: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1247: }
1248: }
1249: DMPlexSetCone(sdm, hybface, coneNew);
1250: for (e = 0, qf = 0; e < supportSize; ++e) {
1251: PetscInt val, face;
1253: DMLabelGetValue(label, support[e], &val);
1254: if (val == dim-1) {
1255: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1256: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1257: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1258: }
1259: }
1260: DMPlexSetSupport(sdm, hybface, supportNew);
1261: }
1262: }
1263: }
1264: for (dep = 0; dep <= depth; ++dep) {
1265: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1266: const PetscInt oldp = unsplitPoints[dep][p];
1267: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1268: const PetscInt *cone, *support, *ornt;
1269: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1271: DMPlexGetConeSize(dm, oldp, &coneSize);
1272: DMPlexGetCone(dm, oldp, &cone);
1273: DMPlexGetConeOrientation(dm, oldp, &ornt);
1274: DMPlexGetSupportSize(dm, oldp, &supportSize);
1275: DMPlexGetSupport(dm, oldp, &support);
1276: if (dep == 0) {
1277: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1279: /* Unsplit vertex */
1280: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1281: for (s = 0, q = 0; s < supportSize; ++s) {
1282: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1283: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1284: if (e >= 0) {
1285: supportNew[q++] = e + pMaxNew[dep+1];
1286: }
1287: }
1288: supportNew[q++] = hybedge;
1289: supportNew[q++] = hybedge;
1290: if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1291: DMPlexSetSupport(sdm, newp, supportNew);
1292: /* Hybrid edge */
1293: coneNew[0] = newp;
1294: coneNew[1] = newp;
1295: DMPlexSetCone(sdm, hybedge, coneNew);
1296: for (e = 0, qf = 0; e < supportSize; ++e) {
1297: PetscInt val, edge;
1299: DMLabelGetValue(label, support[e], &val);
1300: if (val == 1) {
1301: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1302: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1303: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1304: } else if (val == (shift2 + 1)) {
1305: PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1306: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1307: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1308: }
1309: }
1310: DMPlexSetSupport(sdm, hybedge, supportNew);
1311: } else if (dep == dim-2) {
1312: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1314: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1315: for (f = 0, qf = 0; f < supportSize; ++f) {
1316: PetscInt val, face;
1318: DMLabelGetValue(label, support[f], &val);
1319: if (val == dim-1) {
1320: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1321: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1322: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1323: supportNew[qf++] = face + pMaxNew[dep+1];
1324: } else {
1325: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1326: }
1327: }
1328: supportNew[qf++] = hybface;
1329: supportNew[qf++] = hybface;
1330: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1331: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1332: DMPlexSetSupport(sdm, newp, supportNew);
1333: /* Add hybrid face */
1334: coneNew[0] = newp;
1335: coneNew[1] = newp;
1336: PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1337: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1338: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1339: PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1340: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1341: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1342: DMPlexSetCone(sdm, hybface, coneNew);
1343: for (f = 0, qf = 0; f < supportSize; ++f) {
1344: PetscInt val, face;
1346: DMLabelGetValue(label, support[f], &val);
1347: if (val == dim-1) {
1348: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1349: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1350: }
1351: }
1352: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1353: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1354: DMPlexSetSupport(sdm, hybface, supportNew);
1355: }
1356: }
1357: }
1358: /* Step 6b: Replace split points in negative side cones */
1359: for (sp = 0; sp < numSP; ++sp) {
1360: PetscInt dep = values[sp];
1361: IS pIS;
1362: PetscInt numPoints;
1363: const PetscInt *points;
1365: if (dep >= 0) continue;
1366: DMLabelGetStratumIS(label, dep, &pIS);
1367: if (!pIS) continue;
1368: dep = -dep - shift;
1369: ISGetLocalSize(pIS, &numPoints);
1370: ISGetIndices(pIS, &points);
1371: for (p = 0; p < numPoints; ++p) {
1372: const PetscInt oldp = points[p];
1373: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1374: const PetscInt *cone;
1375: PetscInt coneSize, c;
1376: /* PetscBool replaced = PETSC_FALSE; */
1378: /* Negative edge: replace split vertex */
1379: /* Negative cell: replace split face */
1380: DMPlexGetConeSize(sdm, newp, &coneSize);
1381: DMPlexGetCone(sdm, newp, &cone);
1382: for (c = 0; c < coneSize; ++c) {
1383: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1384: PetscInt csplitp, cp, val;
1386: DMLabelGetValue(label, coldp, &val);
1387: if (val == dep-1) {
1388: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1389: if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1390: csplitp = pMaxNew[dep-1] + cp;
1391: DMPlexInsertCone(sdm, newp, c, csplitp);
1392: /* replaced = PETSC_TRUE; */
1393: }
1394: }
1395: /* Cells with only a vertex or edge on the submesh have no replacement */
1396: /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1397: }
1398: ISRestoreIndices(pIS, &points);
1399: ISDestroy(&pIS);
1400: }
1401: /* Step 7: Stratify */
1402: DMPlexStratify(sdm);
1403: /* Step 8: Coordinates */
1404: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1405: DMGetCoordinateSection(sdm, &coordSection);
1406: DMGetCoordinatesLocal(sdm, &coordinates);
1407: VecGetArray(coordinates, &coords);
1408: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1409: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1410: const PetscInt splitp = pMaxNew[0] + v;
1411: PetscInt dof, off, soff, d;
1413: PetscSectionGetDof(coordSection, newp, &dof);
1414: PetscSectionGetOffset(coordSection, newp, &off);
1415: PetscSectionGetOffset(coordSection, splitp, &soff);
1416: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1417: }
1418: VecRestoreArray(coordinates, &coords);
1419: /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1420: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1421: /* Step 10: Labels */
1422: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1423: DMGetNumLabels(sdm, &numLabels);
1424: for (dep = 0; dep <= depth; ++dep) {
1425: for (p = 0; p < numSplitPoints[dep]; ++p) {
1426: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1427: const PetscInt splitp = pMaxNew[dep] + p;
1428: PetscInt l;
1430: for (l = 0; l < numLabels; ++l) {
1431: DMLabel mlabel;
1432: const char *lname;
1433: PetscInt val;
1434: PetscBool isDepth;
1436: DMGetLabelName(sdm, l, &lname);
1437: PetscStrcmp(lname, "depth", &isDepth);
1438: if (isDepth) continue;
1439: DMGetLabel(sdm, lname, &mlabel);
1440: DMLabelGetValue(mlabel, newp, &val);
1441: if (val >= 0) {
1442: DMLabelSetValue(mlabel, splitp, val);
1443: }
1444: }
1445: }
1446: }
1447: for (sp = 0; sp < numSP; ++sp) {
1448: const PetscInt dep = values[sp];
1450: if ((dep < 0) || (dep > depth)) continue;
1451: if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1452: ISDestroy(&splitIS[dep]);
1453: if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1454: ISDestroy(&unsplitIS[dep]);
1455: }
1456: if (label) {
1457: ISRestoreIndices(valueIS, &values);
1458: ISDestroy(&valueIS);
1459: }
1460: for (d = 0; d <= depth; ++d) {
1461: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1462: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1463: }
1464: DMPlexSetHybridBounds(sdm, depth >= 0 ? pMaxNew[depth] : PETSC_DETERMINE, depth>1 ? pMaxNew[depth-1] : PETSC_DETERMINE, depth>2 ? pMaxNew[1] : PETSC_DETERMINE, depth >= 0 ? pMaxNew[0] : PETSC_DETERMINE);
1465: PetscFree3(coneNew, coneONew, supportNew);
1466: PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1467: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1468: return(0);
1469: }
1471: /*@C
1472: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1474: Collective on dm
1476: Input Parameters:
1477: + dm - The original DM
1478: - label - The label specifying the boundary faces (this could be auto-generated)
1480: Output Parameters:
1481: - dmSplit - The new DM
1483: Level: developer
1485: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1486: @*/
1487: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1488: {
1489: DM sdm;
1490: PetscInt dim;
1496: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1497: DMSetType(sdm, DMPLEX);
1498: DMGetDimension(dm, &dim);
1499: DMSetDimension(sdm, dim);
1500: switch (dim) {
1501: case 2:
1502: case 3:
1503: DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1504: break;
1505: default:
1506: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1507: }
1508: *dmSplit = sdm;
1509: return(0);
1510: }
1512: /* Returns the side of the surface for a given cell with a face on the surface */
1513: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1514: {
1515: const PetscInt *cone, *ornt;
1516: PetscInt dim, coneSize, c;
1517: PetscErrorCode ierr;
1520: *pos = PETSC_TRUE;
1521: DMGetDimension(dm, &dim);
1522: DMPlexGetConeSize(dm, cell, &coneSize);
1523: DMPlexGetCone(dm, cell, &cone);
1524: DMPlexGetConeOrientation(dm, cell, &ornt);
1525: for (c = 0; c < coneSize; ++c) {
1526: if (cone[c] == face) {
1527: PetscInt o = ornt[c];
1529: if (subdm) {
1530: const PetscInt *subcone, *subornt;
1531: PetscInt subpoint, subface, subconeSize, sc;
1533: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1534: PetscFindInt(face, numSubpoints, subpoints, &subface);
1535: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1536: DMPlexGetCone(subdm, subpoint, &subcone);
1537: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1538: for (sc = 0; sc < subconeSize; ++sc) {
1539: if (subcone[sc] == subface) {
1540: o = subornt[0];
1541: break;
1542: }
1543: }
1544: 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);
1545: }
1546: if (o >= 0) *pos = PETSC_TRUE;
1547: else *pos = PETSC_FALSE;
1548: break;
1549: }
1550: }
1551: 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);
1552: return(0);
1553: }
1555: /*@
1556: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1557: to complete the surface
1559: Input Parameters:
1560: + dm - The DM
1561: . label - A DMLabel marking the surface
1562: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1563: . flip - Flag to flip the submesh normal and replace points on the other side
1564: - subdm - The subDM associated with the label, or NULL
1566: Output Parameter:
1567: . label - A DMLabel marking all surface points
1569: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1571: Level: developer
1573: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1574: @*/
1575: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1576: {
1577: DMLabel depthLabel;
1578: IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1579: const PetscInt *points, *subpoints;
1580: const PetscInt rev = flip ? -1 : 1;
1581: PetscInt *pMax;
1582: PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1583: PetscErrorCode ierr;
1586: DMPlexGetDepth(dm, &depth);
1587: DMGetDimension(dm, &dim);
1588: pSize = PetscMax(depth, dim) + 1;
1589: PetscMalloc1(pSize,&pMax);
1590: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1591: DMPlexGetDepthLabel(dm, &depthLabel);
1592: DMGetDimension(dm, &dim);
1593: if (subdm) {
1594: DMPlexCreateSubpointIS(subdm, &subpointIS);
1595: if (subpointIS) {
1596: ISGetLocalSize(subpointIS, &numSubpoints);
1597: ISGetIndices(subpointIS, &subpoints);
1598: }
1599: }
1600: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1601: DMLabelGetStratumIS(label, dim-1, &dimIS);
1602: if (!dimIS) {
1603: PetscFree(pMax);
1604: ISDestroy(&subpointIS);
1605: return(0);
1606: }
1607: ISGetLocalSize(dimIS, &numPoints);
1608: ISGetIndices(dimIS, &points);
1609: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1610: const PetscInt *support;
1611: PetscInt supportSize, s;
1613: DMPlexGetSupportSize(dm, points[p], &supportSize);
1614: if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1615: DMPlexGetSupport(dm, points[p], &support);
1616: for (s = 0; s < supportSize; ++s) {
1617: const PetscInt *cone;
1618: PetscInt coneSize, c;
1619: PetscBool pos;
1621: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1622: if (pos) {DMLabelSetValue(label, support[s], rev*(shift+dim));}
1623: else {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1624: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1625: /* Put faces touching the fault in the label */
1626: DMPlexGetConeSize(dm, support[s], &coneSize);
1627: DMPlexGetCone(dm, support[s], &cone);
1628: for (c = 0; c < coneSize; ++c) {
1629: const PetscInt point = cone[c];
1631: DMLabelGetValue(label, point, &val);
1632: if (val == -1) {
1633: PetscInt *closure = NULL;
1634: PetscInt closureSize, cl;
1636: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1637: for (cl = 0; cl < closureSize*2; cl += 2) {
1638: const PetscInt clp = closure[cl];
1639: PetscInt bval = -1;
1641: DMLabelGetValue(label, clp, &val);
1642: if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1643: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1644: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1645: break;
1646: }
1647: }
1648: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1649: }
1650: }
1651: }
1652: }
1653: ISRestoreIndices(dimIS, &points);
1654: ISDestroy(&dimIS);
1655: if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1656: ISDestroy(&subpointIS);
1657: /* Mark boundary points as unsplit */
1658: if (blabel) {
1659: DMLabelGetStratumIS(blabel, 1, &dimIS);
1660: ISGetLocalSize(dimIS, &numPoints);
1661: ISGetIndices(dimIS, &points);
1662: for (p = 0; p < numPoints; ++p) {
1663: const PetscInt point = points[p];
1664: PetscInt val, bval;
1666: DMLabelGetValue(blabel, point, &bval);
1667: if (bval >= 0) {
1668: DMLabelGetValue(label, point, &val);
1669: if ((val < 0) || (val > dim)) {
1670: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1671: DMLabelClearValue(blabel, point, bval);
1672: }
1673: }
1674: }
1675: for (p = 0; p < numPoints; ++p) {
1676: const PetscInt point = points[p];
1677: PetscInt val, bval;
1679: DMLabelGetValue(blabel, point, &bval);
1680: if (bval >= 0) {
1681: const PetscInt *cone, *support;
1682: PetscInt coneSize, supportSize, s, valA, valB, valE;
1684: /* Mark as unsplit */
1685: DMLabelGetValue(label, point, &val);
1686: 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);
1687: DMLabelClearValue(label, point, val);
1688: DMLabelSetValue(label, point, shift2+val);
1689: /* Check for cross-edge
1690: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1691: if (val != 0) continue;
1692: DMPlexGetSupport(dm, point, &support);
1693: DMPlexGetSupportSize(dm, point, &supportSize);
1694: for (s = 0; s < supportSize; ++s) {
1695: DMPlexGetCone(dm, support[s], &cone);
1696: DMPlexGetConeSize(dm, support[s], &coneSize);
1697: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1698: DMLabelGetValue(blabel, cone[0], &valA);
1699: DMLabelGetValue(blabel, cone[1], &valB);
1700: DMLabelGetValue(blabel, support[s], &valE);
1701: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1702: }
1703: }
1704: }
1705: ISRestoreIndices(dimIS, &points);
1706: ISDestroy(&dimIS);
1707: }
1708: /* Search for other cells/faces/edges connected to the fault by a vertex */
1709: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1710: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1711: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1712: cMax = cMax < 0 ? cEnd : cMax;
1713: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1714: DMLabelGetStratumIS(label, 0, &dimIS);
1715: if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1716: if (dimIS && crossEdgeIS) {
1717: IS vertIS = dimIS;
1719: ISExpand(vertIS, crossEdgeIS, &dimIS);
1720: ISDestroy(&crossEdgeIS);
1721: ISDestroy(&vertIS);
1722: }
1723: if (!dimIS) {
1724: PetscFree(pMax);
1725: return(0);
1726: }
1727: ISGetLocalSize(dimIS, &numPoints);
1728: ISGetIndices(dimIS, &points);
1729: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1730: PetscInt *star = NULL;
1731: PetscInt starSize, s;
1732: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1734: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1735: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1736: while (again) {
1737: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1738: again = 0;
1739: for (s = 0; s < starSize*2; s += 2) {
1740: const PetscInt point = star[s];
1741: const PetscInt *cone;
1742: PetscInt coneSize, c;
1744: if ((point < cStart) || (point >= cMax)) continue;
1745: DMLabelGetValue(label, point, &val);
1746: if (val != -1) continue;
1747: again = again == 1 ? 1 : 2;
1748: DMPlexGetConeSize(dm, point, &coneSize);
1749: DMPlexGetCone(dm, point, &cone);
1750: for (c = 0; c < coneSize; ++c) {
1751: DMLabelGetValue(label, cone[c], &val);
1752: if (val != -1) {
1753: const PetscInt *ccone;
1754: PetscInt cconeSize, cc, side;
1756: 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);
1757: if (val > 0) side = 1;
1758: else side = -1;
1759: DMLabelSetValue(label, point, side*(shift+dim));
1760: /* Mark cell faces which touch the fault */
1761: DMPlexGetConeSize(dm, point, &cconeSize);
1762: DMPlexGetCone(dm, point, &ccone);
1763: for (cc = 0; cc < cconeSize; ++cc) {
1764: PetscInt *closure = NULL;
1765: PetscInt closureSize, cl;
1767: DMLabelGetValue(label, ccone[cc], &val);
1768: if (val != -1) continue;
1769: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1770: for (cl = 0; cl < closureSize*2; cl += 2) {
1771: const PetscInt clp = closure[cl];
1773: DMLabelGetValue(label, clp, &val);
1774: if (val == -1) continue;
1775: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1776: break;
1777: }
1778: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1779: }
1780: again = 1;
1781: break;
1782: }
1783: }
1784: }
1785: }
1786: /* Classify the rest by cell membership */
1787: for (s = 0; s < starSize*2; s += 2) {
1788: const PetscInt point = star[s];
1790: DMLabelGetValue(label, point, &val);
1791: if (val == -1) {
1792: PetscInt *sstar = NULL;
1793: PetscInt sstarSize, ss;
1794: PetscBool marked = PETSC_FALSE;
1796: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1797: for (ss = 0; ss < sstarSize*2; ss += 2) {
1798: const PetscInt spoint = sstar[ss];
1800: if ((spoint < cStart) || (spoint >= cMax)) continue;
1801: DMLabelGetValue(label, spoint, &val);
1802: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1803: DMLabelGetValue(depthLabel, point, &dep);
1804: if (val > 0) {
1805: DMLabelSetValue(label, point, shift+dep);
1806: } else {
1807: DMLabelSetValue(label, point, -(shift+dep));
1808: }
1809: marked = PETSC_TRUE;
1810: break;
1811: }
1812: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1813: DMLabelGetValue(depthLabel, point, &dep);
1814: if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1815: }
1816: }
1817: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1818: }
1819: ISRestoreIndices(dimIS, &points);
1820: ISDestroy(&dimIS);
1821: /* If any faces touching the fault divide cells on either side, split them */
1822: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1823: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1824: ISExpand(facePosIS, faceNegIS, &dimIS);
1825: ISDestroy(&facePosIS);
1826: ISDestroy(&faceNegIS);
1827: ISGetLocalSize(dimIS, &numPoints);
1828: ISGetIndices(dimIS, &points);
1829: for (p = 0; p < numPoints; ++p) {
1830: const PetscInt point = points[p];
1831: const PetscInt *support;
1832: PetscInt supportSize, valA, valB;
1834: DMPlexGetSupportSize(dm, point, &supportSize);
1835: if (supportSize != 2) continue;
1836: DMPlexGetSupport(dm, point, &support);
1837: DMLabelGetValue(label, support[0], &valA);
1838: DMLabelGetValue(label, support[1], &valB);
1839: if ((valA == -1) || (valB == -1)) continue;
1840: if (valA*valB > 0) continue;
1841: /* Split the face */
1842: DMLabelGetValue(label, point, &valA);
1843: DMLabelClearValue(label, point, valA);
1844: DMLabelSetValue(label, point, dim-1);
1845: /* Label its closure:
1846: unmarked: label as unsplit
1847: incident: relabel as split
1848: split: do nothing
1849: */
1850: {
1851: PetscInt *closure = NULL;
1852: PetscInt closureSize, cl;
1854: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1855: for (cl = 0; cl < closureSize*2; cl += 2) {
1856: DMLabelGetValue(label, closure[cl], &valA);
1857: if (valA == -1) { /* Mark as unsplit */
1858: DMLabelGetValue(depthLabel, closure[cl], &dep);
1859: DMLabelSetValue(label, closure[cl], shift2+dep);
1860: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1861: DMLabelGetValue(depthLabel, closure[cl], &dep);
1862: DMLabelClearValue(label, closure[cl], valA);
1863: DMLabelSetValue(label, closure[cl], dep);
1864: }
1865: }
1866: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1867: }
1868: }
1869: ISRestoreIndices(dimIS, &points);
1870: ISDestroy(&dimIS);
1871: PetscFree(pMax);
1872: return(0);
1873: }
1875: /*@
1876: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
1878: Collective on dm
1880: Input Parameters:
1881: + dm - The original DM
1882: - labelName - The label specifying the interface vertices
1884: Output Parameters:
1885: + hybridLabel - The label fully marking the interface
1886: - dmHybrid - The new DM
1888: Level: developer
1890: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1891: @*/
1892: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1893: {
1894: DM idm;
1895: DMLabel subpointMap, hlabel;
1896: PetscInt dim;
1903: DMGetDimension(dm, &dim);
1904: DMPlexCreateSubmesh(dm, label, 1, &idm);
1905: DMPlexOrient(idm);
1906: DMPlexGetSubpointMap(idm, &subpointMap);
1907: DMLabelDuplicate(subpointMap, &hlabel);
1908: DMLabelClearStratum(hlabel, dim);
1909: DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1910: DMDestroy(&idm);
1911: DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1912: if (hybridLabel) *hybridLabel = hlabel;
1913: else {DMLabelDestroy(&hlabel);}
1914: return(0);
1915: }
1917: /* Here we need the explicit assumption that:
1919: For any marked cell, the marked vertices constitute a single face
1920: */
1921: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1922: {
1923: IS subvertexIS = NULL;
1924: const PetscInt *subvertices;
1925: PetscInt *pStart, *pEnd, *pMax, pSize;
1926: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
1927: PetscErrorCode ierr;
1930: *numFaces = 0;
1931: *nFV = 0;
1932: DMPlexGetDepth(dm, &depth);
1933: DMGetDimension(dm, &dim);
1934: pSize = PetscMax(depth, dim) + 1;
1935: PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1936: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1937: for (d = 0; d <= depth; ++d) {
1938: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1939: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1940: }
1941: /* Loop over initial vertices and mark all faces in the collective star() */
1942: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1943: if (subvertexIS) {
1944: ISGetSize(subvertexIS, &numSubVerticesInitial);
1945: ISGetIndices(subvertexIS, &subvertices);
1946: }
1947: for (v = 0; v < numSubVerticesInitial; ++v) {
1948: const PetscInt vertex = subvertices[v];
1949: PetscInt *star = NULL;
1950: PetscInt starSize, s, numCells = 0, c;
1952: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1953: for (s = 0; s < starSize*2; s += 2) {
1954: const PetscInt point = star[s];
1955: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1956: }
1957: for (c = 0; c < numCells; ++c) {
1958: const PetscInt cell = star[c];
1959: PetscInt *closure = NULL;
1960: PetscInt closureSize, cl;
1961: PetscInt cellLoc, numCorners = 0, faceSize = 0;
1963: DMLabelGetValue(subpointMap, cell, &cellLoc);
1964: if (cellLoc == 2) continue;
1965: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1966: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1967: for (cl = 0; cl < closureSize*2; cl += 2) {
1968: const PetscInt point = closure[cl];
1969: PetscInt vertexLoc;
1971: if ((point >= pStart[0]) && (point < pEnd[0])) {
1972: ++numCorners;
1973: DMLabelGetValue(vertexLabel, point, &vertexLoc);
1974: if (vertexLoc == value) closure[faceSize++] = point;
1975: }
1976: }
1977: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1978: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1979: if (faceSize == *nFV) {
1980: const PetscInt *cells = NULL;
1981: PetscInt numCells, nc;
1983: ++(*numFaces);
1984: for (cl = 0; cl < faceSize; ++cl) {
1985: DMLabelSetValue(subpointMap, closure[cl], 0);
1986: }
1987: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1988: for (nc = 0; nc < numCells; ++nc) {
1989: DMLabelSetValue(subpointMap, cells[nc], 2);
1990: }
1991: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1992: }
1993: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1994: }
1995: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1996: }
1997: if (subvertexIS) {
1998: ISRestoreIndices(subvertexIS, &subvertices);
1999: }
2000: ISDestroy(&subvertexIS);
2001: PetscFree3(pStart,pEnd,pMax);
2002: return(0);
2003: }
2005: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
2006: {
2007: IS subvertexIS = NULL;
2008: const PetscInt *subvertices;
2009: PetscInt *pStart, *pEnd, *pMax;
2010: PetscInt dim, d, numSubVerticesInitial = 0, v;
2011: PetscErrorCode ierr;
2014: DMGetDimension(dm, &dim);
2015: PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2016: DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2017: for (d = 0; d <= dim; ++d) {
2018: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2019: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2020: }
2021: /* Loop over initial vertices and mark all faces in the collective star() */
2022: if (vertexLabel) {
2023: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2024: if (subvertexIS) {
2025: ISGetSize(subvertexIS, &numSubVerticesInitial);
2026: ISGetIndices(subvertexIS, &subvertices);
2027: }
2028: }
2029: for (v = 0; v < numSubVerticesInitial; ++v) {
2030: const PetscInt vertex = subvertices[v];
2031: PetscInt *star = NULL;
2032: PetscInt starSize, s, numFaces = 0, f;
2034: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2035: for (s = 0; s < starSize*2; s += 2) {
2036: const PetscInt point = star[s];
2037: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
2038: }
2039: for (f = 0; f < numFaces; ++f) {
2040: const PetscInt face = star[f];
2041: PetscInt *closure = NULL;
2042: PetscInt closureSize, c;
2043: PetscInt faceLoc;
2045: DMLabelGetValue(subpointMap, face, &faceLoc);
2046: if (faceLoc == dim-1) continue;
2047: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2048: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2049: for (c = 0; c < closureSize*2; c += 2) {
2050: const PetscInt point = closure[c];
2051: PetscInt vertexLoc;
2053: if ((point >= pStart[0]) && (point < pEnd[0])) {
2054: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2055: if (vertexLoc != value) break;
2056: }
2057: }
2058: if (c == closureSize*2) {
2059: const PetscInt *support;
2060: PetscInt supportSize, s;
2062: for (c = 0; c < closureSize*2; c += 2) {
2063: const PetscInt point = closure[c];
2065: for (d = 0; d < dim; ++d) {
2066: if ((point >= pStart[d]) && (point < pEnd[d])) {
2067: DMLabelSetValue(subpointMap, point, d);
2068: break;
2069: }
2070: }
2071: }
2072: DMPlexGetSupportSize(dm, face, &supportSize);
2073: DMPlexGetSupport(dm, face, &support);
2074: for (s = 0; s < supportSize; ++s) {
2075: DMLabelSetValue(subpointMap, support[s], dim);
2076: }
2077: }
2078: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2079: }
2080: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2081: }
2082: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2083: ISDestroy(&subvertexIS);
2084: PetscFree3(pStart,pEnd,pMax);
2085: return(0);
2086: }
2088: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2089: {
2090: DMLabel label = NULL;
2091: const PetscInt *cone;
2092: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2093: PetscErrorCode ierr;
2096: *numFaces = 0;
2097: *nFV = 0;
2098: if (labelname) {DMGetLabel(dm, labelname, &label);}
2099: *subCells = NULL;
2100: DMGetDimension(dm, &dim);
2101: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2102: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2103: if (cMax < 0) return(0);
2104: if (label) {
2105: for (c = cMax; c < cEnd; ++c) {
2106: PetscInt val;
2108: DMLabelGetValue(label, c, &val);
2109: if (val == value) {
2110: ++(*numFaces);
2111: DMPlexGetConeSize(dm, c, &coneSize);
2112: }
2113: }
2114: } else {
2115: *numFaces = cEnd - cMax;
2116: DMPlexGetConeSize(dm, cMax, &coneSize);
2117: }
2118: PetscMalloc1(*numFaces *2, subCells);
2119: if (!(*numFaces)) return(0);
2120: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2121: for (c = cMax; c < cEnd; ++c) {
2122: const PetscInt *cells;
2123: PetscInt numCells;
2125: if (label) {
2126: PetscInt val;
2128: DMLabelGetValue(label, c, &val);
2129: if (val != value) continue;
2130: }
2131: DMPlexGetCone(dm, c, &cone);
2132: for (p = 0; p < *nFV; ++p) {
2133: DMLabelSetValue(subpointMap, cone[p], 0);
2134: }
2135: /* Negative face */
2136: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2137: /* Not true in parallel
2138: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2139: for (p = 0; p < numCells; ++p) {
2140: DMLabelSetValue(subpointMap, cells[p], 2);
2141: (*subCells)[subc++] = cells[p];
2142: }
2143: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2144: /* Positive face is not included */
2145: }
2146: return(0);
2147: }
2149: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2150: {
2151: PetscInt *pStart, *pEnd;
2152: PetscInt dim, cMax, cEnd, c, d;
2156: DMGetDimension(dm, &dim);
2157: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2158: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2159: if (cMax < 0) return(0);
2160: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2161: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2162: for (c = cMax; c < cEnd; ++c) {
2163: const PetscInt *cone;
2164: PetscInt *closure = NULL;
2165: PetscInt fconeSize, coneSize, closureSize, cl, val;
2167: if (label) {
2168: DMLabelGetValue(label, c, &val);
2169: if (val != value) continue;
2170: }
2171: DMPlexGetConeSize(dm, c, &coneSize);
2172: DMPlexGetCone(dm, c, &cone);
2173: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2174: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2175: /* Negative face */
2176: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2177: for (cl = 0; cl < closureSize*2; cl += 2) {
2178: const PetscInt point = closure[cl];
2180: for (d = 0; d <= dim; ++d) {
2181: if ((point >= pStart[d]) && (point < pEnd[d])) {
2182: DMLabelSetValue(subpointMap, point, d);
2183: break;
2184: }
2185: }
2186: }
2187: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2188: /* Cells -- positive face is not included */
2189: for (cl = 0; cl < 1; ++cl) {
2190: const PetscInt *support;
2191: PetscInt supportSize, s;
2193: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2194: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2195: DMPlexGetSupport(dm, cone[cl], &support);
2196: for (s = 0; s < supportSize; ++s) {
2197: DMLabelSetValue(subpointMap, support[s], dim);
2198: }
2199: }
2200: }
2201: PetscFree2(pStart, pEnd);
2202: return(0);
2203: }
2205: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2206: {
2207: MPI_Comm comm;
2208: PetscBool posOrient = PETSC_FALSE;
2209: const PetscInt debug = 0;
2210: PetscInt cellDim, faceSize, f;
2214: PetscObjectGetComm((PetscObject)dm,&comm);
2215: DMGetDimension(dm, &cellDim);
2216: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2218: if (cellDim == 1 && numCorners == 2) {
2219: /* Triangle */
2220: faceSize = numCorners-1;
2221: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2222: } else if (cellDim == 2 && numCorners == 3) {
2223: /* Triangle */
2224: faceSize = numCorners-1;
2225: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2226: } else if (cellDim == 3 && numCorners == 4) {
2227: /* Tetrahedron */
2228: faceSize = numCorners-1;
2229: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2230: } else if (cellDim == 1 && numCorners == 3) {
2231: /* Quadratic line */
2232: faceSize = 1;
2233: posOrient = PETSC_TRUE;
2234: } else if (cellDim == 2 && numCorners == 4) {
2235: /* Quads */
2236: faceSize = 2;
2237: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2238: posOrient = PETSC_TRUE;
2239: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2240: posOrient = PETSC_TRUE;
2241: } else {
2242: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2243: posOrient = PETSC_FALSE;
2244: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2245: }
2246: } else if (cellDim == 2 && numCorners == 6) {
2247: /* Quadratic triangle (I hate this) */
2248: /* Edges are determined by the first 2 vertices (corners of edges) */
2249: const PetscInt faceSizeTri = 3;
2250: PetscInt sortedIndices[3], i, iFace;
2251: PetscBool found = PETSC_FALSE;
2252: PetscInt faceVerticesTriSorted[9] = {
2253: 0, 3, 4, /* bottom */
2254: 1, 4, 5, /* right */
2255: 2, 3, 5, /* left */
2256: };
2257: PetscInt faceVerticesTri[9] = {
2258: 0, 3, 4, /* bottom */
2259: 1, 4, 5, /* right */
2260: 2, 5, 3, /* left */
2261: };
2263: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2264: PetscSortInt(faceSizeTri, sortedIndices);
2265: for (iFace = 0; iFace < 3; ++iFace) {
2266: const PetscInt ii = iFace*faceSizeTri;
2267: PetscInt fVertex, cVertex;
2269: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2270: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2271: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2272: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2273: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2274: faceVertices[fVertex] = origVertices[cVertex];
2275: break;
2276: }
2277: }
2278: }
2279: found = PETSC_TRUE;
2280: break;
2281: }
2282: }
2283: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2284: if (posOriented) *posOriented = PETSC_TRUE;
2285: return(0);
2286: } else if (cellDim == 2 && numCorners == 9) {
2287: /* Quadratic quad (I hate this) */
2288: /* Edges are determined by the first 2 vertices (corners of edges) */
2289: const PetscInt faceSizeQuad = 3;
2290: PetscInt sortedIndices[3], i, iFace;
2291: PetscBool found = PETSC_FALSE;
2292: PetscInt faceVerticesQuadSorted[12] = {
2293: 0, 1, 4, /* bottom */
2294: 1, 2, 5, /* right */
2295: 2, 3, 6, /* top */
2296: 0, 3, 7, /* left */
2297: };
2298: PetscInt faceVerticesQuad[12] = {
2299: 0, 1, 4, /* bottom */
2300: 1, 2, 5, /* right */
2301: 2, 3, 6, /* top */
2302: 3, 0, 7, /* left */
2303: };
2305: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2306: PetscSortInt(faceSizeQuad, sortedIndices);
2307: for (iFace = 0; iFace < 4; ++iFace) {
2308: const PetscInt ii = iFace*faceSizeQuad;
2309: PetscInt fVertex, cVertex;
2311: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2312: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2313: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2314: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2315: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2316: faceVertices[fVertex] = origVertices[cVertex];
2317: break;
2318: }
2319: }
2320: }
2321: found = PETSC_TRUE;
2322: break;
2323: }
2324: }
2325: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2326: if (posOriented) *posOriented = PETSC_TRUE;
2327: return(0);
2328: } else if (cellDim == 3 && numCorners == 8) {
2329: /* Hexes
2330: A hex is two oriented quads with the normal of the first
2331: pointing up at the second.
2333: 7---6
2334: /| /|
2335: 4---5 |
2336: | 1-|-2
2337: |/ |/
2338: 0---3
2340: Faces are determined by the first 4 vertices (corners of faces) */
2341: const PetscInt faceSizeHex = 4;
2342: PetscInt sortedIndices[4], i, iFace;
2343: PetscBool found = PETSC_FALSE;
2344: PetscInt faceVerticesHexSorted[24] = {
2345: 0, 1, 2, 3, /* bottom */
2346: 4, 5, 6, 7, /* top */
2347: 0, 3, 4, 5, /* front */
2348: 2, 3, 5, 6, /* right */
2349: 1, 2, 6, 7, /* back */
2350: 0, 1, 4, 7, /* left */
2351: };
2352: PetscInt faceVerticesHex[24] = {
2353: 1, 2, 3, 0, /* bottom */
2354: 4, 5, 6, 7, /* top */
2355: 0, 3, 5, 4, /* front */
2356: 3, 2, 6, 5, /* right */
2357: 2, 1, 7, 6, /* back */
2358: 1, 0, 4, 7, /* left */
2359: };
2361: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2362: PetscSortInt(faceSizeHex, sortedIndices);
2363: for (iFace = 0; iFace < 6; ++iFace) {
2364: const PetscInt ii = iFace*faceSizeHex;
2365: PetscInt fVertex, cVertex;
2367: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2368: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2369: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2370: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2371: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2372: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2373: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2374: faceVertices[fVertex] = origVertices[cVertex];
2375: break;
2376: }
2377: }
2378: }
2379: found = PETSC_TRUE;
2380: break;
2381: }
2382: }
2383: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2384: if (posOriented) *posOriented = PETSC_TRUE;
2385: return(0);
2386: } else if (cellDim == 3 && numCorners == 10) {
2387: /* Quadratic tet */
2388: /* Faces are determined by the first 3 vertices (corners of faces) */
2389: const PetscInt faceSizeTet = 6;
2390: PetscInt sortedIndices[6], i, iFace;
2391: PetscBool found = PETSC_FALSE;
2392: PetscInt faceVerticesTetSorted[24] = {
2393: 0, 1, 2, 6, 7, 8, /* bottom */
2394: 0, 3, 4, 6, 7, 9, /* front */
2395: 1, 4, 5, 7, 8, 9, /* right */
2396: 2, 3, 5, 6, 8, 9, /* left */
2397: };
2398: PetscInt faceVerticesTet[24] = {
2399: 0, 1, 2, 6, 7, 8, /* bottom */
2400: 0, 4, 3, 6, 7, 9, /* front */
2401: 1, 5, 4, 7, 8, 9, /* right */
2402: 2, 3, 5, 8, 6, 9, /* left */
2403: };
2405: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2406: PetscSortInt(faceSizeTet, sortedIndices);
2407: for (iFace=0; iFace < 4; ++iFace) {
2408: const PetscInt ii = iFace*faceSizeTet;
2409: PetscInt fVertex, cVertex;
2411: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2412: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2413: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2414: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2415: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2416: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2417: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2418: faceVertices[fVertex] = origVertices[cVertex];
2419: break;
2420: }
2421: }
2422: }
2423: found = PETSC_TRUE;
2424: break;
2425: }
2426: }
2427: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2428: if (posOriented) *posOriented = PETSC_TRUE;
2429: return(0);
2430: } else if (cellDim == 3 && numCorners == 27) {
2431: /* Quadratic hexes (I hate this)
2432: A hex is two oriented quads with the normal of the first
2433: pointing up at the second.
2435: 7---6
2436: /| /|
2437: 4---5 |
2438: | 3-|-2
2439: |/ |/
2440: 0---1
2442: Faces are determined by the first 4 vertices (corners of faces) */
2443: const PetscInt faceSizeQuadHex = 9;
2444: PetscInt sortedIndices[9], i, iFace;
2445: PetscBool found = PETSC_FALSE;
2446: PetscInt faceVerticesQuadHexSorted[54] = {
2447: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2448: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2449: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2450: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2451: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2452: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2453: };
2454: PetscInt faceVerticesQuadHex[54] = {
2455: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2456: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2457: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2458: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2459: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2460: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2461: };
2463: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2464: PetscSortInt(faceSizeQuadHex, sortedIndices);
2465: for (iFace = 0; iFace < 6; ++iFace) {
2466: const PetscInt ii = iFace*faceSizeQuadHex;
2467: PetscInt fVertex, cVertex;
2469: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2470: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2471: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2472: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2473: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2474: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2475: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2476: faceVertices[fVertex] = origVertices[cVertex];
2477: break;
2478: }
2479: }
2480: }
2481: found = PETSC_TRUE;
2482: break;
2483: }
2484: }
2485: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2486: if (posOriented) *posOriented = PETSC_TRUE;
2487: return(0);
2488: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2489: if (!posOrient) {
2490: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2491: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2492: } else {
2493: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2494: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2495: }
2496: if (posOriented) *posOriented = posOrient;
2497: return(0);
2498: }
2500: /*@
2501: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2502: in faceVertices. The orientation is such that the face normal points out of the cell
2504: Not collective
2506: Input Parameters:
2507: + dm - The original mesh
2508: . cell - The cell mesh point
2509: . faceSize - The number of vertices on the face
2510: . face - The face vertices
2511: . numCorners - The number of vertices on the cell
2512: . indices - Local numbering of face vertices in cell cone
2513: - origVertices - Original face vertices
2515: Output Parameter:
2516: + faceVertices - The face vertices properly oriented
2517: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2519: Level: developer
2521: .seealso: DMPlexGetCone()
2522: @*/
2523: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2524: {
2525: const PetscInt *cone = NULL;
2526: PetscInt coneSize, v, f, v2;
2527: PetscInt oppositeVertex = -1;
2528: PetscErrorCode ierr;
2531: DMPlexGetConeSize(dm, cell, &coneSize);
2532: DMPlexGetCone(dm, cell, &cone);
2533: for (v = 0, v2 = 0; v < coneSize; ++v) {
2534: PetscBool found = PETSC_FALSE;
2536: for (f = 0; f < faceSize; ++f) {
2537: if (face[f] == cone[v]) {
2538: found = PETSC_TRUE; break;
2539: }
2540: }
2541: if (found) {
2542: indices[v2] = v;
2543: origVertices[v2] = cone[v];
2544: ++v2;
2545: } else {
2546: oppositeVertex = v;
2547: }
2548: }
2549: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2550: return(0);
2551: }
2553: /*
2554: DMPlexInsertFace_Internal - Puts a face into the mesh
2556: Not collective
2558: Input Parameters:
2559: + dm - The DMPlex
2560: . numFaceVertex - The number of vertices in the face
2561: . faceVertices - The vertices in the face for dm
2562: . subfaceVertices - The vertices in the face for subdm
2563: . numCorners - The number of vertices in the cell
2564: . cell - A cell in dm containing the face
2565: . subcell - A cell in subdm containing the face
2566: . firstFace - First face in the mesh
2567: - newFacePoint - Next face in the mesh
2569: Output Parameters:
2570: . newFacePoint - Contains next face point number on input, updated on output
2572: Level: developer
2573: */
2574: 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)
2575: {
2576: MPI_Comm comm;
2577: DM_Plex *submesh = (DM_Plex*) subdm->data;
2578: const PetscInt *faces;
2579: PetscInt numFaces, coneSize;
2580: PetscErrorCode ierr;
2583: PetscObjectGetComm((PetscObject)dm,&comm);
2584: DMPlexGetConeSize(subdm, subcell, &coneSize);
2585: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2586: #if 0
2587: /* Cannot use this because support() has not been constructed yet */
2588: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2589: #else
2590: {
2591: PetscInt f;
2593: numFaces = 0;
2594: DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2595: for (f = firstFace; f < *newFacePoint; ++f) {
2596: PetscInt dof, off, d;
2598: PetscSectionGetDof(submesh->coneSection, f, &dof);
2599: PetscSectionGetOffset(submesh->coneSection, f, &off);
2600: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2601: for (d = 0; d < dof; ++d) {
2602: const PetscInt p = submesh->cones[off+d];
2603: PetscInt v;
2605: for (v = 0; v < numFaceVertices; ++v) {
2606: if (subfaceVertices[v] == p) break;
2607: }
2608: if (v == numFaceVertices) break;
2609: }
2610: if (d == dof) {
2611: numFaces = 1;
2612: ((PetscInt*) faces)[0] = f;
2613: }
2614: }
2615: }
2616: #endif
2617: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2618: else if (numFaces == 1) {
2619: /* Add the other cell neighbor for this face */
2620: DMPlexSetCone(subdm, subcell, faces);
2621: } else {
2622: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2623: PetscBool posOriented;
2625: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2626: origVertices = &orientedVertices[numFaceVertices];
2627: indices = &orientedVertices[numFaceVertices*2];
2628: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2629: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2630: /* TODO: I know that routine should return a permutation, not the indices */
2631: for (v = 0; v < numFaceVertices; ++v) {
2632: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2633: for (ov = 0; ov < numFaceVertices; ++ov) {
2634: if (orientedVertices[ov] == vertex) {
2635: orientedSubVertices[ov] = subvertex;
2636: break;
2637: }
2638: }
2639: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2640: }
2641: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2642: DMPlexSetCone(subdm, subcell, newFacePoint);
2643: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2644: ++(*newFacePoint);
2645: }
2646: #if 0
2647: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2648: #else
2649: DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2650: #endif
2651: return(0);
2652: }
2654: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2655: {
2656: MPI_Comm comm;
2657: DMLabel subpointMap;
2658: IS subvertexIS, subcellIS;
2659: const PetscInt *subVertices, *subCells;
2660: PetscInt numSubVertices, firstSubVertex, numSubCells;
2661: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2662: PetscInt vStart, vEnd, c, f;
2663: PetscErrorCode ierr;
2666: PetscObjectGetComm((PetscObject)dm,&comm);
2667: /* Create subpointMap which marks the submesh */
2668: DMLabelCreate("subpoint_map", &subpointMap);
2669: DMPlexSetSubpointMap(subdm, subpointMap);
2670: DMLabelDestroy(&subpointMap);
2671: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2672: /* Setup chart */
2673: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2674: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2675: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2676: DMPlexSetVTKCellHeight(subdm, 1);
2677: /* Set cone sizes */
2678: firstSubVertex = numSubCells;
2679: firstSubFace = numSubCells+numSubVertices;
2680: newFacePoint = firstSubFace;
2681: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2682: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2683: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2684: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2685: for (c = 0; c < numSubCells; ++c) {
2686: DMPlexSetConeSize(subdm, c, 1);
2687: }
2688: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2689: DMPlexSetConeSize(subdm, f, nFV);
2690: }
2691: DMSetUp(subdm);
2692: /* Create face cones */
2693: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2694: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2695: DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2696: for (c = 0; c < numSubCells; ++c) {
2697: const PetscInt cell = subCells[c];
2698: const PetscInt subcell = c;
2699: PetscInt *closure = NULL;
2700: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2702: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2703: for (cl = 0; cl < closureSize*2; cl += 2) {
2704: const PetscInt point = closure[cl];
2705: PetscInt subVertex;
2707: if ((point >= vStart) && (point < vEnd)) {
2708: ++numCorners;
2709: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2710: if (subVertex >= 0) {
2711: closure[faceSize] = point;
2712: subface[faceSize] = firstSubVertex+subVertex;
2713: ++faceSize;
2714: }
2715: }
2716: }
2717: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2718: if (faceSize == nFV) {
2719: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2720: }
2721: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2722: }
2723: DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2724: DMPlexSymmetrize(subdm);
2725: DMPlexStratify(subdm);
2726: /* Build coordinates */
2727: {
2728: PetscSection coordSection, subCoordSection;
2729: Vec coordinates, subCoordinates;
2730: PetscScalar *coords, *subCoords;
2731: PetscInt numComp, coordSize, v;
2732: const char *name;
2734: DMGetCoordinateSection(dm, &coordSection);
2735: DMGetCoordinatesLocal(dm, &coordinates);
2736: DMGetCoordinateSection(subdm, &subCoordSection);
2737: PetscSectionSetNumFields(subCoordSection, 1);
2738: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2739: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2740: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2741: for (v = 0; v < numSubVertices; ++v) {
2742: const PetscInt vertex = subVertices[v];
2743: const PetscInt subvertex = firstSubVertex+v;
2744: PetscInt dof;
2746: PetscSectionGetDof(coordSection, vertex, &dof);
2747: PetscSectionSetDof(subCoordSection, subvertex, dof);
2748: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2749: }
2750: PetscSectionSetUp(subCoordSection);
2751: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2752: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2753: PetscObjectGetName((PetscObject)coordinates,&name);
2754: PetscObjectSetName((PetscObject)subCoordinates,name);
2755: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2756: VecSetType(subCoordinates,VECSTANDARD);
2757: if (coordSize) {
2758: VecGetArray(coordinates, &coords);
2759: VecGetArray(subCoordinates, &subCoords);
2760: for (v = 0; v < numSubVertices; ++v) {
2761: const PetscInt vertex = subVertices[v];
2762: const PetscInt subvertex = firstSubVertex+v;
2763: PetscInt dof, off, sdof, soff, d;
2765: PetscSectionGetDof(coordSection, vertex, &dof);
2766: PetscSectionGetOffset(coordSection, vertex, &off);
2767: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2768: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2769: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2770: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2771: }
2772: VecRestoreArray(coordinates, &coords);
2773: VecRestoreArray(subCoordinates, &subCoords);
2774: }
2775: DMSetCoordinatesLocal(subdm, subCoordinates);
2776: VecDestroy(&subCoordinates);
2777: }
2778: /* Cleanup */
2779: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2780: ISDestroy(&subvertexIS);
2781: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2782: ISDestroy(&subcellIS);
2783: return(0);
2784: }
2786: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2787: {
2788: PetscInt subPoint;
2791: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2792: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2793: }
2795: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2796: {
2797: MPI_Comm comm;
2798: DMLabel subpointMap;
2799: IS *subpointIS;
2800: const PetscInt **subpoints;
2801: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2802: PetscInt totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2803: PetscMPIInt rank;
2804: PetscErrorCode ierr;
2807: PetscObjectGetComm((PetscObject)dm,&comm);
2808: MPI_Comm_rank(comm, &rank);
2809: /* Create subpointMap which marks the submesh */
2810: DMLabelCreate("subpoint_map", &subpointMap);
2811: DMPlexSetSubpointMap(subdm, subpointMap);
2812: if (cellHeight) {
2813: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2814: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2815: } else {
2816: DMLabel depth;
2817: IS pointIS;
2818: const PetscInt *points;
2819: PetscInt numPoints;
2821: DMPlexGetDepthLabel(dm, &depth);
2822: DMLabelGetStratumSize(label, value, &numPoints);
2823: DMLabelGetStratumIS(label, value, &pointIS);
2824: ISGetIndices(pointIS, &points);
2825: for (p = 0; p < numPoints; ++p) {
2826: PetscInt *closure = NULL;
2827: PetscInt closureSize, c, pdim;
2829: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2830: for (c = 0; c < closureSize*2; c += 2) {
2831: DMLabelGetValue(depth, closure[c], &pdim);
2832: DMLabelSetValue(subpointMap, closure[c], pdim);
2833: }
2834: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2835: }
2836: ISRestoreIndices(pointIS, &points);
2837: ISDestroy(&pointIS);
2838: }
2839: DMLabelDestroy(&subpointMap);
2840: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2841: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2842: cMax = (cMax < 0) ? cEnd : cMax;
2843: /* Setup chart */
2844: DMGetDimension(dm, &dim);
2845: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2846: for (d = 0; d <= dim; ++d) {
2847: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2848: totSubPoints += numSubPoints[d];
2849: }
2850: DMPlexSetChart(subdm, 0, totSubPoints);
2851: DMPlexSetVTKCellHeight(subdm, cellHeight);
2852: /* Set cone sizes */
2853: firstSubPoint[dim] = 0;
2854: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
2855: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
2856: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2857: for (d = 0; d <= dim; ++d) {
2858: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2859: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2860: }
2861: for (d = 0; d <= dim; ++d) {
2862: for (p = 0; p < numSubPoints[d]; ++p) {
2863: const PetscInt point = subpoints[d][p];
2864: const PetscInt subpoint = firstSubPoint[d] + p;
2865: const PetscInt *cone;
2866: PetscInt coneSize, coneSizeNew, c, val;
2868: DMPlexGetConeSize(dm, point, &coneSize);
2869: DMPlexSetConeSize(subdm, subpoint, coneSize);
2870: if (cellHeight && (d == dim)) {
2871: DMPlexGetCone(dm, point, &cone);
2872: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2873: DMLabelGetValue(subpointMap, cone[c], &val);
2874: if (val >= 0) coneSizeNew++;
2875: }
2876: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2877: }
2878: }
2879: }
2880: DMSetUp(subdm);
2881: /* Set cones */
2882: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2883: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2884: for (d = 0; d <= dim; ++d) {
2885: for (p = 0; p < numSubPoints[d]; ++p) {
2886: const PetscInt point = subpoints[d][p];
2887: const PetscInt subpoint = firstSubPoint[d] + p;
2888: const PetscInt *cone, *ornt;
2889: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
2891: if (d == dim-1) {
2892: const PetscInt *support, *cone, *ornt;
2893: PetscInt supportSize, coneSize, s, subc;
2895: DMPlexGetSupport(dm, point, &support);
2896: DMPlexGetSupportSize(dm, point, &supportSize);
2897: for (s = 0; s < supportSize; ++s) {
2898: if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
2899: PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
2900: if (subc >= 0) {
2901: const PetscInt ccell = subpoints[d+1][subc];
2903: DMPlexGetCone(dm, ccell, &cone);
2904: DMPlexGetConeSize(dm, ccell, &coneSize);
2905: DMPlexGetConeOrientation(dm, ccell, &ornt);
2906: for (c = 0; c < coneSize; ++c) {
2907: if (cone[c] == point) {
2908: fornt = ornt[c];
2909: break;
2910: }
2911: }
2912: break;
2913: }
2914: }
2915: }
2916: DMPlexGetConeSize(dm, point, &coneSize);
2917: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2918: DMPlexGetCone(dm, point, &cone);
2919: DMPlexGetConeOrientation(dm, point, &ornt);
2920: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2921: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2922: if (subc >= 0) {
2923: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
2924: orntNew[coneSizeNew] = ornt[c];
2925: ++coneSizeNew;
2926: }
2927: }
2928: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2929: if (fornt < 0) {
2930: /* This should be replaced by a call to DMPlexReverseCell() */
2931: #if 0
2932: DMPlexReverseCell(subdm, subpoint);
2933: #else
2934: for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
2935: PetscInt faceSize, tmp;
2937: tmp = coneNew[c];
2938: coneNew[c] = coneNew[coneSizeNew-1-c];
2939: coneNew[coneSizeNew-1-c] = tmp;
2940: DMPlexGetConeSize(dm, cone[c], &faceSize);
2941: tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
2942: orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
2943: orntNew[coneSizeNew-1-c] = tmp;
2944: }
2945: }
2946: DMPlexSetCone(subdm, subpoint, coneNew);
2947: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2948: #endif
2949: }
2950: }
2951: PetscFree2(coneNew,orntNew);
2952: DMPlexSymmetrize(subdm);
2953: DMPlexStratify(subdm);
2954: /* Build coordinates */
2955: {
2956: PetscSection coordSection, subCoordSection;
2957: Vec coordinates, subCoordinates;
2958: PetscScalar *coords, *subCoords;
2959: PetscInt cdim, numComp, coordSize;
2960: const char *name;
2962: DMGetCoordinateDim(dm, &cdim);
2963: DMGetCoordinateSection(dm, &coordSection);
2964: DMGetCoordinatesLocal(dm, &coordinates);
2965: DMGetCoordinateSection(subdm, &subCoordSection);
2966: PetscSectionSetNumFields(subCoordSection, 1);
2967: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2968: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2969: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2970: for (v = 0; v < numSubPoints[0]; ++v) {
2971: const PetscInt vertex = subpoints[0][v];
2972: const PetscInt subvertex = firstSubPoint[0]+v;
2973: PetscInt dof;
2975: PetscSectionGetDof(coordSection, vertex, &dof);
2976: PetscSectionSetDof(subCoordSection, subvertex, dof);
2977: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2978: }
2979: PetscSectionSetUp(subCoordSection);
2980: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2981: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2982: PetscObjectGetName((PetscObject)coordinates,&name);
2983: PetscObjectSetName((PetscObject)subCoordinates,name);
2984: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2985: VecSetBlockSize(subCoordinates, cdim);
2986: VecSetType(subCoordinates,VECSTANDARD);
2987: VecGetArray(coordinates, &coords);
2988: VecGetArray(subCoordinates, &subCoords);
2989: for (v = 0; v < numSubPoints[0]; ++v) {
2990: const PetscInt vertex = subpoints[0][v];
2991: const PetscInt subvertex = firstSubPoint[0]+v;
2992: PetscInt dof, off, sdof, soff, d;
2994: PetscSectionGetDof(coordSection, vertex, &dof);
2995: PetscSectionGetOffset(coordSection, vertex, &off);
2996: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2997: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2998: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2999: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3000: }
3001: VecRestoreArray(coordinates, &coords);
3002: VecRestoreArray(subCoordinates, &subCoords);
3003: DMSetCoordinatesLocal(subdm, subCoordinates);
3004: VecDestroy(&subCoordinates);
3005: }
3006: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3007: {
3008: PetscSF sfPoint, sfPointSub;
3009: IS subpIS;
3010: const PetscSFNode *remotePoints;
3011: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3012: const PetscInt *localPoints, *subpoints;
3013: PetscInt *slocalPoints;
3014: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3015: PetscMPIInt rank;
3017: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3018: DMGetPointSF(dm, &sfPoint);
3019: DMGetPointSF(subdm, &sfPointSub);
3020: DMPlexGetChart(dm, &pStart, &pEnd);
3021: DMPlexGetChart(subdm, NULL, &numSubroots);
3022: DMPlexCreateSubpointIS(subdm, &subpIS);
3023: if (subpIS) {
3024: ISGetIndices(subpIS, &subpoints);
3025: ISGetLocalSize(subpIS, &numSubpoints);
3026: }
3027: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3028: if (numRoots >= 0) {
3029: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3030: for (p = 0; p < pEnd-pStart; ++p) {
3031: newLocalPoints[p].rank = -2;
3032: newLocalPoints[p].index = -2;
3033: }
3034: /* Set subleaves */
3035: for (l = 0; l < numLeaves; ++l) {
3036: const PetscInt point = localPoints[l];
3037: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3039: if (subpoint < 0) continue;
3040: newLocalPoints[point-pStart].rank = rank;
3041: newLocalPoints[point-pStart].index = subpoint;
3042: ++numSubleaves;
3043: }
3044: /* Must put in owned subpoints */
3045: for (p = pStart; p < pEnd; ++p) {
3046: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3048: if (subpoint < 0) {
3049: newOwners[p-pStart].rank = -3;
3050: newOwners[p-pStart].index = -3;
3051: } else {
3052: newOwners[p-pStart].rank = rank;
3053: newOwners[p-pStart].index = subpoint;
3054: }
3055: }
3056: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3057: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3058: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3059: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3060: PetscMalloc1(numSubleaves, &slocalPoints);
3061: PetscMalloc1(numSubleaves, &sremotePoints);
3062: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3063: const PetscInt point = localPoints[l];
3064: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3066: if (subpoint < 0) continue;
3067: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3068: slocalPoints[sl] = subpoint;
3069: sremotePoints[sl].rank = newLocalPoints[point].rank;
3070: sremotePoints[sl].index = newLocalPoints[point].index;
3071: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3072: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3073: ++sl;
3074: }
3075: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3076: PetscFree2(newLocalPoints,newOwners);
3077: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3078: }
3079: if (subpIS) {
3080: ISRestoreIndices(subpIS, &subpoints);
3081: ISDestroy(&subpIS);
3082: }
3083: }
3084: /* Cleanup */
3085: for (d = 0; d <= dim; ++d) {
3086: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3087: ISDestroy(&subpointIS[d]);
3088: }
3089: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3090: return(0);
3091: }
3093: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
3094: {
3098: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);
3099: return(0);
3100: }
3102: /*@
3103: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3105: Input Parameters:
3106: + dm - The original mesh
3107: . vertexLabel - The DMLabel marking vertices contained in the surface
3108: - value - The label value to use
3110: Output Parameter:
3111: . subdm - The surface mesh
3113: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3115: Level: developer
3117: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3118: @*/
3119: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
3120: {
3121: PetscInt dim, cdim, depth;
3127: DMGetDimension(dm, &dim);
3128: DMPlexGetDepth(dm, &depth);
3129: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3130: DMSetType(*subdm, DMPLEX);
3131: DMSetDimension(*subdm, dim-1);
3132: DMGetCoordinateDim(dm, &cdim);
3133: DMSetCoordinateDim(*subdm, cdim);
3134: if (depth == dim) {
3135: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
3136: } else {
3137: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3138: }
3139: return(0);
3140: }
3142: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3143: {
3144: MPI_Comm comm;
3145: DMLabel subpointMap;
3146: IS subvertexIS;
3147: const PetscInt *subVertices;
3148: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3149: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3150: PetscInt cMax, c, f;
3151: PetscErrorCode ierr;
3154: PetscObjectGetComm((PetscObject)dm, &comm);
3155: /* Create subpointMap which marks the submesh */
3156: DMLabelCreate("subpoint_map", &subpointMap);
3157: DMPlexSetSubpointMap(subdm, subpointMap);
3158: DMLabelDestroy(&subpointMap);
3159: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3160: /* Setup chart */
3161: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3162: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3163: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3164: DMPlexSetVTKCellHeight(subdm, 1);
3165: /* Set cone sizes */
3166: firstSubVertex = numSubCells;
3167: firstSubFace = numSubCells+numSubVertices;
3168: newFacePoint = firstSubFace;
3169: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3170: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3171: for (c = 0; c < numSubCells; ++c) {
3172: DMPlexSetConeSize(subdm, c, 1);
3173: }
3174: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3175: DMPlexSetConeSize(subdm, f, nFV);
3176: }
3177: DMSetUp(subdm);
3178: /* Create face cones */
3179: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3180: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3181: DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
3182: for (c = 0; c < numSubCells; ++c) {
3183: const PetscInt cell = subCells[c];
3184: const PetscInt subcell = c;
3185: const PetscInt *cone, *cells;
3186: PetscInt numCells, subVertex, p, v;
3188: if (cell < cMax) continue;
3189: DMPlexGetCone(dm, cell, &cone);
3190: for (v = 0; v < nFV; ++v) {
3191: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3192: subface[v] = firstSubVertex+subVertex;
3193: }
3194: DMPlexSetCone(subdm, newFacePoint, subface);
3195: DMPlexSetCone(subdm, subcell, &newFacePoint);
3196: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3197: /* Not true in parallel
3198: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3199: for (p = 0; p < numCells; ++p) {
3200: PetscInt negsubcell;
3202: if (cells[p] >= cMax) continue;
3203: /* I know this is a crap search */
3204: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3205: if (subCells[negsubcell] == cells[p]) break;
3206: }
3207: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3208: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3209: }
3210: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3211: ++newFacePoint;
3212: }
3213: DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
3214: DMPlexSymmetrize(subdm);
3215: DMPlexStratify(subdm);
3216: /* Build coordinates */
3217: {
3218: PetscSection coordSection, subCoordSection;
3219: Vec coordinates, subCoordinates;
3220: PetscScalar *coords, *subCoords;
3221: PetscInt cdim, numComp, coordSize, v;
3222: const char *name;
3224: DMGetCoordinateDim(dm, &cdim);
3225: DMGetCoordinateSection(dm, &coordSection);
3226: DMGetCoordinatesLocal(dm, &coordinates);
3227: DMGetCoordinateSection(subdm, &subCoordSection);
3228: PetscSectionSetNumFields(subCoordSection, 1);
3229: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3230: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3231: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3232: for (v = 0; v < numSubVertices; ++v) {
3233: const PetscInt vertex = subVertices[v];
3234: const PetscInt subvertex = firstSubVertex+v;
3235: PetscInt dof;
3237: PetscSectionGetDof(coordSection, vertex, &dof);
3238: PetscSectionSetDof(subCoordSection, subvertex, dof);
3239: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3240: }
3241: PetscSectionSetUp(subCoordSection);
3242: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3243: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3244: PetscObjectGetName((PetscObject)coordinates,&name);
3245: PetscObjectSetName((PetscObject)subCoordinates,name);
3246: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3247: VecSetBlockSize(subCoordinates, cdim);
3248: VecSetType(subCoordinates,VECSTANDARD);
3249: VecGetArray(coordinates, &coords);
3250: VecGetArray(subCoordinates, &subCoords);
3251: for (v = 0; v < numSubVertices; ++v) {
3252: const PetscInt vertex = subVertices[v];
3253: const PetscInt subvertex = firstSubVertex+v;
3254: PetscInt dof, off, sdof, soff, d;
3256: PetscSectionGetDof(coordSection, vertex, &dof);
3257: PetscSectionGetOffset(coordSection, vertex, &off);
3258: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3259: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3260: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3261: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3262: }
3263: VecRestoreArray(coordinates, &coords);
3264: VecRestoreArray(subCoordinates, &subCoords);
3265: DMSetCoordinatesLocal(subdm, subCoordinates);
3266: VecDestroy(&subCoordinates);
3267: }
3268: /* Build SF */
3269: CHKMEMQ;
3270: {
3271: PetscSF sfPoint, sfPointSub;
3272: const PetscSFNode *remotePoints;
3273: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3274: const PetscInt *localPoints;
3275: PetscInt *slocalPoints;
3276: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3277: PetscMPIInt rank;
3279: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3280: DMGetPointSF(dm, &sfPoint);
3281: DMGetPointSF(subdm, &sfPointSub);
3282: DMPlexGetChart(dm, &pStart, &pEnd);
3283: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3284: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3285: if (numRoots >= 0) {
3286: /* Only vertices should be shared */
3287: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3288: for (p = 0; p < pEnd-pStart; ++p) {
3289: newLocalPoints[p].rank = -2;
3290: newLocalPoints[p].index = -2;
3291: }
3292: /* Set subleaves */
3293: for (l = 0; l < numLeaves; ++l) {
3294: const PetscInt point = localPoints[l];
3295: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3297: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3298: if (subPoint < 0) continue;
3299: newLocalPoints[point-pStart].rank = rank;
3300: newLocalPoints[point-pStart].index = subPoint;
3301: ++numSubLeaves;
3302: }
3303: /* Must put in owned subpoints */
3304: for (p = pStart; p < pEnd; ++p) {
3305: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3307: if (subPoint < 0) {
3308: newOwners[p-pStart].rank = -3;
3309: newOwners[p-pStart].index = -3;
3310: } else {
3311: newOwners[p-pStart].rank = rank;
3312: newOwners[p-pStart].index = subPoint;
3313: }
3314: }
3315: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3316: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3317: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3318: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3319: PetscMalloc1(numSubLeaves, &slocalPoints);
3320: PetscMalloc1(numSubLeaves, &sremotePoints);
3321: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3322: const PetscInt point = localPoints[l];
3323: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3325: if (subPoint < 0) continue;
3326: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3327: slocalPoints[sl] = subPoint;
3328: sremotePoints[sl].rank = newLocalPoints[point].rank;
3329: sremotePoints[sl].index = newLocalPoints[point].index;
3330: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3331: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3332: ++sl;
3333: }
3334: PetscFree2(newLocalPoints,newOwners);
3335: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3336: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3337: }
3338: }
3339: CHKMEMQ;
3340: /* Cleanup */
3341: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3342: ISDestroy(&subvertexIS);
3343: PetscFree(subCells);
3344: return(0);
3345: }
3347: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3348: {
3349: DMLabel label = NULL;
3353: if (labelname) {DMGetLabel(dm, labelname, &label);}
3354: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);
3355: return(0);
3356: }
3358: /*@C
3359: 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.
3361: Input Parameters:
3362: + dm - The original mesh
3363: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3364: . label - A label name, or NULL
3365: - value - A label value
3367: Output Parameter:
3368: . subdm - The surface mesh
3370: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3372: Level: developer
3374: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3375: @*/
3376: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3377: {
3378: PetscInt dim, cdim, depth;
3384: DMGetDimension(dm, &dim);
3385: DMPlexGetDepth(dm, &depth);
3386: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3387: DMSetType(*subdm, DMPLEX);
3388: DMSetDimension(*subdm, dim-1);
3389: DMGetCoordinateDim(dm, &cdim);
3390: DMSetCoordinateDim(*subdm, cdim);
3391: if (depth == dim) {
3392: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3393: } else {
3394: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3395: }
3396: return(0);
3397: }
3399: /*@
3400: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3402: Input Parameters:
3403: + dm - The original mesh
3404: . cellLabel - The DMLabel marking cells contained in the new mesh
3405: - value - The label value to use
3407: Output Parameter:
3408: . subdm - The new mesh
3410: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3412: Level: developer
3414: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3415: @*/
3416: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3417: {
3418: PetscInt dim;
3424: DMGetDimension(dm, &dim);
3425: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3426: DMSetType(*subdm, DMPLEX);
3427: DMSetDimension(*subdm, dim);
3428: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3429: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);
3430: return(0);
3431: }
3433: /*@
3434: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3436: Input Parameter:
3437: . dm - The submesh DM
3439: Output Parameter:
3440: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3442: Level: developer
3444: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3445: @*/
3446: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3447: {
3451: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3452: return(0);
3453: }
3455: /*@
3456: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3458: Input Parameters:
3459: + dm - The submesh DM
3460: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3462: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3464: Level: developer
3466: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3467: @*/
3468: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3469: {
3470: DM_Plex *mesh = (DM_Plex *) dm->data;
3471: DMLabel tmp;
3476: tmp = mesh->subpointMap;
3477: mesh->subpointMap = subpointMap;
3478: ++mesh->subpointMap->refct;
3479: DMLabelDestroy(&tmp);
3480: return(0);
3481: }
3483: /*@
3484: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
3486: Input Parameter:
3487: . dm - The submesh DM
3489: Output Parameter:
3490: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3492: Note: This IS is guaranteed to be sorted by the construction of the submesh
3494: Level: developer
3496: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3497: @*/
3498: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3499: {
3500: MPI_Comm comm;
3501: DMLabel subpointMap;
3502: IS is;
3503: const PetscInt *opoints;
3504: PetscInt *points, *depths;
3505: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3506: PetscErrorCode ierr;
3511: PetscObjectGetComm((PetscObject)dm,&comm);
3512: *subpointIS = NULL;
3513: DMPlexGetSubpointMap(dm, &subpointMap);
3514: DMPlexGetDepth(dm, &depth);
3515: if (subpointMap && depth >= 0) {
3516: DMPlexGetChart(dm, &pStart, &pEnd);
3517: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3518: DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
3519: depths[0] = depth;
3520: depths[1] = 0;
3521: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3522: PetscMalloc1(pEnd, &points);
3523: for(d = 0, off = 0; d <= depth; ++d) {
3524: const PetscInt dep = depths[d];
3526: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3527: DMLabelGetStratumSize(subpointMap, dep, &n);
3528: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3529: 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);
3530: } else {
3531: if (!n) {
3532: if (d == 0) {
3533: /* Missing cells */
3534: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3535: } else {
3536: /* Missing faces */
3537: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3538: }
3539: }
3540: }
3541: if (n) {
3542: DMLabelGetStratumIS(subpointMap, dep, &is);
3543: ISGetIndices(is, &opoints);
3544: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3545: ISRestoreIndices(is, &opoints);
3546: ISDestroy(&is);
3547: }
3548: }
3549: DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
3550: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3551: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3552: }
3553: return(0);
3554: }