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:   DMPlexGetDepth(dm, &depth);
 11:   DMPlexGetChart(dm, &pStart, &pEnd);
 12:   PetscMalloc1(pEnd - pStart, &perm);
 13:   PetscMalloc1(pEnd - pStart, &iperm);
 14:   for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
 15:   for (d = depth; d > 0; --d) {
 16:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
 17:     DMPlexGetDepthStratum(dm, d - 1, &fStart, &fEnd);
 18:     fMax = fStart;
 19:     for (p = pStart; p < pEnd; ++p) {
 20:       const PetscInt *cone;
 21:       PetscInt        point, coneSize, c;

 23:       if (d == depth) {
 24:         perm[p]         = pperm[p];
 25:         iperm[pperm[p]] = p;
 26:       }
 27:       point = perm[p];
 28:       DMPlexGetConeSize(dm, point, &coneSize);
 29:       DMPlexGetCone(dm, point, &cone);
 30:       for (c = 0; c < coneSize; ++c) {
 31:         const PetscInt oldc = cone[c];
 32:         const PetscInt newc = iperm[oldc];

 34:         if (newc < 0) {
 35:           perm[fMax]  = oldc;
 36:           iperm[oldc] = fMax++;
 37:         }
 38:       }
 39:     }
 41:   }
 42:   *clperm    = perm;
 43:   *invclperm = iperm;
 44:   return 0;
 45: }

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

 50:   Collective on dm

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

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

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

 68:   Level: intermediate

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

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

 99:     // Can't directly sort the valueIS, since it is a view into the DMLabel
100:     DMLabelGetValueIS(label, &valueIS);
101:     ISGetLocalSize(valueIS, &numValues);
102:     ISGetIndices(valueIS, &valuesTmp);
103:     PetscCalloc4(numCells, &sperm, numValues, &values, numValues, &vsize, numValues + 1, &voff);
104:     PetscArraycpy(values, valuesTmp, numValues);
105:     PetscSortInt(numValues, values);
106:     ISRestoreIndices(valueIS, &valuesTmp);
107:     ISDestroy(&valueIS);
108:     for (v = 0; v < numValues; ++v) {
109:       DMLabelGetStratumSize(label, values[v], &vsize[v]);
110:       if (v < numValues - 1) voff[v + 2] += vsize[v] + voff[v + 1];
111:       numPoints += vsize[v];
112:     }
114:     for (c = 0; c < numCells; ++c) {
115:       const PetscInt oldc = cperm[c];
116:       PetscInt       val, vloc;

118:       DMLabelGetValue(label, oldc, &val);
120:       PetscFindInt(val, numValues, values, &vloc);
122:       sperm[voff[vloc + 1]++] = oldc;
123:     }
125:     PetscArraycpy(cperm, sperm, numCells);
126:     PetscFree4(sperm, values, vsize, voff);
127:   }
128:   /* Construct closure */
129:   DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
130:   PetscFree3(cperm, mask, xls);
131:   PetscFree(clperm);
132:   /* Invert permutation */
133:   DMPlexGetChart(dm, &pStart, &pEnd);
134:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, invclperm, PETSC_OWN_POINTER, perm);
135:   return 0;
136: }

138: /*@
139:   DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line

141:   Collective on dm

143:   Input Parameter:
144: . dm - The DMPlex object

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

149:   Level: intermediate

151: .seealso: `DMPlexGetOrdering()`, `DMPlexPermute()`, `MatGetOrdering()`
152: @*/
153: PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm)
154: {
155:   PetscInt       *points;
156:   const PetscInt *support, *cone;
157:   PetscInt        dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0;

159:   DMGetDimension(dm, &dim);
161:   DMPlexGetChart(dm, &pStart, &pEnd);
162:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
163:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
164:   PetscMalloc1(pEnd - pStart, &points);
165:   for (c = cStart; c < cEnd; ++c) points[c] = c;
166:   for (v = vStart; v < vEnd; ++v) points[v] = v;
167:   for (v = vStart; v < vEnd; ++v) {
168:     DMPlexGetSupportSize(dm, v, &suppSize);
169:     DMPlexGetSupport(dm, v, &support);
170:     if (suppSize == 1) {
171:       lastCell = support[0];
172:       break;
173:     }
174:   }
175:   if (v < vEnd) {
176:     PetscInt pos = cEnd;

178:     points[v] = pos++;
179:     while (lastCell >= cStart) {
180:       DMPlexGetCone(dm, lastCell, &cone);
181:       if (cone[0] == v) v = cone[1];
182:       else v = cone[0];
183:       DMPlexGetSupport(dm, v, &support);
184:       DMPlexGetSupportSize(dm, v, &suppSize);
185:       if (suppSize == 1) {
186:         lastCell = -1;
187:       } else {
188:         if (support[0] == lastCell) lastCell = support[1];
189:         else lastCell = support[0];
190:       }
191:       points[v] = pos++;
192:     }
194:   }
195:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, points, PETSC_OWN_POINTER, perm);
196:   return 0;
197: }

