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