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, &section));
267:     PetscCall(PetscSectionPermute(section, perm, &sectionNew));
268:     PetscCall(DMSetLocalSection(*pdm, sectionNew));
269:     PetscCall(PetscSectionDestroy(&sectionNew));
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(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: PetscErrorCode DMCreateSectionPermutation_Plex(DM dm, IS *perm, PetscBT *blockStarts)
545: {
546:   DMReorderDefaultFlag reorder;
547:   MatOrderingType      otype;
548:   PetscBool            iscohesive, isreverse;

550:   PetscFunctionBegin;
551:   PetscCall(DMReorderSectionGetDefault(dm, &reorder));
552:   if (reorder != DM_REORDER_DEFAULT_TRUE) PetscFunctionReturn(PETSC_SUCCESS);
553:   PetscCall(DMReorderSectionGetType(dm, &otype));
554:   if (!otype) PetscFunctionReturn(PETSC_SUCCESS);
555:   PetscCall(PetscStrncmp(otype, "cohesive", 1024, &iscohesive));
556:   PetscCall(PetscStrncmp(otype, "reverse", 1024, &isreverse));
557:   if (iscohesive) {
558:     PetscCall(DMCreateSectionPermutation_Plex_Cohesive(dm, perm, blockStarts));
559:   } else if (isreverse) {
560:     PetscCall(DMCreateSectionPermutation_Plex_Reverse(dm, perm, blockStarts));
561:   }
562:   PetscFunctionReturn(PETSC_SUCCESS);
563: }

565: PetscErrorCode DMReorderSectionSetDefault_Plex(DM dm, DMReorderDefaultFlag reorder)
566: {
567:   PetscFunctionBegin;
568:   dm->reorderSection = reorder;
569:   PetscFunctionReturn(PETSC_SUCCESS);
570: }

572: PetscErrorCode DMReorderSectionGetDefault_Plex(DM dm, DMReorderDefaultFlag *reorder)
573: {
574:   PetscFunctionBegin;
575:   *reorder = dm->reorderSection;
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

579: PetscErrorCode DMReorderSectionSetType_Plex(DM dm, MatOrderingType reorder)
580: {
581:   PetscFunctionBegin;
582:   PetscCall(PetscFree(dm->reorderSectionType));
583:   PetscCall(PetscStrallocpy(reorder, (char **)&dm->reorderSectionType));
584:   PetscFunctionReturn(PETSC_SUCCESS);
585: }

587: PetscErrorCode DMReorderSectionGetType_Plex(DM dm, MatOrderingType *reorder)
588: {
589:   PetscFunctionBegin;
590:   *reorder = dm->reorderSectionType;
591:   PetscFunctionReturn(PETSC_SUCCESS);
592: }