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, §ion));
267: PetscCall(PetscSectionPermute(section, perm, §ionNew));
268: PetscCall(DMSetLocalSection(*pdm, sectionNew));
269: PetscCall(PetscSectionDestroy(§ionNew));
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, DMPlexReorderDefaultFlag 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, DMPlexReorderDefaultFlag reorder)
392: {
393: PetscFunctionBegin;
395: PetscTryMethod(dm, "DMPlexReorderSetDefault_C", (DM, DMPlexReorderDefaultFlag), (dm, reorder));
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: PetscErrorCode DMPlexReorderGetDefault_Plex(DM dm, DMPlexReorderDefaultFlag *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, DMPlexReorderDefaultFlag *reorder)
424: {
425: PetscFunctionBegin;
427: PetscAssertPointer(reorder, 2);
428: PetscUseMethod(dm, "DMPlexReorderGetDefault_C", (DM, DMPlexReorderDefaultFlag *), (dm, reorder));
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }