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 Parameters:
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
63: Output Parameter:
64: . perm - The point permutation as an IS, perm[old point number] = new point number
66: Note: The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which
67: has different types of cells, and then loop over each set of reordered cells for assembly.
69: Level: intermediate
71: .seealso: MatGetOrdering()
72: @*/
73: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
74: {
75: PetscInt numCells = 0;
76: PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i;
82: DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);
83: PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);
84: if (numCells) {
85: /* Shift for Fortran numbering */
86: for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
87: for (i = 0; i <= numCells; ++i) ++start[i];
88: SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);
89: }
90: PetscFree(start);
91: PetscFree(adjacency);
92: /* Shift for Fortran numbering */
93: for (c = 0; c < numCells; ++c) --cperm[c];
94: /* Segregate */
95: if (label) {
96: IS valueIS;
97: const PetscInt *values;
98: PetscInt numValues, numPoints = 0;
99: PetscInt *sperm, *vsize, *voff, v;
101: DMLabelGetValueIS(label, &valueIS);
102: ISSort(valueIS);
103: ISGetLocalSize(valueIS, &numValues);
104: ISGetIndices(valueIS, &values);
105: PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff);
106: for (v = 0; v < numValues; ++v) {
107: 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: if (numPoints != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %D cells < %D total", numPoints, numCells);
112: for (c = 0; c < numCells; ++c) {
113: const PetscInt oldc = cperm[c];
114: PetscInt val, vloc;
116: DMLabelGetValue(label, oldc, &val);
117: if (val == -1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %D not present in label", oldc);
118: PetscFindInt(val, numValues, values, &vloc);
119: if (vloc < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %D not present label", val);
120: sperm[voff[vloc+1]++] = oldc;
121: }
122: for (v = 0; v < numValues; ++v) {
123: 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]);
124: }
125: ISRestoreIndices(valueIS, &values);
126: ISDestroy(&valueIS);
127: PetscArraycpy(cperm, sperm, numCells);
128: PetscFree3(sperm, vsize, voff);
129: }
130: /* Construct closure */
131: DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
132: PetscFree3(cperm,mask,xls);
133: PetscFree(clperm);
134: /* Invert permutation */
135: DMPlexGetChart(dm, &pStart, &pEnd);
136: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);
137: return(0);
138: }
140: /*@
141: DMPlexPermute - Reorder the mesh according to the input permutation
143: Collective on dm
145: Input Parameters:
146: + dm - The DMPlex object
147: - perm - The point permutation, perm[old point number] = new point number
149: Output Parameter:
150: . pdm - The permuted DM
152: Level: intermediate
154: .seealso: MatPermute()
155: @*/
156: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
157: {
158: DM_Plex *plex = (DM_Plex *) dm->data, *plexNew;
159: PetscSection section, sectionNew;
160: PetscInt dim, cdim;
167: DMCreate(PetscObjectComm((PetscObject) dm), pdm);
168: DMSetType(*pdm, DMPLEX);
169: DMGetDimension(dm, &dim);
170: DMSetDimension(*pdm, dim);
171: DMGetCoordinateDim(dm, &cdim);
172: DMSetCoordinateDim(*pdm, cdim);
173: DMCopyDisc(dm, *pdm);
174: DMGetLocalSection(dm, §ion);
175: if (section) {
176: PetscSectionPermute(section, perm, §ionNew);
177: DMSetLocalSection(*pdm, sectionNew);
178: PetscSectionDestroy(§ionNew);
179: }
180: plexNew = (DM_Plex *) (*pdm)->data;
181: /* Ignore ltogmap, ltogmapb */
182: /* Ignore sf, sectionSF */
183: /* Ignore globalVertexNumbers, globalCellNumbers */
184: /* Reorder labels */
185: {
186: PetscInt numLabels, l;
187: DMLabel label, labelNew;
189: DMGetNumLabels(dm, &numLabels);
190: for (l = 0; l < numLabels; ++l) {
191: DMGetLabelByNum(dm, l, &label);
192: DMLabelPermute(label, perm, &labelNew);
193: DMAddLabel(*pdm, labelNew);
194: DMLabelDestroy(&labelNew);
195: }
196: DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel);
197: if (plex->subpointMap) {DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);}
198: }
199: /* Reorder topology */
200: {
201: const PetscInt *pperm;
202: PetscInt maxConeSize, maxSupportSize, n, pStart, pEnd, p;
204: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
205: plexNew->maxConeSize = maxConeSize;
206: plexNew->maxSupportSize = maxSupportSize;
207: PetscSectionDestroy(&plexNew->coneSection);
208: PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
209: PetscSectionGetStorageSize(plexNew->coneSection, &n);
210: PetscMalloc1(n, &plexNew->cones);
211: PetscMalloc1(n, &plexNew->coneOrientations);
212: ISGetIndices(perm, &pperm);
213: PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
214: for (p = pStart; p < pEnd; ++p) {
215: PetscInt dof, off, offNew, d;
217: PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
218: PetscSectionGetOffset(plex->coneSection, p, &off);
219: PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
220: for (d = 0; d < dof; ++d) {
221: plexNew->cones[offNew+d] = pperm[plex->cones[off+d]];
222: plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
223: }
224: }
225: PetscSectionDestroy(&plexNew->supportSection);
226: PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
227: PetscSectionGetStorageSize(plexNew->supportSection, &n);
228: PetscMalloc1(n, &plexNew->supports);
229: PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
230: for (p = pStart; p < pEnd; ++p) {
231: PetscInt dof, off, offNew, d;
233: PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
234: PetscSectionGetOffset(plex->supportSection, p, &off);
235: PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
236: for (d = 0; d < dof; ++d) {
237: plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
238: }
239: }
240: ISRestoreIndices(perm, &pperm);
241: }
242: /* Remap coordinates */
243: {
244: DM cdm, cdmNew;
245: PetscSection csection, csectionNew;
246: Vec coordinates, coordinatesNew;
247: PetscScalar *coords, *coordsNew;
248: const PetscInt *pperm;
249: PetscInt pStart, pEnd, p;
250: const char *name;
252: DMGetCoordinateDM(dm, &cdm);
253: DMGetLocalSection(cdm, &csection);
254: PetscSectionPermute(csection, perm, &csectionNew);
255: DMGetCoordinatesLocal(dm, &coordinates);
256: VecDuplicate(coordinates, &coordinatesNew);
257: PetscObjectGetName((PetscObject)coordinates,&name);
258: PetscObjectSetName((PetscObject)coordinatesNew,name);
259: VecGetArray(coordinates, &coords);
260: VecGetArray(coordinatesNew, &coordsNew);
261: PetscSectionGetChart(csectionNew, &pStart, &pEnd);
262: ISGetIndices(perm, &pperm);
263: for (p = pStart; p < pEnd; ++p) {
264: PetscInt dof, off, offNew, d;
266: PetscSectionGetDof(csectionNew, p, &dof);
267: PetscSectionGetOffset(csection, p, &off);
268: PetscSectionGetOffset(csectionNew, pperm[p], &offNew);
269: for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
270: }
271: ISRestoreIndices(perm, &pperm);
272: VecRestoreArray(coordinates, &coords);
273: VecRestoreArray(coordinatesNew, &coordsNew);
274: DMGetCoordinateDM(*pdm, &cdmNew);
275: DMSetLocalSection(cdmNew, csectionNew);
276: DMSetCoordinatesLocal(*pdm, coordinatesNew);
277: PetscSectionDestroy(&csectionNew);
278: VecDestroy(&coordinatesNew);
279: }
280: (*pdm)->setupcalled = PETSC_TRUE;
281: return(0);
282: }