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;
9: DMPlexGetDepth(dm, &depth);
10: DMPlexGetChart(dm, &pStart, &pEnd);
11: PetscMalloc1(pEnd-pStart,&perm);
12: PetscMalloc1(pEnd-pStart,&iperm);
13: for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
14: for (d = depth; d > 0; --d) {
15: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
16: DMPlexGetDepthStratum(dm, d-1, &fStart, &fEnd);
17: fMax = fStart;
18: for (p = pStart; p < pEnd; ++p) {
19: const PetscInt *cone;
20: PetscInt point, coneSize, c;
22: if (d == depth) {
23: perm[p] = pperm[p];
24: iperm[pperm[p]] = p;
25: }
26: point = perm[p];
27: DMPlexGetConeSize(dm, point, &coneSize);
28: DMPlexGetCone(dm, point, &cone);
29: for (c = 0; c < coneSize; ++c) {
30: const PetscInt oldc = cone[c];
31: const PetscInt newc = iperm[oldc];
33: if (newc < 0) {
34: perm[fMax] = oldc;
35: iperm[oldc] = fMax++;
36: }
37: }
38: }
40: }
41: *clperm = perm;
42: *invclperm = iperm;
43: return 0;
44: }
46: /*@
47: DMPlexGetOrdering - Calculate a reordering of the mesh
49: Collective on dm
51: Input Parameters:
52: + dm - The DMPlex object
53: . otype - type of reordering, one of the following:
54: $ MATORDERINGNATURAL - Natural
55: $ MATORDERINGND - Nested Dissection
56: $ MATORDERING1WD - One-way Dissection
57: $ MATORDERINGRCM - Reverse Cuthill-McKee
58: $ MATORDERINGQMD - Quotient Minimum Degree
59: - label - [Optional] Label used to segregate ordering into sets, or NULL
61: Output Parameter:
62: . perm - The point permutation as an IS, perm[old point number] = new point number
64: Note: 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: Level: intermediate
69: .seealso: DMPlexPermute(), MatGetOrdering()
70: @*/
71: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
72: {
73: PetscInt numCells = 0;
74: PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i;
78: DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);
79: PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);
80: if (numCells) {
81: /* Shift for Fortran numbering */
82: for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
83: for (i = 0; i <= numCells; ++i) ++start[i];
84: SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);
85: }
86: PetscFree(start);
87: PetscFree(adjacency);
88: /* Shift for Fortran numbering */
89: for (c = 0; c < numCells; ++c) --cperm[c];
90: /* Segregate */
91: if (label) {
92: IS valueIS;
93: const PetscInt *values;
94: PetscInt numValues, numPoints = 0;
95: PetscInt *sperm, *vsize, *voff, v;
97: DMLabelGetValueIS(label, &valueIS);
98: ISSort(valueIS);
99: ISGetLocalSize(valueIS, &numValues);
100: ISGetIndices(valueIS, &values);
101: PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff);
102: for (v = 0; v < numValues; ++v) {
103: DMLabelGetStratumSize(label, values[v], &vsize[v]);
104: if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1];
105: numPoints += vsize[v];
106: }
108: for (c = 0; c < numCells; ++c) {
109: const PetscInt oldc = cperm[c];
110: PetscInt val, vloc;
112: DMLabelGetValue(label, oldc, &val);
114: PetscFindInt(val, numValues, values, &vloc);
116: sperm[voff[vloc+1]++] = oldc;
117: }
118: for (v = 0; v < numValues; ++v) {
120: }
121: ISRestoreIndices(valueIS, &values);
122: ISDestroy(&valueIS);
123: PetscArraycpy(cperm, sperm, numCells);
124: PetscFree3(sperm, vsize, voff);
125: }
126: /* Construct closure */
127: DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
128: PetscFree3(cperm,mask,xls);
129: PetscFree(clperm);
130: /* Invert permutation */
131: DMPlexGetChart(dm, &pStart, &pEnd);
132: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);
133: return 0;
134: }
136: /*@
137: DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line
139: Collective on dm
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: 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: DMGetDimension(dm, &dim);
159: DMPlexGetChart(dm, &pStart, &pEnd);
160: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
161: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
162: PetscMalloc1(pEnd-pStart, &points);
163: for (c = cStart; c < cEnd; ++c) points[c] = c;
164: for (v = vStart; v < vEnd; ++v) points[v] = v;
165: for (v = vStart; v < vEnd; ++v) {
166: DMPlexGetSupportSize(dm, v, &suppSize);
167: DMPlexGetSupport(dm, v, &support);
168: if (suppSize == 1) {lastCell = support[0]; break;}
169: }
170: if (v < vEnd) {
171: PetscInt pos = cEnd;
173: points[v] = pos++;
174: while (lastCell >= cStart) {
175: DMPlexGetCone(dm, lastCell, &cone);
176: if (cone[0] == v) v = cone[1];
177: else v = cone[0];
178: DMPlexGetSupport(dm, v, &support);
179: DMPlexGetSupportSize(dm, v, &suppSize);
180: if (suppSize == 1) {lastCell = -1;}
181: else {
182: if (support[0] == lastCell) lastCell = support[1];
183: else lastCell = support[0];
184: }
185: points[v] = pos++;
186: }
188: }
189: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, points, PETSC_OWN_POINTER, perm);
190: return 0;
191: }
193: /*@
194: DMPlexPermute - Reorder the mesh according to the input permutation
196: Collective on dm
198: Input Parameters:
199: + dm - The DMPlex object
200: - perm - The point permutation, perm[old point number] = new point number
202: Output Parameter:
203: . pdm - The permuted DM
205: Level: intermediate
207: .seealso: MatPermute()
208: @*/
209: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
210: {
211: DM_Plex *plex = (DM_Plex *) dm->data, *plexNew;
212: PetscInt dim, cdim;
213: const char *name;
218: DMCreate(PetscObjectComm((PetscObject) dm), pdm);
219: DMSetType(*pdm, DMPLEX);
220: PetscObjectGetName((PetscObject) dm, &name);
221: PetscObjectSetName((PetscObject) *pdm, name);
222: DMGetDimension(dm, &dim);
223: DMSetDimension(*pdm, dim);
224: DMGetCoordinateDim(dm, &cdim);
225: DMSetCoordinateDim(*pdm, cdim);
226: DMCopyDisc(dm, *pdm);
227: if (dm->localSection) {
228: PetscSection section, sectionNew;
230: DMGetLocalSection(dm, §ion);
231: PetscSectionPermute(section, perm, §ionNew);
232: DMSetLocalSection(*pdm, sectionNew);
233: PetscSectionDestroy(§ionNew);
234: }
235: plexNew = (DM_Plex *) (*pdm)->data;
236: /* Ignore ltogmap, ltogmapb */
237: /* Ignore sf, sectionSF */
238: /* Ignore globalVertexNumbers, globalCellNumbers */
239: /* Reorder labels */
240: {
241: PetscInt numLabels, l;
242: DMLabel label, labelNew;
244: DMGetNumLabels(dm, &numLabels);
245: for (l = 0; l < numLabels; ++l) {
246: DMGetLabelByNum(dm, l, &label);
247: DMLabelPermute(label, perm, &labelNew);
248: DMAddLabel(*pdm, labelNew);
249: DMLabelDestroy(&labelNew);
250: }
251: DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel);
252: if (plex->subpointMap) DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);
253: }
254: /* Reorder topology */
255: {
256: const PetscInt *pperm;
257: PetscInt n, pStart, pEnd, p;
259: PetscSectionDestroy(&plexNew->coneSection);
260: PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
261: PetscSectionGetStorageSize(plexNew->coneSection, &n);
262: PetscMalloc1(n, &plexNew->cones);
263: PetscMalloc1(n, &plexNew->coneOrientations);
264: ISGetIndices(perm, &pperm);
265: PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
266: for (p = pStart; p < pEnd; ++p) {
267: PetscInt dof, off, offNew, d;
269: PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
270: PetscSectionGetOffset(plex->coneSection, p, &off);
271: PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
272: for (d = 0; d < dof; ++d) {
273: plexNew->cones[offNew+d] = pperm[plex->cones[off+d]];
274: plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
275: }
276: }
277: PetscSectionDestroy(&plexNew->supportSection);
278: PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
279: PetscSectionGetStorageSize(plexNew->supportSection, &n);
280: PetscMalloc1(n, &plexNew->supports);
281: PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
282: for (p = pStart; p < pEnd; ++p) {
283: PetscInt dof, off, offNew, d;
285: PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
286: PetscSectionGetOffset(plex->supportSection, p, &off);
287: PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
288: for (d = 0; d < dof; ++d) {
289: plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
290: }
291: }
292: ISRestoreIndices(perm, &pperm);
293: }
294: /* Remap coordinates */
295: {
296: DM cdm, cdmNew;
297: PetscSection csection, csectionNew;
298: Vec coordinates, coordinatesNew;
299: PetscScalar *coords, *coordsNew;
300: const PetscInt *pperm;
301: PetscInt pStart, pEnd, p;
302: const char *name;
304: DMGetCoordinateDM(dm, &cdm);
305: DMGetLocalSection(cdm, &csection);
306: PetscSectionPermute(csection, perm, &csectionNew);
307: DMGetCoordinatesLocal(dm, &coordinates);
308: VecDuplicate(coordinates, &coordinatesNew);
309: PetscObjectGetName((PetscObject)coordinates,&name);
310: PetscObjectSetName((PetscObject)coordinatesNew,name);
311: VecGetArray(coordinates, &coords);
312: VecGetArray(coordinatesNew, &coordsNew);
313: PetscSectionGetChart(csectionNew, &pStart, &pEnd);
314: ISGetIndices(perm, &pperm);
315: for (p = pStart; p < pEnd; ++p) {
316: PetscInt dof, off, offNew, d;
318: PetscSectionGetDof(csectionNew, p, &dof);
319: PetscSectionGetOffset(csection, p, &off);
320: PetscSectionGetOffset(csectionNew, pperm[p], &offNew);
321: for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
322: }
323: ISRestoreIndices(perm, &pperm);
324: VecRestoreArray(coordinates, &coords);
325: VecRestoreArray(coordinatesNew, &coordsNew);
326: DMGetCoordinateDM(*pdm, &cdmNew);
327: DMSetLocalSection(cdmNew, csectionNew);
328: DMSetCoordinatesLocal(*pdm, coordinatesNew);
329: PetscSectionDestroy(&csectionNew);
330: VecDestroy(&coordinatesNew);
331: }
332: DMPlexCopy_Internal(dm, PETSC_TRUE, *pdm);
333: (*pdm)->setupcalled = PETSC_TRUE;
334: return 0;
335: }