Actual source code: plexreorder.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/matorderimpl.h>
3: #include <petsc/private/dmlabelimpl.h>
5: static PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm)
6: {
7: PetscInt *perm, *iperm;
8: PetscInt depth, d, pStart, pEnd, fStart, fMax, fEnd, p;
10: PetscFunctionBegin;
11: PetscCall(DMPlexGetDepth(dm, &depth));
12: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
13: PetscCall(PetscMalloc1(pEnd - pStart, &perm));
14: PetscCall(PetscMalloc1(pEnd - pStart, &iperm));
15: for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
16: for (d = depth; d > 0; --d) {
17: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
18: PetscCall(DMPlexGetDepthStratum(dm, d - 1, &fStart, &fEnd));
19: fMax = fStart;
20: for (p = pStart; p < pEnd; ++p) {
21: const PetscInt *cone;
22: PetscInt point, coneSize, c;
24: if (d == depth) {
25: perm[p] = pperm[p];
26: iperm[pperm[p]] = p;
27: }
28: point = perm[p];
29: PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
30: PetscCall(DMPlexGetCone(dm, point, &cone));
31: for (c = 0; c < coneSize; ++c) {
32: const PetscInt oldc = cone[c];
33: const PetscInt newc = iperm[oldc];
35: if (newc < 0) {
36: perm[fMax] = oldc;
37: iperm[oldc] = fMax++;
38: }
39: }
40: }
41: PetscCheck(fMax == fEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %" PetscInt_FMT " faces %" PetscInt_FMT " does not match permuted number %" PetscInt_FMT, d, fEnd - fStart, fMax - fStart);
42: }
43: *clperm = perm;
44: *invclperm = iperm;
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
48: /*@
49: DMPlexGetOrdering - Calculate a reordering of the mesh
51: Collective
53: Input Parameters:
54: + dm - The `DMPLEX` object
55: . otype - type of reordering, see `MatOrderingType`; `DMPLEXCURVEMORTON` selects a space-filling curve
56: - label - [Optional] Label used to segregate ordering into sets, or `NULL`
58: Output Parameter:
59: . perm - The point permutation as an `IS`, `perm`[old point number] = new point number
61: Level: intermediate
63: Note:
64: The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which
65: has different types of cells, and then loop over each set of reordered cells for assembly.
67: Passing `DMPLEXCURVEMORTON` orders the cells along a Morton (Z-order) curve computed from the cell
68: centroids. This requires the `DMPLEX` to have coordinates. It needs no adjacency graph. Every other
69: value currently gives reverse Cuthill-McKee.
71: .seealso: `DMPLEX`, `DMPlexPermute()`, `MatOrderingType`, `MatGetOrdering()`
72: @*/
73: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
74: {
75: PetscInt numCells = 0;
76: PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, pStart, pEnd, c, i;
77: PetscBool iscurve;
79: PetscFunctionBegin;
81: PetscAssertPointer(perm, 4);
82: PetscCall(DMPlexCurveTypeResolve_Internal(PetscObjectComm((PetscObject)dm), otype, &iscurve));
83: if (iscurve) {
84: /* A space-filling curve orders the cells from their coordinates, so it needs no adjacency graph */
85: Vec coordinates;
86: PetscInt cStart, cEnd, cdim;
88: /* The curve interleaves at most three axes. Check that before anything is allocated, so a mesh
89: in a higher-dimensional space reports the limit instead of leaking the work arrays. */
90: PetscCall(DMGetCoordinateDim(dm, &cdim));
91: PetscCheck(cdim >= 1 && cdim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension %" PetscInt_FMT " must be in [1, 3] for a space-filling curve", cdim);
92: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
93: PetscCheck(coordinates, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "A space-filling curve orders cells by centroid, but this DMPLEX has no coordinates");
94: PetscCall(DMGetCellCoordinatesLocalSetUp(dm));
95: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
96: /* DMPlexCreateOrderingClosure_Static() indexes the permutation by the cell point number */
97: PetscCheck(cStart == 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cells must start at point 0, not %" PetscInt_FMT, cStart);
98: numCells = cEnd - cStart;
99: PetscCall(PetscMalloc1(numCells, &cperm));
100: PetscCall(DMPlexGetCellOrderingByCurve_Internal(dm, otype, cStart, cEnd, cperm));
101: } else {
102: PetscInt *mask, *xls;
104: PetscCall(DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency));
105: PetscCall(PetscMalloc1(numCells, &cperm));
106: PetscCall(PetscMalloc2(numCells, &mask, numCells * 2, &xls));
107: if (numCells) {
108: /* Shift for Fortran numbering */
109: for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
110: for (i = 0; i <= numCells; ++i) ++start[i];
111: PetscCall(SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls));
112: }
113: PetscCall(PetscFree2(mask, xls));
114: PetscCall(PetscFree(start));
115: PetscCall(PetscFree(adjacency));
116: /* Shift for Fortran numbering */
117: for (c = 0; c < numCells; ++c) --cperm[c];
118: }
119: /* Segregate */
120: if (label) {
121: IS valueIS;
122: const PetscInt *valuesTmp;
123: PetscInt *values;
124: PetscInt numValues, numPoints = 0;
125: PetscInt *sperm, *vsize, *voff, v;
127: // Can't directly sort the valueIS, since it is a view into the DMLabel
128: PetscCall(DMLabelGetValueIS(label, &valueIS));
129: PetscCall(ISGetLocalSize(valueIS, &numValues));
130: PetscCall(ISGetIndices(valueIS, &valuesTmp));
131: PetscCall(PetscCalloc4(numCells, &sperm, numValues, &values, numValues, &vsize, numValues + 1, &voff));
132: PetscCall(PetscArraycpy(values, valuesTmp, numValues));
133: PetscCall(PetscSortInt(numValues, values));
134: PetscCall(ISRestoreIndices(valueIS, &valuesTmp));
135: PetscCall(ISDestroy(&valueIS));
136: for (v = 0; v < numValues; ++v) {
137: PetscCall(DMLabelGetStratumSize(label, values[v], &vsize[v]));
138: if (v < numValues - 1) voff[v + 2] += vsize[v] + voff[v + 1];
139: numPoints += vsize[v];
140: }
141: PetscCheck(numPoints == numCells, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %" PetscInt_FMT " cells != %" PetscInt_FMT " total cells", numPoints, numCells);
142: for (c = 0; c < numCells; ++c) {
143: const PetscInt oldc = cperm[c];
144: PetscInt val, vloc;
146: PetscCall(DMLabelGetValue(label, oldc, &val));
147: PetscCheck(val != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " not present in label", oldc);
148: PetscCall(PetscFindInt(val, numValues, values, &vloc));
149: PetscCheck(vloc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " not present label", val);
150: sperm[voff[vloc + 1]++] = oldc;
151: }
152: for (v = 0; v < numValues; ++v) PetscCheck(voff[v + 1] - voff[v] == vsize[v], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of %" PetscInt_FMT " values found is %" PetscInt_FMT " != %" PetscInt_FMT, values[v], voff[v + 1] - voff[v], vsize[v]);
153: PetscCall(PetscArraycpy(cperm, sperm, numCells));
154: PetscCall(PetscFree4(sperm, values, vsize, voff));
155: }
156: /* Construct closure */
157: PetscCall(DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm));
158: PetscCall(PetscFree(cperm));
159: PetscCall(PetscFree(clperm));
160: /* Invert permutation */
161: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
162: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, invclperm, PETSC_OWN_POINTER, perm));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: /*@
167: DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line
169: Collective
171: Input Parameter:
172: . dm - The `DMPLEX` object
174: Output Parameter:
175: . perm - The point permutation as an `IS`, `perm`[old point number] = new point number
177: Level: intermediate
179: .seealso: `DMPLEX`, `DMPlexGetOrdering()`, `DMPlexPermute()`, `MatGetOrdering()`
180: @*/
181: PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm)
182: {
183: PetscInt *points;
184: const PetscInt *support, *cone;
185: PetscInt dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0;
187: PetscFunctionBegin;
188: PetscCall(DMGetDimension(dm, &dim));
189: PetscCheck(dim == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Input mesh must be one dimensional, not %" PetscInt_FMT, dim);
190: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
191: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
192: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
193: PetscCall(PetscMalloc1(pEnd - pStart, &points));
194: for (c = cStart; c < cEnd; ++c) points[c] = c;
195: for (v = vStart; v < vEnd; ++v) points[v] = v;
196: for (v = vStart; v < vEnd; ++v) {
197: PetscCall(DMPlexGetSupportSize(dm, v, &suppSize));
198: PetscCall(DMPlexGetSupport(dm, v, &support));
199: if (suppSize == 1) {
200: lastCell = support[0];
201: break;
202: }
203: }
204: if (v < vEnd) {
205: PetscInt pos = cEnd;
207: points[v] = pos++;
208: while (lastCell >= cStart) {
209: PetscCall(DMPlexGetCone(dm, lastCell, &cone));
210: if (cone[0] == v) v = cone[1];
211: else v = cone[0];
212: PetscCall(DMPlexGetSupport(dm, v, &support));
213: PetscCall(DMPlexGetSupportSize(dm, v, &suppSize));
214: if (suppSize == 1) {
215: lastCell = -1;
216: } else {
217: if (support[0] == lastCell) lastCell = support[1];
218: else lastCell = support[0];
219: }
220: points[v] = pos++;
221: }
222: PetscCheck(pos == pEnd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Last vertex was %" PetscInt_FMT ", not %" PetscInt_FMT, pos, pEnd);
223: }
224: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, points, PETSC_OWN_POINTER, perm));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static PetscErrorCode DMPlexRemapCoordinates_Private(IS perm, PetscSection cs, Vec coordinates, PetscSection *csNew, Vec *coordinatesNew)
229: {
230: PetscScalar *coords, *coordsNew;
231: const PetscInt *pperm;
232: PetscInt pStart, pEnd, p;
233: const char *name;
235: PetscFunctionBegin;
236: PetscCall(PetscSectionPermute(cs, perm, csNew));
237: PetscCall(VecDuplicate(coordinates, coordinatesNew));
238: PetscCall(PetscObjectGetName((PetscObject)coordinates, &name));
239: PetscCall(PetscObjectSetName((PetscObject)*coordinatesNew, name));
240: PetscCall(VecGetArray(coordinates, &coords));
241: PetscCall(VecGetArray(*coordinatesNew, &coordsNew));
242: PetscCall(PetscSectionGetChart(*csNew, &pStart, &pEnd));
243: PetscCall(ISGetIndices(perm, &pperm));
244: for (p = pStart; p < pEnd; ++p) {
245: PetscInt dof, off, offNew, d;
247: PetscCall(PetscSectionGetDof(*csNew, p, &dof));
248: PetscCall(PetscSectionGetOffset(cs, p, &off));
249: PetscCall(PetscSectionGetOffset(*csNew, pperm[p], &offNew));
250: for (d = 0; d < dof; ++d) coordsNew[offNew + d] = coords[off + d];
251: }
252: PetscCall(ISRestoreIndices(perm, &pperm));
253: PetscCall(VecRestoreArray(coordinates, &coords));
254: PetscCall(VecRestoreArray(*coordinatesNew, &coordsNew));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*@
259: DMPlexPermute - Reorder the mesh according to the input permutation
261: Collective
263: Input Parameters:
264: + dm - The `DMPLEX` object
265: - perm - The point permutation, `perm`[old point number] = new point number
267: Output Parameter:
268: . pdm - The permuted `DM`
270: Level: intermediate
272: .seealso: `DMPLEX`, `MatPermute()`
273: @*/
274: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
275: {
276: DM_Plex *plex = (DM_Plex *)dm->data, *plexNew;
277: PetscInt dim, cdim;
278: const char *name;
280: PetscFunctionBegin;
283: PetscAssertPointer(pdm, 3);
284: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), pdm));
285: PetscCall(DMSetType(*pdm, DMPLEX));
286: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
287: PetscCall(PetscObjectSetName((PetscObject)*pdm, name));
288: PetscCall(DMGetDimension(dm, &dim));
289: PetscCall(DMSetDimension(*pdm, dim));
290: PetscCall(DMGetCoordinateDim(dm, &cdim));
291: PetscCall(DMSetCoordinateDim(*pdm, cdim));
292: PetscCall(DMCopyDisc(dm, *pdm));
293: if (dm->localSection) {
294: PetscSection section, sectionNew;
296: PetscCall(DMGetLocalSection(dm, §ion));
297: PetscCall(PetscSectionPermute(section, perm, §ionNew));
298: PetscCall(DMSetLocalSection(*pdm, sectionNew));
299: PetscCall(PetscSectionDestroy(§ionNew));
300: }
301: plexNew = (DM_Plex *)(*pdm)->data;
302: /* Ignore ltogmap, ltogmapb */
303: /* Ignore sf, sectionSF */
304: /* Ignore globalVertexNumbers, globalCellNumbers */
305: /* Reorder labels */
306: {
307: PetscInt numLabels;
308: DMLabel label, labelNew;
310: PetscCall(DMGetNumLabels(dm, &numLabels));
311: for (PetscInt l = 0; l < numLabels; ++l) {
312: PetscCall(DMGetLabelByNum(dm, l, &label));
313: PetscCall(DMLabelPermute(label, perm, &labelNew));
314: PetscCall(DMAddLabel(*pdm, labelNew));
315: PetscCall(DMLabelDestroy(&labelNew));
316: }
317: PetscCall(DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel));
318: if (plex->subpointMap) PetscCall(DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap));
319: }
320: if ((*pdm)->celltypeLabel) {
321: DMLabel ctLabel;
323: // Reset label for fast lookup
324: PetscCall(DMPlexGetCellTypeLabel(*pdm, &ctLabel));
325: PetscCall(DMLabelMakeAllInvalid_Internal(ctLabel));
326: }
327: /* Reorder topology */
328: {
329: const PetscInt *pperm;
330: PetscInt n, pStart, pEnd, p;
332: PetscCall(PetscSectionDestroy(&plexNew->coneSection));
333: PetscCall(PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection));
334: PetscCall(PetscSectionGetStorageSize(plexNew->coneSection, &n));
335: PetscCall(PetscMalloc1(n, &plexNew->cones));
336: PetscCall(PetscMalloc1(n, &plexNew->coneOrientations));
337: PetscCall(ISGetIndices(perm, &pperm));
338: PetscCall(PetscSectionGetChart(plex->coneSection, &pStart, &pEnd));
339: for (p = pStart; p < pEnd; ++p) {
340: PetscInt dof, off, offNew, d;
342: PetscCall(PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof));
343: PetscCall(PetscSectionGetOffset(plex->coneSection, p, &off));
344: PetscCall(PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew));
345: for (d = 0; d < dof; ++d) {
346: plexNew->cones[offNew + d] = pperm[plex->cones[off + d]];
347: plexNew->coneOrientations[offNew + d] = plex->coneOrientations[off + d];
348: }
349: }
350: PetscCall(PetscSectionDestroy(&plexNew->supportSection));
351: PetscCall(PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection));
352: PetscCall(PetscSectionGetStorageSize(plexNew->supportSection, &n));
353: PetscCall(PetscMalloc1(n, &plexNew->supports));
354: PetscCall(PetscSectionGetChart(plex->supportSection, &pStart, &pEnd));
355: for (p = pStart; p < pEnd; ++p) {
356: PetscInt dof, off, offNew, d;
358: PetscCall(PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof));
359: PetscCall(PetscSectionGetOffset(plex->supportSection, p, &off));
360: PetscCall(PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew));
361: for (d = 0; d < dof; ++d) plexNew->supports[offNew + d] = pperm[plex->supports[off + d]];
362: }
363: PetscCall(ISRestoreIndices(perm, &pperm));
364: }
365: /* Remap coordinates */
366: {
367: DM cdm, cdmNew;
368: PetscSection cs, csNew;
369: Vec coordinates, coordinatesNew;
371: PetscCall(DMGetCoordinateSection(dm, &cs));
372: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
373: PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew));
374: PetscCall(DMSetCoordinateSection(*pdm, PETSC_DETERMINE, csNew));
375: PetscCall(DMSetCoordinatesLocal(*pdm, coordinatesNew));
376: PetscCall(PetscSectionDestroy(&csNew));
377: PetscCall(VecDestroy(&coordinatesNew));
379: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
380: if (cdm) {
381: PetscCall(DMGetCoordinateDM(*pdm, &cdm));
382: PetscCall(DMClone(cdm, &cdmNew));
383: PetscCall(DMSetCellCoordinateDM(*pdm, cdmNew));
384: PetscCall(DMDestroy(&cdmNew));
385: PetscCall(DMGetCellCoordinateSection(dm, &cs));
386: PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates));
387: PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew));
388: PetscCall(DMSetCellCoordinateSection(*pdm, PETSC_DETERMINE, csNew));
389: PetscCall(DMSetCellCoordinatesLocal(*pdm, coordinatesNew));
390: PetscCall(PetscSectionDestroy(&csNew));
391: PetscCall(VecDestroy(&coordinatesNew));
392: }
393: }
394: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *pdm));
395: (*pdm)->setupcalled = PETSC_TRUE;
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: PetscErrorCode DMPlexReorderSetDefault_Plex(DM dm, DMReorderDefaultFlag reorder)
400: {
401: DM_Plex *mesh = (DM_Plex *)dm->data;
403: PetscFunctionBegin;
404: mesh->reorderDefault = reorder;
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: DMPlexReorderSetDefault - Set flag indicating whether the DM should be reordered by default
411: Logically Collective
413: Input Parameters:
414: + dm - The `DM`
415: - reorder - Flag for reordering
417: Level: intermediate
419: .seealso: `DMPlexReorderGetDefault()`
420: @*/
421: PetscErrorCode DMPlexReorderSetDefault(DM dm, DMReorderDefaultFlag reorder)
422: {
423: PetscFunctionBegin;
425: PetscTryMethod(dm, "DMPlexReorderSetDefault_C", (DM, DMReorderDefaultFlag), (dm, reorder));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: PetscErrorCode DMPlexReorderGetDefault_Plex(DM dm, DMReorderDefaultFlag *reorder)
430: {
431: DM_Plex *mesh = (DM_Plex *)dm->data;
433: PetscFunctionBegin;
434: *reorder = mesh->reorderDefault;
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@
439: DMPlexReorderGetDefault - Get flag indicating whether the DM should be reordered by default
441: Not Collective
443: Input Parameter:
444: . dm - The `DM`
446: Output Parameter:
447: . reorder - Flag for reordering
449: Level: intermediate
451: .seealso: `DMPlexReorderSetDefault()`
452: @*/
453: PetscErrorCode DMPlexReorderGetDefault(DM dm, DMReorderDefaultFlag *reorder)
454: {
455: PetscFunctionBegin;
457: PetscAssertPointer(reorder, 2);
458: PetscUseMethod(dm, "DMPlexReorderGetDefault_C", (DM, DMReorderDefaultFlag *), (dm, reorder));
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: static PetscErrorCode DMCreateSectionPermutation_Plex_Reverse(DM dm, IS *permutation, PetscBT *blockStarts)
463: {
464: IS permIS;
465: PetscInt *perm;
466: PetscInt pStart, pEnd;
468: PetscFunctionBegin;
469: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
470: PetscCall(PetscMalloc1(pEnd - pStart, &perm));
471: for (PetscInt p = pStart; p < pEnd; ++p) perm[pEnd - 1 - p] = p;
472: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
473: PetscCall(ISSetPermutation(permIS));
474: *permutation = permIS;
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: // Reorder to group split nodes
479: static PetscErrorCode DMCreateSectionPermutation_Plex_Cohesive_Old(DM dm, IS *permutation, PetscBT *blockStarts)
480: {
481: IS permIS;
482: PetscBT bt, blst;
483: PetscInt *perm;
484: PetscInt pStart, pEnd, i = 0;
486: PetscFunctionBegin;
487: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
488: PetscCall(PetscMalloc1(pEnd - pStart, &perm));
489: PetscCall(PetscBTCreate(pEnd - pStart, &bt));
490: PetscCall(PetscBTCreate(pEnd - pStart, &blst));
491: for (PetscInt p = pStart; p < pEnd; ++p) {
492: const PetscInt *supp, *cone;
493: PetscInt suppSize;
494: DMPolytopeType ct;
496: PetscCall(DMPlexGetCellType(dm, p, &ct));
497: // Do not order tensor cells until they appear below
498: if (ct == DM_POLYTOPE_POINT_PRISM_TENSOR || ct == DM_POLYTOPE_SEG_PRISM_TENSOR || ct == DM_POLYTOPE_TRI_PRISM_TENSOR || ct == DM_POLYTOPE_QUAD_PRISM_TENSOR) continue;
499: if (PetscBTLookupSet(bt, p)) continue;
500: PetscCall(PetscBTSet(blst, p));
501: perm[i++] = p;
502: // Check for tensor cells in the support
503: PetscCall(DMPlexGetSupport(dm, p, &supp));
504: PetscCall(DMPlexGetSupportSize(dm, p, &suppSize));
505: for (PetscInt s = 0; s < suppSize; ++s) {
506: DMPolytopeType sct;
507: PetscInt q, qq;
509: PetscCall(DMPlexGetCellType(dm, supp[s], &sct));
510: switch (sct) {
511: case DM_POLYTOPE_POINT_PRISM_TENSOR:
512: case DM_POLYTOPE_SEG_PRISM_TENSOR:
513: case DM_POLYTOPE_TRI_PRISM_TENSOR:
514: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
515: // If found, move up the split partner of the tensor cell, and the cell itself
516: PetscCall(DMPlexGetCone(dm, supp[s], &cone));
517: qq = supp[s];
518: q = (cone[0] == p) ? cone[1] : cone[0];
519: if (!PetscBTLookupSet(bt, q)) {
520: perm[i++] = q;
521: s = suppSize;
522: // At T-junctions, we can have an unsplit point at the other end, so also order that loop
523: {
524: const PetscInt *qsupp, *qcone;
525: PetscInt qsuppSize;
527: PetscCall(DMPlexGetSupport(dm, q, &qsupp));
528: PetscCall(DMPlexGetSupportSize(dm, q, &qsuppSize));
529: for (PetscInt qs = 0; qs < qsuppSize; ++qs) {
530: DMPolytopeType qsct;
532: PetscCall(DMPlexGetCellType(dm, qsupp[qs], &qsct));
533: switch (qsct) {
534: case DM_POLYTOPE_POINT_PRISM_TENSOR:
535: case DM_POLYTOPE_SEG_PRISM_TENSOR:
536: case DM_POLYTOPE_TRI_PRISM_TENSOR:
537: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
538: PetscCall(DMPlexGetCone(dm, qsupp[qs], &qcone));
539: if (qcone[0] == qcone[1]) {
540: if (!PetscBTLookupSet(bt, qsupp[qs])) {
541: perm[i++] = qsupp[qs];
542: qs = qsuppSize;
543: }
544: }
545: break;
546: default:
547: break;
548: }
549: }
550: }
551: }
552: if (!PetscBTLookupSet(bt, qq)) {
553: perm[i++] = qq;
554: s = suppSize;
555: }
556: break;
557: default:
558: break;
559: }
560: }
561: }
562: if (PetscDefined(USE_DEBUG)) {
563: for (PetscInt p = pStart; p < pEnd; ++p) PetscCheck(PetscBTLookup(bt, p), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " missed in permutation of [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, pStart, pEnd);
564: }
565: PetscCall(PetscBTDestroy(&bt));
566: PetscCheck(i == pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of points in permutation %" PetscInt_FMT " does not match chart size %" PetscInt_FMT, i, pEnd - pStart);
567: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
568: PetscCall(ISSetPermutation(permIS));
569: *permutation = permIS;
570: *blockStarts = blst;
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: // Mark the block associated with a cohesive cell p
575: static PetscErrorCode InsertCohesiveBlock_Private(DM dm, PetscBT bt, PetscBT blst, PetscInt p, PetscInt *idx, PetscInt perm[])
576: {
577: const PetscInt *cone;
578: PetscInt cS;
580: PetscFunctionBegin;
581: if (PetscBTLookupSet(bt, p)) PetscFunctionReturn(PETSC_SUCCESS);
582: // Order the endcaps
583: PetscCall(DMPlexGetCone(dm, p, &cone));
584: PetscCall(DMPlexGetConeSize(dm, p, &cS));
585: if (blst) PetscCall(PetscBTSet(blst, cone[0]));
586: if (!PetscBTLookupSet(bt, cone[0])) perm[(*idx)++] = cone[0];
587: if (!PetscBTLookupSet(bt, cone[1])) perm[(*idx)++] = cone[1];
588: // Order sides
589: for (PetscInt c = 2; c < cS; ++c) PetscCall(InsertCohesiveBlock_Private(dm, bt, NULL, cone[c], idx, perm));
590: // Order cell
591: perm[(*idx)++] = p;
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: // Reorder to group split nodes
596: static PetscErrorCode DMCreateSectionPermutation_Plex_Cohesive(DM dm, IS *permutation, PetscBT *blockStarts)
597: {
598: IS permIS;
599: PetscBT bt, blst;
600: PetscInt *perm;
601: PetscInt dim, pStart, pEnd, i = 0;
603: PetscFunctionBegin;
604: PetscCall(DMGetDimension(dm, &dim));
605: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
606: PetscCall(PetscMalloc1(pEnd - pStart, &perm));
607: PetscCall(PetscBTCreate(pEnd - pStart, &bt));
608: PetscCall(PetscBTCreate(pEnd - pStart, &blst));
609: // Add cohesive blocks
610: for (PetscInt p = pStart; p < pEnd; ++p) {
611: DMPolytopeType ct;
613: PetscCall(DMPlexGetCellType(dm, p, &ct));
614: switch (dim) {
615: case 2:
616: if (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) PetscCall(InsertCohesiveBlock_Private(dm, bt, blst, p, &i, perm));
617: break;
618: case 3:
619: if (ct == DM_POLYTOPE_TRI_PRISM_TENSOR || ct == DM_POLYTOPE_QUAD_PRISM_TENSOR) PetscCall(InsertCohesiveBlock_Private(dm, bt, blst, p, &i, perm));
620: break;
621: default:
622: break;
623: }
624: }
625: // Add normal blocks
626: for (PetscInt p = pStart; p < pEnd; ++p) {
627: if (PetscBTLookupSet(bt, p)) continue;
628: PetscCall(PetscBTSet(blst, p));
629: perm[i++] = p;
630: }
631: if (PetscDefined(USE_DEBUG)) {
632: for (PetscInt p = pStart; p < pEnd; ++p) PetscCheck(PetscBTLookup(bt, p), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " missed in permutation of [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, pStart, pEnd);
633: }
634: PetscCall(PetscBTDestroy(&bt));
635: PetscCheck(i == pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of points in permutation %" PetscInt_FMT " does not match chart size %" PetscInt_FMT, i, pEnd - pStart);
636: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
637: PetscCall(ISSetPermutation(permIS));
638: *permutation = permIS;
639: *blockStarts = blst;
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: PetscErrorCode DMCreateSectionPermutation_Plex(DM dm, IS *perm, PetscBT *blockStarts)
644: {
645: DMReorderDefaultFlag reorder;
646: MatOrderingType otype;
647: PetscBool iscohesive, iscohesiveOld, isreverse;
649: PetscFunctionBegin;
650: PetscCall(DMReorderSectionGetDefault(dm, &reorder));
651: if (reorder != DM_REORDER_DEFAULT_TRUE) PetscFunctionReturn(PETSC_SUCCESS);
652: PetscCall(DMReorderSectionGetType(dm, &otype));
653: if (!otype) PetscFunctionReturn(PETSC_SUCCESS);
654: PetscCall(PetscStrncmp(otype, "cohesive_old", 1024, &iscohesiveOld));
655: PetscCall(PetscStrncmp(otype, "cohesive", 1024, &iscohesive));
656: PetscCall(PetscStrncmp(otype, "reverse", 1024, &isreverse));
657: if (iscohesive) PetscCall(DMCreateSectionPermutation_Plex_Cohesive(dm, perm, blockStarts));
658: else if (iscohesiveOld) PetscCall(DMCreateSectionPermutation_Plex_Cohesive_Old(dm, perm, blockStarts));
659: else if (isreverse) PetscCall(DMCreateSectionPermutation_Plex_Reverse(dm, perm, blockStarts));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: PetscErrorCode DMReorderSectionSetDefault_Plex(DM dm, DMReorderDefaultFlag reorder)
664: {
665: PetscFunctionBegin;
666: dm->reorderSection = reorder;
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: PetscErrorCode DMReorderSectionGetDefault_Plex(DM dm, DMReorderDefaultFlag *reorder)
671: {
672: PetscFunctionBegin;
673: *reorder = dm->reorderSection;
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: PetscErrorCode DMReorderSectionSetType_Plex(DM dm, MatOrderingType reorder)
678: {
679: PetscFunctionBegin;
680: PetscCall(PetscFree(dm->reorderSectionType));
681: PetscCall(PetscStrallocpy(reorder, (char **)&dm->reorderSectionType));
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: PetscErrorCode DMReorderSectionGetType_Plex(DM dm, MatOrderingType *reorder)
686: {
687: PetscFunctionBegin;
688: *reorder = dm->reorderSectionType;
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }