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: /*@C
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`
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: .seealso: `DMPLEX`, `DMPlexPermute()`, `MatOrderingType`, `MatGetOrdering()`
68: @*/
69: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
70: {
71: PetscInt numCells = 0;
72: PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i;
74: PetscFunctionBegin;
76: PetscAssertPointer(perm, 4);
77: PetscCall(DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency));
78: PetscCall(PetscMalloc3(numCells, &cperm, numCells, &mask, numCells * 2, &xls));
79: if (numCells) {
80: /* Shift for Fortran numbering */
81: for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
82: for (i = 0; i <= numCells; ++i) ++start[i];
83: PetscCall(SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls));
84: }
85: PetscCall(PetscFree(start));
86: PetscCall(PetscFree(adjacency));
87: /* Shift for Fortran numbering */
88: for (c = 0; c < numCells; ++c) --cperm[c];
89: /* Segregate */
90: if (label) {
91: IS valueIS;
92: const PetscInt *valuesTmp;
93: PetscInt *values;
94: PetscInt numValues, numPoints = 0;
95: PetscInt *sperm, *vsize, *voff, v;
97: // Can't directly sort the valueIS, since it is a view into the DMLabel
98: PetscCall(DMLabelGetValueIS(label, &valueIS));
99: PetscCall(ISGetLocalSize(valueIS, &numValues));
100: PetscCall(ISGetIndices(valueIS, &valuesTmp));
101: PetscCall(PetscCalloc4(numCells, &sperm, numValues, &values, numValues, &vsize, numValues + 1, &voff));
102: PetscCall(PetscArraycpy(values, valuesTmp, numValues));
103: PetscCall(PetscSortInt(numValues, values));
104: PetscCall(ISRestoreIndices(valueIS, &valuesTmp));
105: PetscCall(ISDestroy(&valueIS));
106: for (v = 0; v < numValues; ++v) {
107: PetscCall(DMLabelGetStratumSize(label, values[v], &vsize[v]));
108: if (v < numValues - 1) voff[v + 2] += vsize[v] + voff[v + 1];
109: numPoints += vsize[v];
110: }
111: PetscCheck(numPoints == numCells, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %" PetscInt_FMT " cells != %" PetscInt_FMT " total cells", numPoints, numCells);
112: for (c = 0; c < numCells; ++c) {
113: const PetscInt oldc = cperm[c];
114: PetscInt val, vloc;
116: PetscCall(DMLabelGetValue(label, oldc, &val));
117: PetscCheck(val != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " not present in label", oldc);
118: PetscCall(PetscFindInt(val, numValues, values, &vloc));
119: PetscCheck(vloc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " not present label", val);
120: sperm[voff[vloc + 1]++] = oldc;
121: }
122: 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]);
123: PetscCall(PetscArraycpy(cperm, sperm, numCells));
124: PetscCall(PetscFree4(sperm, values, vsize, voff));
125: }
126: /* Construct closure */
127: PetscCall(DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm));
128: PetscCall(PetscFree3(cperm, mask, xls));
129: PetscCall(PetscFree(clperm));
130: /* Invert permutation */
131: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
132: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, invclperm, PETSC_OWN_POINTER, perm));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: /*@
137: DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line
139: Collective
141: Input Parameter:
142: . dm - The `DMPLEX` object
144: Output Parameter:
145: . perm - The point permutation as an `IS`, `perm`[old point number] = new point number
147: Level: intermediate
149: .seealso: `DMPLEX`, `DMPlexGetOrdering()`, `DMPlexPermute()`, `MatGetOrdering()`
150: @*/
151: PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm)
152: {
153: PetscInt *points;
154: const PetscInt *support, *cone;
155: PetscInt dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0;
157: PetscFunctionBegin;
158: PetscCall(DMGetDimension(dm, &dim));
159: PetscCheck(dim == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Input mesh must be one dimensional, not %" PetscInt_FMT, dim);
160: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
161: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
162: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
163: PetscCall(PetscMalloc1(pEnd - pStart, &points));
164: for (c = cStart; c < cEnd; ++c) points[c] = c;
165: for (v = vStart; v < vEnd; ++v) points[v] = v;
166: for (v = vStart; v < vEnd; ++v) {
167: PetscCall(DMPlexGetSupportSize(dm, v, &suppSize));
168: PetscCall(DMPlexGetSupport(dm, v, &support));
169: if (suppSize == 1) {
170: lastCell = support[0];
171: break;
172: }
173: }
174: if (v < vEnd) {
175: PetscInt pos = cEnd;
177: points[v] = pos++;
178: while (lastCell >= cStart) {
179: PetscCall(DMPlexGetCone(dm, lastCell, &cone));
180: if (cone[0] == v) v = cone[1];
181: else v = cone[0];
182: PetscCall(DMPlexGetSupport(dm, v, &support));
183: PetscCall(DMPlexGetSupportSize(dm, v, &suppSize));
184: if (suppSize == 1) {
185: lastCell = -1;
186: } else {
187: if (support[0] == lastCell) lastCell = support[1];
188: else lastCell = support[0];
189: }
190: points[v] = pos++;
191: }
192: PetscCheck(pos == pEnd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Last vertex was %" PetscInt_FMT ", not %" PetscInt_FMT, pos, pEnd);
193: }
194: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, points, PETSC_OWN_POINTER, perm));
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: static PetscErrorCode DMPlexRemapCoordinates_Private(IS perm, PetscSection cs, Vec coordinates, PetscSection *csNew, Vec *coordinatesNew)
199: {
200: PetscScalar *coords, *coordsNew;
201: const PetscInt *pperm;
202: PetscInt pStart, pEnd, p;
203: const char *name;
205: PetscFunctionBegin;
206: PetscCall(PetscSectionPermute(cs, perm, csNew));
207: PetscCall(VecDuplicate(coordinates, coordinatesNew));
208: PetscCall(PetscObjectGetName((PetscObject)coordinates, &name));
209: PetscCall(PetscObjectSetName((PetscObject)*coordinatesNew, name));
210: PetscCall(VecGetArray(coordinates, &coords));
211: PetscCall(VecGetArray(*coordinatesNew, &coordsNew));
212: PetscCall(PetscSectionGetChart(*csNew, &pStart, &pEnd));
213: PetscCall(ISGetIndices(perm, &pperm));
214: for (p = pStart; p < pEnd; ++p) {
215: PetscInt dof, off, offNew, d;
217: PetscCall(PetscSectionGetDof(*csNew, p, &dof));
218: PetscCall(PetscSectionGetOffset(cs, p, &off));
219: PetscCall(PetscSectionGetOffset(*csNew, pperm[p], &offNew));
220: for (d = 0; d < dof; ++d) coordsNew[offNew + d] = coords[off + d];
221: }
222: PetscCall(ISRestoreIndices(perm, &pperm));
223: PetscCall(VecRestoreArray(coordinates, &coords));
224: PetscCall(VecRestoreArray(*coordinatesNew, &coordsNew));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*@
229: DMPlexPermute - Reorder the mesh according to the input permutation
231: Collective
233: Input Parameters:
234: + dm - The `DMPLEX` object
235: - perm - The point permutation, `perm`[old point number] = new point number
237: Output Parameter:
238: . pdm - The permuted `DM`
240: Level: intermediate
242: .seealso: `DMPLEX`, `MatPermute()`
243: @*/
244: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
245: {
246: DM_Plex *plex = (DM_Plex *)dm->data, *plexNew;
247: PetscInt dim, cdim;
248: const char *name;
250: PetscFunctionBegin;
253: PetscAssertPointer(pdm, 3);
254: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), pdm));
255: PetscCall(DMSetType(*pdm, DMPLEX));
256: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
257: PetscCall(PetscObjectSetName((PetscObject)*pdm, name));
258: PetscCall(DMGetDimension(dm, &dim));
259: PetscCall(DMSetDimension(*pdm, dim));
260: PetscCall(DMGetCoordinateDim(dm, &cdim));
261: PetscCall(DMSetCoordinateDim(*pdm, cdim));
262: PetscCall(DMCopyDisc(dm, *pdm));
263: if (dm->localSection) {
264: PetscSection section, sectionNew;
266: PetscCall(DMGetLocalSection(dm, §ion));
267: PetscCall(PetscSectionPermute(section, perm, §ionNew));
268: PetscCall(DMSetLocalSection(*pdm, sectionNew));
269: PetscCall(PetscSectionDestroy(§ionNew));
270: }
271: plexNew = (DM_Plex *)(*pdm)->data;
272: /* Ignore ltogmap, ltogmapb */
273: /* Ignore sf, sectionSF */
274: /* Ignore globalVertexNumbers, globalCellNumbers */
275: /* Reorder labels */
276: {
277: PetscInt numLabels, l;
278: DMLabel label, labelNew;
280: PetscCall(DMGetNumLabels(dm, &numLabels));
281: for (l = 0; l < numLabels; ++l) {
282: PetscCall(DMGetLabelByNum(dm, l, &label));
283: PetscCall(DMLabelPermute(label, perm, &labelNew));
284: PetscCall(DMAddLabel(*pdm, labelNew));
285: PetscCall(DMLabelDestroy(&labelNew));
286: }
287: PetscCall(DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel));
288: if (plex->subpointMap) PetscCall(DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap));
289: }
290: if ((*pdm)->celltypeLabel) {
291: DMLabel ctLabel;
293: // Reset label for fast lookup
294: PetscCall(DMPlexGetCellTypeLabel(*pdm, &ctLabel));
295: PetscCall(DMLabelMakeAllInvalid_Internal(ctLabel));
296: }
297: /* Reorder topology */
298: {
299: const PetscInt *pperm;
300: PetscInt n, pStart, pEnd, p;
302: PetscCall(PetscSectionDestroy(&plexNew->coneSection));
303: PetscCall(PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection));
304: PetscCall(PetscSectionGetStorageSize(plexNew->coneSection, &n));
305: PetscCall(PetscMalloc1(n, &plexNew->cones));
306: PetscCall(PetscMalloc1(n, &plexNew->coneOrientations));
307: PetscCall(ISGetIndices(perm, &pperm));
308: PetscCall(PetscSectionGetChart(plex->coneSection, &pStart, &pEnd));
309: for (p = pStart; p < pEnd; ++p) {
310: PetscInt dof, off, offNew, d;
312: PetscCall(PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof));
313: PetscCall(PetscSectionGetOffset(plex->coneSection, p, &off));
314: PetscCall(PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew));
315: for (d = 0; d < dof; ++d) {
316: plexNew->cones[offNew + d] = pperm[plex->cones[off + d]];
317: plexNew->coneOrientations[offNew + d] = plex->coneOrientations[off + d];
318: }
319: }
320: PetscCall(PetscSectionDestroy(&plexNew->supportSection));
321: PetscCall(PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection));
322: PetscCall(PetscSectionGetStorageSize(plexNew->supportSection, &n));
323: PetscCall(PetscMalloc1(n, &plexNew->supports));
324: PetscCall(PetscSectionGetChart(plex->supportSection, &pStart, &pEnd));
325: for (p = pStart; p < pEnd; ++p) {
326: PetscInt dof, off, offNew, d;
328: PetscCall(PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof));
329: PetscCall(PetscSectionGetOffset(plex->supportSection, p, &off));
330: PetscCall(PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew));
331: for (d = 0; d < dof; ++d) plexNew->supports[offNew + d] = pperm[plex->supports[off + d]];
332: }
333: PetscCall(ISRestoreIndices(perm, &pperm));
334: }
335: /* Remap coordinates */
336: {
337: DM cdm, cdmNew;
338: PetscSection cs, csNew;
339: Vec coordinates, coordinatesNew;
341: PetscCall(DMGetCoordinateSection(dm, &cs));
342: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
343: PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew));
344: PetscCall(DMSetCoordinateSection(*pdm, PETSC_DETERMINE, csNew));
345: PetscCall(DMSetCoordinatesLocal(*pdm, coordinatesNew));
346: PetscCall(PetscSectionDestroy(&csNew));
347: PetscCall(VecDestroy(&coordinatesNew));
349: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
350: if (cdm) {
351: PetscCall(DMGetCoordinateDM(*pdm, &cdm));
352: PetscCall(DMClone(cdm, &cdmNew));
353: PetscCall(DMSetCellCoordinateDM(*pdm, cdmNew));
354: PetscCall(DMDestroy(&cdmNew));
355: PetscCall(DMGetCellCoordinateSection(dm, &cs));
356: PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates));
357: PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew));
358: PetscCall(DMSetCellCoordinateSection(*pdm, PETSC_DETERMINE, csNew));
359: PetscCall(DMSetCellCoordinatesLocal(*pdm, coordinatesNew));
360: PetscCall(PetscSectionDestroy(&csNew));
361: PetscCall(VecDestroy(&coordinatesNew));
362: }
363: }
364: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *pdm));
365: (*pdm)->setupcalled = PETSC_TRUE;
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: PetscErrorCode DMPlexReorderSetDefault_Plex(DM dm, DMReorderDefaultFlag reorder)
370: {
371: DM_Plex *mesh = (DM_Plex *)dm->data;
373: PetscFunctionBegin;
374: mesh->reorderDefault = reorder;
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*@
379: DMPlexReorderSetDefault - Set flag indicating whether the DM should be reordered by default
381: Logically Collective
383: Input Parameters:
384: + dm - The `DM`
385: - reorder - Flag for reordering
387: Level: intermediate
389: .seealso: `DMPlexReorderGetDefault()`
390: @*/
391: PetscErrorCode DMPlexReorderSetDefault(DM dm, DMReorderDefaultFlag reorder)
392: {
393: PetscFunctionBegin;
395: PetscTryMethod(dm, "DMPlexReorderSetDefault_C", (DM, DMReorderDefaultFlag), (dm, reorder));
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: PetscErrorCode DMPlexReorderGetDefault_Plex(DM dm, DMReorderDefaultFlag *reorder)
400: {
401: DM_Plex *mesh = (DM_Plex *)dm->data;
403: PetscFunctionBegin;
404: *reorder = mesh->reorderDefault;
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: DMPlexReorderGetDefault - Get flag indicating whether the DM should be reordered by default
411: Not Collective
413: Input Parameter:
414: . dm - The `DM`
416: Output Parameter:
417: . reorder - Flag for reordering
419: Level: intermediate
421: .seealso: `DMPlexReorderSetDefault()`
422: @*/
423: PetscErrorCode DMPlexReorderGetDefault(DM dm, DMReorderDefaultFlag *reorder)
424: {
425: PetscFunctionBegin;
427: PetscAssertPointer(reorder, 2);
428: PetscUseMethod(dm, "DMPlexReorderGetDefault_C", (DM, DMReorderDefaultFlag *), (dm, reorder));
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
432: static PetscErrorCode DMCreateSectionPermutation_Plex_Reverse(DM dm, IS *permutation, PetscBT *blockStarts)
433: {
434: IS permIS;
435: PetscInt *perm;
436: PetscInt pStart, pEnd;
438: PetscFunctionBegin;
439: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
440: PetscCall(PetscMalloc1(pEnd - pStart, &perm));
441: for (PetscInt p = pStart; p < pEnd; ++p) perm[pEnd - 1 - p] = p;
442: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
443: PetscCall(ISSetPermutation(permIS));
444: *permutation = permIS;
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: // Reorder to group split nodes
449: static PetscErrorCode DMCreateSectionPermutation_Plex_Cohesive_Old(DM dm, IS *permutation, PetscBT *blockStarts)
450: {
451: IS permIS;
452: PetscBT bt, blst;
453: PetscInt *perm;
454: PetscInt pStart, pEnd, i = 0;
456: PetscFunctionBegin;
457: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
458: PetscCall(PetscMalloc1(pEnd - pStart, &perm));
459: PetscCall(PetscBTCreate(pEnd - pStart, &bt));
460: PetscCall(PetscBTCreate(pEnd - pStart, &blst));
461: for (PetscInt p = pStart; p < pEnd; ++p) {
462: const PetscInt *supp, *cone;
463: PetscInt suppSize;
464: DMPolytopeType ct;
466: PetscCall(DMPlexGetCellType(dm, p, &ct));
467: // Do not order tensor cells until they appear below
468: 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;
469: if (PetscBTLookupSet(bt, p)) continue;
470: PetscCall(PetscBTSet(blst, p));
471: perm[i++] = p;
472: // Check for tensor cells in the support
473: PetscCall(DMPlexGetSupport(dm, p, &supp));
474: PetscCall(DMPlexGetSupportSize(dm, p, &suppSize));
475: for (PetscInt s = 0; s < suppSize; ++s) {
476: DMPolytopeType sct;
477: PetscInt q, qq;
479: PetscCall(DMPlexGetCellType(dm, supp[s], &sct));
480: switch (sct) {
481: case DM_POLYTOPE_POINT_PRISM_TENSOR:
482: case DM_POLYTOPE_SEG_PRISM_TENSOR:
483: case DM_POLYTOPE_TRI_PRISM_TENSOR:
484: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
485: // If found, move up the split partner of the tensor cell, and the cell itself
486: PetscCall(DMPlexGetCone(dm, supp[s], &cone));
487: qq = supp[s];
488: q = (cone[0] == p) ? cone[1] : cone[0];
489: if (!PetscBTLookupSet(bt, q)) {
490: perm[i++] = q;
491: s = suppSize;
492: // At T-junctions, we can have an unsplit point at the other end, so also order that loop
493: {
494: const PetscInt *qsupp, *qcone;
495: PetscInt qsuppSize;
497: PetscCall(DMPlexGetSupport(dm, q, &qsupp));
498: PetscCall(DMPlexGetSupportSize(dm, q, &qsuppSize));
499: for (PetscInt qs = 0; qs < qsuppSize; ++qs) {
500: DMPolytopeType qsct;
502: PetscCall(DMPlexGetCellType(dm, qsupp[qs], &qsct));
503: switch (qsct) {
504: case DM_POLYTOPE_POINT_PRISM_TENSOR:
505: case DM_POLYTOPE_SEG_PRISM_TENSOR:
506: case DM_POLYTOPE_TRI_PRISM_TENSOR:
507: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
508: PetscCall(DMPlexGetCone(dm, qsupp[qs], &qcone));
509: if (qcone[0] == qcone[1]) {
510: if (!PetscBTLookupSet(bt, qsupp[qs])) {
511: perm[i++] = qsupp[qs];
512: qs = qsuppSize;
513: }
514: }
515: break;
516: default:
517: break;
518: }
519: }
520: }
521: }
522: if (!PetscBTLookupSet(bt, qq)) {
523: perm[i++] = qq;
524: s = suppSize;
525: }
526: break;
527: default:
528: break;
529: }
530: }
531: }
532: if (PetscDefined(USE_DEBUG)) {
533: 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);
534: }
535: PetscCall(PetscBTDestroy(&bt));
536: 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);
537: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
538: PetscCall(ISSetPermutation(permIS));
539: *permutation = permIS;
540: *blockStarts = blst;
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: // Mark the block associated with a cohesive cell p
545: static PetscErrorCode InsertCohesiveBlock_Private(DM dm, PetscBT bt, PetscBT blst, PetscInt p, PetscInt *idx, PetscInt perm[])
546: {
547: const PetscInt *cone;
548: PetscInt cS;
550: PetscFunctionBegin;
551: if (PetscBTLookupSet(bt, p)) PetscFunctionReturn(PETSC_SUCCESS);
552: // Order the endcaps
553: PetscCall(DMPlexGetCone(dm, p, &cone));
554: PetscCall(DMPlexGetConeSize(dm, p, &cS));
555: if (blst) PetscCall(PetscBTSet(blst, cone[0]));
556: if (!PetscBTLookupSet(bt, cone[0])) perm[(*idx)++] = cone[0];
557: if (!PetscBTLookupSet(bt, cone[1])) perm[(*idx)++] = cone[1];
558: // Order sides
559: for (PetscInt c = 2; c < cS; ++c) PetscCall(InsertCohesiveBlock_Private(dm, bt, NULL, cone[c], idx, perm));
560: // Order cell
561: perm[(*idx)++] = p;
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: // Reorder to group split nodes
566: static PetscErrorCode DMCreateSectionPermutation_Plex_Cohesive(DM dm, IS *permutation, PetscBT *blockStarts)
567: {
568: IS permIS;
569: PetscBT bt, blst;
570: PetscInt *perm;
571: PetscInt dim, pStart, pEnd, i = 0;
573: PetscFunctionBegin;
574: PetscCall(DMGetDimension(dm, &dim));
575: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
576: PetscCall(PetscMalloc1(pEnd - pStart, &perm));
577: PetscCall(PetscBTCreate(pEnd - pStart, &bt));
578: PetscCall(PetscBTCreate(pEnd - pStart, &blst));
579: // Add cohesive blocks
580: for (PetscInt p = pStart; p < pEnd; ++p) {
581: DMPolytopeType ct;
583: PetscCall(DMPlexGetCellType(dm, p, &ct));
584: switch (dim) {
585: case 2:
586: if (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) PetscCall(InsertCohesiveBlock_Private(dm, bt, blst, p, &i, perm));
587: break;
588: case 3:
589: if (ct == DM_POLYTOPE_TRI_PRISM_TENSOR || ct == DM_POLYTOPE_QUAD_PRISM_TENSOR) PetscCall(InsertCohesiveBlock_Private(dm, bt, blst, p, &i, perm));
590: break;
591: default:
592: break;
593: }
594: }
595: // Add normal blocks
596: for (PetscInt p = pStart; p < pEnd; ++p) {
597: if (PetscBTLookupSet(bt, p)) continue;
598: PetscCall(PetscBTSet(blst, p));
599: perm[i++] = p;
600: }
601: if (PetscDefined(USE_DEBUG)) {
602: 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);
603: }
604: PetscCall(PetscBTDestroy(&bt));
605: 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);
606: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
607: PetscCall(ISSetPermutation(permIS));
608: *permutation = permIS;
609: *blockStarts = blst;
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: PetscErrorCode DMCreateSectionPermutation_Plex(DM dm, IS *perm, PetscBT *blockStarts)
614: {
615: DMReorderDefaultFlag reorder;
616: MatOrderingType otype;
617: PetscBool iscohesive, iscohesiveOld, isreverse;
619: PetscFunctionBegin;
620: PetscCall(DMReorderSectionGetDefault(dm, &reorder));
621: if (reorder != DM_REORDER_DEFAULT_TRUE) PetscFunctionReturn(PETSC_SUCCESS);
622: PetscCall(DMReorderSectionGetType(dm, &otype));
623: if (!otype) PetscFunctionReturn(PETSC_SUCCESS);
624: PetscCall(PetscStrncmp(otype, "cohesive_old", 1024, &iscohesiveOld));
625: PetscCall(PetscStrncmp(otype, "cohesive", 1024, &iscohesive));
626: PetscCall(PetscStrncmp(otype, "reverse", 1024, &isreverse));
627: if (iscohesive) {
628: PetscCall(DMCreateSectionPermutation_Plex_Cohesive(dm, perm, blockStarts));
629: } else if (iscohesiveOld) {
630: PetscCall(DMCreateSectionPermutation_Plex_Cohesive_Old(dm, perm, blockStarts));
631: } else if (isreverse) {
632: PetscCall(DMCreateSectionPermutation_Plex_Reverse(dm, perm, blockStarts));
633: }
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: PetscErrorCode DMReorderSectionSetDefault_Plex(DM dm, DMReorderDefaultFlag reorder)
638: {
639: PetscFunctionBegin;
640: dm->reorderSection = reorder;
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: PetscErrorCode DMReorderSectionGetDefault_Plex(DM dm, DMReorderDefaultFlag *reorder)
645: {
646: PetscFunctionBegin;
647: *reorder = dm->reorderSection;
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
651: PetscErrorCode DMReorderSectionSetType_Plex(DM dm, MatOrderingType reorder)
652: {
653: PetscFunctionBegin;
654: PetscCall(PetscFree(dm->reorderSectionType));
655: PetscCall(PetscStrallocpy(reorder, (char **)&dm->reorderSectionType));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: PetscErrorCode DMReorderSectionGetType_Plex(DM dm, MatOrderingType *reorder)
660: {
661: PetscFunctionBegin;
662: *reorder = dm->reorderSectionType;
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }