Actual source code: plexreorder.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/matorderimpl.h>

      4: static PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm)
      5: {
      6:   PetscInt      *perm, *iperm;
      7:   PetscInt       depth, d, pStart, pEnd, fStart, fMax, fEnd, p;

     11:   DMPlexGetDepth(dm, &depth);
     12:   DMPlexGetChart(dm, &pStart, &pEnd);
     13:   PetscMalloc1(pEnd-pStart,&perm);
     14:   PetscMalloc1(pEnd-pStart,&iperm);
     15:   for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
     16:   for (d = depth; d > 0; --d) {
     17:     DMPlexGetDepthStratum(dm, d,   &pStart, &pEnd);
     18:     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:       DMPlexGetConeSize(dm, point, &coneSize);
     30:       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:     if (fMax != fEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %d faces %d does not match permuted number %d", d, fEnd-fStart, fMax-fStart);
     42:   }
     43:   *clperm    = perm;
     44:   *invclperm = iperm;
     45:   return(0);
     46: }

     48: /*@
     49:   DMPlexGetOrdering - Calculate a reordering of the mesh

     51:   Collective on dm

     53:   Input Parameter:
     54: + dm - The DMPlex object
     55: . otype - type of reordering, one of the following:
     56: $     MATORDERINGNATURAL - Natural
     57: $     MATORDERINGND - Nested Dissection
     58: $     MATORDERING1WD - One-way Dissection
     59: $     MATORDERINGRCM - Reverse Cuthill-McKee
     60: $     MATORDERINGQMD - Quotient Minimum Degree
     61: - label - [Optional] Label used to segregate ordering into sets, or NULL

     64:   Output Parameter:
     65: . perm - The point permutation as an IS, perm[old point number] = new point number

     67:   Note: The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which
     68:   has different types of cells, and then loop over each set of reordered cells for assembly.

     70:   Level: intermediate

     72: .seealso: MatGetOrdering()
     73: @*/
     74: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
     75: {
     76:   PetscInt       numCells = 0;
     77:   PetscInt      *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i;

     83:   DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);
     84:   PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);
     85:   if (numCells) {
     86:     /* Shift for Fortran numbering */
     87:     for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
     88:     for (i = 0; i <= numCells; ++i)       ++start[i];
     89:     SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);
     90:   }
     91:   PetscFree(start);
     92:   PetscFree(adjacency);
     93:   /* Shift for Fortran numbering */
     94:   for (c = 0; c < numCells; ++c) --cperm[c];
     95:   /* Segregate */
     96:   if (label) {
     97:     IS              valueIS;
     98:     const PetscInt *values;
     99:     PetscInt        numValues, numPoints = 0;
    100:     PetscInt       *sperm, *vsize, *voff, v;

    102:     DMLabelGetValueIS(label, &valueIS);
    103:     ISSort(valueIS);
    104:     ISGetLocalSize(valueIS, &numValues);
    105:     ISGetIndices(valueIS, &values);
    106:     PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff);
    107:     for (v = 0; v < numValues; ++v) {
    108:       DMLabelGetStratumSize(label, values[v], &vsize[v]);
    109:       if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1];
    110:       numPoints += vsize[v];
    111:     }
    112:     if (numPoints != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %D cells < %D total", numPoints, numCells);
    113:     for (c = 0; c < numCells; ++c) {
    114:       const PetscInt oldc = cperm[c];
    115:       PetscInt       val, vloc;

    117:       DMLabelGetValue(label, oldc, &val);
    118:       if (val == -1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %D not present in label", oldc);
    119:       PetscFindInt(val, numValues, values, &vloc);
    120:       if (vloc < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %D not present label", val);
    121:       sperm[voff[vloc+1]++] = oldc;
    122:     }
    123:     for (v = 0; v < numValues; ++v) {
    124:       if (voff[v+1] - voff[v] != vsize[v]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of %D values found is %D != %D", values[v], voff[v+1] - voff[v], vsize[v]);
    125:     }
    126:     ISRestoreIndices(valueIS, &values);
    127:     ISDestroy(&valueIS);
    128:     PetscArraycpy(cperm, sperm, numCells);
    129:     PetscFree3(sperm, vsize, voff);
    130:   }
    131:   /* Construct closure */
    132:   DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
    133:   PetscFree3(cperm,mask,xls);
    134:   PetscFree(clperm);
    135:   /* Invert permutation */
    136:   DMPlexGetChart(dm, &pStart, &pEnd);
    137:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);
    138:   return(0);
    139: }

    141: /*@
    142:   DMPlexPermute - Reorder the mesh according to the input permutation

    144:   Collective on dm

    146:   Input Parameter:
    147: + dm - The DMPlex object
    148: - perm - The point permutation, perm[old point number] = new point number

    150:   Output Parameter:
    151: . pdm - The permuted DM

    153:   Level: intermediate

    155: .seealso: MatPermute()
    156: @*/
    157: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
    158: {
    159:   DM_Plex       *plex = (DM_Plex *) dm->data, *plexNew;
    160:   PetscSection   section, sectionNew;
    161:   PetscInt       dim;

    168:   DMCreate(PetscObjectComm((PetscObject) dm), pdm);
    169:   DMSetType(*pdm, DMPLEX);
    170:   DMGetDimension(dm, &dim);
    171:   DMSetDimension(*pdm, dim);
    172:   DMCopyDisc(dm, *pdm);
    173:   DMGetLocalSection(dm, &section);
    174:   if (section) {
    175:     PetscSectionPermute(section, perm, &sectionNew);
    176:     DMSetLocalSection(*pdm, sectionNew);
    177:     PetscSectionDestroy(&sectionNew);
    178:   }
    179:   plexNew = (DM_Plex *) (*pdm)->data;
    180:   /* Ignore ltogmap, ltogmapb */
    181:   /* Ignore sf, sectionSF */
    182:   /* Ignore globalVertexNumbers, globalCellNumbers */
    183:   /* Reorder labels */
    184:   {
    185:     PetscInt numLabels, l;
    186:     DMLabel  label, labelNew;

    188:     DMGetNumLabels(dm, &numLabels);
    189:     for (l = 0; l < numLabels; ++l) {
    190:       DMGetLabelByNum(dm, l, &label);
    191:       DMLabelPermute(label, perm, &labelNew);
    192:       DMAddLabel(*pdm, labelNew);
    193:       DMLabelDestroy(&labelNew);
    194:     }
    195:     DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel);
    196:     if (plex->subpointMap) {DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);}
    197:   }
    198:   /* Reorder topology */
    199:   {
    200:     const PetscInt *pperm;
    201:     PetscInt        maxConeSize, maxSupportSize, n, pStart, pEnd, p;

    203:     DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
    204:     plexNew->maxConeSize    = maxConeSize;
    205:     plexNew->maxSupportSize = maxSupportSize;
    206:     PetscSectionDestroy(&plexNew->coneSection);
    207:     PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
    208:     PetscSectionGetStorageSize(plexNew->coneSection, &n);
    209:     PetscMalloc1(n, &plexNew->cones);
    210:     PetscMalloc1(n, &plexNew->coneOrientations);
    211:     ISGetIndices(perm, &pperm);
    212:     PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
    213:     for (p = pStart; p < pEnd; ++p) {
    214:       PetscInt dof, off, offNew, d;

    216:       PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
    217:       PetscSectionGetOffset(plex->coneSection, p, &off);
    218:       PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
    219:       for (d = 0; d < dof; ++d) {
    220:         plexNew->cones[offNew+d]            = pperm[plex->cones[off+d]];
    221:         plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
    222:       }
    223:     }
    224:     PetscSectionDestroy(&plexNew->supportSection);
    225:     PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
    226:     PetscSectionGetStorageSize(plexNew->supportSection, &n);
    227:     PetscMalloc1(n, &plexNew->supports);
    228:     PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
    229:     for (p = pStart; p < pEnd; ++p) {
    230:       PetscInt dof, off, offNew, d;

    232:       PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
    233:       PetscSectionGetOffset(plex->supportSection, p, &off);
    234:       PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
    235:       for (d = 0; d < dof; ++d) {
    236:         plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
    237:       }
    238:     }
    239:     ISRestoreIndices(perm, &pperm);
    240:   }
    241:   /* Remap coordinates */
    242:   {
    243:     DM              cdm, cdmNew;
    244:     PetscSection    csection, csectionNew;
    245:     Vec             coordinates, coordinatesNew;
    246:     PetscScalar    *coords, *coordsNew;
    247:     const PetscInt *pperm;
    248:     PetscInt        pStart, pEnd, p;
    249:     const char     *name;

    251:     DMGetCoordinateDM(dm, &cdm);
    252:     DMGetLocalSection(cdm, &csection);
    253:     PetscSectionPermute(csection, perm, &csectionNew);
    254:     DMGetCoordinatesLocal(dm, &coordinates);
    255:     VecDuplicate(coordinates, &coordinatesNew);
    256:     PetscObjectGetName((PetscObject)coordinates,&name);
    257:     PetscObjectSetName((PetscObject)coordinatesNew,name);
    258:     VecGetArray(coordinates, &coords);
    259:     VecGetArray(coordinatesNew, &coordsNew);
    260:     PetscSectionGetChart(csectionNew, &pStart, &pEnd);
    261:     ISGetIndices(perm, &pperm);
    262:     for (p = pStart; p < pEnd; ++p) {
    263:       PetscInt dof, off, offNew, d;

    265:       PetscSectionGetDof(csectionNew, p, &dof);
    266:       PetscSectionGetOffset(csection, p, &off);
    267:       PetscSectionGetOffset(csectionNew, pperm[p], &offNew);
    268:       for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
    269:     }
    270:     ISRestoreIndices(perm, &pperm);
    271:     VecRestoreArray(coordinates, &coords);
    272:     VecRestoreArray(coordinatesNew, &coordsNew);
    273:     DMGetCoordinateDM(*pdm, &cdmNew);
    274:     DMSetLocalSection(cdmNew, csectionNew);
    275:     DMSetCoordinatesLocal(*pdm, coordinatesNew);
    276:     PetscSectionDestroy(&csectionNew);
    277:     VecDestroy(&coordinatesNew);
    278:   }
    279:   (*pdm)->setupcalled = PETSC_TRUE;
    280:   return(0);
    281: }