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, §ion);
266: PetscSectionPermute(section, perm, §ionNew);
267: DMSetLocalSection(*pdm, sectionNew);
268: PetscSectionDestroy(§ionNew);
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: }