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_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) PetscCall(DMCreateSectionPermutation_Plex_Cohesive(dm, perm, blockStarts));
    628:   else if (iscohesiveOld) PetscCall(DMCreateSectionPermutation_Plex_Cohesive_Old(dm, perm, blockStarts));
    629:   else if (isreverse) PetscCall(DMCreateSectionPermutation_Plex_Reverse(dm, perm, blockStarts));
    630:   PetscFunctionReturn(PETSC_SUCCESS);
    631: }

    633: PetscErrorCode DMReorderSectionSetDefault_Plex(DM dm, DMReorderDefaultFlag reorder)
    634: {
    635:   PetscFunctionBegin;
    636:   dm->reorderSection = reorder;
    637:   PetscFunctionReturn(PETSC_SUCCESS);
    638: }

    640: PetscErrorCode DMReorderSectionGetDefault_Plex(DM dm, DMReorderDefaultFlag *reorder)
    641: {
    642:   PetscFunctionBegin;
    643:   *reorder = dm->reorderSection;
    644:   PetscFunctionReturn(PETSC_SUCCESS);
    645: }

    647: PetscErrorCode DMReorderSectionSetType_Plex(DM dm, MatOrderingType reorder)
    648: {
    649:   PetscFunctionBegin;
    650:   PetscCall(PetscFree(dm->reorderSectionType));
    651:   PetscCall(PetscStrallocpy(reorder, (char **)&dm->reorderSectionType));
    652:   PetscFunctionReturn(PETSC_SUCCESS);
    653: }

    655: PetscErrorCode DMReorderSectionGetType_Plex(DM dm, MatOrderingType *reorder)
    656: {
    657:   PetscFunctionBegin;
    658:   *reorder = dm->reorderSectionType;
    659:   PetscFunctionReturn(PETSC_SUCCESS);
    660: }