199: static PetscErrorCode DMPlexRemapCoordinates_Private(IS perm, PetscSection cs, Vec coordinates, PetscSection *csNew, Vec *coordinatesNew)
200: {
201:   PetscScalar    *coords, *coordsNew;
202:   const PetscInt *pperm;
203:   PetscInt        pStart, pEnd, p;
204:   const char     *name;

206:   PetscSectionPermute(cs, perm, csNew);
207:   VecDuplicate(coordinates, coordinatesNew);
208:   PetscObjectGetName((PetscObject)coordinates, &name);
209:   PetscObjectSetName((PetscObject)*coordinatesNew, name);
210:   VecGetArray(coordinates, &coords);
211:   VecGetArray(*coordinatesNew, &coordsNew);
212:   PetscSectionGetChart(*csNew, &pStart, &pEnd);
213:   ISGetIndices(perm, &pperm);
214:   for (p = pStart; p < pEnd; ++p) {
215:     PetscInt dof, off, offNew, d;

217:     PetscSectionGetDof(*csNew, p, &dof);
218:     PetscSectionGetOffset(cs, p, &off);
219:     PetscSectionGetOffset(*csNew, pperm[p], &offNew);
220:     for (d = 0; d < dof; ++d) coordsNew[offNew + d] = coords[off + d];
221:   }
222:   ISRestoreIndices(perm, &pperm);
223:   VecRestoreArray(coordinates, &coords);
224:   VecRestoreArray(*coordinatesNew, &coordsNew);
225:   return 0;
226: }

228: /*@
229:   DMPlexPermute - Reorder the mesh according to the input permutation

231:   Collective on dm

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: `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;

253:   DMCreate(PetscObjectComm((PetscObject)dm), pdm);
254:   DMSetType(*pdm, DMPLEX);
255:   PetscObjectGetName((PetscObject)dm, &name);
256:   PetscObjectSetName((PetscObject)*pdm, name);
257:   DMGetDimension(dm, &dim);
258:   DMSetDimension(*pdm, dim);
259:   DMGetCoordinateDim(dm, &cdim);
260:   DMSetCoordinateDim(*pdm, cdim);
261:   DMCopyDisc(dm, *pdm);
262:   if (dm->localSection) {
263:     PetscSection section, sectionNew;

265:     DMGetLocalSection(dm, &section);
266:     PetscSectionPermute(section, perm, &sectionNew);
267:     DMSetLocalSection(*pdm, sectionNew);
268:     PetscSectionDestroy(&sectionNew);
269:   }
270:   plexNew = (DM_Plex *)(*pdm)->data;
271:   /* Ignore ltogmap, ltogmapb */
272:   /* Ignore sf, sectionSF */
273:   /* Ignore globalVertexNumbers, globalCellNumbers */
274:   /* Reorder labels */
275:   {
276:     PetscInt numLabels, l;
277:     DMLabel  label, labelNew;

279:     DMGetNumLabels(dm, &numLabels);
280:     for (l = 0; l < numLabels; ++l) {
281:       DMGetLabelByNum(dm, l, &label);
282:       DMLabelPermute(label, perm, &labelNew);
283:       DMAddLabel(*pdm, labelNew);
284:       DMLabelDestroy(&labelNew);
285:     }
286:     DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel);
287:     if (plex->subpointMap) DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);
288:   }
289:   if ((*pdm)->celltypeLabel) {
290:     DMLabel ctLabel;

292:     // Reset label for fast lookup
293:     DMPlexGetCellTypeLabel(*pdm, &ctLabel);
294:     DMLabelMakeAllInvalid_Internal(ctLabel);
295:   }
296:   /* Reorder topology */
297:   {
298:     const PetscInt *pperm;
299:     PetscInt        n, pStart, pEnd, p;

301:     PetscSectionDestroy(&plexNew->coneSection);
302:     PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
303:     PetscSectionGetStorageSize(plexNew->coneSection, &n);
304:     PetscMalloc1(n, &plexNew->cones);
305:     PetscMalloc1(n, &plexNew->coneOrientations);
306:     ISGetIndices(perm, &pperm);
307:     PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
308:     for (p = pStart; p < pEnd; ++p) {
309:       PetscInt dof, off, offNew, d;

311:       PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
312:       PetscSectionGetOffset(plex->coneSection, p, &off);
313:       PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
314:       for (d = 0; d < dof; ++d) {
315:         plexNew->cones[offNew + d]            = pperm[plex->cones[off + d]];
316:         plexNew->coneOrientations[offNew + d] = plex->coneOrientations[off + d];
317:       }
318:     }
319:     PetscSectionDestroy(&plexNew->supportSection);
320:     PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
321:     PetscSectionGetStorageSize(plexNew->supportSection, &n);
322:     PetscMalloc1(n, &plexNew->supports);
323:     PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
324:     for (p = pStart; p < pEnd; ++p) {
325:       PetscInt dof, off, offNew, d;

327:       PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
328:       PetscSectionGetOffset(plex->supportSection, p, &off);
329:       PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
330:       for (d = 0; d < dof; ++d) plexNew->supports[offNew + d] = pperm[plex->supports[off + d]];
331:     }
332:     ISRestoreIndices(perm, &pperm);
333:   }
334:   /* Remap coordinates */
335:   {
336:     DM           cdm, cdmNew;
337:     PetscSection cs, csNew;
338:     Vec          coordinates, coordinatesNew;

340:     DMGetCoordinateSection(dm, &cs);
341:     DMGetCoordinatesLocal(dm, &coordinates);
342:     DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew);
343:     DMSetCoordinateSection(*pdm, PETSC_DETERMINE, csNew);
344:     DMSetCoordinatesLocal(*pdm, coordinatesNew);
345:     PetscSectionDestroy(&csNew);
346:     VecDestroy(&coordinatesNew);

