Actual source code: plexsubmesh.c
petsc-3.6.4 2016-04-12
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petscsf.h>
6: PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt cellHeight, DMLabel label)
7: {
8: PetscInt fStart, fEnd, f;
13: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
14: for (f = fStart; f < fEnd; ++f) {
15: PetscInt supportSize;
17: DMPlexGetSupportSize(dm, f, &supportSize);
18: if (supportSize == 1) {DMLabelSetValue(label, f, 1);}
19: }
20: return(0);
21: }
25: /*@
26: DMPlexMarkBoundaryFaces - Mark all faces on the boundary
28: Not Collective
30: Input Parameter:
31: . dm - The original DM
33: Output Parameter:
34: . label - The DMLabel marking boundary faces with value 1
36: Level: developer
38: .seealso: DMLabelCreate(), DMPlexCreateLabel()
39: @*/
40: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
41: {
45: DMPlexMarkBoundaryFaces_Internal(dm, 0, label);
46: return(0);
47: }
51: /*@
52: DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface
54: Input Parameters:
55: + dm - The DM
56: - label - A DMLabel marking the surface points
58: Output Parameter:
59: . label - A DMLabel marking all surface points in the transitive closure
61: Level: developer
63: .seealso: DMPlexLabelCohesiveComplete()
64: @*/
65: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
66: {
67: IS valueIS;
68: const PetscInt *values;
69: PetscInt numValues, v;
70: PetscErrorCode ierr;
73: DMLabelGetNumValues(label, &numValues);
74: DMLabelGetValueIS(label, &valueIS);
75: ISGetIndices(valueIS, &values);
76: for (v = 0; v < numValues; ++v) {
77: IS pointIS;
78: const PetscInt *points;
79: PetscInt numPoints, p;
81: DMLabelGetStratumSize(label, values[v], &numPoints);
82: DMLabelGetStratumIS(label, values[v], &pointIS);
83: ISGetIndices(pointIS, &points);
84: for (p = 0; p < numPoints; ++p) {
85: PetscInt *closure = NULL;
86: PetscInt closureSize, c;
88: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
89: for (c = 0; c < closureSize*2; c += 2) {
90: DMLabelSetValue(label, closure[c], values[v]);
91: }
92: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
93: }
94: ISRestoreIndices(pointIS, &points);
95: ISDestroy(&pointIS);
96: }
97: ISRestoreIndices(valueIS, &values);
98: ISDestroy(&valueIS);
99: return(0);
100: }
104: /*@
105: DMPlexLabelAddCells - Starting with a label marking faces on a surface, we add a cell for each face
107: Input Parameters:
108: + dm - The DM
109: - label - A DMLabel marking the surface points
111: Output Parameter:
112: . label - A DMLabel incorporating cells
114: Level: developer
116: Note: The cells allow FEM boundary conditions to be applied using the cell geometry
118: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
119: @*/
120: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
121: {
122: IS valueIS;
123: const PetscInt *values;
124: PetscInt numValues, v, cStart, cEnd, cEndInterior;
125: PetscErrorCode ierr;
128: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
129: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
130: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
131: DMLabelGetNumValues(label, &numValues);
132: DMLabelGetValueIS(label, &valueIS);
133: ISGetIndices(valueIS, &values);
134: for (v = 0; v < numValues; ++v) {
135: IS pointIS;
136: const PetscInt *points;
137: PetscInt numPoints, p;
139: DMLabelGetStratumSize(label, values[v], &numPoints);
140: DMLabelGetStratumIS(label, values[v], &pointIS);
141: ISGetIndices(pointIS, &points);
142: for (p = 0; p < numPoints; ++p) {
143: PetscInt *closure = NULL;
144: PetscInt closureSize, point, cl;
146: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
147: for (cl = closureSize-1; cl > 0; --cl) {
148: point = closure[cl*2];
149: if ((point >= cStart) && (point < cEnd)) {DMLabelSetValue(label, point, values[v]); break;}
150: }
151: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
152: }
153: ISRestoreIndices(pointIS, &points);
154: ISDestroy(&pointIS);
155: }
156: ISRestoreIndices(valueIS, &values);
157: ISDestroy(&valueIS);
158: return(0);
159: }
163: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthMax[], PetscInt depthEnd[], PetscInt depthShift[])
164: {
165: if (depth < 0) return p;
166: /* Normal Cells */ if (p < depthMax[depth]) return p;
167: /* Hybrid Cells+Normal Vertices */ if (p < depthMax[0]) return p + depthShift[depth];
168: /* Hybrid Vertices+Normal Faces */ if (depth < 2 || p < depthMax[depth-1]) return p + depthShift[depth] + depthShift[0];
169: /* Hybrid Faces+Normal Edges */ if (depth < 3 || p < depthMax[depth-2]) return p + depthShift[depth] + depthShift[0] + depthShift[depth-1];
170: /* Hybrid Edges */ return p + depthShift[depth] + depthShift[0] + depthShift[depth-1] + depthShift[depth-2];
171: }
175: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
176: {
177: PetscInt *depthMax, *depthEnd;
178: PetscInt depth = 0, d, pStart, pEnd, p;
182: DMPlexGetDepth(dm, &depth);
183: if (depth < 0) return(0);
184: PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
185: /* Step 1: Expand chart */
186: DMPlexGetChart(dm, &pStart, &pEnd);
187: DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);
188: for (d = 0; d <= depth; ++d) {
189: pEnd += depthShift[d];
190: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
191: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
192: }
193: DMPlexSetChart(dmNew, pStart, pEnd);
194: /* Step 2: Set cone and support sizes */
195: for (d = 0; d <= depth; ++d) {
196: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
197: for (p = pStart; p < pEnd; ++p) {
198: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift);
199: PetscInt size;
201: DMPlexGetConeSize(dm, p, &size);
202: DMPlexSetConeSize(dmNew, newp, size);
203: DMPlexGetSupportSize(dm, p, &size);
204: DMPlexSetSupportSize(dmNew, newp, size);
205: }
206: }
207: PetscFree2(depthMax,depthEnd);
208: return(0);
209: }
213: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
214: {
215: PetscInt *depthEnd, *depthMax, *newpoints;
216: PetscInt depth = 0, d, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;
220: DMPlexGetDepth(dm, &depth);
221: if (depth < 0) return(0);
222: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
223: DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
224: PetscMalloc3(depth+1,&depthEnd,depth+1,&depthMax,PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
225: DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);
226: for (d = 0; d <= depth; ++d) {
227: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
228: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
229: }
230: /* Step 5: Set cones and supports */
231: DMPlexGetChart(dm, &pStart, &pEnd);
232: for (p = pStart; p < pEnd; ++p) {
233: const PetscInt *points = NULL, *orientations = NULL;
234: PetscInt size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift);
236: DMPlexGetConeSize(dm, p, &size);
237: DMPlexGetCone(dm, p, &points);
238: DMPlexGetConeOrientation(dm, p, &orientations);
239: for (i = 0; i < size; ++i) {
240: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift);
241: }
242: DMPlexSetCone(dmNew, newp, newpoints);
243: DMPlexSetConeOrientation(dmNew, newp, orientations);
244: DMPlexGetSupportSize(dm, p, &size);
245: DMPlexGetSupportSize(dmNew, newp, &sizeNew);
246: DMPlexGetSupport(dm, p, &points);
247: for (i = 0; i < size; ++i) {
248: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift);
249: }
250: for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
251: DMPlexSetSupport(dmNew, newp, newpoints);
252: }
253: PetscFree3(depthEnd,depthMax,newpoints);
254: return(0);
255: }
259: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
260: {
261: PetscSection coordSection, newCoordSection;
262: Vec coordinates, newCoordinates;
263: PetscScalar *coords, *newCoords;
264: PetscInt *depthEnd, coordSize;
265: PetscInt dim, depth = 0, d, vStart, vEnd, vStartNew, vEndNew, v;
269: DMGetDimension(dm, &dim);
270: DMPlexGetDepth(dm, &depth);
271: PetscMalloc1(depth+1, &depthEnd);
272: for (d = 0; d <= depth; ++d) {
273: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
274: }
275: /* Step 8: Convert coordinates */
276: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
277: DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
278: DMGetCoordinateSection(dm, &coordSection);
279: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
280: PetscSectionSetNumFields(newCoordSection, 1);
281: PetscSectionSetFieldComponents(newCoordSection, 0, dim);
282: PetscSectionSetChart(newCoordSection, vStartNew, vEndNew);
283: for (v = vStartNew; v < vEndNew; ++v) {
284: PetscSectionSetDof(newCoordSection, v, dim);
285: PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
286: }
287: PetscSectionSetUp(newCoordSection);
288: DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
289: PetscSectionGetStorageSize(newCoordSection, &coordSize);
290: VecCreate(PetscObjectComm((PetscObject)dm), &newCoordinates);
291: PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
292: VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
293: VecSetBlockSize(newCoordinates, dim);
294: VecSetType(newCoordinates,VECSTANDARD);
295: DMSetCoordinatesLocal(dmNew, newCoordinates);
296: DMGetCoordinatesLocal(dm, &coordinates);
297: VecGetArray(coordinates, &coords);
298: VecGetArray(newCoordinates, &newCoords);
299: for (v = vStart; v < vEnd; ++v) {
300: PetscInt dof, off, noff, d;
302: PetscSectionGetDof(coordSection, v, &dof);
303: PetscSectionGetOffset(coordSection, v, &off);
304: PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthEnd, depthEnd, depthShift), &noff);
305: for (d = 0; d < dof; ++d) {
306: newCoords[noff+d] = coords[off+d];
307: }
308: }
309: VecRestoreArray(coordinates, &coords);
310: VecRestoreArray(newCoordinates, &newCoords);
311: VecDestroy(&newCoordinates);
312: PetscSectionDestroy(&newCoordSection);
313: PetscFree(depthEnd);
314: return(0);
315: }
319: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
320: {
321: PetscInt *depthMax, *depthEnd;
322: PetscInt depth = 0, d;
323: PetscSF sfPoint, sfPointNew;
324: const PetscSFNode *remotePoints;
325: PetscSFNode *gremotePoints;
326: const PetscInt *localPoints;
327: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
328: PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
329: PetscErrorCode ierr;
332: DMPlexGetDepth(dm, &depth);
333: PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
334: DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
335: for (d = 0; d <= depth; ++d) {
336: totShift += depthShift[d];
337: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
338: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
339: }
340: /* Step 9: Convert pointSF */
341: DMGetPointSF(dm, &sfPoint);
342: DMGetPointSF(dmNew, &sfPointNew);
343: DMPlexGetChart(dm, &pStart, &pEnd);
344: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
345: if (numRoots >= 0) {
346: PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
347: for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthMax, depthEnd, depthShift);
348: PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
349: PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
350: PetscMalloc1(numLeaves, &glocalPoints);
351: PetscMalloc1(numLeaves, &gremotePoints);
352: for (l = 0; l < numLeaves; ++l) {
353: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthMax, depthEnd, depthShift);
354: gremotePoints[l].rank = remotePoints[l].rank;
355: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
356: }
357: PetscFree2(newLocation,newRemoteLocation);
358: PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
359: }
360: PetscFree2(depthMax,depthEnd);
361: return(0);
362: }
366: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
367: {
368: PetscSF sfPoint;
369: DMLabel vtkLabel, ghostLabel;
370: PetscInt *depthMax, *depthEnd;
371: const PetscSFNode *leafRemote;
372: const PetscInt *leafLocal;
373: PetscInt depth = 0, d, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
374: PetscMPIInt rank;
375: PetscErrorCode ierr;
378: DMPlexGetDepth(dm, &depth);
379: PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
380: DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
381: for (d = 0; d <= depth; ++d) {
382: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
383: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
384: }
385: /* Step 10: Convert labels */
386: DMPlexGetNumLabels(dm, &numLabels);
387: for (l = 0; l < numLabels; ++l) {
388: DMLabel label, newlabel;
389: const char *lname;
390: PetscBool isDepth;
391: IS valueIS;
392: const PetscInt *values;
393: PetscInt numValues, val;
395: DMPlexGetLabelName(dm, l, &lname);
396: PetscStrcmp(lname, "depth", &isDepth);
397: if (isDepth) continue;
398: DMPlexCreateLabel(dmNew, lname);
399: DMPlexGetLabel(dm, lname, &label);
400: DMPlexGetLabel(dmNew, lname, &newlabel);
401: DMLabelGetValueIS(label, &valueIS);
402: ISGetLocalSize(valueIS, &numValues);
403: ISGetIndices(valueIS, &values);
404: for (val = 0; val < numValues; ++val) {
405: IS pointIS;
406: const PetscInt *points;
407: PetscInt numPoints, p;
409: DMLabelGetStratumIS(label, values[val], &pointIS);
410: ISGetLocalSize(pointIS, &numPoints);
411: ISGetIndices(pointIS, &points);
412: for (p = 0; p < numPoints; ++p) {
413: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthMax, depthEnd, depthShift);
415: DMLabelSetValue(newlabel, newpoint, values[val]);
416: }
417: ISRestoreIndices(pointIS, &points);
418: ISDestroy(&pointIS);
419: }
420: ISRestoreIndices(valueIS, &values);
421: ISDestroy(&valueIS);
422: }
423: PetscFree2(depthMax,depthEnd);
424: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
425: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
426: DMGetPointSF(dm, &sfPoint);
427: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
428: PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
429: DMPlexCreateLabel(dmNew, "vtk");
430: DMPlexCreateLabel(dmNew, "ghost");
431: DMPlexGetLabel(dmNew, "vtk", &vtkLabel);
432: DMPlexGetLabel(dmNew, "ghost", &ghostLabel);
433: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
434: for (; c < leafLocal[l] && c < cEnd; ++c) {
435: DMLabelSetValue(vtkLabel, c, 1);
436: }
437: if (leafLocal[l] >= cEnd) break;
438: if (leafRemote[l].rank == rank) {
439: DMLabelSetValue(vtkLabel, c, 1);
440: } else {
441: DMLabelSetValue(ghostLabel, c, 2);
442: }
443: }
444: for (; c < cEnd; ++c) {
445: DMLabelSetValue(vtkLabel, c, 1);
446: }
447: if (0) {
448: PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
449: DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
450: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
451: }
452: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
453: for (f = fStart; f < fEnd; ++f) {
454: PetscInt numCells;
456: DMPlexGetSupportSize(dmNew, f, &numCells);
457: if (numCells < 2) {
458: DMLabelSetValue(ghostLabel, f, 1);
459: } else {
460: const PetscInt *cells = NULL;
461: PetscInt vA, vB;
463: DMPlexGetSupport(dmNew, f, &cells);
464: DMLabelGetValue(vtkLabel, cells[0], &vA);
465: DMLabelGetValue(vtkLabel, cells[1], &vB);
466: if (!vA && !vB) {DMLabelSetValue(ghostLabel, f, 1);}
467: }
468: }
469: if (0) {
470: PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
471: DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
472: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
473: }
474: return(0);
475: }
479: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
480: {
481: PetscSF sf;
482: IS valueIS;
483: const PetscInt *values, *leaves;
484: PetscInt *depthShift;
485: PetscInt depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
486: PetscErrorCode ierr;
489: DMGetPointSF(dm, &sf);
490: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
491: nleaves = PetscMax(0, nleaves);
492: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
493: /* Count ghost cells */
494: DMLabelGetValueIS(label, &valueIS);
495: ISGetLocalSize(valueIS, &numFS);
496: ISGetIndices(valueIS, &values);
497: Ng = 0;
498: for (fs = 0; fs < numFS; ++fs) {
499: IS faceIS;
500: const PetscInt *faces;
501: PetscInt numFaces, f, numBdFaces = 0;
503: DMLabelGetStratumIS(label, values[fs], &faceIS);
504: ISGetLocalSize(faceIS, &numFaces);
505: ISGetIndices(faceIS, &faces);
506: for (f = 0; f < numFaces; ++f) {
507: PetscFindInt(faces[f], nleaves, leaves, &loc);
508: if (loc >= 0) continue;
509: if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
510: }
511: Ng += numBdFaces;
512: ISDestroy(&faceIS);
513: }
514: DMPlexGetDepth(dm, &depth);
515: PetscMalloc1(depth+1, &depthShift);
516: PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));
517: if (depth >= 0) depthShift[depth] = Ng;
518: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
519: /* Step 3: Set cone/support sizes for new points */
520: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
521: DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
522: for (c = cEnd; c < cEnd + Ng; ++c) {
523: DMPlexSetConeSize(gdm, c, 1);
524: }
525: for (fs = 0; fs < numFS; ++fs) {
526: IS faceIS;
527: const PetscInt *faces;
528: PetscInt numFaces, f;
530: DMLabelGetStratumIS(label, values[fs], &faceIS);
531: ISGetLocalSize(faceIS, &numFaces);
532: ISGetIndices(faceIS, &faces);
533: for (f = 0; f < numFaces; ++f) {
534: PetscInt size;
536: PetscFindInt(faces[f], nleaves, leaves, &loc);
537: if (loc >= 0) continue;
538: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
539: DMPlexGetSupportSize(dm, faces[f], &size);
540: if (size != 1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
541: DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
542: }
543: ISRestoreIndices(faceIS, &faces);
544: ISDestroy(&faceIS);
545: }
546: /* Step 4: Setup ghosted DM */
547: DMSetUp(gdm);
548: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
549: /* Step 6: Set cones and supports for new points */
550: ghostCell = cEnd;
551: for (fs = 0; fs < numFS; ++fs) {
552: IS faceIS;
553: const PetscInt *faces;
554: PetscInt numFaces, f;
556: DMLabelGetStratumIS(label, values[fs], &faceIS);
557: ISGetLocalSize(faceIS, &numFaces);
558: ISGetIndices(faceIS, &faces);
559: for (f = 0; f < numFaces; ++f) {
560: PetscInt newFace = faces[f] + Ng;
562: PetscFindInt(faces[f], nleaves, leaves, &loc);
563: if (loc >= 0) continue;
564: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
565: DMPlexSetCone(gdm, ghostCell, &newFace);
566: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
567: ++ghostCell;
568: }
569: ISRestoreIndices(faceIS, &faces);
570: ISDestroy(&faceIS);
571: }
572: ISRestoreIndices(valueIS, &values);
573: ISDestroy(&valueIS);
574: /* Step 7: Stratify */
575: DMPlexStratify(gdm);
576: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
577: DMPlexShiftSF_Internal(dm, depthShift, gdm);
578: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
579: PetscFree(depthShift);
580: /* Step 7: Periodicity */
581: if (dm->maxCell) {
582: const PetscReal *maxCell, *L;
583: const DMBoundaryType *bd;
584: DMGetPeriodicity(dm, &maxCell, &L, &bd);
585: DMSetPeriodicity(gdm, maxCell, L, bd);
586: }
587: if (numGhostCells) *numGhostCells = Ng;
588: return(0);
589: }
593: /*@C
594: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
596: Collective on dm
598: Input Parameters:
599: + dm - The original DM
600: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
602: Output Parameters:
603: + numGhostCells - The number of ghost cells added to the DM
604: - dmGhosted - The new DM
606: Note: If no label exists of that name, one will be created marking all boundary faces
608: Level: developer
610: .seealso: DMCreate()
611: @*/
612: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
613: {
614: DM gdm;
615: DMLabel label;
616: const char *name = labelName ? labelName : "Face Sets";
617: PetscInt dim;
618: PetscBool flag;
625: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
626: DMSetType(gdm, DMPLEX);
627: DMGetDimension(dm, &dim);
628: DMSetDimension(gdm, dim);
629: DMPlexGetAdjacencyUseCone(dm, &flag);
630: DMPlexSetAdjacencyUseCone(gdm, flag);
631: DMPlexGetAdjacencyUseClosure(dm, &flag);
632: DMPlexSetAdjacencyUseClosure(gdm, flag);
633: DMPlexGetLabel(dm, name, &label);
634: if (!label) {
635: /* Get label for boundary faces */
636: DMPlexCreateLabel(dm, name);
637: DMPlexGetLabel(dm, name, &label);
638: DMPlexMarkBoundaryFaces(dm, label);
639: }
640: DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
641: DMPlexCopyBoundary(dm, gdm);
642: *dmGhosted = gdm;
643: return(0);
644: }
648: /*
649: We are adding three kinds of points here:
650: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
651: Non-replicated: Points which exist on the fault, but are not replicated
652: Hybrid: Entirely new points, such as cohesive cells
654: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
655: each depth so that the new split/hybrid points can be inserted as a block.
656: */
657: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
658: {
659: MPI_Comm comm;
660: IS valueIS;
661: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
662: const PetscInt *values; /* List of depths for which we have replicated points */
663: IS *splitIS;
664: IS *unsplitIS;
665: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
666: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
667: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
668: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
669: const PetscInt **splitPoints; /* Replicated points for each depth */
670: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
671: PetscSection coordSection;
672: Vec coordinates;
673: PetscScalar *coords;
674: PetscInt depths[4]; /* Depths in the order that plex points are numbered */
675: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
676: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
677: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
678: PetscInt *depthOffset; /* Prefix sums of depthShift */
679: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
680: PetscInt *coneNew, *coneONew, *supportNew;
681: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
682: PetscErrorCode ierr;
685: PetscObjectGetComm((PetscObject)dm,&comm);
686: DMGetDimension(dm, &dim);
687: DMPlexGetDepth(dm, &depth);
688: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
689: depths[0] = depth;
690: depths[1] = 0;
691: depths[2] = depth-1;
692: depths[3] = 1;
693: /* Count split points and add cohesive cells */
694: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
695: PetscMalloc6(depth+1,&depthMax,depth+1,&depthEnd,depth+1,&depthShift,depth+1,&depthOffset,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
696: PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
697: PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));
698: PetscMemzero(depthOffset, (depth+1) * sizeof(PetscInt));
699: DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
700: for (d = 0; d <= depth; ++d) {
701: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
702: depthEnd[d] = pMaxNew[d];
703: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
704: numSplitPoints[d] = 0;
705: numUnsplitPoints[d] = 0;
706: numHybridPoints[d] = 0;
707: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
708: splitPoints[d] = NULL;
709: unsplitPoints[d] = NULL;
710: splitIS[d] = NULL;
711: unsplitIS[d] = NULL;
712: }
713: if (label) {
714: DMLabelGetValueIS(label, &valueIS);
715: ISGetLocalSize(valueIS, &numSP);
716: ISGetIndices(valueIS, &values);
717: }
718: for (sp = 0; sp < numSP; ++sp) {
719: const PetscInt dep = values[sp];
721: if ((dep < 0) || (dep > depth)) continue;
722: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
723: if (splitIS[dep]) {
724: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
725: ISGetIndices(splitIS[dep], &splitPoints[dep]);
726: }
727: DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
728: if (unsplitIS[dep]) {
729: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
730: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
731: }
732: }
733: /* Calculate number of hybrid points */
734: 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 */
735: for (d = 0; d <= depth; ++d) depthShift[d] = numSplitPoints[d] + numHybridPoints[d];
736: for (d = 1; d <= depth; ++d) depthOffset[depths[d]] = depthOffset[depths[d-1]] + depthShift[depths[d-1]];
737: for (d = 0; d <= depth; ++d) pMaxNew[d] += depthOffset[d] - numHybridPointsOld[d];
738: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
739: /* Step 3: Set cone/support sizes for new points */
740: for (dep = 0; dep <= depth; ++dep) {
741: for (p = 0; p < numSplitPoints[dep]; ++p) {
742: const PetscInt oldp = splitPoints[dep][p];
743: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
744: const PetscInt splitp = p + pMaxNew[dep];
745: const PetscInt *support;
746: PetscInt coneSize, supportSize, qf, qn, qp, e;
748: DMPlexGetConeSize(dm, oldp, &coneSize);
749: DMPlexSetConeSize(sdm, splitp, coneSize);
750: DMPlexGetSupportSize(dm, oldp, &supportSize);
751: DMPlexSetSupportSize(sdm, splitp, supportSize);
752: if (dep == depth-1) {
753: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
755: /* Add cohesive cells, they are prisms */
756: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
757: } else if (dep == 0) {
758: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
760: DMPlexGetSupport(dm, oldp, &support);
761: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
762: PetscInt val;
764: DMLabelGetValue(label, support[e], &val);
765: if (val == 1) ++qf;
766: if ((val == 1) || (val == (shift + 1))) ++qn;
767: if ((val == 1) || (val == -(shift + 1))) ++qp;
768: }
769: /* Split old vertex: Edges into original vertex and new cohesive edge */
770: DMPlexSetSupportSize(sdm, newp, qn+1);
771: /* Split new vertex: Edges into split vertex and new cohesive edge */
772: DMPlexSetSupportSize(sdm, splitp, qp+1);
773: /* Add hybrid edge */
774: DMPlexSetConeSize(sdm, hybedge, 2);
775: DMPlexSetSupportSize(sdm, hybedge, qf);
776: } else if (dep == dim-2) {
777: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
779: DMPlexGetSupport(dm, oldp, &support);
780: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
781: PetscInt val;
783: DMLabelGetValue(label, support[e], &val);
784: if (val == dim-1) ++qf;
785: if ((val == dim-1) || (val == (shift + dim-1))) ++qn;
786: if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
787: }
788: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
789: DMPlexSetSupportSize(sdm, newp, qn+1);
790: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
791: DMPlexSetSupportSize(sdm, splitp, qp+1);
792: /* Add hybrid face */
793: DMPlexSetConeSize(sdm, hybface, 4);
794: DMPlexSetSupportSize(sdm, hybface, qf);
795: }
796: }
797: }
798: for (dep = 0; dep <= depth; ++dep) {
799: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
800: const PetscInt oldp = unsplitPoints[dep][p];
801: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
802: const PetscInt *support;
803: PetscInt coneSize, supportSize, qf, e, s;
805: DMPlexGetConeSize(dm, oldp, &coneSize);
806: DMPlexGetSupportSize(dm, oldp, &supportSize);
807: DMPlexGetSupport(dm, oldp, &support);
808: if (dep == 0) {
809: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
811: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
812: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
813: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
814: if (e >= 0) ++qf;
815: }
816: DMPlexSetSupportSize(sdm, newp, qf+2);
817: /* Add hybrid edge */
818: DMPlexSetConeSize(sdm, hybedge, 2);
819: for (e = 0, qf = 0; e < supportSize; ++e) {
820: PetscInt val;
822: DMLabelGetValue(label, support[e], &val);
823: /* Split and unsplit edges produce hybrid faces */
824: if (val == 1) ++qf;
825: if (val == (shift2 + 1)) ++qf;
826: }
827: DMPlexSetSupportSize(sdm, hybedge, qf);
828: } else if (dep == dim-2) {
829: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
830: PetscInt val;
832: for (e = 0, qf = 0; e < supportSize; ++e) {
833: DMLabelGetValue(label, support[e], &val);
834: if (val == dim-1) qf += 2;
835: else ++qf;
836: }
837: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
838: DMPlexSetSupportSize(sdm, newp, qf+2);
839: /* Add hybrid face */
840: for (e = 0, qf = 0; e < supportSize; ++e) {
841: DMLabelGetValue(label, support[e], &val);
842: if (val == dim-1) ++qf;
843: }
844: DMPlexSetConeSize(sdm, hybface, 4);
845: DMPlexSetSupportSize(sdm, hybface, qf);
846: }
847: }
848: }
849: /* Step 4: Setup split DM */
850: DMSetUp(sdm);
851: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
852: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
853: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
854: /* Step 6: Set cones and supports for new points */
855: for (dep = 0; dep <= depth; ++dep) {
856: for (p = 0; p < numSplitPoints[dep]; ++p) {
857: const PetscInt oldp = splitPoints[dep][p];
858: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
859: const PetscInt splitp = p + pMaxNew[dep];
860: const PetscInt *cone, *support, *ornt;
861: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
863: DMPlexGetConeSize(dm, oldp, &coneSize);
864: DMPlexGetCone(dm, oldp, &cone);
865: DMPlexGetConeOrientation(dm, oldp, &ornt);
866: DMPlexGetSupportSize(dm, oldp, &supportSize);
867: DMPlexGetSupport(dm, oldp, &support);
868: if (dep == depth-1) {
869: PetscBool hasUnsplit = PETSC_FALSE;
870: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
871: const PetscInt *supportF;
873: /* Split face: copy in old face to new face to start */
874: DMPlexGetSupport(sdm, newp, &supportF);
875: DMPlexSetSupport(sdm, splitp, supportF);
876: /* Split old face: old vertices/edges in cone so no change */
877: /* Split new face: new vertices/edges in cone */
878: for (q = 0; q < coneSize; ++q) {
879: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
880: if (v < 0) {
881: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
882: 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);
883: coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/;
884: hasUnsplit = PETSC_TRUE;
885: } else {
886: coneNew[2+q] = v + pMaxNew[dep-1];
887: if (dep > 1) {
888: const PetscInt *econe;
889: PetscInt econeSize, r, vs, vu;
891: DMPlexGetConeSize(dm, cone[q], &econeSize);
892: DMPlexGetCone(dm, cone[q], &econe);
893: for (r = 0; r < econeSize; ++r) {
894: PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);
895: PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
896: if (vs >= 0) continue;
897: 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);
898: hasUnsplit = PETSC_TRUE;
899: }
900: }
901: }
902: }
903: DMPlexSetCone(sdm, splitp, &coneNew[2]);
904: DMPlexSetConeOrientation(sdm, splitp, ornt);
905: /* Face support */
906: for (s = 0; s < supportSize; ++s) {
907: PetscInt val;
909: DMLabelGetValue(label, support[s], &val);
910: if (val < 0) {
911: /* Split old face: Replace negative side cell with cohesive cell */
912: DMPlexInsertSupport(sdm, newp, s, hybcell);
913: } else {
914: /* Split new face: Replace positive side cell with cohesive cell */
915: DMPlexInsertSupport(sdm, splitp, s, hybcell);
916: /* Get orientation for cohesive face */
917: {
918: const PetscInt *ncone, *nconeO;
919: PetscInt nconeSize, nc;
921: DMPlexGetConeSize(dm, support[s], &nconeSize);
922: DMPlexGetCone(dm, support[s], &ncone);
923: DMPlexGetConeOrientation(dm, support[s], &nconeO);
924: for (nc = 0; nc < nconeSize; ++nc) {
925: if (ncone[nc] == oldp) {
926: coneONew[0] = nconeO[nc];
927: break;
928: }
929: }
930: if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
931: }
932: }
933: }
934: /* Cohesive cell: Old and new split face, then new cohesive faces */
935: coneNew[0] = newp; /* Extracted negative side orientation above */
936: coneNew[1] = splitp;
937: coneONew[1] = coneONew[0];
938: for (q = 0; q < coneSize; ++q) {
939: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
940: if (v < 0) {
941: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
942: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
943: coneONew[2+q] = 0;
944: } else {
945: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep];
946: }
947: coneONew[2+q] = 0;
948: }
949: DMPlexSetCone(sdm, hybcell, coneNew);
950: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
951: /* Label the hybrid cells on the boundary of the split */
952: if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
953: } else if (dep == 0) {
954: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
956: /* Split old vertex: Edges in old split faces and new cohesive edge */
957: for (e = 0, qn = 0; e < supportSize; ++e) {
958: PetscInt val;
960: DMLabelGetValue(label, support[e], &val);
961: if ((val == 1) || (val == (shift + 1))) {
962: supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
963: }
964: }
965: supportNew[qn] = hybedge;
966: DMPlexSetSupport(sdm, newp, supportNew);
967: /* Split new vertex: Edges in new split faces and new cohesive edge */
968: for (e = 0, qp = 0; e < supportSize; ++e) {
969: PetscInt val, edge;
971: DMLabelGetValue(label, support[e], &val);
972: if (val == 1) {
973: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
974: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
975: supportNew[qp++] = edge + pMaxNew[dep+1];
976: } else if (val == -(shift + 1)) {
977: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
978: }
979: }
980: supportNew[qp] = hybedge;
981: DMPlexSetSupport(sdm, splitp, supportNew);
982: /* Hybrid edge: Old and new split vertex */
983: coneNew[0] = newp;
984: coneNew[1] = splitp;
985: DMPlexSetCone(sdm, hybedge, coneNew);
986: for (e = 0, qf = 0; e < supportSize; ++e) {
987: PetscInt val, edge;
989: DMLabelGetValue(label, support[e], &val);
990: if (val == 1) {
991: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
992: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
993: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
994: }
995: }
996: DMPlexSetSupport(sdm, hybedge, supportNew);
997: } else if (dep == dim-2) {
998: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1000: /* Split old edge: old vertices in cone so no change */
1001: /* Split new edge: new vertices in cone */
1002: for (q = 0; q < coneSize; ++q) {
1003: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1004: if (v < 0) {
1005: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1006: 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);
1007: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1008: } else {
1009: coneNew[q] = v + pMaxNew[dep-1];
1010: }
1011: }
1012: DMPlexSetCone(sdm, splitp, coneNew);
1013: /* Split old edge: Faces in positive side cells and old split faces */
1014: for (e = 0, q = 0; e < supportSize; ++e) {
1015: PetscInt val;
1017: DMLabelGetValue(label, support[e], &val);
1018: if (val == dim-1) {
1019: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
1020: } else if (val == (shift + dim-1)) {
1021: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
1022: }
1023: }
1024: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1025: DMPlexSetSupport(sdm, newp, supportNew);
1026: /* Split new edge: Faces in negative side cells and new split faces */
1027: for (e = 0, q = 0; e < supportSize; ++e) {
1028: PetscInt val, face;
1030: DMLabelGetValue(label, support[e], &val);
1031: if (val == dim-1) {
1032: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1033: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1034: supportNew[q++] = face + pMaxNew[dep+1];
1035: } else if (val == -(shift + dim-1)) {
1036: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
1037: }
1038: }
1039: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1040: DMPlexSetSupport(sdm, splitp, supportNew);
1041: /* Hybrid face */
1042: coneNew[0] = newp;
1043: coneNew[1] = splitp;
1044: for (v = 0; v < coneSize; ++v) {
1045: PetscInt vertex;
1046: PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1047: if (vertex < 0) {
1048: PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1049: 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);
1050: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1051: } else {
1052: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1053: }
1054: }
1055: DMPlexSetCone(sdm, hybface, coneNew);
1056: for (e = 0, qf = 0; e < supportSize; ++e) {
1057: PetscInt val, face;
1059: DMLabelGetValue(label, support[e], &val);
1060: if (val == dim-1) {
1061: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1062: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1063: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1064: }
1065: }
1066: DMPlexSetSupport(sdm, hybface, supportNew);
1067: }
1068: }
1069: }
1070: for (dep = 0; dep <= depth; ++dep) {
1071: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1072: const PetscInt oldp = unsplitPoints[dep][p];
1073: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
1074: const PetscInt *cone, *support, *ornt;
1075: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1077: DMPlexGetConeSize(dm, oldp, &coneSize);
1078: DMPlexGetCone(dm, oldp, &cone);
1079: DMPlexGetConeOrientation(dm, oldp, &ornt);
1080: DMPlexGetSupportSize(dm, oldp, &supportSize);
1081: DMPlexGetSupport(dm, oldp, &support);
1082: if (dep == 0) {
1083: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1085: /* Unsplit vertex */
1086: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1087: for (s = 0, q = 0; s < supportSize; ++s) {
1088: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthMax, depthEnd, depthShift) /*support[s] + depthOffset[dep+1]*/;
1089: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1090: if (e >= 0) {
1091: supportNew[q++] = e + pMaxNew[dep+1];
1092: }
1093: }
1094: supportNew[q++] = hybedge;
1095: supportNew[q++] = hybedge;
1096: if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1097: DMPlexSetSupport(sdm, newp, supportNew);
1098: /* Hybrid edge */
1099: coneNew[0] = newp;
1100: coneNew[1] = newp;
1101: DMPlexSetCone(sdm, hybedge, coneNew);
1102: for (e = 0, qf = 0; e < supportSize; ++e) {
1103: PetscInt val, edge;
1105: DMLabelGetValue(label, support[e], &val);
1106: if (val == 1) {
1107: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1108: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1109: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1110: } else if (val == (shift2 + 1)) {
1111: PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1112: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1113: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1114: }
1115: }
1116: DMPlexSetSupport(sdm, hybedge, supportNew);
1117: } else if (dep == dim-2) {
1118: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1120: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1121: for (f = 0, qf = 0; f < supportSize; ++f) {
1122: PetscInt val, face;
1124: DMLabelGetValue(label, support[f], &val);
1125: if (val == dim-1) {
1126: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1127: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1128: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthMax, depthEnd, depthShift) /*support[f] + depthOffset[dep+1]*/;
1129: supportNew[qf++] = face + pMaxNew[dep+1];
1130: } else {
1131: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthMax, depthEnd, depthShift) /*support[f] + depthOffset[dep+1]*/;
1132: }
1133: }
1134: supportNew[qf++] = hybface;
1135: supportNew[qf++] = hybface;
1136: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1137: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1138: DMPlexSetSupport(sdm, newp, supportNew);
1139: /* Add hybrid face */
1140: coneNew[0] = newp;
1141: coneNew[1] = newp;
1142: PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1143: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1144: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1145: PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1146: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1147: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1148: DMPlexSetCone(sdm, hybface, coneNew);
1149: for (f = 0, qf = 0; f < supportSize; ++f) {
1150: PetscInt val, face;
1152: DMLabelGetValue(label, support[f], &val);
1153: if (val == dim-1) {
1154: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1155: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1156: }
1157: }
1158: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1159: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1160: DMPlexSetSupport(sdm, hybface, supportNew);
1161: }
1162: }
1163: }
1164: /* Step 6b: Replace split points in negative side cones */
1165: for (sp = 0; sp < numSP; ++sp) {
1166: PetscInt dep = values[sp];
1167: IS pIS;
1168: PetscInt numPoints;
1169: const PetscInt *points;
1171: if (dep >= 0) continue;
1172: DMLabelGetStratumIS(label, dep, &pIS);
1173: if (!pIS) continue;
1174: dep = -dep - shift;
1175: ISGetLocalSize(pIS, &numPoints);
1176: ISGetIndices(pIS, &points);
1177: for (p = 0; p < numPoints; ++p) {
1178: const PetscInt oldp = points[p];
1179: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + oldp*/;
1180: const PetscInt *cone;
1181: PetscInt coneSize, c;
1182: /* PetscBool replaced = PETSC_FALSE; */
1184: /* Negative edge: replace split vertex */
1185: /* Negative cell: replace split face */
1186: DMPlexGetConeSize(sdm, newp, &coneSize);
1187: DMPlexGetCone(sdm, newp, &cone);
1188: for (c = 0; c < coneSize; ++c) {
1189: const PetscInt coldp = cone[c] - depthOffset[dep-1];
1190: PetscInt csplitp, cp, val;
1192: DMLabelGetValue(label, coldp, &val);
1193: if (val == dep-1) {
1194: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1195: if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1196: csplitp = pMaxNew[dep-1] + cp;
1197: DMPlexInsertCone(sdm, newp, c, csplitp);
1198: /* replaced = PETSC_TRUE; */
1199: }
1200: }
1201: /* Cells with only a vertex or edge on the submesh have no replacement */
1202: /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1203: }
1204: ISRestoreIndices(pIS, &points);
1205: ISDestroy(&pIS);
1206: }
1207: /* Step 7: Stratify */
1208: DMPlexStratify(sdm);
1209: /* Step 8: Coordinates */
1210: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1211: DMGetCoordinateSection(sdm, &coordSection);
1212: DMGetCoordinatesLocal(sdm, &coordinates);
1213: VecGetArray(coordinates, &coords);
1214: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1215: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthMax, depthEnd, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1216: const PetscInt splitp = pMaxNew[0] + v;
1217: PetscInt dof, off, soff, d;
1219: PetscSectionGetDof(coordSection, newp, &dof);
1220: PetscSectionGetOffset(coordSection, newp, &off);
1221: PetscSectionGetOffset(coordSection, splitp, &soff);
1222: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1223: }
1224: VecRestoreArray(coordinates, &coords);
1225: /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1226: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1227: /* Step 10: Labels */
1228: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1229: DMPlexGetNumLabels(sdm, &numLabels);
1230: for (dep = 0; dep <= depth; ++dep) {
1231: for (p = 0; p < numSplitPoints[dep]; ++p) {
1232: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1233: const PetscInt splitp = pMaxNew[dep] + p;
1234: PetscInt l;
1236: for (l = 0; l < numLabels; ++l) {
1237: DMLabel mlabel;
1238: const char *lname;
1239: PetscInt val;
1240: PetscBool isDepth;
1242: DMPlexGetLabelName(sdm, l, &lname);
1243: PetscStrcmp(lname, "depth", &isDepth);
1244: if (isDepth) continue;
1245: DMPlexGetLabel(sdm, lname, &mlabel);
1246: DMLabelGetValue(mlabel, newp, &val);
1247: if (val >= 0) {
1248: DMLabelSetValue(mlabel, splitp, val);
1249: #if 0
1250: /* Do not put cohesive edges into the label */
1251: if (dep == 0) {
1252: const PetscInt cedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1253: DMLabelSetValue(mlabel, cedge, val);
1254: } else if (dep == dim-2) {
1255: const PetscInt cface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1256: DMLabelSetValue(mlabel, cface, val);
1257: }
1258: /* Do not put cohesive faces into the label */
1259: #endif
1260: }
1261: }
1262: }
1263: }
1264: for (sp = 0; sp < numSP; ++sp) {
1265: const PetscInt dep = values[sp];
1267: if ((dep < 0) || (dep > depth)) continue;
1268: if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1269: ISDestroy(&splitIS[dep]);
1270: if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1271: ISDestroy(&unsplitIS[dep]);
1272: }
1273: if (label) {
1274: ISRestoreIndices(valueIS, &values);
1275: ISDestroy(&valueIS);
1276: }
1277: for (d = 0; d <= depth; ++d) {
1278: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1279: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1280: }
1281: 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);
1282: PetscFree3(coneNew, coneONew, supportNew);
1283: PetscFree6(depthMax, depthEnd, depthShift, depthOffset, pMaxNew, numHybridPointsOld);
1284: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1285: return(0);
1286: }
1290: /*@C
1291: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1293: Collective on dm
1295: Input Parameters:
1296: + dm - The original DM
1297: - label - The label specifying the boundary faces (this could be auto-generated)
1299: Output Parameters:
1300: - dmSplit - The new DM
1302: Level: developer
1304: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1305: @*/
1306: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1307: {
1308: DM sdm;
1309: PetscInt dim;
1315: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1316: DMSetType(sdm, DMPLEX);
1317: DMGetDimension(dm, &dim);
1318: DMSetDimension(sdm, dim);
1319: switch (dim) {
1320: case 2:
1321: case 3:
1322: DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1323: break;
1324: default:
1325: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1326: }
1327: *dmSplit = sdm;
1328: return(0);
1329: }
1333: /* Returns the side of the surface for a given cell with a face on the surface */
1334: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1335: {
1336: const PetscInt *cone, *ornt;
1337: PetscInt dim, coneSize, c;
1338: PetscErrorCode ierr;
1341: *pos = PETSC_TRUE;
1342: DMGetDimension(dm, &dim);
1343: DMPlexGetConeSize(dm, cell, &coneSize);
1344: DMPlexGetCone(dm, cell, &cone);
1345: DMPlexGetConeOrientation(dm, cell, &ornt);
1346: for (c = 0; c < coneSize; ++c) {
1347: if (cone[c] == face) {
1348: PetscInt o = ornt[c];
1350: if (subdm) {
1351: const PetscInt *subcone, *subornt;
1352: PetscInt subpoint, subface, subconeSize, sc;
1354: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1355: PetscFindInt(face, numSubpoints, subpoints, &subface);
1356: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1357: DMPlexGetCone(subdm, subpoint, &subcone);
1358: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1359: for (sc = 0; sc < subconeSize; ++sc) {
1360: if (subcone[sc] == subface) {
1361: o = subornt[0];
1362: break;
1363: }
1364: }
1365: 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);
1366: }
1367: if (o >= 0) *pos = PETSC_TRUE;
1368: else *pos = PETSC_FALSE;
1369: break;
1370: }
1371: }
1372: 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);
1373: return(0);
1374: }
1378: /*@
1379: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1380: to complete the surface
1382: Input Parameters:
1383: + dm - The DM
1384: . label - A DMLabel marking the surface
1385: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1386: . flip - Flag to flip the submesh normal and replace points on the other side
1387: - subdm - The subDM associated with the label, or NULL
1389: Output Parameter:
1390: . label - A DMLabel marking all surface points
1392: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1394: Level: developer
1396: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1397: @*/
1398: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1399: {
1400: DMLabel depthLabel;
1401: IS dimIS, subpointIS, facePosIS, faceNegIS, crossEdgeIS = NULL;
1402: const PetscInt *points, *subpoints;
1403: const PetscInt rev = flip ? -1 : 1;
1404: PetscInt *pMax;
1405: PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1406: PetscErrorCode ierr;
1409: DMPlexGetDepth(dm, &depth);
1410: DMGetDimension(dm, &dim);
1411: pSize = PetscMax(depth, dim) + 1;
1412: PetscMalloc1(pSize,&pMax);
1413: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1414: DMPlexGetDepthLabel(dm, &depthLabel);
1415: DMGetDimension(dm, &dim);
1416: if (subdm) {
1417: DMPlexCreateSubpointIS(subdm, &subpointIS);
1418: if (subpointIS) {
1419: ISGetLocalSize(subpointIS, &numSubpoints);
1420: ISGetIndices(subpointIS, &subpoints);
1421: }
1422: }
1423: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1424: DMLabelGetStratumIS(label, dim-1, &dimIS);
1425: if (!dimIS) {
1426: PetscFree(pMax);
1427: return(0);
1428: }
1429: ISGetLocalSize(dimIS, &numPoints);
1430: ISGetIndices(dimIS, &points);
1431: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1432: const PetscInt *support;
1433: PetscInt supportSize, s;
1435: DMPlexGetSupportSize(dm, points[p], &supportSize);
1436: if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1437: DMPlexGetSupport(dm, points[p], &support);
1438: for (s = 0; s < supportSize; ++s) {
1439: const PetscInt *cone;
1440: PetscInt coneSize, c;
1441: PetscBool pos;
1443: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1444: if (pos) {DMLabelSetValue(label, support[s], rev*(shift+dim));}
1445: else {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1446: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1447: /* Put faces touching the fault in the label */
1448: DMPlexGetConeSize(dm, support[s], &coneSize);
1449: DMPlexGetCone(dm, support[s], &cone);
1450: for (c = 0; c < coneSize; ++c) {
1451: const PetscInt point = cone[c];
1453: DMLabelGetValue(label, point, &val);
1454: if (val == -1) {
1455: PetscInt *closure = NULL;
1456: PetscInt closureSize, cl;
1458: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1459: for (cl = 0; cl < closureSize*2; cl += 2) {
1460: const PetscInt clp = closure[cl];
1461: PetscInt bval = -1;
1463: DMLabelGetValue(label, clp, &val);
1464: if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1465: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1466: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1467: break;
1468: }
1469: }
1470: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1471: }
1472: }
1473: }
1474: }
1475: ISRestoreIndices(dimIS, &points);
1476: ISDestroy(&dimIS);
1477: if (subdm) {
1478: if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1479: ISDestroy(&subpointIS);
1480: }
1481: /* Mark boundary points as unsplit */
1482: if (blabel) {
1483: DMLabelGetStratumIS(blabel, 1, &dimIS);
1484: ISGetLocalSize(dimIS, &numPoints);
1485: ISGetIndices(dimIS, &points);
1486: for (p = 0; p < numPoints; ++p) {
1487: const PetscInt point = points[p];
1488: PetscInt val, bval;
1490: DMLabelGetValue(blabel, point, &bval);
1491: if (bval >= 0) {
1492: DMLabelGetValue(label, point, &val);
1493: if ((val < 0) || (val > dim)) {
1494: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1495: DMLabelClearValue(blabel, point, bval);
1496: }
1497: }
1498: }
1499: for (p = 0; p < numPoints; ++p) {
1500: const PetscInt point = points[p];
1501: PetscInt val, bval;
1503: DMLabelGetValue(blabel, point, &bval);
1504: if (bval >= 0) {
1505: const PetscInt *cone, *support;
1506: PetscInt coneSize, supportSize, s, valA, valB, valE;
1508: /* Mark as unsplit */
1509: DMLabelGetValue(label, point, &val);
1510: 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);
1511: DMLabelClearValue(label, point, val);
1512: DMLabelSetValue(label, point, shift2+val);
1513: /* Check for cross-edge
1514: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1515: if (val != 0) continue;
1516: DMPlexGetSupport(dm, point, &support);
1517: DMPlexGetSupportSize(dm, point, &supportSize);
1518: for (s = 0; s < supportSize; ++s) {
1519: DMPlexGetCone(dm, support[s], &cone);
1520: DMPlexGetConeSize(dm, support[s], &coneSize);
1521: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1522: DMLabelGetValue(blabel, cone[0], &valA);
1523: DMLabelGetValue(blabel, cone[1], &valB);
1524: DMLabelGetValue(blabel, support[s], &valE);
1525: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1526: }
1527: }
1528: }
1529: ISRestoreIndices(dimIS, &points);
1530: ISDestroy(&dimIS);
1531: }
1532: /* Search for other cells/faces/edges connected to the fault by a vertex */
1533: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1534: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1535: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1536: cMax = cMax < 0 ? cEnd : cMax;
1537: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1538: DMLabelGetStratumIS(label, 0, &dimIS);
1539: if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1540: if (dimIS && crossEdgeIS) {
1541: IS vertIS = dimIS;
1543: ISExpand(vertIS, crossEdgeIS, &dimIS);
1544: ISDestroy(&crossEdgeIS);
1545: ISDestroy(&vertIS);
1546: }
1547: if (!dimIS) {
1548: PetscFree(pMax);
1549: return(0);
1550: }
1551: ISGetLocalSize(dimIS, &numPoints);
1552: ISGetIndices(dimIS, &points);
1553: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1554: PetscInt *star = NULL;
1555: PetscInt starSize, s;
1556: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1558: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1559: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1560: while (again) {
1561: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1562: again = 0;
1563: for (s = 0; s < starSize*2; s += 2) {
1564: const PetscInt point = star[s];
1565: const PetscInt *cone;
1566: PetscInt coneSize, c;
1568: if ((point < cStart) || (point >= cMax)) continue;
1569: DMLabelGetValue(label, point, &val);
1570: if (val != -1) continue;
1571: again = again == 1 ? 1 : 2;
1572: DMPlexGetConeSize(dm, point, &coneSize);
1573: DMPlexGetCone(dm, point, &cone);
1574: for (c = 0; c < coneSize; ++c) {
1575: DMLabelGetValue(label, cone[c], &val);
1576: if (val != -1) {
1577: const PetscInt *ccone;
1578: PetscInt cconeSize, cc, side;
1580: 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);
1581: if (val > 0) side = 1;
1582: else side = -1;
1583: DMLabelSetValue(label, point, side*(shift+dim));
1584: /* Mark cell faces which touch the fault */
1585: DMPlexGetConeSize(dm, point, &cconeSize);
1586: DMPlexGetCone(dm, point, &ccone);
1587: for (cc = 0; cc < cconeSize; ++cc) {
1588: PetscInt *closure = NULL;
1589: PetscInt closureSize, cl;
1591: DMLabelGetValue(label, ccone[cc], &val);
1592: if (val != -1) continue;
1593: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1594: for (cl = 0; cl < closureSize*2; cl += 2) {
1595: const PetscInt clp = closure[cl];
1597: DMLabelGetValue(label, clp, &val);
1598: if (val == -1) continue;
1599: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1600: break;
1601: }
1602: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1603: }
1604: again = 1;
1605: break;
1606: }
1607: }
1608: }
1609: }
1610: /* Classify the rest by cell membership */
1611: for (s = 0; s < starSize*2; s += 2) {
1612: const PetscInt point = star[s];
1614: DMLabelGetValue(label, point, &val);
1615: if (val == -1) {
1616: PetscInt *sstar = NULL;
1617: PetscInt sstarSize, ss;
1618: PetscBool marked = PETSC_FALSE;
1620: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1621: for (ss = 0; ss < sstarSize*2; ss += 2) {
1622: const PetscInt spoint = sstar[ss];
1624: if ((spoint < cStart) || (spoint >= cMax)) continue;
1625: DMLabelGetValue(label, spoint, &val);
1626: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1627: DMLabelGetValue(depthLabel, point, &dep);
1628: if (val > 0) {
1629: DMLabelSetValue(label, point, shift+dep);
1630: } else {
1631: DMLabelSetValue(label, point, -(shift+dep));
1632: }
1633: marked = PETSC_TRUE;
1634: break;
1635: }
1636: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1637: DMLabelGetValue(depthLabel, point, &dep);
1638: if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1639: }
1640: }
1641: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1642: }
1643: ISRestoreIndices(dimIS, &points);
1644: ISDestroy(&dimIS);
1645: /* If any faces touching the fault divide cells on either side, split them */
1646: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1647: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1648: ISExpand(facePosIS, faceNegIS, &dimIS);
1649: ISDestroy(&facePosIS);
1650: ISDestroy(&faceNegIS);
1651: ISGetLocalSize(dimIS, &numPoints);
1652: ISGetIndices(dimIS, &points);
1653: for (p = 0; p < numPoints; ++p) {
1654: const PetscInt point = points[p];
1655: const PetscInt *support;
1656: PetscInt supportSize, valA, valB;
1658: DMPlexGetSupportSize(dm, point, &supportSize);
1659: if (supportSize != 2) continue;
1660: DMPlexGetSupport(dm, point, &support);
1661: DMLabelGetValue(label, support[0], &valA);
1662: DMLabelGetValue(label, support[1], &valB);
1663: if ((valA == -1) || (valB == -1)) continue;
1664: if (valA*valB > 0) continue;
1665: /* Split the face */
1666: DMLabelGetValue(label, point, &valA);
1667: DMLabelClearValue(label, point, valA);
1668: DMLabelSetValue(label, point, dim-1);
1669: /* Label its closure:
1670: unmarked: label as unsplit
1671: incident: relabel as split
1672: split: do nothing
1673: */
1674: {
1675: PetscInt *closure = NULL;
1676: PetscInt closureSize, cl;
1678: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1679: for (cl = 0; cl < closureSize*2; cl += 2) {
1680: DMLabelGetValue(label, closure[cl], &valA);
1681: if (valA == -1) { /* Mark as unsplit */
1682: DMLabelGetValue(depthLabel, closure[cl], &dep);
1683: DMLabelSetValue(label, closure[cl], shift2+dep);
1684: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1685: DMLabelGetValue(depthLabel, closure[cl], &dep);
1686: DMLabelClearValue(label, closure[cl], valA);
1687: DMLabelSetValue(label, closure[cl], dep);
1688: }
1689: }
1690: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1691: }
1692: }
1693: ISRestoreIndices(dimIS, &points);
1694: ISDestroy(&dimIS);
1695: PetscFree(pMax);
1696: return(0);
1697: }
1701: /*@C
1702: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
1704: Collective on dm
1706: Input Parameters:
1707: + dm - The original DM
1708: - labelName - The label specifying the interface vertices
1710: Output Parameters:
1711: + hybridLabel - The label fully marking the interface
1712: - dmHybrid - The new DM
1714: Level: developer
1716: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1717: @*/
1718: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1719: {
1720: DM idm;
1721: DMLabel subpointMap, hlabel;
1722: PetscInt dim;
1729: DMGetDimension(dm, &dim);
1730: DMPlexCreateSubmesh(dm, label, 1, &idm);
1731: DMPlexOrient(idm);
1732: DMPlexGetSubpointMap(idm, &subpointMap);
1733: DMLabelDuplicate(subpointMap, &hlabel);
1734: DMLabelClearStratum(hlabel, dim);
1735: DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1736: DMDestroy(&idm);
1737: DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1738: if (hybridLabel) *hybridLabel = hlabel;
1739: else {DMLabelDestroy(&hlabel);}
1740: return(0);
1741: }
1745: /* Here we need the explicit assumption that:
1747: For any marked cell, the marked vertices constitute a single face
1748: */
1749: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1750: {
1751: IS subvertexIS = NULL;
1752: const PetscInt *subvertices;
1753: PetscInt *pStart, *pEnd, *pMax, pSize;
1754: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
1755: PetscErrorCode ierr;
1758: *numFaces = 0;
1759: *nFV = 0;
1760: DMPlexGetDepth(dm, &depth);
1761: DMGetDimension(dm, &dim);
1762: pSize = PetscMax(depth, dim) + 1;
1763: PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1764: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1765: for (d = 0; d <= depth; ++d) {
1766: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1767: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1768: }
1769: /* Loop over initial vertices and mark all faces in the collective star() */
1770: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1771: if (subvertexIS) {
1772: ISGetSize(subvertexIS, &numSubVerticesInitial);
1773: ISGetIndices(subvertexIS, &subvertices);
1774: }
1775: for (v = 0; v < numSubVerticesInitial; ++v) {
1776: const PetscInt vertex = subvertices[v];
1777: PetscInt *star = NULL;
1778: PetscInt starSize, s, numCells = 0, c;
1780: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1781: for (s = 0; s < starSize*2; s += 2) {
1782: const PetscInt point = star[s];
1783: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1784: }
1785: for (c = 0; c < numCells; ++c) {
1786: const PetscInt cell = star[c];
1787: PetscInt *closure = NULL;
1788: PetscInt closureSize, cl;
1789: PetscInt cellLoc, numCorners = 0, faceSize = 0;
1791: DMLabelGetValue(subpointMap, cell, &cellLoc);
1792: if (cellLoc == 2) continue;
1793: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1794: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1795: for (cl = 0; cl < closureSize*2; cl += 2) {
1796: const PetscInt point = closure[cl];
1797: PetscInt vertexLoc;
1799: if ((point >= pStart[0]) && (point < pEnd[0])) {
1800: ++numCorners;
1801: DMLabelGetValue(vertexLabel, point, &vertexLoc);
1802: if (vertexLoc == value) closure[faceSize++] = point;
1803: }
1804: }
1805: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1806: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1807: if (faceSize == *nFV) {
1808: const PetscInt *cells = NULL;
1809: PetscInt numCells, nc;
1811: ++(*numFaces);
1812: for (cl = 0; cl < faceSize; ++cl) {
1813: DMLabelSetValue(subpointMap, closure[cl], 0);
1814: }
1815: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1816: for (nc = 0; nc < numCells; ++nc) {
1817: DMLabelSetValue(subpointMap, cells[nc], 2);
1818: }
1819: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1820: }
1821: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1822: }
1823: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1824: }
1825: if (subvertexIS) {
1826: ISRestoreIndices(subvertexIS, &subvertices);
1827: }
1828: ISDestroy(&subvertexIS);
1829: PetscFree3(pStart,pEnd,pMax);
1830: return(0);
1831: }
1835: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
1836: {
1837: IS subvertexIS = NULL;
1838: const PetscInt *subvertices;
1839: PetscInt *pStart, *pEnd, *pMax;
1840: PetscInt dim, d, numSubVerticesInitial = 0, v;
1841: PetscErrorCode ierr;
1844: DMGetDimension(dm, &dim);
1845: PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
1846: DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
1847: for (d = 0; d <= dim; ++d) {
1848: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1849: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1850: }
1851: /* Loop over initial vertices and mark all faces in the collective star() */
1852: if (vertexLabel) {
1853: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
1854: if (subvertexIS) {
1855: ISGetSize(subvertexIS, &numSubVerticesInitial);
1856: ISGetIndices(subvertexIS, &subvertices);
1857: }
1858: }
1859: for (v = 0; v < numSubVerticesInitial; ++v) {
1860: const PetscInt vertex = subvertices[v];
1861: PetscInt *star = NULL;
1862: PetscInt starSize, s, numFaces = 0, f;
1864: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1865: for (s = 0; s < starSize*2; s += 2) {
1866: const PetscInt point = star[s];
1867: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
1868: }
1869: for (f = 0; f < numFaces; ++f) {
1870: const PetscInt face = star[f];
1871: PetscInt *closure = NULL;
1872: PetscInt closureSize, c;
1873: PetscInt faceLoc;
1875: DMLabelGetValue(subpointMap, face, &faceLoc);
1876: if (faceLoc == dim-1) continue;
1877: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
1878: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1879: for (c = 0; c < closureSize*2; c += 2) {
1880: const PetscInt point = closure[c];
1881: PetscInt vertexLoc;
1883: if ((point >= pStart[0]) && (point < pEnd[0])) {
1884: DMLabelGetValue(vertexLabel, point, &vertexLoc);
1885: if (vertexLoc != value) break;
1886: }
1887: }
1888: if (c == closureSize*2) {
1889: const PetscInt *support;
1890: PetscInt supportSize, s;
1892: for (c = 0; c < closureSize*2; c += 2) {
1893: const PetscInt point = closure[c];
1895: for (d = 0; d < dim; ++d) {
1896: if ((point >= pStart[d]) && (point < pEnd[d])) {
1897: DMLabelSetValue(subpointMap, point, d);
1898: break;
1899: }
1900: }
1901: }
1902: DMPlexGetSupportSize(dm, face, &supportSize);
1903: DMPlexGetSupport(dm, face, &support);
1904: for (s = 0; s < supportSize; ++s) {
1905: DMLabelSetValue(subpointMap, support[s], dim);
1906: }
1907: }
1908: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1909: }
1910: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1911: }
1912: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
1913: ISDestroy(&subvertexIS);
1914: PetscFree3(pStart,pEnd,pMax);
1915: return(0);
1916: }
1920: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
1921: {
1922: DMLabel label = NULL;
1923: const PetscInt *cone;
1924: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize;
1925: PetscErrorCode ierr;
1928: *numFaces = 0;
1929: *nFV = 0;
1930: if (labelname) {DMPlexGetLabel(dm, labelname, &label);}
1931: *subCells = NULL;
1932: DMGetDimension(dm, &dim);
1933: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1934: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1935: if (cMax < 0) return(0);
1936: if (label) {
1937: for (c = cMax; c < cEnd; ++c) {
1938: PetscInt val;
1940: DMLabelGetValue(label, c, &val);
1941: if (val == value) {
1942: ++(*numFaces);
1943: DMPlexGetConeSize(dm, c, &coneSize);
1944: }
1945: }
1946: } else {
1947: *numFaces = cEnd - cMax;
1948: DMPlexGetConeSize(dm, cMax, &coneSize);
1949: }
1950: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
1951: PetscMalloc1(*numFaces *2, subCells);
1952: for (c = cMax; c < cEnd; ++c) {
1953: const PetscInt *cells;
1954: PetscInt numCells;
1956: if (label) {
1957: PetscInt val;
1959: DMLabelGetValue(label, c, &val);
1960: if (val != value) continue;
1961: }
1962: DMPlexGetCone(dm, c, &cone);
1963: for (p = 0; p < *nFV; ++p) {
1964: DMLabelSetValue(subpointMap, cone[p], 0);
1965: }
1966: /* Negative face */
1967: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
1968: /* Not true in parallel
1969: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
1970: for (p = 0; p < numCells; ++p) {
1971: DMLabelSetValue(subpointMap, cells[p], 2);
1972: (*subCells)[subc++] = cells[p];
1973: }
1974: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
1975: /* Positive face is not included */
1976: }
1977: return(0);
1978: }
1982: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
1983: {
1984: PetscInt *pStart, *pEnd;
1985: PetscInt dim, cMax, cEnd, c, d;
1989: DMGetDimension(dm, &dim);
1990: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1991: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1992: if (cMax < 0) return(0);
1993: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
1994: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1995: for (c = cMax; c < cEnd; ++c) {
1996: const PetscInt *cone;
1997: PetscInt *closure = NULL;
1998: PetscInt fconeSize, coneSize, closureSize, cl, val;
2000: if (label) {
2001: DMLabelGetValue(label, c, &val);
2002: if (val != value) continue;
2003: }
2004: DMPlexGetConeSize(dm, c, &coneSize);
2005: DMPlexGetCone(dm, c, &cone);
2006: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2007: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2008: /* Negative face */
2009: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2010: for (cl = 0; cl < closureSize*2; cl += 2) {
2011: const PetscInt point = closure[cl];
2013: for (d = 0; d <= dim; ++d) {
2014: if ((point >= pStart[d]) && (point < pEnd[d])) {
2015: DMLabelSetValue(subpointMap, point, d);
2016: break;
2017: }
2018: }
2019: }
2020: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2021: /* Cells -- positive face is not included */
2022: for (cl = 0; cl < 1; ++cl) {
2023: const PetscInt *support;
2024: PetscInt supportSize, s;
2026: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2027: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2028: DMPlexGetSupport(dm, cone[cl], &support);
2029: for (s = 0; s < supportSize; ++s) {
2030: DMLabelSetValue(subpointMap, support[s], dim);
2031: }
2032: }
2033: }
2034: PetscFree2(pStart, pEnd);
2035: return(0);
2036: }
2040: PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2041: {
2042: MPI_Comm comm;
2043: PetscBool posOrient = PETSC_FALSE;
2044: const PetscInt debug = 0;
2045: PetscInt cellDim, faceSize, f;
2049: PetscObjectGetComm((PetscObject)dm,&comm);
2050: DMGetDimension(dm, &cellDim);
2051: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2053: if (cellDim == 1 && numCorners == 2) {
2054: /* Triangle */
2055: faceSize = numCorners-1;
2056: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2057: } else if (cellDim == 2 && numCorners == 3) {
2058: /* Triangle */
2059: faceSize = numCorners-1;
2060: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2061: } else if (cellDim == 3 && numCorners == 4) {
2062: /* Tetrahedron */
2063: faceSize = numCorners-1;
2064: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2065: } else if (cellDim == 1 && numCorners == 3) {
2066: /* Quadratic line */
2067: faceSize = 1;
2068: posOrient = PETSC_TRUE;
2069: } else if (cellDim == 2 && numCorners == 4) {
2070: /* Quads */
2071: faceSize = 2;
2072: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2073: posOrient = PETSC_TRUE;
2074: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2075: posOrient = PETSC_TRUE;
2076: } else {
2077: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2078: posOrient = PETSC_FALSE;
2079: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2080: }
2081: } else if (cellDim == 2 && numCorners == 6) {
2082: /* Quadratic triangle (I hate this) */
2083: /* Edges are determined by the first 2 vertices (corners of edges) */
2084: const PetscInt faceSizeTri = 3;
2085: PetscInt sortedIndices[3], i, iFace;
2086: PetscBool found = PETSC_FALSE;
2087: PetscInt faceVerticesTriSorted[9] = {
2088: 0, 3, 4, /* bottom */
2089: 1, 4, 5, /* right */
2090: 2, 3, 5, /* left */
2091: };
2092: PetscInt faceVerticesTri[9] = {
2093: 0, 3, 4, /* bottom */
2094: 1, 4, 5, /* right */
2095: 2, 5, 3, /* left */
2096: };
2098: faceSize = faceSizeTri;
2099: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2100: PetscSortInt(faceSizeTri, sortedIndices);
2101: for (iFace = 0; iFace < 3; ++iFace) {
2102: const PetscInt ii = iFace*faceSizeTri;
2103: PetscInt fVertex, cVertex;
2105: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2106: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2107: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2108: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2109: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2110: faceVertices[fVertex] = origVertices[cVertex];
2111: break;
2112: }
2113: }
2114: }
2115: found = PETSC_TRUE;
2116: break;
2117: }
2118: }
2119: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2120: if (posOriented) *posOriented = PETSC_TRUE;
2121: return(0);
2122: } else if (cellDim == 2 && numCorners == 9) {
2123: /* Quadratic quad (I hate this) */
2124: /* Edges are determined by the first 2 vertices (corners of edges) */
2125: const PetscInt faceSizeQuad = 3;
2126: PetscInt sortedIndices[3], i, iFace;
2127: PetscBool found = PETSC_FALSE;
2128: PetscInt faceVerticesQuadSorted[12] = {
2129: 0, 1, 4, /* bottom */
2130: 1, 2, 5, /* right */
2131: 2, 3, 6, /* top */
2132: 0, 3, 7, /* left */
2133: };
2134: PetscInt faceVerticesQuad[12] = {
2135: 0, 1, 4, /* bottom */
2136: 1, 2, 5, /* right */
2137: 2, 3, 6, /* top */
2138: 3, 0, 7, /* left */
2139: };
2141: faceSize = faceSizeQuad;
2142: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2143: PetscSortInt(faceSizeQuad, sortedIndices);
2144: for (iFace = 0; iFace < 4; ++iFace) {
2145: const PetscInt ii = iFace*faceSizeQuad;
2146: PetscInt fVertex, cVertex;
2148: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2149: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2150: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2151: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2152: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2153: faceVertices[fVertex] = origVertices[cVertex];
2154: break;
2155: }
2156: }
2157: }
2158: found = PETSC_TRUE;
2159: break;
2160: }
2161: }
2162: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2163: if (posOriented) *posOriented = PETSC_TRUE;
2164: return(0);
2165: } else if (cellDim == 3 && numCorners == 8) {
2166: /* Hexes
2167: A hex is two oriented quads with the normal of the first
2168: pointing up at the second.
2170: 7---6
2171: /| /|
2172: 4---5 |
2173: | 1-|-2
2174: |/ |/
2175: 0---3
2177: Faces are determined by the first 4 vertices (corners of faces) */
2178: const PetscInt faceSizeHex = 4;
2179: PetscInt sortedIndices[4], i, iFace;
2180: PetscBool found = PETSC_FALSE;
2181: PetscInt faceVerticesHexSorted[24] = {
2182: 0, 1, 2, 3, /* bottom */
2183: 4, 5, 6, 7, /* top */
2184: 0, 3, 4, 5, /* front */
2185: 2, 3, 5, 6, /* right */
2186: 1, 2, 6, 7, /* back */
2187: 0, 1, 4, 7, /* left */
2188: };
2189: PetscInt faceVerticesHex[24] = {
2190: 1, 2, 3, 0, /* bottom */
2191: 4, 5, 6, 7, /* top */
2192: 0, 3, 5, 4, /* front */
2193: 3, 2, 6, 5, /* right */
2194: 2, 1, 7, 6, /* back */
2195: 1, 0, 4, 7, /* left */
2196: };
2198: faceSize = faceSizeHex;
2199: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2200: PetscSortInt(faceSizeHex, sortedIndices);
2201: for (iFace = 0; iFace < 6; ++iFace) {
2202: const PetscInt ii = iFace*faceSizeHex;
2203: PetscInt fVertex, cVertex;
2205: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2206: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2207: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2208: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2209: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2210: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2211: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2212: faceVertices[fVertex] = origVertices[cVertex];
2213: break;
2214: }
2215: }
2216: }
2217: found = PETSC_TRUE;
2218: break;
2219: }
2220: }
2221: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2222: if (posOriented) *posOriented = PETSC_TRUE;
2223: return(0);
2224: } else if (cellDim == 3 && numCorners == 10) {
2225: /* Quadratic tet */
2226: /* Faces are determined by the first 3 vertices (corners of faces) */
2227: const PetscInt faceSizeTet = 6;
2228: PetscInt sortedIndices[6], i, iFace;
2229: PetscBool found = PETSC_FALSE;
2230: PetscInt faceVerticesTetSorted[24] = {
2231: 0, 1, 2, 6, 7, 8, /* bottom */
2232: 0, 3, 4, 6, 7, 9, /* front */
2233: 1, 4, 5, 7, 8, 9, /* right */
2234: 2, 3, 5, 6, 8, 9, /* left */
2235: };
2236: PetscInt faceVerticesTet[24] = {
2237: 0, 1, 2, 6, 7, 8, /* bottom */
2238: 0, 4, 3, 6, 7, 9, /* front */
2239: 1, 5, 4, 7, 8, 9, /* right */
2240: 2, 3, 5, 8, 6, 9, /* left */
2241: };
2243: faceSize = faceSizeTet;
2244: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2245: PetscSortInt(faceSizeTet, sortedIndices);
2246: for (iFace=0; iFace < 4; ++iFace) {
2247: const PetscInt ii = iFace*faceSizeTet;
2248: PetscInt fVertex, cVertex;
2250: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2251: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2252: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2253: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2254: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2255: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2256: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2257: faceVertices[fVertex] = origVertices[cVertex];
2258: break;
2259: }
2260: }
2261: }
2262: found = PETSC_TRUE;
2263: break;
2264: }
2265: }
2266: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2267: if (posOriented) *posOriented = PETSC_TRUE;
2268: return(0);
2269: } else if (cellDim == 3 && numCorners == 27) {
2270: /* Quadratic hexes (I hate this)
2271: A hex is two oriented quads with the normal of the first
2272: pointing up at the second.
2274: 7---6
2275: /| /|
2276: 4---5 |
2277: | 3-|-2
2278: |/ |/
2279: 0---1
2281: Faces are determined by the first 4 vertices (corners of faces) */
2282: const PetscInt faceSizeQuadHex = 9;
2283: PetscInt sortedIndices[9], i, iFace;
2284: PetscBool found = PETSC_FALSE;
2285: PetscInt faceVerticesQuadHexSorted[54] = {
2286: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2287: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2288: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2289: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2290: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2291: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2292: };
2293: PetscInt faceVerticesQuadHex[54] = {
2294: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2295: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2296: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2297: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2298: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2299: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2300: };
2302: faceSize = faceSizeQuadHex;
2303: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2304: PetscSortInt(faceSizeQuadHex, sortedIndices);
2305: for (iFace = 0; iFace < 6; ++iFace) {
2306: const PetscInt ii = iFace*faceSizeQuadHex;
2307: PetscInt fVertex, cVertex;
2309: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2310: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2311: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2312: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2313: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2314: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2315: if (indices[cVertex] == faceVerticesQuadHex[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 hex crossface");
2326: if (posOriented) *posOriented = PETSC_TRUE;
2327: return(0);
2328: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2329: if (!posOrient) {
2330: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2331: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2332: } else {
2333: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2334: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2335: }
2336: if (posOriented) *posOriented = posOrient;
2337: return(0);
2338: }
2342: /*
2343: Given a cell and a face, as a set of vertices,
2344: return the oriented face, as a set of vertices, in faceVertices
2345: The orientation is such that the face normal points out of the cell
2346: */
2347: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2348: {
2349: const PetscInt *cone = NULL;
2350: PetscInt coneSize, v, f, v2;
2351: PetscInt oppositeVertex = -1;
2352: PetscErrorCode ierr;
2355: DMPlexGetConeSize(dm, cell, &coneSize);
2356: DMPlexGetCone(dm, cell, &cone);
2357: for (v = 0, v2 = 0; v < coneSize; ++v) {
2358: PetscBool found = PETSC_FALSE;
2360: for (f = 0; f < faceSize; ++f) {
2361: if (face[f] == cone[v]) {
2362: found = PETSC_TRUE; break;
2363: }
2364: }
2365: if (found) {
2366: indices[v2] = v;
2367: origVertices[v2] = cone[v];
2368: ++v2;
2369: } else {
2370: oppositeVertex = v;
2371: }
2372: }
2373: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2374: return(0);
2375: }
2379: /*
2380: DMPlexInsertFace_Internal - Puts a face into the mesh
2382: Not collective
2384: Input Parameters:
2385: + dm - The DMPlex
2386: . numFaceVertex - The number of vertices in the face
2387: . faceVertices - The vertices in the face for dm
2388: . subfaceVertices - The vertices in the face for subdm
2389: . numCorners - The number of vertices in the cell
2390: . cell - A cell in dm containing the face
2391: . subcell - A cell in subdm containing the face
2392: . firstFace - First face in the mesh
2393: - newFacePoint - Next face in the mesh
2395: Output Parameters:
2396: . newFacePoint - Contains next face point number on input, updated on output
2398: Level: developer
2399: */
2400: 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)
2401: {
2402: MPI_Comm comm;
2403: DM_Plex *submesh = (DM_Plex*) subdm->data;
2404: const PetscInt *faces;
2405: PetscInt numFaces, coneSize;
2406: PetscErrorCode ierr;
2409: PetscObjectGetComm((PetscObject)dm,&comm);
2410: DMPlexGetConeSize(subdm, subcell, &coneSize);
2411: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2412: #if 0
2413: /* Cannot use this because support() has not been constructed yet */
2414: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2415: #else
2416: {
2417: PetscInt f;
2419: numFaces = 0;
2420: DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2421: for (f = firstFace; f < *newFacePoint; ++f) {
2422: PetscInt dof, off, d;
2424: PetscSectionGetDof(submesh->coneSection, f, &dof);
2425: PetscSectionGetOffset(submesh->coneSection, f, &off);
2426: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2427: for (d = 0; d < dof; ++d) {
2428: const PetscInt p = submesh->cones[off+d];
2429: PetscInt v;
2431: for (v = 0; v < numFaceVertices; ++v) {
2432: if (subfaceVertices[v] == p) break;
2433: }
2434: if (v == numFaceVertices) break;
2435: }
2436: if (d == dof) {
2437: numFaces = 1;
2438: ((PetscInt*) faces)[0] = f;
2439: }
2440: }
2441: }
2442: #endif
2443: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2444: else if (numFaces == 1) {
2445: /* Add the other cell neighbor for this face */
2446: DMPlexSetCone(subdm, subcell, faces);
2447: } else {
2448: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2449: PetscBool posOriented;
2451: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2452: origVertices = &orientedVertices[numFaceVertices];
2453: indices = &orientedVertices[numFaceVertices*2];
2454: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2455: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2456: /* TODO: I know that routine should return a permutation, not the indices */
2457: for (v = 0; v < numFaceVertices; ++v) {
2458: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2459: for (ov = 0; ov < numFaceVertices; ++ov) {
2460: if (orientedVertices[ov] == vertex) {
2461: orientedSubVertices[ov] = subvertex;
2462: break;
2463: }
2464: }
2465: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2466: }
2467: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2468: DMPlexSetCone(subdm, subcell, newFacePoint);
2469: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2470: ++(*newFacePoint);
2471: }
2472: #if 0
2473: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2474: #else
2475: DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2476: #endif
2477: return(0);
2478: }
2482: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2483: {
2484: MPI_Comm comm;
2485: DMLabel subpointMap;
2486: IS subvertexIS, subcellIS;
2487: const PetscInt *subVertices, *subCells;
2488: PetscInt numSubVertices, firstSubVertex, numSubCells;
2489: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2490: PetscInt vStart, vEnd, c, f;
2491: PetscErrorCode ierr;
2494: PetscObjectGetComm((PetscObject)dm,&comm);
2495: /* Create subpointMap which marks the submesh */
2496: DMLabelCreate("subpoint_map", &subpointMap);
2497: DMPlexSetSubpointMap(subdm, subpointMap);
2498: DMLabelDestroy(&subpointMap);
2499: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2500: /* Setup chart */
2501: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2502: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2503: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2504: DMPlexSetVTKCellHeight(subdm, 1);
2505: /* Set cone sizes */
2506: firstSubVertex = numSubCells;
2507: firstSubFace = numSubCells+numSubVertices;
2508: newFacePoint = firstSubFace;
2509: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2510: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2511: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2512: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2513: for (c = 0; c < numSubCells; ++c) {
2514: DMPlexSetConeSize(subdm, c, 1);
2515: }
2516: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2517: DMPlexSetConeSize(subdm, f, nFV);
2518: }
2519: DMSetUp(subdm);
2520: /* Create face cones */
2521: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2522: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2523: DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2524: for (c = 0; c < numSubCells; ++c) {
2525: const PetscInt cell = subCells[c];
2526: const PetscInt subcell = c;
2527: PetscInt *closure = NULL;
2528: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2530: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2531: for (cl = 0; cl < closureSize*2; cl += 2) {
2532: const PetscInt point = closure[cl];
2533: PetscInt subVertex;
2535: if ((point >= vStart) && (point < vEnd)) {
2536: ++numCorners;
2537: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2538: if (subVertex >= 0) {
2539: closure[faceSize] = point;
2540: subface[faceSize] = firstSubVertex+subVertex;
2541: ++faceSize;
2542: }
2543: }
2544: }
2545: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2546: if (faceSize == nFV) {
2547: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2548: }
2549: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2550: }
2551: DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2552: DMPlexSymmetrize(subdm);
2553: DMPlexStratify(subdm);
2554: /* Build coordinates */
2555: {
2556: PetscSection coordSection, subCoordSection;
2557: Vec coordinates, subCoordinates;
2558: PetscScalar *coords, *subCoords;
2559: PetscInt numComp, coordSize, v;
2561: DMGetCoordinateSection(dm, &coordSection);
2562: DMGetCoordinatesLocal(dm, &coordinates);
2563: DMGetCoordinateSection(subdm, &subCoordSection);
2564: PetscSectionSetNumFields(subCoordSection, 1);
2565: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2566: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2567: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2568: for (v = 0; v < numSubVertices; ++v) {
2569: const PetscInt vertex = subVertices[v];
2570: const PetscInt subvertex = firstSubVertex+v;
2571: PetscInt dof;
2573: PetscSectionGetDof(coordSection, vertex, &dof);
2574: PetscSectionSetDof(subCoordSection, subvertex, dof);
2575: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2576: }
2577: PetscSectionSetUp(subCoordSection);
2578: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2579: VecCreate(comm, &subCoordinates);
2580: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2581: VecSetType(subCoordinates,VECSTANDARD);
2582: if (coordSize) {
2583: VecGetArray(coordinates, &coords);
2584: VecGetArray(subCoordinates, &subCoords);
2585: for (v = 0; v < numSubVertices; ++v) {
2586: const PetscInt vertex = subVertices[v];
2587: const PetscInt subvertex = firstSubVertex+v;
2588: PetscInt dof, off, sdof, soff, d;
2590: PetscSectionGetDof(coordSection, vertex, &dof);
2591: PetscSectionGetOffset(coordSection, vertex, &off);
2592: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2593: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2594: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2595: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2596: }
2597: VecRestoreArray(coordinates, &coords);
2598: VecRestoreArray(subCoordinates, &subCoords);
2599: }
2600: DMSetCoordinatesLocal(subdm, subCoordinates);
2601: VecDestroy(&subCoordinates);
2602: }
2603: /* Cleanup */
2604: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2605: ISDestroy(&subvertexIS);
2606: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2607: ISDestroy(&subcellIS);
2608: return(0);
2609: }
2613: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2614: {
2615: PetscInt subPoint;
2618: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2619: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2620: }
2624: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, DM subdm)
2625: {
2626: MPI_Comm comm;
2627: DMLabel subpointMap;
2628: IS *subpointIS;
2629: const PetscInt **subpoints;
2630: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2631: PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v;
2632: PetscErrorCode ierr;
2635: PetscObjectGetComm((PetscObject)dm,&comm);
2636: /* Create subpointMap which marks the submesh */
2637: DMLabelCreate("subpoint_map", &subpointMap);
2638: DMPlexSetSubpointMap(subdm, subpointMap);
2639: DMLabelDestroy(&subpointMap);
2640: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2641: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2642: /* Setup chart */
2643: DMGetDimension(dm, &dim);
2644: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2645: for (d = 0; d <= dim; ++d) {
2646: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2647: totSubPoints += numSubPoints[d];
2648: }
2649: DMPlexSetChart(subdm, 0, totSubPoints);
2650: DMPlexSetVTKCellHeight(subdm, 1);
2651: /* Set cone sizes */
2652: firstSubPoint[dim] = 0;
2653: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
2654: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
2655: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2656: for (d = 0; d <= dim; ++d) {
2657: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2658: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2659: }
2660: for (d = 0; d <= dim; ++d) {
2661: for (p = 0; p < numSubPoints[d]; ++p) {
2662: const PetscInt point = subpoints[d][p];
2663: const PetscInt subpoint = firstSubPoint[d] + p;
2664: const PetscInt *cone;
2665: PetscInt coneSize, coneSizeNew, c, val;
2667: DMPlexGetConeSize(dm, point, &coneSize);
2668: DMPlexSetConeSize(subdm, subpoint, coneSize);
2669: if (d == dim) {
2670: DMPlexGetCone(dm, point, &cone);
2671: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2672: DMLabelGetValue(subpointMap, cone[c], &val);
2673: if (val >= 0) coneSizeNew++;
2674: }
2675: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2676: }
2677: }
2678: }
2679: DMSetUp(subdm);
2680: /* Set cones */
2681: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2682: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2683: for (d = 0; d <= dim; ++d) {
2684: for (p = 0; p < numSubPoints[d]; ++p) {
2685: const PetscInt point = subpoints[d][p];
2686: const PetscInt subpoint = firstSubPoint[d] + p;
2687: const PetscInt *cone, *ornt;
2688: PetscInt coneSize, subconeSize, coneSizeNew, c, subc;
2690: DMPlexGetConeSize(dm, point, &coneSize);
2691: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2692: DMPlexGetCone(dm, point, &cone);
2693: DMPlexGetConeOrientation(dm, point, &ornt);
2694: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2695: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2696: if (subc >= 0) {
2697: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
2698: orntNew[coneSizeNew] = ornt[c];
2699: ++coneSizeNew;
2700: }
2701: }
2702: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2703: DMPlexSetCone(subdm, subpoint, coneNew);
2704: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2705: }
2706: }
2707: PetscFree2(coneNew,orntNew);
2708: DMPlexSymmetrize(subdm);
2709: DMPlexStratify(subdm);
2710: /* Build coordinates */
2711: {
2712: PetscSection coordSection, subCoordSection;
2713: Vec coordinates, subCoordinates;
2714: PetscScalar *coords, *subCoords;
2715: PetscInt numComp, coordSize;
2717: DMGetCoordinateSection(dm, &coordSection);
2718: DMGetCoordinatesLocal(dm, &coordinates);
2719: DMGetCoordinateSection(subdm, &subCoordSection);
2720: PetscSectionSetNumFields(subCoordSection, 1);
2721: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2722: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2723: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2724: for (v = 0; v < numSubPoints[0]; ++v) {
2725: const PetscInt vertex = subpoints[0][v];
2726: const PetscInt subvertex = firstSubPoint[0]+v;
2727: PetscInt dof;
2729: PetscSectionGetDof(coordSection, vertex, &dof);
2730: PetscSectionSetDof(subCoordSection, subvertex, dof);
2731: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2732: }
2733: PetscSectionSetUp(subCoordSection);
2734: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2735: VecCreate(comm, &subCoordinates);
2736: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2737: VecSetType(subCoordinates,VECSTANDARD);
2738: VecGetArray(coordinates, &coords);
2739: VecGetArray(subCoordinates, &subCoords);
2740: for (v = 0; v < numSubPoints[0]; ++v) {
2741: const PetscInt vertex = subpoints[0][v];
2742: const PetscInt subvertex = firstSubPoint[0]+v;
2743: PetscInt dof, off, sdof, soff, d;
2745: PetscSectionGetDof(coordSection, vertex, &dof);
2746: PetscSectionGetOffset(coordSection, vertex, &off);
2747: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2748: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2749: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2750: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2751: }
2752: VecRestoreArray(coordinates, &coords);
2753: VecRestoreArray(subCoordinates, &subCoords);
2754: DMSetCoordinatesLocal(subdm, subCoordinates);
2755: VecDestroy(&subCoordinates);
2756: }
2757: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
2758: {
2759: PetscSF sfPoint, sfPointSub;
2760: IS subpIS;
2761: const PetscSFNode *remotePoints;
2762: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
2763: const PetscInt *localPoints, *subpoints;
2764: PetscInt *slocalPoints;
2765: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
2766: PetscMPIInt rank;
2768: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2769: DMGetPointSF(dm, &sfPoint);
2770: DMGetPointSF(subdm, &sfPointSub);
2771: DMPlexGetChart(dm, &pStart, &pEnd);
2772: DMPlexGetChart(subdm, NULL, &numSubroots);
2773: DMPlexCreateSubpointIS(subdm, &subpIS);
2774: if (subpIS) {
2775: ISGetIndices(subpIS, &subpoints);
2776: ISGetLocalSize(subpIS, &numSubpoints);
2777: }
2778: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2779: if (numRoots >= 0) {
2780: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
2781: for (p = 0; p < pEnd-pStart; ++p) {
2782: newLocalPoints[p].rank = -2;
2783: newLocalPoints[p].index = -2;
2784: }
2785: /* Set subleaves */
2786: for (l = 0; l < numLeaves; ++l) {
2787: const PetscInt point = localPoints[l];
2788: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
2790: if (subpoint < 0) continue;
2791: newLocalPoints[point-pStart].rank = rank;
2792: newLocalPoints[point-pStart].index = subpoint;
2793: ++numSubleaves;
2794: }
2795: /* Must put in owned subpoints */
2796: for (p = pStart; p < pEnd; ++p) {
2797: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
2799: if (subpoint < 0) {
2800: newOwners[p-pStart].rank = -3;
2801: newOwners[p-pStart].index = -3;
2802: } else {
2803: newOwners[p-pStart].rank = rank;
2804: newOwners[p-pStart].index = subpoint;
2805: }
2806: }
2807: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2808: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2809: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2810: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2811: PetscMalloc1(numSubleaves, &slocalPoints);
2812: PetscMalloc1(numSubleaves, &sremotePoints);
2813: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
2814: const PetscInt point = localPoints[l];
2815: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
2817: if (subpoint < 0) continue;
2818: if (newLocalPoints[point].rank == rank) {++ll; continue;}
2819: slocalPoints[sl] = subpoint;
2820: sremotePoints[sl].rank = newLocalPoints[point].rank;
2821: sremotePoints[sl].index = newLocalPoints[point].index;
2822: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
2823: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
2824: ++sl;
2825: }
2826: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
2827: PetscFree2(newLocalPoints,newOwners);
2828: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
2829: }
2830: if (subpIS) {
2831: ISRestoreIndices(subpIS, &subpoints);
2832: ISDestroy(&subpIS);
2833: }
2834: }
2835: /* Cleanup */
2836: for (d = 0; d <= dim; ++d) {
2837: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
2838: ISDestroy(&subpointIS[d]);
2839: }
2840: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
2841: return(0);
2842: }
2846: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2847: {
2851: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, subdm);
2852: return(0);
2853: }
2857: /*@
2858: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
2860: Input Parameters:
2861: + dm - The original mesh
2862: . vertexLabel - The DMLabel marking vertices contained in the surface
2863: - value - The label value to use
2865: Output Parameter:
2866: . subdm - The surface mesh
2868: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
2870: Level: developer
2872: .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue()
2873: @*/
2874: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
2875: {
2876: PetscInt dim, depth;
2882: DMGetDimension(dm, &dim);
2883: DMPlexGetDepth(dm, &depth);
2884: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
2885: DMSetType(*subdm, DMPLEX);
2886: DMSetDimension(*subdm, dim-1);
2887: if (depth == dim) {
2888: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
2889: } else {
2890: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
2891: }
2892: return(0);
2893: }
2897: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
2898: {
2899: MPI_Comm comm;
2900: DMLabel subpointMap;
2901: IS subvertexIS;
2902: const PetscInt *subVertices;
2903: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
2904: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
2905: PetscInt cMax, c, f;
2906: PetscErrorCode ierr;
2909: PetscObjectGetComm((PetscObject)dm, &comm);
2910: /* Create subpointMap which marks the submesh */
2911: DMLabelCreate("subpoint_map", &subpointMap);
2912: DMPlexSetSubpointMap(subdm, subpointMap);
2913: DMLabelDestroy(&subpointMap);
2914: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
2915: /* Setup chart */
2916: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2917: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2918: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2919: DMPlexSetVTKCellHeight(subdm, 1);
2920: /* Set cone sizes */
2921: firstSubVertex = numSubCells;
2922: firstSubFace = numSubCells+numSubVertices;
2923: newFacePoint = firstSubFace;
2924: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2925: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2926: for (c = 0; c < numSubCells; ++c) {
2927: DMPlexSetConeSize(subdm, c, 1);
2928: }
2929: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2930: DMPlexSetConeSize(subdm, f, nFV);
2931: }
2932: DMSetUp(subdm);
2933: /* Create face cones */
2934: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2935: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2936: DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2937: for (c = 0; c < numSubCells; ++c) {
2938: const PetscInt cell = subCells[c];
2939: const PetscInt subcell = c;
2940: const PetscInt *cone, *cells;
2941: PetscInt numCells, subVertex, p, v;
2943: if (cell < cMax) continue;
2944: DMPlexGetCone(dm, cell, &cone);
2945: for (v = 0; v < nFV; ++v) {
2946: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
2947: subface[v] = firstSubVertex+subVertex;
2948: }
2949: DMPlexSetCone(subdm, newFacePoint, subface);
2950: DMPlexSetCone(subdm, subcell, &newFacePoint);
2951: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
2952: /* Not true in parallel
2953: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2954: for (p = 0; p < numCells; ++p) {
2955: PetscInt negsubcell;
2957: if (cells[p] >= cMax) continue;
2958: /* I know this is a crap search */
2959: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
2960: if (subCells[negsubcell] == cells[p]) break;
2961: }
2962: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
2963: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
2964: }
2965: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
2966: ++newFacePoint;
2967: }
2968: DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2969: DMPlexSymmetrize(subdm);
2970: DMPlexStratify(subdm);
2971: /* Build coordinates */
2972: {
2973: PetscSection coordSection, subCoordSection;
2974: Vec coordinates, subCoordinates;
2975: PetscScalar *coords, *subCoords;
2976: PetscInt numComp, coordSize, v;
2978: DMGetCoordinateSection(dm, &coordSection);
2979: DMGetCoordinatesLocal(dm, &coordinates);
2980: DMGetCoordinateSection(subdm, &subCoordSection);
2981: PetscSectionSetNumFields(subCoordSection, 1);
2982: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2983: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2984: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2985: for (v = 0; v < numSubVertices; ++v) {
2986: const PetscInt vertex = subVertices[v];
2987: const PetscInt subvertex = firstSubVertex+v;
2988: PetscInt dof;
2990: PetscSectionGetDof(coordSection, vertex, &dof);
2991: PetscSectionSetDof(subCoordSection, subvertex, dof);
2992: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2993: }
2994: PetscSectionSetUp(subCoordSection);
2995: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2996: VecCreate(comm, &subCoordinates);
2997: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2998: VecSetType(subCoordinates,VECSTANDARD);
2999: VecGetArray(coordinates, &coords);
3000: VecGetArray(subCoordinates, &subCoords);
3001: for (v = 0; v < numSubVertices; ++v) {
3002: const PetscInt vertex = subVertices[v];
3003: const PetscInt subvertex = firstSubVertex+v;
3004: PetscInt dof, off, sdof, soff, d;
3006: PetscSectionGetDof(coordSection, vertex, &dof);
3007: PetscSectionGetOffset(coordSection, vertex, &off);
3008: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3009: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3010: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3011: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3012: }
3013: VecRestoreArray(coordinates, &coords);
3014: VecRestoreArray(subCoordinates, &subCoords);
3015: DMSetCoordinatesLocal(subdm, subCoordinates);
3016: VecDestroy(&subCoordinates);
3017: }
3018: /* Build SF */
3019: CHKMEMQ;
3020: {
3021: PetscSF sfPoint, sfPointSub;
3022: const PetscSFNode *remotePoints;
3023: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3024: const PetscInt *localPoints;
3025: PetscInt *slocalPoints;
3026: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3027: PetscMPIInt rank;
3029: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3030: DMGetPointSF(dm, &sfPoint);
3031: DMGetPointSF(subdm, &sfPointSub);
3032: DMPlexGetChart(dm, &pStart, &pEnd);
3033: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3034: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3035: if (numRoots >= 0) {
3036: /* Only vertices should be shared */
3037: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3038: for (p = 0; p < pEnd-pStart; ++p) {
3039: newLocalPoints[p].rank = -2;
3040: newLocalPoints[p].index = -2;
3041: }
3042: /* Set subleaves */
3043: for (l = 0; l < numLeaves; ++l) {
3044: const PetscInt point = localPoints[l];
3045: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3047: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3048: if (subPoint < 0) continue;
3049: newLocalPoints[point-pStart].rank = rank;
3050: newLocalPoints[point-pStart].index = subPoint;
3051: ++numSubLeaves;
3052: }
3053: /* Must put in owned subpoints */
3054: for (p = pStart; p < pEnd; ++p) {
3055: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3057: if (subPoint < 0) {
3058: newOwners[p-pStart].rank = -3;
3059: newOwners[p-pStart].index = -3;
3060: } else {
3061: newOwners[p-pStart].rank = rank;
3062: newOwners[p-pStart].index = subPoint;
3063: }
3064: }
3065: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3066: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3067: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3068: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3069: PetscMalloc1(numSubLeaves, &slocalPoints);
3070: PetscMalloc1(numSubLeaves, &sremotePoints);
3071: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3072: const PetscInt point = localPoints[l];
3073: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3075: if (subPoint < 0) continue;
3076: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3077: slocalPoints[sl] = subPoint;
3078: sremotePoints[sl].rank = newLocalPoints[point].rank;
3079: sremotePoints[sl].index = newLocalPoints[point].index;
3080: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3081: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3082: ++sl;
3083: }
3084: PetscFree2(newLocalPoints,newOwners);
3085: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3086: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3087: }
3088: }
3089: CHKMEMQ;
3090: /* Cleanup */
3091: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3092: ISDestroy(&subvertexIS);
3093: PetscFree(subCells);
3094: return(0);
3095: }
3099: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3100: {
3101: DMLabel label = NULL;
3105: if (labelname) {DMPlexGetLabel(dm, labelname, &label);}
3106: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, subdm);
3107: return(0);
3108: }
3112: /*
3113: 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.
3115: Input Parameters:
3116: + dm - The original mesh
3117: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3118: . label - A label name, or NULL
3119: - value - A label value
3121: Output Parameter:
3122: . subdm - The surface mesh
3124: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3126: Level: developer
3128: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3129: */
3130: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3131: {
3132: PetscInt dim, depth;
3138: DMGetDimension(dm, &dim);
3139: DMPlexGetDepth(dm, &depth);
3140: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3141: DMSetType(*subdm, DMPLEX);
3142: DMSetDimension(*subdm, dim-1);
3143: if (depth == dim) {
3144: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3145: } else {
3146: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3147: }
3148: return(0);
3149: }
3153: /*@
3154: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3156: Input Parameter:
3157: . dm - The submesh DM
3159: Output Parameter:
3160: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3162: Level: developer
3164: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3165: @*/
3166: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3167: {
3168: DM_Plex *mesh = (DM_Plex*) dm->data;
3173: *subpointMap = mesh->subpointMap;
3174: return(0);
3175: }
3179: /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */
3180: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3181: {
3182: DM_Plex *mesh = (DM_Plex *) dm->data;
3183: DMLabel tmp;
3188: tmp = mesh->subpointMap;
3189: mesh->subpointMap = subpointMap;
3190: ++mesh->subpointMap->refct;
3191: DMLabelDestroy(&tmp);
3192: return(0);
3193: }
3197: /*@
3198: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
3200: Input Parameter:
3201: . dm - The submesh DM
3203: Output Parameter:
3204: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3206: Note: This IS is guaranteed to be sorted by the construction of the submesh
3208: Level: developer
3210: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3211: @*/
3212: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3213: {
3214: MPI_Comm comm;
3215: DMLabel subpointMap;
3216: IS is;
3217: const PetscInt *opoints;
3218: PetscInt *points, *depths;
3219: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3220: PetscErrorCode ierr;
3225: PetscObjectGetComm((PetscObject)dm,&comm);
3226: *subpointIS = NULL;
3227: DMPlexGetSubpointMap(dm, &subpointMap);
3228: DMPlexGetDepth(dm, &depth);
3229: if (subpointMap && depth >= 0) {
3230: DMPlexGetChart(dm, &pStart, &pEnd);
3231: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3232: DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
3233: depths[0] = depth;
3234: depths[1] = 0;
3235: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3236: PetscMalloc1(pEnd, &points);
3237: for(d = 0, off = 0; d <= depth; ++d) {
3238: const PetscInt dep = depths[d];
3240: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3241: DMLabelGetStratumSize(subpointMap, dep, &n);
3242: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3243: 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);
3244: } else {
3245: if (!n) {
3246: if (d == 0) {
3247: /* Missing cells */
3248: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3249: } else {
3250: /* Missing faces */
3251: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3252: }
3253: }
3254: }
3255: if (n) {
3256: DMLabelGetStratumIS(subpointMap, dep, &is);
3257: ISGetIndices(is, &opoints);
3258: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3259: ISRestoreIndices(is, &opoints);
3260: ISDestroy(&is);
3261: }
3262: }
3263: DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
3264: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3265: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3266: }
3267: return(0);
3268: }