Actual source code: plexsubmesh.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petscsf.h>
6: /*@
7: DMPlexMarkBoundaryFaces - Mark all faces on the boundary
9: Not Collective
11: Input Parameter:
12: . dm - The original DM
14: Output Parameter:
15: . label - The DMLabel marking boundary faces with value 1
17: Level: developer
19: .seealso: DMLabelCreate(), DMPlexCreateLabel()
20: @*/
21: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
22: {
23: PetscInt fStart, fEnd, f;
28: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
29: for (f = fStart; f < fEnd; ++f) {
30: PetscInt supportSize;
32: DMPlexGetSupportSize(dm, f, &supportSize);
33: if (supportSize == 1) {
34: DMLabelSetValue(label, f, 1);
35: }
36: }
37: return(0);
38: }
42: /*@
43: DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface
45: Input Parameters:
46: + dm - The DM
47: - label - A DMLabel marking the surface points
49: Output Parameter:
50: . label - A DMLabel marking all surface points in the transitive closure
52: Level: developer
54: .seealso: DMPlexLabelCohesiveComplete()
55: @*/
56: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
57: {
58: IS valueIS;
59: const PetscInt *values;
60: PetscInt numValues, v;
61: PetscErrorCode ierr;
64: DMLabelGetNumValues(label, &numValues);
65: DMLabelGetValueIS(label, &valueIS);
66: ISGetIndices(valueIS, &values);
67: for (v = 0; v < numValues; ++v) {
68: IS pointIS;
69: const PetscInt *points;
70: PetscInt numPoints, p;
72: DMLabelGetStratumSize(label, values[v], &numPoints);
73: DMLabelGetStratumIS(label, values[v], &pointIS);
74: ISGetIndices(pointIS, &points);
75: for (p = 0; p < numPoints; ++p) {
76: PetscInt *closure = NULL;
77: PetscInt closureSize, c;
79: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
80: for (c = 0; c < closureSize*2; c += 2) {
81: DMLabelSetValue(label, closure[c], values[v]);
82: }
83: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
84: }
85: ISRestoreIndices(pointIS, &points);
86: ISDestroy(&pointIS);
87: }
88: ISRestoreIndices(valueIS, &values);
89: ISDestroy(&valueIS);
90: return(0);
91: }
95: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthEnd[], PetscInt depthShift[])
96: {
97: if (depth < 0) return p;
98: /* Cells */ if (p < depthEnd[depth]) return p;
99: /* Vertices */ if (p < depthEnd[0]) return p + depthShift[depth];
100: /* Faces */ if (p < depthEnd[depth-1]) return p + depthShift[depth] + depthShift[0];
101: /* Edges */ return p + depthShift[depth] + depthShift[0] + depthShift[depth-1];
102: }
106: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
107: {
108: PetscInt *depthEnd;
109: PetscInt depth = 0, d, pStart, pEnd, p;
113: DMPlexGetDepth(dm, &depth);
114: if (depth < 0) return(0);
115: PetscMalloc((depth+1) * sizeof(PetscInt), &depthEnd);
116: /* Step 1: Expand chart */
117: DMPlexGetChart(dm, &pStart, &pEnd);
118: for (d = 0; d <= depth; ++d) {
119: pEnd += depthShift[d];
120: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
121: }
122: DMPlexSetChart(dmNew, pStart, pEnd);
123: /* Step 2: Set cone and support sizes */
124: for (d = 0; d <= depth; ++d) {
125: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
126: for (p = pStart; p < pEnd; ++p) {
127: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthEnd, depthShift);
128: PetscInt size;
130: DMPlexGetConeSize(dm, p, &size);
131: DMPlexSetConeSize(dmNew, newp, size);
132: DMPlexGetSupportSize(dm, p, &size);
133: DMPlexSetSupportSize(dmNew, newp, size);
134: }
135: }
136: PetscFree(depthEnd);
137: return(0);
138: }
142: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
143: {
144: PetscInt *depthEnd, *newpoints;
145: PetscInt depth = 0, d, maxConeSize, maxSupportSize, pStart, pEnd, p;
149: DMPlexGetDepth(dm, &depth);
150: if (depth < 0) return(0);
151: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
152: PetscMalloc2(depth+1,PetscInt,&depthEnd,PetscMax(maxConeSize, maxSupportSize),PetscInt,&newpoints);
153: for (d = 0; d <= depth; ++d) {
154: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
155: }
156: /* Step 5: Set cones and supports */
157: DMPlexGetChart(dm, &pStart, &pEnd);
158: for (p = pStart; p < pEnd; ++p) {
159: const PetscInt *points = NULL, *orientations = NULL;
160: PetscInt size, i, newp = DMPlexShiftPoint_Internal(p, depth, depthEnd, depthShift);
162: DMPlexGetConeSize(dm, p, &size);
163: DMPlexGetCone(dm, p, &points);
164: DMPlexGetConeOrientation(dm, p, &orientations);
165: for (i = 0; i < size; ++i) {
166: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthEnd, depthShift);
167: }
168: DMPlexSetCone(dmNew, newp, newpoints);
169: DMPlexSetConeOrientation(dmNew, newp, orientations);
170: DMPlexGetSupportSize(dm, p, &size);
171: DMPlexGetSupport(dm, p, &points);
172: for (i = 0; i < size; ++i) {
173: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthEnd, depthShift);
174: }
175: DMPlexSetSupport(dmNew, newp, newpoints);
176: }
177: PetscFree2(depthEnd,newpoints);
178: return(0);
179: }
183: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
184: {
185: PetscSection coordSection, newCoordSection;
186: Vec coordinates, newCoordinates;
187: PetscScalar *coords, *newCoords;
188: PetscInt *depthEnd, coordSize;
189: PetscInt dim, depth = 0, d, vStart, vEnd, vStartNew, vEndNew, v;
193: DMPlexGetDimension(dm, &dim);
194: DMPlexGetDepth(dm, &depth);
195: PetscMalloc((depth+1) * sizeof(PetscInt), &depthEnd);
196: for (d = 0; d <= depth; ++d) {
197: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
198: }
199: /* Step 8: Convert coordinates */
200: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
201: DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
202: DMPlexGetCoordinateSection(dm, &coordSection);
203: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
204: PetscSectionSetNumFields(newCoordSection, 1);
205: PetscSectionSetFieldComponents(newCoordSection, 0, dim);
206: PetscSectionSetChart(newCoordSection, vStartNew, vEndNew);
207: for (v = vStartNew; v < vEndNew; ++v) {
208: PetscSectionSetDof(newCoordSection, v, dim);
209: PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
210: }
211: PetscSectionSetUp(newCoordSection);
212: DMPlexSetCoordinateSection(dmNew, newCoordSection);
213: PetscSectionGetStorageSize(newCoordSection, &coordSize);
214: VecCreate(PetscObjectComm((PetscObject)dm), &newCoordinates);
215: PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
216: VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
217: VecSetFromOptions(newCoordinates);
218: DMSetCoordinatesLocal(dmNew, newCoordinates);
219: DMGetCoordinatesLocal(dm, &coordinates);
220: VecGetArray(coordinates, &coords);
221: VecGetArray(newCoordinates, &newCoords);
222: for (v = vStart; v < vEnd; ++v) {
223: PetscInt dof, off, noff, d;
225: PetscSectionGetDof(coordSection, v, &dof);
226: PetscSectionGetOffset(coordSection, v, &off);
227: PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthEnd, depthShift), &noff);
228: for (d = 0; d < dof; ++d) {
229: newCoords[noff+d] = coords[off+d];
230: }
231: }
232: VecRestoreArray(coordinates, &coords);
233: VecRestoreArray(newCoordinates, &newCoords);
234: VecDestroy(&newCoordinates);
235: PetscSectionDestroy(&newCoordSection);
236: PetscFree(depthEnd);
237: return(0);
238: }
242: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
243: {
244: PetscInt *depthEnd;
245: PetscInt depth = 0, d;
246: PetscSF sfPoint, sfPointNew;
247: const PetscSFNode *remotePoints;
248: PetscSFNode *gremotePoints;
249: const PetscInt *localPoints;
250: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
251: PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
252: PetscMPIInt numProcs;
253: PetscErrorCode ierr;
256: DMPlexGetDepth(dm, &depth);
257: PetscMalloc((depth+1) * sizeof(PetscInt), &depthEnd);
258: for (d = 0; d <= depth; ++d) {
259: totShift += depthShift[d];
260: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
261: }
262: /* Step 9: Convert pointSF */
263: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
264: DMGetPointSF(dm, &sfPoint);
265: DMGetPointSF(dmNew, &sfPointNew);
266: DMPlexGetChart(dm, &pStart, &pEnd);
267: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
268: if (numRoots >= 0) {
269: PetscMalloc2(numRoots,PetscInt,&newLocation,pEnd-pStart,PetscInt,&newRemoteLocation);
270: for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthEnd, depthShift);
271: PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
272: PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
273: PetscMalloc(numLeaves * sizeof(PetscInt), &glocalPoints);
274: PetscMalloc(numLeaves * sizeof(PetscSFNode), &gremotePoints);
275: for (l = 0; l < numLeaves; ++l) {
276: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthEnd, depthShift);
277: gremotePoints[l].rank = remotePoints[l].rank;
278: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
279: }
280: PetscFree2(newLocation,newRemoteLocation);
281: PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
282: }
283: PetscFree(depthEnd);
284: return(0);
285: }
289: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
290: {
291: PetscSF sfPoint;
292: DMLabel vtkLabel, ghostLabel;
293: PetscInt *depthEnd;
294: const PetscSFNode *leafRemote;
295: const PetscInt *leafLocal;
296: PetscInt depth = 0, d, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
297: PetscMPIInt rank;
298: PetscErrorCode ierr;
301: DMPlexGetDepth(dm, &depth);
302: PetscMalloc((depth+1) * sizeof(PetscInt), &depthEnd);
303: for (d = 0; d <= depth; ++d) {
304: DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
305: }
306: /* Step 10: Convert labels */
307: DMPlexGetNumLabels(dm, &numLabels);
308: for (l = 0; l < numLabels; ++l) {
309: DMLabel label, newlabel;
310: const char *lname;
311: PetscBool isDepth;
312: IS valueIS;
313: const PetscInt *values;
314: PetscInt numValues, val;
316: DMPlexGetLabelName(dm, l, &lname);
317: PetscStrcmp(lname, "depth", &isDepth);
318: if (isDepth) continue;
319: DMPlexCreateLabel(dmNew, lname);
320: DMPlexGetLabel(dm, lname, &label);
321: DMPlexGetLabel(dmNew, lname, &newlabel);
322: DMLabelGetValueIS(label, &valueIS);
323: ISGetLocalSize(valueIS, &numValues);
324: ISGetIndices(valueIS, &values);
325: for (val = 0; val < numValues; ++val) {
326: IS pointIS;
327: const PetscInt *points;
328: PetscInt numPoints, p;
330: DMLabelGetStratumIS(label, values[val], &pointIS);
331: ISGetLocalSize(pointIS, &numPoints);
332: ISGetIndices(pointIS, &points);
333: for (p = 0; p < numPoints; ++p) {
334: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthEnd, depthShift);
336: DMLabelSetValue(newlabel, newpoint, values[val]);
337: }
338: ISRestoreIndices(pointIS, &points);
339: ISDestroy(&pointIS);
340: }
341: ISRestoreIndices(valueIS, &values);
342: ISDestroy(&valueIS);
343: }
344: PetscFree(depthEnd);
345: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
346: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
347: DMGetPointSF(dm, &sfPoint);
348: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
349: PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
350: DMPlexCreateLabel(dmNew, "vtk");
351: DMPlexCreateLabel(dmNew, "ghost");
352: DMPlexGetLabel(dmNew, "vtk", &vtkLabel);
353: DMPlexGetLabel(dmNew, "ghost", &ghostLabel);
354: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
355: for (; c < leafLocal[l] && c < cEnd; ++c) {
356: DMLabelSetValue(vtkLabel, c, 1);
357: }
358: if (leafLocal[l] >= cEnd) break;
359: if (leafRemote[l].rank == rank) {
360: DMLabelSetValue(vtkLabel, c, 1);
361: } else {
362: DMLabelSetValue(ghostLabel, c, 2);
363: }
364: }
365: for (; c < cEnd; ++c) {
366: DMLabelSetValue(vtkLabel, c, 1);
367: }
368: if (0) {
369: PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
370: DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
371: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
372: }
373: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
374: for (f = fStart; f < fEnd; ++f) {
375: PetscInt numCells;
377: DMPlexGetSupportSize(dmNew, f, &numCells);
378: if (numCells < 2) {
379: DMLabelSetValue(ghostLabel, f, 1);
380: } else {
381: const PetscInt *cells = NULL;
382: PetscInt vA, vB;
384: DMPlexGetSupport(dmNew, f, &cells);
385: DMLabelGetValue(vtkLabel, cells[0], &vA);
386: DMLabelGetValue(vtkLabel, cells[1], &vB);
387: if (!vA && !vB) {DMLabelSetValue(ghostLabel, f, 1);}
388: }
389: }
390: if (0) {
391: PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
392: DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
393: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
394: }
395: return(0);
396: }
400: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
401: {
402: IS valueIS;
403: const PetscInt *values;
404: PetscInt *depthShift;
405: PetscInt depth = 0, numFS, fs, ghostCell, cEnd, c;
406: PetscErrorCode ierr;
409: /* Count ghost cells */
410: DMLabelGetValueIS(label, &valueIS);
411: ISGetLocalSize(valueIS, &numFS);
412: ISGetIndices(valueIS, &values);
414: *numGhostCells = 0;
415: for (fs = 0; fs < numFS; ++fs) {
416: PetscInt numBdFaces;
418: DMLabelGetStratumSize(label, values[fs], &numBdFaces);
420: *numGhostCells += numBdFaces;
421: }
422: DMPlexGetDepth(dm, &depth);
423: PetscMalloc((depth+1) * sizeof(PetscInt), &depthShift);
424: PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));
425: if (depth >= 0) depthShift[depth] = *numGhostCells;
426: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
427: /* Step 3: Set cone/support sizes for new points */
428: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
429: for (c = cEnd; c < cEnd + *numGhostCells; ++c) {
430: DMPlexSetConeSize(gdm, c, 1);
431: }
432: for (fs = 0; fs < numFS; ++fs) {
433: IS faceIS;
434: const PetscInt *faces;
435: PetscInt numFaces, f;
437: DMLabelGetStratumIS(label, values[fs], &faceIS);
438: ISGetLocalSize(faceIS, &numFaces);
439: ISGetIndices(faceIS, &faces);
440: for (f = 0; f < numFaces; ++f) {
441: PetscInt size;
443: DMPlexGetSupportSize(dm, faces[f], &size);
444: if (size != 1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
445: DMPlexSetSupportSize(gdm, faces[f] + *numGhostCells, 2);
446: }
447: ISRestoreIndices(faceIS, &faces);
448: ISDestroy(&faceIS);
449: }
450: /* Step 4: Setup ghosted DM */
451: DMSetUp(gdm);
452: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
453: /* Step 6: Set cones and supports for new points */
454: ghostCell = cEnd;
455: for (fs = 0; fs < numFS; ++fs) {
456: IS faceIS;
457: const PetscInt *faces;
458: PetscInt numFaces, f;
460: DMLabelGetStratumIS(label, values[fs], &faceIS);
461: ISGetLocalSize(faceIS, &numFaces);
462: ISGetIndices(faceIS, &faces);
463: for (f = 0; f < numFaces; ++f, ++ghostCell) {
464: PetscInt newFace = faces[f] + *numGhostCells;
466: DMPlexSetCone(gdm, ghostCell, &newFace);
467: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
468: }
469: ISRestoreIndices(faceIS, &faces);
470: ISDestroy(&faceIS);
471: }
472: ISRestoreIndices(valueIS, &values);
473: ISDestroy(&valueIS);
474: /* Step 7: Stratify */
475: DMPlexStratify(gdm);
476: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
477: DMPlexShiftSF_Internal(dm, depthShift, gdm);
478: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
479: PetscFree(depthShift);
480: return(0);
481: }
485: /*@C
486: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
488: Collective on dm
490: Input Parameters:
491: + dm - The original DM
492: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
494: Output Parameters:
495: + numGhostCells - The number of ghost cells added to the DM
496: - dmGhosted - The new DM
498: Note: If no label exists of that name, one will be created marking all boundary faces
500: Level: developer
502: .seealso: DMCreate()
503: @*/
504: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
505: {
506: DM gdm;
507: DMLabel label;
508: const char *name = labelName ? labelName : "Face Sets";
509: PetscInt dim;
516: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
517: DMSetType(gdm, DMPLEX);
518: DMPlexGetDimension(dm, &dim);
519: DMPlexSetDimension(gdm, dim);
520: DMPlexGetLabel(dm, name, &label);
521: if (!label) {
522: /* Get label for boundary faces */
523: DMPlexCreateLabel(dm, name);
524: DMPlexGetLabel(dm, name, &label);
525: DMPlexMarkBoundaryFaces(dm, label);
526: }
527: DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
528: DMSetFromOptions(gdm);
529: *dmGhosted = gdm;
530: return(0);
531: }
535: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
536: {
537: MPI_Comm comm;
538: IS valueIS, *pointIS;
539: const PetscInt *values, **splitPoints;
540: PetscSection coordSection;
541: Vec coordinates;
542: PetscScalar *coords;
543: PetscInt *depthShift, *depthOffset, *pMaxNew, *numSplitPoints, *coneNew, *supportNew;
544: PetscInt shift = 100, depth = 0, dep, dim, d, numSP = 0, sp, maxConeSize, maxSupportSize, numLabels, p, v;
545: PetscErrorCode ierr;
548: PetscObjectGetComm((PetscObject)dm,&comm);
549: DMPlexGetDimension(dm, &dim);
550: /* Count split points and add cohesive cells */
551: if (label) {
552: DMLabelGetValueIS(label, &valueIS);
553: ISGetLocalSize(valueIS, &numSP);
554: ISGetIndices(valueIS, &values);
555: }
556: DMPlexGetDepth(dm, &depth);
557: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
558: PetscMalloc5(depth+1,PetscInt,&depthShift,depth+1,PetscInt,&depthOffset,depth+1,PetscInt,&pMaxNew,maxConeSize*3,PetscInt,&coneNew,maxSupportSize,PetscInt,&supportNew);
559: PetscMalloc3(depth+1,IS,&pointIS,depth+1,PetscInt,&numSplitPoints,depth+1,const PetscInt*,&splitPoints);
560: PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));
561: for (d = 0; d <= depth; ++d) {
562: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
563: numSplitPoints[d] = 0;
564: splitPoints[d] = NULL;
565: pointIS[d] = NULL;
566: }
567: for (sp = 0; sp < numSP; ++sp) {
568: const PetscInt dep = values[sp];
570: if ((dep < 0) || (dep > depth)) continue;
571: DMLabelGetStratumSize(label, dep, &depthShift[dep]);
572: DMLabelGetStratumIS(label, dep, &pointIS[dep]);
573: if (pointIS[dep]) {
574: ISGetLocalSize(pointIS[dep], &numSplitPoints[dep]);
575: ISGetIndices(pointIS[dep], &splitPoints[dep]);
576: }
577: }
578: if (depth >= 0) {
579: /* Calculate number of additional points */
580: depthShift[depth] = depthShift[depth-1]; /* There is a cohesive cell for every split face */
581: depthShift[1] += depthShift[0]; /* There is a cohesive edge for every split vertex */
582: /* Calculate hybrid bound for each dimension */
583: pMaxNew[0] += depthShift[depth];
584: if (depth > 1) pMaxNew[dim-1] += depthShift[depth] + depthShift[0];
585: if (depth > 2) pMaxNew[1] += depthShift[depth] + depthShift[0] + depthShift[dim-1];
587: /* Calculate point offset for each dimension */
588: depthOffset[depth] = 0;
589: depthOffset[0] = depthOffset[depth] + depthShift[depth];
590: if (depth > 1) depthOffset[dim-1] = depthOffset[0] + depthShift[0];
591: if (depth > 2) depthOffset[1] = depthOffset[dim-1] + depthShift[dim-1];
592: }
593: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
594: /* Step 3: Set cone/support sizes for new points */
595: for (dep = 0; dep <= depth; ++dep) {
596: for (p = 0; p < numSplitPoints[dep]; ++p) {
597: const PetscInt oldp = splitPoints[dep][p];
598: const PetscInt newp = depthOffset[dep] + oldp;
599: const PetscInt splitp = pMaxNew[dep] + p;
600: const PetscInt *support;
601: PetscInt coneSize, supportSize, q, e;
603: DMPlexGetConeSize(dm, oldp, &coneSize);
604: DMPlexSetConeSize(sdm, splitp, coneSize);
605: DMPlexGetSupportSize(dm, oldp, &supportSize);
606: DMPlexSetSupportSize(sdm, splitp, supportSize);
607: if (dep == depth-1) {
608: const PetscInt ccell = pMaxNew[depth] + p;
609: /* Add cohesive cells, they are prisms */
610: DMPlexSetConeSize(sdm, ccell, 2 + coneSize);
611: } else if (dep == 0) {
612: const PetscInt cedge = pMaxNew[1] + (depthShift[1] - depthShift[0]) + p;
614: DMPlexGetSupport(dm, oldp, &support);
615: /* Split old vertex: Edges in old split faces and new cohesive edge */
616: for (e = 0, q = 0; e < supportSize; ++e) {
617: PetscInt val;
619: DMLabelGetValue(label, support[e], &val);
620: if ((val == 1) || (val == (shift + 1))) ++q;
621: }
622: DMPlexSetSupportSize(sdm, newp, q+1);
623: /* Split new vertex: Edges in new split faces and new cohesive edge */
624: for (e = 0, q = 0; e < supportSize; ++e) {
625: PetscInt val;
627: DMLabelGetValue(label, support[e], &val);
628: if ((val == 1) || (val == -(shift + 1))) ++q;
629: }
630: DMPlexSetSupportSize(sdm, splitp, q+1);
631: /* Add cohesive edges */
632: DMPlexSetConeSize(sdm, cedge, 2);
633: /* Punt for now on support, you loop over closure, extract faces, check which ones are in the label */
634: } else if (dep == dim-2) {
635: DMPlexGetSupport(dm, oldp, &support);
636: /* Split old edge: Faces in positive side cells and old split faces */
637: for (e = 0, q = 0; e < supportSize; ++e) {
638: PetscInt val;
640: DMLabelGetValue(label, support[e], &val);
641: if ((val == dim-1) || (val == (shift + dim-1))) ++q;
642: }
643: DMPlexSetSupportSize(sdm, newp, q);
644: /* Split new edge: Faces in negative side cells and new split faces */
645: for (e = 0, q = 0; e < supportSize; ++e) {
646: PetscInt val;
648: DMLabelGetValue(label, support[e], &val);
649: if ((val == dim-1) || (val == -(shift + dim-1))) ++q;
650: }
651: DMPlexSetSupportSize(sdm, splitp, q);
652: }
653: }
654: }
655: /* Step 4: Setup split DM */
656: DMSetUp(sdm);
657: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
658: /* Step 6: Set cones and supports for new points */
659: for (dep = 0; dep <= depth; ++dep) {
660: for (p = 0; p < numSplitPoints[dep]; ++p) {
661: const PetscInt oldp = splitPoints[dep][p];
662: const PetscInt newp = depthOffset[dep] + oldp;
663: const PetscInt splitp = pMaxNew[dep] + p;
664: const PetscInt *cone, *support, *ornt;
665: PetscInt coneSize, supportSize, q, v, e, s;
667: DMPlexGetConeSize(dm, oldp, &coneSize);
668: DMPlexGetCone(dm, oldp, &cone);
669: DMPlexGetConeOrientation(dm, oldp, &ornt);
670: DMPlexGetSupportSize(dm, oldp, &supportSize);
671: DMPlexGetSupport(dm, oldp, &support);
672: if (dep == depth-1) {
673: const PetscInt ccell = pMaxNew[depth] + p;
674: const PetscInt *supportF;
676: /* Split face: copy in old face to new face to start */
677: DMPlexGetSupport(sdm, newp, &supportF);
678: DMPlexSetSupport(sdm, splitp, supportF);
679: /* Split old face: old vertices/edges in cone so no change */
680: /* Split new face: new vertices/edges in cone */
681: for (q = 0; q < coneSize; ++q) {
682: PetscFindInt(cone[q], numSplitPoints[dim-2], splitPoints[dim-2], &v);
684: coneNew[2+q] = pMaxNew[dim-2] + v;
685: }
686: DMPlexSetCone(sdm, splitp, &coneNew[2]);
687: DMPlexSetConeOrientation(sdm, splitp, ornt);
688: /* Cohesive cell: Old and new split face, then new cohesive edges */
689: coneNew[0] = newp;
690: coneNew[1] = splitp;
691: for (q = 0; q < coneSize; ++q) {
692: coneNew[2+q] = (pMaxNew[1] - pMaxNew[dim-2]) + (depthShift[1] - depthShift[0]) + coneNew[2+q];
693: }
694: DMPlexSetCone(sdm, ccell, coneNew);
697: for (s = 0; s < supportSize; ++s) {
698: PetscInt val;
700: DMLabelGetValue(label, support[s], &val);
701: if (val < 0) {
702: /* Split old face: Replace negative side cell with cohesive cell */
703: DMPlexInsertSupport(sdm, newp, s, ccell);
704: } else {
705: /* Split new face: Replace positive side cell with cohesive cell */
706: DMPlexInsertSupport(sdm, splitp, s, ccell);
707: }
708: }
709: } else if (dep == 0) {
710: const PetscInt cedge = pMaxNew[1] + (depthShift[1] - depthShift[0]) + p;
712: /* Split old vertex: Edges in old split faces and new cohesive edge */
713: for (e = 0, q = 0; e < supportSize; ++e) {
714: PetscInt val;
716: DMLabelGetValue(label, support[e], &val);
717: if ((val == 1) || (val == (shift + 1))) {
718: supportNew[q++] = depthOffset[1] + support[e];
719: }
720: }
721: supportNew[q] = cedge;
723: DMPlexSetSupport(sdm, newp, supportNew);
724: /* Split new vertex: Edges in new split faces and new cohesive edge */
725: for (e = 0, q = 0; e < supportSize; ++e) {
726: PetscInt val, edge;
728: DMLabelGetValue(label, support[e], &val);
729: if (val == 1) {
730: PetscFindInt(support[e], numSplitPoints[1], splitPoints[1], &edge);
731: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
732: supportNew[q++] = pMaxNew[1] + edge;
733: } else if (val == -(shift + 1)) {
734: supportNew[q++] = depthOffset[1] + support[e];
735: }
736: }
737: supportNew[q] = cedge;
738: DMPlexSetSupport(sdm, splitp, supportNew);
739: /* Cohesive edge: Old and new split vertex, punting on support */
740: coneNew[0] = newp;
741: coneNew[1] = splitp;
742: DMPlexSetCone(sdm, cedge, coneNew);
743: } else if (dep == dim-2) {
744: /* Split old edge: old vertices in cone so no change */
745: /* Split new edge: new vertices in cone */
746: for (q = 0; q < coneSize; ++q) {
747: PetscFindInt(cone[q], numSplitPoints[dim-3], splitPoints[dim-3], &v);
749: coneNew[q] = pMaxNew[dim-3] + v;
750: }
751: DMPlexSetCone(sdm, splitp, coneNew);
752: /* Split old edge: Faces in positive side cells and old split faces */
753: for (e = 0, q = 0; e < supportSize; ++e) {
754: PetscInt val;
756: DMLabelGetValue(label, support[e], &val);
757: if ((val == dim-1) || (val == (shift + dim-1))) {
758: supportNew[q++] = depthOffset[dim-1] + support[e];
759: }
760: }
761: DMPlexSetSupport(sdm, newp, supportNew);
762: /* Split new edge: Faces in negative side cells and new split faces */
763: for (e = 0, q = 0; e < supportSize; ++e) {
764: PetscInt val, face;
766: DMLabelGetValue(label, support[e], &val);
767: if (val == dim-1) {
768: PetscFindInt(support[e], numSplitPoints[dim-1], splitPoints[dim-1], &face);
769: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
770: supportNew[q++] = pMaxNew[dim-1] + face;
771: } else if (val == -(shift + dim-1)) {
772: supportNew[q++] = depthOffset[dim-1] + support[e];
773: }
774: }
775: DMPlexSetSupport(sdm, splitp, supportNew);
776: }
777: }
778: }
779: /* Step 6b: Replace split points in negative side cones */
780: for (sp = 0; sp < numSP; ++sp) {
781: PetscInt dep = values[sp];
782: IS pIS;
783: PetscInt numPoints;
784: const PetscInt *points;
786: if (dep >= 0) continue;
787: DMLabelGetStratumIS(label, dep, &pIS);
788: if (!pIS) continue;
789: dep = -dep - shift;
790: ISGetLocalSize(pIS, &numPoints);
791: ISGetIndices(pIS, &points);
792: for (p = 0; p < numPoints; ++p) {
793: const PetscInt oldp = points[p];
794: const PetscInt newp = depthOffset[dep] + oldp;
795: const PetscInt *cone;
796: PetscInt coneSize, c;
797: PetscBool replaced = PETSC_FALSE;
799: /* Negative edge: replace split vertex */
800: /* Negative cell: replace split face */
801: DMPlexGetConeSize(sdm, newp, &coneSize);
802: DMPlexGetCone(sdm, newp, &cone);
803: for (c = 0; c < coneSize; ++c) {
804: const PetscInt coldp = cone[c] - depthOffset[dep-1];
805: PetscInt csplitp, cp, val;
807: DMLabelGetValue(label, coldp, &val);
808: if (val == dep-1) {
809: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
810: if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
811: csplitp = pMaxNew[dep-1] + cp;
812: DMPlexInsertCone(sdm, newp, c, csplitp);
813: replaced = PETSC_TRUE;
814: }
815: }
816: if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp);
817: }
818: ISRestoreIndices(pIS, &points);
819: ISDestroy(&pIS);
820: }
821: /* Step 7: Stratify */
822: DMPlexStratify(sdm);
823: /* Step 8: Coordinates */
824: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
825: DMPlexGetCoordinateSection(sdm, &coordSection);
826: DMGetCoordinatesLocal(sdm, &coordinates);
827: VecGetArray(coordinates, &coords);
828: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
829: const PetscInt newp = depthOffset[0] + splitPoints[0][v];
830: const PetscInt splitp = pMaxNew[0] + v;
831: PetscInt dof, off, soff, d;
833: PetscSectionGetDof(coordSection, newp, &dof);
834: PetscSectionGetOffset(coordSection, newp, &off);
835: PetscSectionGetOffset(coordSection, splitp, &soff);
836: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
837: }
838: VecRestoreArray(coordinates, &coords);
839: /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
840: DMPlexShiftSF_Internal(dm, depthShift, sdm);
841: /* Step 10: Labels */
842: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
843: DMPlexGetNumLabels(sdm, &numLabels);
844: for (dep = 0; dep <= depth; ++dep) {
845: for (p = 0; p < numSplitPoints[dep]; ++p) {
846: const PetscInt newp = depthOffset[dep] + splitPoints[dep][p];
847: const PetscInt splitp = pMaxNew[dep] + p;
848: PetscInt l;
850: for (l = 0; l < numLabels; ++l) {
851: DMLabel mlabel;
852: const char *lname;
853: PetscInt val;
855: DMPlexGetLabelName(sdm, l, &lname);
856: DMPlexGetLabel(sdm, lname, &mlabel);
857: DMLabelGetValue(mlabel, newp, &val);
858: if (val >= 0) {
859: DMLabelSetValue(mlabel, splitp, val);
860: if (dep == 0) {
861: const PetscInt cedge = pMaxNew[1] + (depthShift[1] - depthShift[0]) + p;
862: DMLabelSetValue(mlabel, cedge, val);
863: }
864: }
865: }
866: }
867: }
868: for (sp = 0; sp < numSP; ++sp) {
869: const PetscInt dep = values[sp];
871: if ((dep < 0) || (dep > depth)) continue;
872: if (pointIS[dep]) {ISRestoreIndices(pointIS[dep], &splitPoints[dep]);}
873: ISDestroy(&pointIS[dep]);
874: }
875: if (label) {
876: ISRestoreIndices(valueIS, &values);
877: ISDestroy(&valueIS);
878: }
879: PetscFree5(depthShift, depthOffset, pMaxNew, coneNew, supportNew);
880: PetscFree3(pointIS, numSplitPoints, splitPoints);
881: return(0);
882: }
886: /*@C
887: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
889: Collective on dm
891: Input Parameters:
892: + dm - The original DM
893: - labelName - The label specifying the boundary faces (this could be auto-generated)
895: Output Parameters:
896: - dmSplit - The new DM
898: Level: developer
900: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
901: @*/
902: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
903: {
904: DM sdm;
905: PetscInt dim;
911: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
912: DMSetType(sdm, DMPLEX);
913: DMPlexGetDimension(dm, &dim);
914: DMPlexSetDimension(sdm, dim);
915: switch (dim) {
916: case 2:
917: case 3:
918: DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
919: break;
920: default:
921: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
922: }
923: *dmSplit = sdm;
924: return(0);
925: }
929: /*@
930: DMPlexLabelCohesiveComplete - Starting with a label marking vertices on an internal surface, we add all other mesh pieces
931: to complete the surface
933: Input Parameters:
934: + dm - The DM
935: - label - A DMLabel marking the surface vertices
937: Output Parameter:
938: . label - A DMLabel marking all surface points
940: Level: developer
942: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
943: @*/
944: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label)
945: {
946: IS dimIS;
947: const PetscInt *points;
948: PetscInt shift = 100, dim, dep, cStart, cEnd, numPoints, p, val;
949: PetscErrorCode ierr;
952: DMPlexGetDimension(dm, &dim);
953: /* Cell orientation for face gives the side of the fault */
954: DMLabelGetStratumIS(label, dim-1, &dimIS);
955: if (!dimIS) return(0);
956: ISGetLocalSize(dimIS, &numPoints);
957: ISGetIndices(dimIS, &points);
958: for (p = 0; p < numPoints; ++p) {
959: const PetscInt *support;
960: PetscInt supportSize, s;
962: DMPlexGetSupportSize(dm, points[p], &supportSize);
963: if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
964: DMPlexGetSupport(dm, points[p], &support);
965: for (s = 0; s < supportSize; ++s) {
966: const PetscInt *cone, *ornt;
967: PetscInt coneSize, c;
968: PetscBool pos = PETSC_TRUE;
970: DMPlexGetConeSize(dm, support[s], &coneSize);
971: DMPlexGetCone(dm, support[s], &cone);
972: DMPlexGetConeOrientation(dm, support[s], &ornt);
973: for (c = 0; c < coneSize; ++c) {
974: if (cone[c] == points[p]) {
975: if (ornt[c] >= 0) {
976: DMLabelSetValue(label, support[s], shift+dim);
977: } else {
978: DMLabelSetValue(label, support[s], -(shift+dim));
979: pos = PETSC_FALSE;
980: }
981: break;
982: }
983: }
984: if (c == coneSize) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cell split face %d support does not have it in the cone", points[p]);
985: /* Put faces touching the fault in the label */
986: for (c = 0; c < coneSize; ++c) {
987: const PetscInt point = cone[c];
989: DMLabelGetValue(label, point, &val);
990: if (val == -1) {
991: PetscInt *closure = NULL;
992: PetscInt closureSize, cl;
994: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
995: for (cl = 0; cl < closureSize*2; cl += 2) {
996: const PetscInt clp = closure[cl];
998: DMLabelGetValue(label, clp, &val);
999: if ((val >= 0) && (val < dim-1)) {
1000: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1001: break;
1002: }
1003: }
1004: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1005: }
1006: }
1007: }
1008: }
1009: ISRestoreIndices(dimIS, &points);
1010: ISDestroy(&dimIS);
1011: /* Search for other cells/faces/edges connected to the fault by a vertex */
1012: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1013: DMLabelGetStratumIS(label, 0, &dimIS);
1014: if (!dimIS) return(0);
1015: ISGetLocalSize(dimIS, &numPoints);
1016: ISGetIndices(dimIS, &points);
1017: for (p = 0; p < numPoints; ++p) {
1018: PetscInt *star = NULL;
1019: PetscInt starSize, s;
1020: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1022: /* First mark cells connected to the fault */
1023: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1024: while (again) {
1025: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1026: again = 0;
1027: for (s = 0; s < starSize*2; s += 2) {
1028: const PetscInt point = star[s];
1029: const PetscInt *cone;
1030: PetscInt coneSize, c;
1032: if ((point < cStart) || (point >= cEnd)) continue;
1033: DMLabelGetValue(label, point, &val);
1034: if (val != -1) continue;
1035: again = 2;
1036: DMPlexGetConeSize(dm, point, &coneSize);
1037: DMPlexGetCone(dm, point, &cone);
1038: for (c = 0; c < coneSize; ++c) {
1039: DMLabelGetValue(label, cone[c], &val);
1040: if (val != -1) {
1041: if (abs(val) < shift) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d on cell %d has an invalid label %d", cone[c], point, val);
1042: if (val > 0) {
1043: DMLabelSetValue(label, point, shift+dim);
1044: } else {
1045: DMLabelSetValue(label, point, -(shift+dim));
1046: }
1047: again = 1;
1048: break;
1049: }
1050: }
1051: }
1052: }
1053: /* Classify the rest by cell membership */
1054: for (s = 0; s < starSize*2; s += 2) {
1055: const PetscInt point = star[s];
1057: DMLabelGetValue(label, point, &val);
1058: if (val == -1) {
1059: PetscInt *sstar = NULL;
1060: PetscInt sstarSize, ss;
1061: PetscBool marked = PETSC_FALSE;
1063: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1064: for (ss = 0; ss < sstarSize*2; ss += 2) {
1065: const PetscInt spoint = sstar[ss];
1067: if ((spoint < cStart) || (spoint >= cEnd)) continue;
1068: DMLabelGetValue(label, spoint, &val);
1069: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1070: DMPlexGetLabelValue(dm, "depth", point, &dep);
1071: if (val > 0) {
1072: DMLabelSetValue(label, point, shift+dep);
1073: } else {
1074: DMLabelSetValue(label, point, -(shift+dep));
1075: }
1076: marked = PETSC_TRUE;
1077: break;
1078: }
1079: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1080: if (!marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1081: }
1082: }
1083: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1084: }
1085: ISRestoreIndices(dimIS, &points);
1086: ISDestroy(&dimIS);
1087: return(0);
1088: }
1092: /* Here we need the explicit assumption that:
1094: For any marked cell, the marked vertices constitute a single face
1095: */
1096: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1097: {
1098: IS subvertexIS;
1099: const PetscInt *subvertices;
1100: PetscInt *pStart, *pEnd, *pMax, pSize;
1101: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
1102: PetscErrorCode ierr;
1105: *numFaces = 0;
1106: *nFV = 0;
1107: DMPlexGetDepth(dm, &depth);
1108: DMPlexGetDimension(dm, &dim);
1109: pSize = PetscMax(depth, dim) + 1;
1110: PetscMalloc3(pSize,PetscInt,&pStart,pSize,PetscInt,&pEnd,pSize,PetscInt,&pMax);
1111: DMPlexGetHybridBounds(dm, &pMax[depth], depth>1 ? &pMax[depth-1] : NULL, depth > 2 ? &pMax[1] : NULL, &pMax[0]);
1112: for (d = 0; d <= depth; ++d) {
1113: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1114: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1115: }
1116: /* Loop over initial vertices and mark all faces in the collective star() */
1117: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
1118: if (subvertexIS) {
1119: ISGetSize(subvertexIS, &numSubVerticesInitial);
1120: ISGetIndices(subvertexIS, &subvertices);
1121: }
1122: for (v = 0; v < numSubVerticesInitial; ++v) {
1123: const PetscInt vertex = subvertices[v];
1124: PetscInt *star = NULL;
1125: PetscInt starSize, s, numCells = 0, c;
1127: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1128: for (s = 0; s < starSize*2; s += 2) {
1129: const PetscInt point = star[s];
1130: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1131: }
1132: for (c = 0; c < numCells; ++c) {
1133: const PetscInt cell = star[c];
1134: PetscInt *closure = NULL;
1135: PetscInt closureSize, cl;
1136: PetscInt cellLoc, numCorners = 0, faceSize = 0;
1138: DMLabelGetValue(subpointMap, cell, &cellLoc);
1139: if (cellLoc == 2) continue;
1140: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1141: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1142: for (cl = 0; cl < closureSize*2; cl += 2) {
1143: const PetscInt point = closure[cl];
1144: PetscInt vertexLoc;
1146: if ((point >= pStart[0]) && (point < pEnd[0])) {
1147: ++numCorners;
1148: DMLabelGetValue(vertexLabel, point, &vertexLoc);
1149: if (vertexLoc == value) closure[faceSize++] = point;
1150: }
1151: }
1152: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1153: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1154: if (faceSize == *nFV) {
1155: const PetscInt *cells = NULL;
1156: PetscInt numCells, nc;
1158: ++(*numFaces);
1159: for (cl = 0; cl < faceSize; ++cl) {
1160: DMLabelSetValue(subpointMap, closure[cl], 0);
1161: }
1162: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1163: for (nc = 0; nc < numCells; ++nc) {
1164: DMLabelSetValue(subpointMap, cells[nc], 2);
1165: }
1166: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1167: }
1168: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1169: }
1170: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1171: }
1172: if (subvertexIS) {
1173: ISRestoreIndices(subvertexIS, &subvertices);
1174: }
1175: ISDestroy(&subvertexIS);
1176: PetscFree3(pStart,pEnd,pMax);
1177: return(0);
1178: }
1182: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
1183: {
1184: IS subvertexIS;
1185: const PetscInt *subvertices;
1186: PetscInt *pStart, *pEnd, *pMax;
1187: PetscInt dim, d, numSubVerticesInitial = 0, v;
1188: PetscErrorCode ierr;
1191: DMPlexGetDimension(dm, &dim);
1192: PetscMalloc3(dim+1,PetscInt,&pStart,dim+1,PetscInt,&pEnd,dim+1,PetscInt,&pMax);
1193: DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
1194: for (d = 0; d <= dim; ++d) {
1195: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1196: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1197: }
1198: /* Loop over initial vertices and mark all faces in the collective star() */
1199: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
1200: if (subvertexIS) {
1201: ISGetSize(subvertexIS, &numSubVerticesInitial);
1202: ISGetIndices(subvertexIS, &subvertices);
1203: }
1204: for (v = 0; v < numSubVerticesInitial; ++v) {
1205: const PetscInt vertex = subvertices[v];
1206: PetscInt *star = NULL;
1207: PetscInt starSize, s, numFaces = 0, f;
1209: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1210: for (s = 0; s < starSize*2; s += 2) {
1211: const PetscInt point = star[s];
1212: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
1213: }
1214: for (f = 0; f < numFaces; ++f) {
1215: const PetscInt face = star[f];
1216: PetscInt *closure = NULL;
1217: PetscInt closureSize, c;
1218: PetscInt faceLoc;
1220: DMLabelGetValue(subpointMap, face, &faceLoc);
1221: if (faceLoc == dim-1) continue;
1222: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
1223: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1224: for (c = 0; c < closureSize*2; c += 2) {
1225: const PetscInt point = closure[c];
1226: PetscInt vertexLoc;
1228: if ((point >= pStart[0]) && (point < pEnd[0])) {
1229: DMLabelGetValue(vertexLabel, point, &vertexLoc);
1230: if (vertexLoc != value) break;
1231: }
1232: }
1233: if (c == closureSize*2) {
1234: const PetscInt *support;
1235: PetscInt supportSize, s;
1237: for (c = 0; c < closureSize*2; c += 2) {
1238: const PetscInt point = closure[c];
1240: for (d = 0; d < dim; ++d) {
1241: if ((point >= pStart[d]) && (point < pEnd[d])) {
1242: DMLabelSetValue(subpointMap, point, d);
1243: break;
1244: }
1245: }
1246: }
1247: DMPlexGetSupportSize(dm, face, &supportSize);
1248: DMPlexGetSupport(dm, face, &support);
1249: for (s = 0; s < supportSize; ++s) {
1250: DMLabelSetValue(subpointMap, support[s], dim);
1251: }
1252: }
1253: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1254: }
1255: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1256: }
1257: if (subvertexIS) {
1258: ISRestoreIndices(subvertexIS, &subvertices);
1259: }
1260: ISDestroy(&subvertexIS);
1261: PetscFree3(pStart,pEnd,pMax);
1262: return(0);
1263: }
1267: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
1268: {
1269: const PetscInt *cone;
1270: PetscInt dim, cMax, cEnd, c, p, coneSize;
1271: PetscErrorCode ierr;
1274: *numFaces = 0;
1275: *nFV = 0;
1276: DMPlexGetDimension(dm, &dim);
1277: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1278: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1279: if (cMax < 0) return(0);
1280: DMPlexGetConeSize(dm, cMax, &coneSize);
1281: *numFaces = cEnd - cMax;
1282: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
1283: PetscMalloc(*numFaces *2 * sizeof(PetscInt), subCells);
1284: for (c = cMax; c < cEnd; ++c) {
1285: const PetscInt *cells;
1286: PetscInt numCells;
1288: DMPlexGetCone(dm, c, &cone);
1289: for (p = 0; p < *nFV; ++p) {
1290: DMLabelSetValue(subpointMap, cone[p], 0);
1291: }
1292: /* Negative face */
1293: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
1294: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
1295: for (p = 0; p < numCells; ++p) {
1296: DMLabelSetValue(subpointMap, cells[p], 2);
1297: (*subCells)[(c-cMax)*2+p] = cells[p];
1298: }
1299: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
1300: /* Positive face is not included */
1301: }
1302: return(0);
1303: }
1307: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, PetscBool hasLagrange, DMLabel subpointMap, DM subdm)
1308: {
1309: PetscInt *pStart, *pEnd;
1310: PetscInt dim, cMax, cEnd, c, d;
1314: DMPlexGetDimension(dm, &dim);
1315: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1316: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1317: if (cMax < 0) return(0);
1318: PetscMalloc2(dim+1,PetscInt,&pStart,dim+1,PetscInt,&pEnd);
1319: for (d = 0; d <= dim; ++d) {
1320: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1321: }
1322: for (c = cMax; c < cEnd; ++c) {
1323: const PetscInt *cone;
1324: PetscInt *closure = NULL;
1325: PetscInt coneSize, closureSize, cl;
1327: DMPlexGetConeSize(dm, c, &coneSize);
1328: if (hasLagrange) {
1329: if (coneSize != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
1330: } else {
1331: if (coneSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
1332: }
1333: /* Negative face */
1334: DMPlexGetCone(dm, c, &cone);
1335: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
1336: for (cl = 0; cl < closureSize*2; cl += 2) {
1337: const PetscInt point = closure[cl];
1339: for (d = 0; d <= dim; ++d) {
1340: if ((point >= pStart[d]) && (point < pEnd[d])) {
1341: DMLabelSetValue(subpointMap, point, d);
1342: break;
1343: }
1344: }
1345: }
1346: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
1347: /* Cells -- positive face is not included */
1348: for (cl = 0; cl < 1; ++cl) {
1349: const PetscInt *support;
1350: PetscInt supportSize, s;
1352: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
1353: if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells");
1354: DMPlexGetSupport(dm, cone[cl], &support);
1355: for (s = 0; s < supportSize; ++s) {
1356: DMLabelSetValue(subpointMap, support[s], dim);
1357: }
1358: }
1359: }
1360: PetscFree2(pStart, pEnd);
1361: return(0);
1362: }
1366: PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
1367: {
1368: MPI_Comm comm;
1369: PetscBool posOrient = PETSC_FALSE;
1370: const PetscInt debug = 0;
1371: PetscInt cellDim, faceSize, f;
1375: PetscObjectGetComm((PetscObject)dm,&comm);
1376: DMPlexGetDimension(dm, &cellDim);
1377: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
1379: if (cellDim == numCorners-1) {
1380: /* Simplices */
1381: faceSize = numCorners-1;
1382: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
1383: } else if (cellDim == 1 && numCorners == 3) {
1384: /* Quadratic line */
1385: faceSize = 1;
1386: posOrient = PETSC_TRUE;
1387: } else if (cellDim == 2 && numCorners == 4) {
1388: /* Quads */
1389: faceSize = 2;
1390: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
1391: posOrient = PETSC_TRUE;
1392: } else if ((indices[0] == 3) && (indices[1] == 0)) {
1393: posOrient = PETSC_TRUE;
1394: } else {
1395: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
1396: posOrient = PETSC_FALSE;
1397: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
1398: }
1399: } else if (cellDim == 2 && numCorners == 6) {
1400: /* Quadratic triangle (I hate this) */
1401: /* Edges are determined by the first 2 vertices (corners of edges) */
1402: const PetscInt faceSizeTri = 3;
1403: PetscInt sortedIndices[3], i, iFace;
1404: PetscBool found = PETSC_FALSE;
1405: PetscInt faceVerticesTriSorted[9] = {
1406: 0, 3, 4, /* bottom */
1407: 1, 4, 5, /* right */
1408: 2, 3, 5, /* left */
1409: };
1410: PetscInt faceVerticesTri[9] = {
1411: 0, 3, 4, /* bottom */
1412: 1, 4, 5, /* right */
1413: 2, 5, 3, /* left */
1414: };
1416: faceSize = faceSizeTri;
1417: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
1418: PetscSortInt(faceSizeTri, sortedIndices);
1419: for (iFace = 0; iFace < 3; ++iFace) {
1420: const PetscInt ii = iFace*faceSizeTri;
1421: PetscInt fVertex, cVertex;
1423: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
1424: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
1425: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
1426: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
1427: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
1428: faceVertices[fVertex] = origVertices[cVertex];
1429: break;
1430: }
1431: }
1432: }
1433: found = PETSC_TRUE;
1434: break;
1435: }
1436: }
1437: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
1438: if (posOriented) *posOriented = PETSC_TRUE;
1439: return(0);
1440: } else if (cellDim == 2 && numCorners == 9) {
1441: /* Quadratic quad (I hate this) */
1442: /* Edges are determined by the first 2 vertices (corners of edges) */
1443: const PetscInt faceSizeQuad = 3;
1444: PetscInt sortedIndices[3], i, iFace;
1445: PetscBool found = PETSC_FALSE;
1446: PetscInt faceVerticesQuadSorted[12] = {
1447: 0, 1, 4, /* bottom */
1448: 1, 2, 5, /* right */
1449: 2, 3, 6, /* top */
1450: 0, 3, 7, /* left */
1451: };
1452: PetscInt faceVerticesQuad[12] = {
1453: 0, 1, 4, /* bottom */
1454: 1, 2, 5, /* right */
1455: 2, 3, 6, /* top */
1456: 3, 0, 7, /* left */
1457: };
1459: faceSize = faceSizeQuad;
1460: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
1461: PetscSortInt(faceSizeQuad, sortedIndices);
1462: for (iFace = 0; iFace < 4; ++iFace) {
1463: const PetscInt ii = iFace*faceSizeQuad;
1464: PetscInt fVertex, cVertex;
1466: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
1467: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
1468: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
1469: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
1470: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
1471: faceVertices[fVertex] = origVertices[cVertex];
1472: break;
1473: }
1474: }
1475: }
1476: found = PETSC_TRUE;
1477: break;
1478: }
1479: }
1480: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
1481: if (posOriented) *posOriented = PETSC_TRUE;
1482: return(0);
1483: } else if (cellDim == 3 && numCorners == 8) {
1484: /* Hexes
1485: A hex is two oriented quads with the normal of the first
1486: pointing up at the second.
1488: 7---6
1489: /| /|
1490: 4---5 |
1491: | 3-|-2
1492: |/ |/
1493: 0---1
1495: Faces are determined by the first 4 vertices (corners of faces) */
1496: const PetscInt faceSizeHex = 4;
1497: PetscInt sortedIndices[4], i, iFace;
1498: PetscBool found = PETSC_FALSE;
1499: PetscInt faceVerticesHexSorted[24] = {
1500: 0, 1, 2, 3, /* bottom */
1501: 4, 5, 6, 7, /* top */
1502: 0, 1, 4, 5, /* front */
1503: 1, 2, 5, 6, /* right */
1504: 2, 3, 6, 7, /* back */
1505: 0, 3, 4, 7, /* left */
1506: };
1507: PetscInt faceVerticesHex[24] = {
1508: 3, 2, 1, 0, /* bottom */
1509: 4, 5, 6, 7, /* top */
1510: 0, 1, 5, 4, /* front */
1511: 1, 2, 6, 5, /* right */
1512: 2, 3, 7, 6, /* back */
1513: 3, 0, 4, 7, /* left */
1514: };
1516: faceSize = faceSizeHex;
1517: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
1518: PetscSortInt(faceSizeHex, sortedIndices);
1519: for (iFace = 0; iFace < 6; ++iFace) {
1520: const PetscInt ii = iFace*faceSizeHex;
1521: PetscInt fVertex, cVertex;
1523: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
1524: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
1525: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
1526: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
1527: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
1528: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
1529: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
1530: faceVertices[fVertex] = origVertices[cVertex];
1531: break;
1532: }
1533: }
1534: }
1535: found = PETSC_TRUE;
1536: break;
1537: }
1538: }
1539: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
1540: if (posOriented) *posOriented = PETSC_TRUE;
1541: return(0);
1542: } else if (cellDim == 3 && numCorners == 10) {
1543: /* Quadratic tet */
1544: /* Faces are determined by the first 3 vertices (corners of faces) */
1545: const PetscInt faceSizeTet = 6;
1546: PetscInt sortedIndices[6], i, iFace;
1547: PetscBool found = PETSC_FALSE;
1548: PetscInt faceVerticesTetSorted[24] = {
1549: 0, 1, 2, 6, 7, 8, /* bottom */
1550: 0, 3, 4, 6, 7, 9, /* front */
1551: 1, 4, 5, 7, 8, 9, /* right */
1552: 2, 3, 5, 6, 8, 9, /* left */
1553: };
1554: PetscInt faceVerticesTet[24] = {
1555: 0, 1, 2, 6, 7, 8, /* bottom */
1556: 0, 4, 3, 6, 7, 9, /* front */
1557: 1, 5, 4, 7, 8, 9, /* right */
1558: 2, 3, 5, 8, 6, 9, /* left */
1559: };
1561: faceSize = faceSizeTet;
1562: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
1563: PetscSortInt(faceSizeTet, sortedIndices);
1564: for (iFace=0; iFace < 4; ++iFace) {
1565: const PetscInt ii = iFace*faceSizeTet;
1566: PetscInt fVertex, cVertex;
1568: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
1569: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
1570: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
1571: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
1572: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
1573: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
1574: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
1575: faceVertices[fVertex] = origVertices[cVertex];
1576: break;
1577: }
1578: }
1579: }
1580: found = PETSC_TRUE;
1581: break;
1582: }
1583: }
1584: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
1585: if (posOriented) *posOriented = PETSC_TRUE;
1586: return(0);
1587: } else if (cellDim == 3 && numCorners == 27) {
1588: /* Quadratic hexes (I hate this)
1589: A hex is two oriented quads with the normal of the first
1590: pointing up at the second.
1592: 7---6
1593: /| /|
1594: 4---5 |
1595: | 3-|-2
1596: |/ |/
1597: 0---1
1599: Faces are determined by the first 4 vertices (corners of faces) */
1600: const PetscInt faceSizeQuadHex = 9;
1601: PetscInt sortedIndices[9], i, iFace;
1602: PetscBool found = PETSC_FALSE;
1603: PetscInt faceVerticesQuadHexSorted[54] = {
1604: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
1605: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
1606: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
1607: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
1608: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
1609: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
1610: };
1611: PetscInt faceVerticesQuadHex[54] = {
1612: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
1613: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
1614: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
1615: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
1616: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
1617: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
1618: };
1620: faceSize = faceSizeQuadHex;
1621: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
1622: PetscSortInt(faceSizeQuadHex, sortedIndices);
1623: for (iFace = 0; iFace < 6; ++iFace) {
1624: const PetscInt ii = iFace*faceSizeQuadHex;
1625: PetscInt fVertex, cVertex;
1627: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
1628: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
1629: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
1630: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
1631: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
1632: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
1633: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
1634: faceVertices[fVertex] = origVertices[cVertex];
1635: break;
1636: }
1637: }
1638: }
1639: found = PETSC_TRUE;
1640: break;
1641: }
1642: }
1643: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
1644: if (posOriented) *posOriented = PETSC_TRUE;
1645: return(0);
1646: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
1647: if (!posOrient) {
1648: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
1649: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
1650: } else {
1651: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
1652: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
1653: }
1654: if (posOriented) *posOriented = posOrient;
1655: return(0);
1656: }
1660: /*
1661: Given a cell and a face, as a set of vertices,
1662: return the oriented face, as a set of vertices, in faceVertices
1663: The orientation is such that the face normal points out of the cell
1664: */
1665: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
1666: {
1667: const PetscInt *cone = NULL;
1668: PetscInt coneSize, v, f, v2;
1669: PetscInt oppositeVertex = -1;
1670: PetscErrorCode ierr;
1673: DMPlexGetConeSize(dm, cell, &coneSize);
1674: DMPlexGetCone(dm, cell, &cone);
1675: for (v = 0, v2 = 0; v < coneSize; ++v) {
1676: PetscBool found = PETSC_FALSE;
1678: for (f = 0; f < faceSize; ++f) {
1679: if (face[f] == cone[v]) {
1680: found = PETSC_TRUE; break;
1681: }
1682: }
1683: if (found) {
1684: indices[v2] = v;
1685: origVertices[v2] = cone[v];
1686: ++v2;
1687: } else {
1688: oppositeVertex = v;
1689: }
1690: }
1691: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
1692: return(0);
1693: }
1697: /*
1698: DMPlexInsertFace_Internal - Puts a face into the mesh
1700: Not collective
1702: Input Parameters:
1703: + dm - The DMPlex
1704: . numFaceVertex - The number of vertices in the face
1705: . faceVertices - The vertices in the face for dm
1706: . subfaceVertices - The vertices in the face for subdm
1707: . numCorners - The number of vertices in the cell
1708: . cell - A cell in dm containing the face
1709: . subcell - A cell in subdm containing the face
1710: . firstFace - First face in the mesh
1711: - newFacePoint - Next face in the mesh
1713: Output Parameters:
1714: . newFacePoint - Contains next face point number on input, updated on output
1716: Level: developer
1717: */
1718: 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)
1719: {
1720: MPI_Comm comm;
1721: DM_Plex *submesh = (DM_Plex*) subdm->data;
1722: const PetscInt *faces;
1723: PetscInt numFaces, coneSize;
1724: PetscErrorCode ierr;
1727: PetscObjectGetComm((PetscObject)dm,&comm);
1728: DMPlexGetConeSize(subdm, subcell, &coneSize);
1729: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
1730: #if 0
1731: /* Cannot use this because support() has not been constructed yet */
1732: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
1733: #else
1734: {
1735: PetscInt f;
1737: numFaces = 0;
1738: DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
1739: for (f = firstFace; f < *newFacePoint; ++f) {
1740: PetscInt dof, off, d;
1742: PetscSectionGetDof(submesh->coneSection, f, &dof);
1743: PetscSectionGetOffset(submesh->coneSection, f, &off);
1744: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
1745: for (d = 0; d < dof; ++d) {
1746: const PetscInt p = submesh->cones[off+d];
1747: PetscInt v;
1749: for (v = 0; v < numFaceVertices; ++v) {
1750: if (subfaceVertices[v] == p) break;
1751: }
1752: if (v == numFaceVertices) break;
1753: }
1754: if (d == dof) {
1755: numFaces = 1;
1756: ((PetscInt*) faces)[0] = f;
1757: }
1758: }
1759: }
1760: #endif
1761: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
1762: else if (numFaces == 1) {
1763: /* Add the other cell neighbor for this face */
1764: DMPlexSetCone(subdm, subcell, faces);
1765: } else {
1766: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
1767: PetscBool posOriented;
1769: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
1770: origVertices = &orientedVertices[numFaceVertices];
1771: indices = &orientedVertices[numFaceVertices*2];
1772: orientedSubVertices = &orientedVertices[numFaceVertices*3];
1773: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
1774: /* TODO: I know that routine should return a permutation, not the indices */
1775: for (v = 0; v < numFaceVertices; ++v) {
1776: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
1777: for (ov = 0; ov < numFaceVertices; ++ov) {
1778: if (orientedVertices[ov] == vertex) {
1779: orientedSubVertices[ov] = subvertex;
1780: break;
1781: }
1782: }
1783: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
1784: }
1785: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
1786: DMPlexSetCone(subdm, subcell, newFacePoint);
1787: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
1788: ++(*newFacePoint);
1789: }
1790: #if 0
1791: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
1792: #else
1793: DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
1794: #endif
1795: return(0);
1796: }
1800: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, const char vertexLabelName[], PetscInt value, DM subdm)
1801: {
1802: MPI_Comm comm;
1803: DMLabel vertexLabel, subpointMap;
1804: IS subvertexIS, subcellIS;
1805: const PetscInt *subVertices, *subCells;
1806: PetscInt numSubVertices, firstSubVertex, numSubCells;
1807: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
1808: PetscInt vStart, vEnd, c, f;
1809: PetscErrorCode ierr;
1812: PetscObjectGetComm((PetscObject)dm,&comm);
1813: /* Create subpointMap which marks the submesh */
1814: DMLabelCreate("subpoint_map", &subpointMap);
1815: DMPlexSetSubpointMap(subdm, subpointMap);
1816: DMLabelDestroy(&subpointMap);
1817: DMPlexGetLabel(dm, vertexLabelName, &vertexLabel);
1818: DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);
1819: /* Setup chart */
1820: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
1821: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
1822: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
1823: DMPlexSetVTKCellHeight(subdm, 1);
1824: /* Set cone sizes */
1825: firstSubVertex = numSubCells;
1826: firstSubFace = numSubCells+numSubVertices;
1827: newFacePoint = firstSubFace;
1828: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
1829: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
1830: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
1831: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
1832: for (c = 0; c < numSubCells; ++c) {
1833: DMPlexSetConeSize(subdm, c, 1);
1834: }
1835: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
1836: DMPlexSetConeSize(subdm, f, nFV);
1837: }
1838: DMSetUp(subdm);
1839: /* Create face cones */
1840: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1841: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
1842: DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
1843: for (c = 0; c < numSubCells; ++c) {
1844: const PetscInt cell = subCells[c];
1845: const PetscInt subcell = c;
1846: PetscInt *closure = NULL;
1847: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
1849: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1850: for (cl = 0; cl < closureSize*2; cl += 2) {
1851: const PetscInt point = closure[cl];
1852: PetscInt subVertex;
1854: if ((point >= vStart) && (point < vEnd)) {
1855: ++numCorners;
1856: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
1857: if (subVertex >= 0) {
1858: closure[faceSize] = point;
1859: subface[faceSize] = firstSubVertex+subVertex;
1860: ++faceSize;
1861: }
1862: }
1863: }
1864: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1865: if (faceSize == nFV) {
1866: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
1867: }
1868: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1869: }
1870: DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
1871: DMPlexSymmetrize(subdm);
1872: DMPlexStratify(subdm);
1873: /* Build coordinates */
1874: {
1875: PetscSection coordSection, subCoordSection;
1876: Vec coordinates, subCoordinates;
1877: PetscScalar *coords, *subCoords;
1878: PetscInt numComp, coordSize, v;
1880: DMPlexGetCoordinateSection(dm, &coordSection);
1881: DMGetCoordinatesLocal(dm, &coordinates);
1882: DMPlexGetCoordinateSection(subdm, &subCoordSection);
1883: PetscSectionSetNumFields(subCoordSection, 1);
1884: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
1885: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
1886: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
1887: for (v = 0; v < numSubVertices; ++v) {
1888: const PetscInt vertex = subVertices[v];
1889: const PetscInt subvertex = firstSubVertex+v;
1890: PetscInt dof;
1892: PetscSectionGetDof(coordSection, vertex, &dof);
1893: PetscSectionSetDof(subCoordSection, subvertex, dof);
1894: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
1895: }
1896: PetscSectionSetUp(subCoordSection);
1897: PetscSectionGetStorageSize(subCoordSection, &coordSize);
1898: VecCreate(comm, &subCoordinates);
1899: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
1900: VecSetFromOptions(subCoordinates);
1901: if (coordSize) {
1902: VecGetArray(coordinates, &coords);
1903: VecGetArray(subCoordinates, &subCoords);
1904: for (v = 0; v < numSubVertices; ++v) {
1905: const PetscInt vertex = subVertices[v];
1906: const PetscInt subvertex = firstSubVertex+v;
1907: PetscInt dof, off, sdof, soff, d;
1909: PetscSectionGetDof(coordSection, vertex, &dof);
1910: PetscSectionGetOffset(coordSection, vertex, &off);
1911: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
1912: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
1913: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
1914: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
1915: }
1916: VecRestoreArray(coordinates, &coords);
1917: VecRestoreArray(subCoordinates, &subCoords);
1918: }
1919: DMSetCoordinatesLocal(subdm, subCoordinates);
1920: VecDestroy(&subCoordinates);
1921: }
1922: /* Cleanup */
1923: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
1924: ISDestroy(&subvertexIS);
1925: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
1926: ISDestroy(&subcellIS);
1927: return(0);
1928: }
1932: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, const char vertexLabelName[], PetscInt value, DM subdm)
1933: {
1934: MPI_Comm comm;
1935: DMLabel subpointMap, vertexLabel;
1936: IS *subpointIS;
1937: const PetscInt **subpoints;
1938: PetscInt *numSubPoints, *firstSubPoint, *coneNew;
1939: PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v;
1940: PetscErrorCode ierr;
1943: PetscObjectGetComm((PetscObject)dm,&comm);
1944: /* Create subpointMap which marks the submesh */
1945: DMLabelCreate("subpoint_map", &subpointMap);
1946: DMPlexSetSubpointMap(subdm, subpointMap);
1947: DMLabelDestroy(&subpointMap);
1948: DMPlexGetLabel(dm, vertexLabelName, &vertexLabel);
1949: DMPlexMarkSubmesh_Interpolated(dm, vertexLabel, value, subpointMap, subdm);
1950: /* Setup chart */
1951: DMPlexGetDimension(dm, &dim);
1952: PetscMalloc4(dim+1,PetscInt,&numSubPoints,dim+1,PetscInt,&firstSubPoint,dim+1,IS,&subpointIS,dim+1,const PetscInt *,&subpoints);
1953: for (d = 0; d <= dim; ++d) {
1954: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
1955: totSubPoints += numSubPoints[d];
1956: }
1957: DMPlexSetChart(subdm, 0, totSubPoints);
1958: DMPlexSetVTKCellHeight(subdm, 1);
1959: /* Set cone sizes */
1960: firstSubPoint[dim] = 0;
1961: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
1962: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
1963: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
1964: for (d = 0; d <= dim; ++d) {
1965: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
1966: ISGetIndices(subpointIS[d], &subpoints[d]);
1967: }
1968: for (d = 0; d <= dim; ++d) {
1969: for (p = 0; p < numSubPoints[d]; ++p) {
1970: const PetscInt point = subpoints[d][p];
1971: const PetscInt subpoint = firstSubPoint[d] + p;
1972: const PetscInt *cone;
1973: PetscInt coneSize, coneSizeNew, c, val;
1975: DMPlexGetConeSize(dm, point, &coneSize);
1976: DMPlexSetConeSize(subdm, subpoint, coneSize);
1977: if (d == dim) {
1978: DMPlexGetCone(dm, point, &cone);
1979: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
1980: DMLabelGetValue(subpointMap, cone[c], &val);
1981: if (val >= 0) coneSizeNew++;
1982: }
1983: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
1984: }
1985: }
1986: }
1987: DMSetUp(subdm);
1988: /* Set cones */
1989: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
1990: PetscMalloc(maxConeSize * sizeof(PetscInt), &coneNew);
1991: for (d = 0; d <= dim; ++d) {
1992: for (p = 0; p < numSubPoints[d]; ++p) {
1993: const PetscInt point = subpoints[d][p];
1994: const PetscInt subpoint = firstSubPoint[d] + p;
1995: const PetscInt *cone;
1996: PetscInt coneSize, subconeSize, coneSizeNew, c, subc;
1998: DMPlexGetConeSize(dm, point, &coneSize);
1999: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2000: DMPlexGetCone(dm, point, &cone);
2001: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2002: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2003: if (subc >= 0) coneNew[coneSizeNew++] = firstSubPoint[d-1] + subc;
2004: }
2005: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2006: DMPlexSetCone(subdm, subpoint, coneNew);
2007: }
2008: }
2009: PetscFree(coneNew);
2010: DMPlexSymmetrize(subdm);
2011: DMPlexStratify(subdm);
2012: /* Build coordinates */
2013: {
2014: PetscSection coordSection, subCoordSection;
2015: Vec coordinates, subCoordinates;
2016: PetscScalar *coords, *subCoords;
2017: PetscInt numComp, coordSize;
2019: DMPlexGetCoordinateSection(dm, &coordSection);
2020: DMGetCoordinatesLocal(dm, &coordinates);
2021: DMPlexGetCoordinateSection(subdm, &subCoordSection);
2022: PetscSectionSetNumFields(subCoordSection, 1);
2023: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2024: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2025: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2026: for (v = 0; v < numSubPoints[0]; ++v) {
2027: const PetscInt vertex = subpoints[0][v];
2028: const PetscInt subvertex = firstSubPoint[0]+v;
2029: PetscInt dof;
2031: PetscSectionGetDof(coordSection, vertex, &dof);
2032: PetscSectionSetDof(subCoordSection, subvertex, dof);
2033: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2034: }
2035: PetscSectionSetUp(subCoordSection);
2036: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2037: VecCreate(comm, &subCoordinates);
2038: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2039: VecSetFromOptions(subCoordinates);
2040: VecGetArray(coordinates, &coords);
2041: VecGetArray(subCoordinates, &subCoords);
2042: for (v = 0; v < numSubPoints[0]; ++v) {
2043: const PetscInt vertex = subpoints[0][v];
2044: const PetscInt subvertex = firstSubPoint[0]+v;
2045: PetscInt dof, off, sdof, soff, d;
2047: PetscSectionGetDof(coordSection, vertex, &dof);
2048: PetscSectionGetOffset(coordSection, vertex, &off);
2049: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2050: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2051: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2052: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2053: }
2054: VecRestoreArray(coordinates, &coords);
2055: VecRestoreArray(subCoordinates, &subCoords);
2056: DMSetCoordinatesLocal(subdm, subCoordinates);
2057: VecDestroy(&subCoordinates);
2058: }
2059: /* Cleanup */
2060: for (d = 0; d <= dim; ++d) {
2061: ISRestoreIndices(subpointIS[d], &subpoints[d]);
2062: ISDestroy(&subpointIS[d]);
2063: }
2064: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
2065: return(0);
2066: }
2070: /*@C
2071: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
2073: Input Parameters:
2074: + dm - The original mesh
2075: . vertexLabel - The DMLabel marking vertices contained in the surface
2076: - value - The label value to use
2078: Output Parameter:
2079: . subdm - The surface mesh
2081: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
2083: Level: developer
2085: .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue()
2086: @*/
2087: PetscErrorCode DMPlexCreateSubmesh(DM dm, const char vertexLabel[], PetscInt value, DM *subdm)
2088: {
2089: PetscInt dim, depth;
2096: DMPlexGetDimension(dm, &dim);
2097: DMPlexGetDepth(dm, &depth);
2098: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
2099: DMSetType(*subdm, DMPLEX);
2100: DMPlexSetDimension(*subdm, dim-1);
2101: if (depth == dim) {
2102: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
2103: } else {
2104: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
2105: }
2106: return(0);
2107: }
2111: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, DM subdm)
2112: {
2113: MPI_Comm comm;
2114: DMLabel subpointMap;
2115: IS subvertexIS;
2116: const PetscInt *subVertices;
2117: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells;
2118: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
2119: PetscInt cMax, c, f;
2120: PetscErrorCode ierr;
2123: PetscObjectGetComm((PetscObject)dm, &comm);
2124: /* Create subpointMap which marks the submesh */
2125: DMLabelCreate("subpoint_map", &subpointMap);
2126: DMPlexSetSubpointMap(subdm, subpointMap);
2127: DMLabelDestroy(&subpointMap);
2128: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
2129: /* Setup chart */
2130: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2131: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2132: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2133: DMPlexSetVTKCellHeight(subdm, 1);
2134: /* Set cone sizes */
2135: firstSubVertex = numSubCells;
2136: firstSubFace = numSubCells+numSubVertices;
2137: newFacePoint = firstSubFace;
2138: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2139: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2140: for (c = 0; c < numSubCells; ++c) {
2141: DMPlexSetConeSize(subdm, c, 1);
2142: }
2143: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2144: DMPlexSetConeSize(subdm, f, nFV);
2145: }
2146: DMSetUp(subdm);
2147: /* Create face cones */
2148: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2149: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2150: DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2151: for (c = 0; c < numSubCells; ++c) {
2152: const PetscInt cell = subCells[c];
2153: const PetscInt subcell = c;
2154: const PetscInt *cone, *cells;
2155: PetscInt numCells, subVertex, p, v;
2157: if (cell < cMax) continue;
2158: DMPlexGetCone(dm, cell, &cone);
2159: for (v = 0; v < nFV; ++v) {
2160: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
2161: subface[v] = firstSubVertex+subVertex;
2162: }
2163: DMPlexSetCone(subdm, newFacePoint, subface);
2164: DMPlexSetCone(subdm, subcell, &newFacePoint);
2165: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
2166: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2167: for (p = 0; p < numCells; ++p) {
2168: PetscInt negsubcell;
2170: if (cells[p] >= cMax) continue;
2171: /* I know this is a crap search */
2172: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
2173: if (subCells[negsubcell] == cells[p]) break;
2174: }
2175: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
2176: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
2177: }
2178: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
2179: ++newFacePoint;
2180: }
2181: DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2182: DMPlexSymmetrize(subdm);
2183: DMPlexStratify(subdm);
2184: /* Build coordinates */
2185: {
2186: PetscSection coordSection, subCoordSection;
2187: Vec coordinates, subCoordinates;
2188: PetscScalar *coords, *subCoords;
2189: PetscInt numComp, coordSize, v;
2191: DMPlexGetCoordinateSection(dm, &coordSection);
2192: DMGetCoordinatesLocal(dm, &coordinates);
2193: DMPlexGetCoordinateSection(subdm, &subCoordSection);
2194: PetscSectionSetNumFields(subCoordSection, 1);
2195: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2196: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2197: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2198: for (v = 0; v < numSubVertices; ++v) {
2199: const PetscInt vertex = subVertices[v];
2200: const PetscInt subvertex = firstSubVertex+v;
2201: PetscInt dof;
2203: PetscSectionGetDof(coordSection, vertex, &dof);
2204: PetscSectionSetDof(subCoordSection, subvertex, dof);
2205: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2206: }
2207: PetscSectionSetUp(subCoordSection);
2208: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2209: VecCreate(comm, &subCoordinates);
2210: if (coordSize) {
2211: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2212: VecSetFromOptions(subCoordinates);
2213: VecGetArray(coordinates, &coords);
2214: VecGetArray(subCoordinates, &subCoords);
2215: for (v = 0; v < numSubVertices; ++v) {
2216: const PetscInt vertex = subVertices[v];
2217: const PetscInt subvertex = firstSubVertex+v;
2218: PetscInt dof, off, sdof, soff, d;
2220: PetscSectionGetDof(coordSection, vertex, &dof);
2221: PetscSectionGetOffset(coordSection, vertex, &off);
2222: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2223: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2224: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2225: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2226: }
2227: VecRestoreArray(coordinates, &coords);
2228: VecRestoreArray(subCoordinates, &subCoords);
2229: }
2230: DMSetCoordinatesLocal(subdm, subCoordinates);
2231: VecDestroy(&subCoordinates);
2232: }
2233: /* Cleanup */
2234: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2235: ISDestroy(&subvertexIS);
2236: PetscFree(subCells);
2237: return(0);
2238: }
2242: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, PetscBool hasLagrange, DM subdm)
2243: {
2244: MPI_Comm comm;
2245: DMLabel subpointMap;
2246: IS *subpointIS;
2247: const PetscInt **subpoints;
2248: PetscInt *numSubPoints, *firstSubPoint, *coneNew;
2249: PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v;
2250: PetscErrorCode ierr;
2253: PetscObjectGetComm((PetscObject)dm,&comm);
2254: /* Create subpointMap which marks the submesh */
2255: DMLabelCreate("subpoint_map", &subpointMap);
2256: DMPlexSetSubpointMap(subdm, subpointMap);
2257: DMLabelDestroy(&subpointMap);
2258: DMPlexMarkCohesiveSubmesh_Interpolated(dm, hasLagrange, subpointMap, subdm);
2259: /* Setup chart */
2260: DMPlexGetDimension(dm, &dim);
2261: PetscMalloc4(dim+1,PetscInt,&numSubPoints,dim+1,PetscInt,&firstSubPoint,dim+1,IS,&subpointIS,dim+1,const PetscInt *,&subpoints);
2262: for (d = 0; d <= dim; ++d) {
2263: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2264: totSubPoints += numSubPoints[d];
2265: }
2266: DMPlexSetChart(subdm, 0, totSubPoints);
2267: DMPlexSetVTKCellHeight(subdm, 1);
2268: /* Set cone sizes */
2269: firstSubPoint[dim] = 0;
2270: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
2271: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
2272: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2273: for (d = 0; d <= dim; ++d) {
2274: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2275: ISGetIndices(subpointIS[d], &subpoints[d]);
2276: }
2277: for (d = 0; d <= dim; ++d) {
2278: for (p = 0; p < numSubPoints[d]; ++p) {
2279: const PetscInt point = subpoints[d][p];
2280: const PetscInt subpoint = firstSubPoint[d] + p;
2281: const PetscInt *cone;
2282: PetscInt coneSize, coneSizeNew, c, val;
2284: DMPlexGetConeSize(dm, point, &coneSize);
2285: DMPlexSetConeSize(subdm, subpoint, coneSize);
2286: if (d == dim) {
2287: DMPlexGetCone(dm, point, &cone);
2288: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2289: DMLabelGetValue(subpointMap, cone[c], &val);
2290: if (val >= 0) coneSizeNew++;
2291: }
2292: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2293: }
2294: }
2295: }
2296: DMSetUp(subdm);
2297: /* Set cones */
2298: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2299: PetscMalloc(maxConeSize * sizeof(PetscInt), &coneNew);
2300: for (d = 0; d <= dim; ++d) {
2301: for (p = 0; p < numSubPoints[d]; ++p) {
2302: const PetscInt point = subpoints[d][p];
2303: const PetscInt subpoint = firstSubPoint[d] + p;
2304: const PetscInt *cone;
2305: PetscInt coneSize, subconeSize, coneSizeNew, c, subc;
2307: DMPlexGetConeSize(dm, point, &coneSize);
2308: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2309: DMPlexGetCone(dm, point, &cone);
2310: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2311: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2312: if (subc >= 0) coneNew[coneSizeNew++] = firstSubPoint[d-1] + subc;
2313: }
2314: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2315: DMPlexSetCone(subdm, subpoint, coneNew);
2316: }
2317: }
2318: PetscFree(coneNew);
2319: DMPlexSymmetrize(subdm);
2320: DMPlexStratify(subdm);
2321: /* Build coordinates */
2322: {
2323: PetscSection coordSection, subCoordSection;
2324: Vec coordinates, subCoordinates;
2325: PetscScalar *coords, *subCoords;
2326: PetscInt numComp, coordSize;
2328: DMPlexGetCoordinateSection(dm, &coordSection);
2329: DMGetCoordinatesLocal(dm, &coordinates);
2330: DMPlexGetCoordinateSection(subdm, &subCoordSection);
2331: PetscSectionSetNumFields(subCoordSection, 1);
2332: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2333: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2334: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2335: for (v = 0; v < numSubPoints[0]; ++v) {
2336: const PetscInt vertex = subpoints[0][v];
2337: const PetscInt subvertex = firstSubPoint[0]+v;
2338: PetscInt dof;
2340: PetscSectionGetDof(coordSection, vertex, &dof);
2341: PetscSectionSetDof(subCoordSection, subvertex, dof);
2342: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2343: }
2344: PetscSectionSetUp(subCoordSection);
2345: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2346: VecCreate(comm, &subCoordinates);
2347: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2348: VecSetFromOptions(subCoordinates);
2349: VecGetArray(coordinates, &coords);
2350: VecGetArray(subCoordinates, &subCoords);
2351: for (v = 0; v < numSubPoints[0]; ++v) {
2352: const PetscInt vertex = subpoints[0][v];
2353: const PetscInt subvertex = firstSubPoint[0]+v;
2354: PetscInt dof, off, sdof, soff, d;
2356: PetscSectionGetDof(coordSection, vertex, &dof);
2357: PetscSectionGetOffset(coordSection, vertex, &off);
2358: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2359: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2360: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2361: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2362: }
2363: VecRestoreArray(coordinates, &coords);
2364: VecRestoreArray(subCoordinates, &subCoords);
2365: DMSetCoordinatesLocal(subdm, subCoordinates);
2366: VecDestroy(&subCoordinates);
2367: }
2368: /* Cleanup */
2369: for (d = 0; d <= dim; ++d) {
2370: ISRestoreIndices(subpointIS[d], &subpoints[d]);
2371: ISDestroy(&subpointIS[d]);
2372: }
2373: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
2374: return(0);
2375: }
2379: /*
2380: DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells.
2382: Input Parameters:
2383: + dm - The original mesh
2384: - hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
2386: Output Parameter:
2387: . subdm - The surface mesh
2389: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
2391: Level: developer
2393: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
2394: */
2395: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, DM *subdm)
2396: {
2397: PetscInt dim, depth;
2403: DMPlexGetDimension(dm, &dim);
2404: DMPlexGetDepth(dm, &depth);
2405: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
2406: DMSetType(*subdm, DMPLEX);
2407: DMPlexSetDimension(*subdm, dim-1);
2408: if (depth == dim) {
2409: DMPlexCreateCohesiveSubmesh_Interpolated(dm, hasLagrange, *subdm);
2410: } else {
2411: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, *subdm);
2412: }
2413: return(0);
2414: }
2418: /*@
2419: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
2421: Input Parameter:
2422: . dm - The submesh DM
2424: Output Parameter:
2425: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
2427: Level: developer
2429: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
2430: @*/
2431: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
2432: {
2433: DM_Plex *mesh = (DM_Plex*) dm->data;
2438: *subpointMap = mesh->subpointMap;
2439: return(0);
2440: }
2444: /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */
2445: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
2446: {
2447: DM_Plex *mesh = (DM_Plex *) dm->data;
2448: DMLabel tmp;
2453: tmp = mesh->subpointMap;
2454: mesh->subpointMap = subpointMap;
2455: ++mesh->subpointMap->refct;
2456: DMLabelDestroy(&tmp);
2457: return(0);
2458: }
2462: /*@
2463: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
2465: Input Parameter:
2466: . dm - The submesh DM
2468: Output Parameter:
2469: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
2471: Note: This is IS is guaranteed to be sorted by the construction of the submesh
2473: Level: developer
2475: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
2476: @*/
2477: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
2478: {
2479: MPI_Comm comm;
2480: DMLabel subpointMap;
2481: IS is;
2482: const PetscInt *opoints;
2483: PetscInt *points, *depths;
2484: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
2485: PetscErrorCode ierr;
2490: PetscObjectGetComm((PetscObject)dm,&comm);
2491: *subpointIS = NULL;
2492: DMPlexGetSubpointMap(dm, &subpointMap);
2493: if (subpointMap) {
2494: DMPlexGetDepth(dm, &depth);
2495: DMPlexGetChart(dm, &pStart, &pEnd);
2496: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
2497: DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
2498: depths[0] = depth;
2499: depths[1] = 0;
2500: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
2501: PetscMalloc(pEnd * sizeof(PetscInt), &points);
2502: for(d = 0, off = 0; d <= depth; ++d) {
2503: const PetscInt dep = depths[d];
2505: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
2506: DMLabelGetStratumSize(subpointMap, dep, &n);
2507: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
2508: 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);
2509: } else {
2510: if (!n) {
2511: if (d == 0) {
2512: /* Missing cells */
2513: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
2514: } else {
2515: /* Missing faces */
2516: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
2517: }
2518: }
2519: }
2520: if (n) {
2521: DMLabelGetStratumIS(subpointMap, dep, &is);
2522: ISGetIndices(is, &opoints);
2523: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
2524: ISRestoreIndices(is, &opoints);
2525: ISDestroy(&is);
2526: }
2527: }
2528: DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
2529: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
2530: ISCreateGeneral(comm, pEnd, points, PETSC_OWN_POINTER, subpointIS);
2531: }
2532: return(0);
2533: }