348:     DMGetCellCoordinateDM(dm, &cdm);
349:     if (cdm) {
350:       DMGetCoordinateDM(*pdm, &cdm);
351:       DMClone(cdm, &cdmNew);
352:       DMSetCellCoordinateDM(*pdm, cdmNew);
353:       DMDestroy(&cdmNew);
354:       DMGetCellCoordinateSection(dm, &cs);
355:       DMGetCellCoordinatesLocal(dm, &coordinates);
356:       DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew);
357:       DMSetCellCoordinateSection(*pdm, PETSC_DETERMINE, csNew);
358:       DMSetCellCoordinatesLocal(*pdm, coordinatesNew);
359:       PetscSectionDestroy(&csNew);
360:       VecDestroy(&coordinatesNew);
361:     }
362:   }
363:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *pdm);
364:   (*pdm)->setupcalled = PETSC_TRUE;
365:   return 0;
366: }

368: PetscErrorCode DMPlexReorderSetDefault_Plex(DM dm, DMPlexReorderDefaultFlag reorder)
369: {
370:   DM_Plex *mesh = (DM_Plex *)dm->data;

372:   mesh->reorderDefault = reorder;
373:   return 0;
374: }

376: /*@
377:   DMPlexReorderSetDefault - Set flag indicating whether the DM should be reordered by default

379:   Logically collective

381:   Input Parameters:
382: + dm        - The DM
383: - reorder   - Flag for reordering

385:   Level: intermediate

387: .seealso: `DMPlexReorderGetDefault()`
388: @*/
389: PetscErrorCode DMPlexReorderSetDefault(DM dm, DMPlexReorderDefaultFlag reorder)
390: {
392:   PetscTryMethod(dm, "DMPlexReorderSetDefault_C", (DM, DMPlexReorderDefaultFlag), (dm, reorder));
393:   return 0;
394: }

396: PetscErrorCode DMPlexReorderGetDefault_Plex(DM dm, DMPlexReorderDefaultFlag *reorder)
397: {
398:   DM_Plex *mesh = (DM_Plex *)dm->data;

400:   *reorder = mesh->reorderDefault;
401:   return 0;
402: }

404: /*@
405:   DMPlexReorderGetDefault - Get flag indicating whether the DM should be reordered by default

407:   Not collective

409:   Input Parameter:
410: . dm      - The DM

412:   Output Parameter:
413: . reorder - Flag for reordering

415:   Level: intermediate

417: .seealso: `DMPlexReorderSetDefault()`
418: @*/
419: PetscErrorCode DMPlexReorderGetDefault(DM dm, DMPlexReorderDefaultFlag *reorder)
420: {
423:   PetscUseMethod(dm, "DMPlexReorderGetDefault_C", (DM, DMPlexReorderDefaultFlag *), (dm, reorder));
424:   return 0;
425: }