Actual source code: plexreorder.c
petsc-3.7.3 2016-08-01
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc/private/matorderimpl.h> /*I "petscmat.h" I*/
6: PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm)
7: {
8: PetscInt *perm, *iperm;
9: PetscInt depth, d, pStart, pEnd, fStart, fMax, fEnd, p;
13: DMPlexGetDepth(dm, &depth);
14: DMPlexGetChart(dm, &pStart, &pEnd);
15: PetscMalloc1(pEnd-pStart,&perm);
16: PetscMalloc1(pEnd-pStart,&iperm);
17: for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
18: for (d = depth; d > 0; --d) {
19: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
20: DMPlexGetDepthStratum(dm, d-1, &fStart, &fEnd);
21: fMax = fStart;
22: for (p = pStart; p < pEnd; ++p) {
23: const PetscInt *cone;
24: PetscInt point, coneSize, c;
26: if (d == depth) {
27: perm[p] = pperm[p];
28: iperm[pperm[p]] = p;
29: }
30: point = perm[p];
31: DMPlexGetConeSize(dm, point, &coneSize);
32: DMPlexGetCone(dm, point, &cone);
33: for (c = 0; c < coneSize; ++c) {
34: const PetscInt oldc = cone[c];
35: const PetscInt newc = iperm[oldc];
37: if (newc < 0) {
38: perm[fMax] = oldc;
39: iperm[oldc] = fMax++;
40: }
41: }
42: }
43: if (fMax != fEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %d faces %d does not match permuted nubmer %d", d, fEnd-fStart, fMax-fStart);
44: }
45: *clperm = perm;
46: *invclperm = iperm;
47: return(0);
48: }
52: /*@
53: DMPlexGetOrdering - Calculate a reordering of the mesh
55: Collective on DM
57: Input Parameter:
58: + dm - The DMPlex object
59: . otype - type of reordering, one of the following:
60: $ MATORDERINGNATURAL - Natural
61: $ MATORDERINGND - Nested Dissection
62: $ MATORDERING1WD - One-way Dissection
63: $ MATORDERINGRCM - Reverse Cuthill-McKee
64: $ MATORDERINGQMD - Quotient Minimum Degree
65: - label - [Optional] Label used to segregate ordering into sets, or NULL
68: Output Parameter:
69: . perm - The point permutation as an IS, perm[old point number] = new point number
71: Note: The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which
72: has different types of cells, and then loop over each set of reordered cells for assembly.
74: Level: intermediate
76: .keywords: mesh
77: .seealso: MatGetOrdering()
78: @*/
79: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
80: {
81: PetscInt numCells = 0;
82: PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm, *invclperm, *mask, *xls, pStart, pEnd, c, i;
88: DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);
89: PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);
90: if (numCells) {
91: /* Shift for Fortran numbering */
92: for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
93: for (i = 0; i <= numCells; ++i) ++start[i];
94: SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);
95: }
96: PetscFree(start);
97: PetscFree(adjacency);
98: /* Shift for Fortran numbering */
99: for (c = 0; c < numCells; ++c) --cperm[c];
100: /* Segregate */
101: if (label) {
102: IS valueIS;
103: const PetscInt *values;
104: PetscInt numValues, numPoints = 0;
105: PetscInt *sperm, *vsize, *voff, v;
107: DMLabelGetValueIS(label, &valueIS);
108: ISGetLocalSize(valueIS, &numValues);
109: ISGetIndices(valueIS, &values);
110: PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff);
111: for (v = 0; v < numValues; ++v) {
112: DMLabelGetStratumSize(label, values[v], &vsize[v]);
113: if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1];
114: numPoints += vsize[v];
115: }
116: if (numPoints != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %D cells < %D total", numPoints, numCells);
117: for (c = 0; c < numCells; ++c) {
118: const PetscInt oldc = cperm[c];
119: PetscInt val, vloc;
121: DMLabelGetValue(label, oldc, &val);
122: if (val == -1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %D not present in label", oldc);
123: PetscFindInt(val, numValues, values, &vloc);
124: if (vloc < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %D not present label", val);
125: sperm[voff[vloc+1]++] = oldc;
126: }
127: for (v = 0; v < numValues; ++v) {
128: 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]);
129: }
130: ISRestoreIndices(valueIS, &values);
131: ISDestroy(&valueIS);
132: PetscMemcpy(cperm, sperm, numCells * sizeof(PetscInt));
133: PetscFree3(sperm, vsize, voff);
134: }
135: /* Construct closure */
136: DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
137: PetscFree3(cperm,mask,xls);
138: PetscFree(clperm);
139: /* Invert permutation */
140: DMPlexGetChart(dm, &pStart, &pEnd);
141: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);
142: return(0);
143: }
147: /*@
148: DMPlexPermute - Reorder the mesh according to the input permutation
150: Collective on DM
152: Input Parameter:
153: + dm - The DMPlex object
154: - perm - The point permutation, perm[old point number] = new point number
156: Output Parameter:
157: . pdm - The permuted DM
159: Level: intermediate
161: .keywords: mesh
162: .seealso: MatPermute()
163: @*/
164: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
165: {
166: DM_Plex *plex = (DM_Plex *) dm->data, *plexNew;
167: PetscSection section, sectionNew;
168: PetscInt dim;
175: DMCreate(PetscObjectComm((PetscObject) dm), pdm);
176: DMSetType(*pdm, DMPLEX);
177: DMGetDimension(dm, &dim);
178: DMSetDimension(*pdm, dim);
179: DMGetDefaultSection(dm, §ion);
180: if (section) {
181: PetscSectionPermute(section, perm, §ionNew);
182: DMSetDefaultSection(*pdm, sectionNew);
183: PetscSectionDestroy(§ionNew);
184: }
185: plexNew = (DM_Plex *) (*pdm)->data;
186: /* Ignore ltogmap, ltogmapb */
187: /* Ignore sf, defaultSF */
188: /* Ignore globalVertexNumbers, globalCellNumbers */
189: /* Remap coordinates */
190: {
191: DM cdm, cdmNew;
192: PetscSection csection, csectionNew;
193: Vec coordinates, coordinatesNew;
194: PetscScalar *coords, *coordsNew;
195: const PetscInt *pperm;
196: PetscInt pStart, pEnd, p;
197: const char *name;
199: DMGetCoordinateDM(dm, &cdm);
200: DMGetDefaultSection(cdm, &csection);
201: PetscSectionPermute(csection, perm, &csectionNew);
202: DMGetCoordinatesLocal(dm, &coordinates);
203: VecDuplicate(coordinates, &coordinatesNew);
204: PetscObjectGetName((PetscObject)coordinates,&name);
205: PetscObjectSetName((PetscObject)coordinatesNew,name);
206: VecGetArray(coordinates, &coords);
207: VecGetArray(coordinatesNew, &coordsNew);
208: PetscSectionGetChart(csectionNew, &pStart, &pEnd);
209: ISGetIndices(perm, &pperm);
210: for (p = pStart; p < pEnd; ++p) {
211: PetscInt dof, off, offNew, d;
213: PetscSectionGetDof(csectionNew, p, &dof);
214: PetscSectionGetOffset(csection, p, &off);
215: PetscSectionGetOffset(csectionNew, pperm[p], &offNew);
216: for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
217: }
218: ISRestoreIndices(perm, &pperm);
219: VecRestoreArray(coordinates, &coords);
220: VecRestoreArray(coordinatesNew, &coordsNew);
221: DMGetCoordinateDM(*pdm, &cdmNew);
222: DMSetDefaultSection(cdmNew, csectionNew);
223: DMSetCoordinatesLocal(*pdm, coordinatesNew);
224: PetscSectionDestroy(&csectionNew);
225: VecDestroy(&coordinatesNew);
226: }
227: /* Reorder labels */
228: {
229: PetscInt numLabels, l;
230: DMLabel label, labelNew;
232: DMGetNumLabels(dm, &numLabels);
233: for (l = numLabels-1; l >= 0; --l) {
234: DMGetLabelByNum(dm, l, &label);
235: DMLabelPermute(label, perm, &labelNew);
236: DMAddLabel(*pdm, labelNew);
237: }
238: if (plex->subpointMap) {DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);}
239: }
240: /* Reorder topology */
241: {
242: const PetscInt *pperm;
243: PetscInt maxConeSize, maxSupportSize, n, pStart, pEnd, p;
245: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
246: plexNew->maxConeSize = maxConeSize;
247: plexNew->maxSupportSize = maxSupportSize;
248: PetscSectionDestroy(&plexNew->coneSection);
249: PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
250: PetscSectionGetStorageSize(plexNew->coneSection, &n);
251: PetscMalloc1(n, &plexNew->cones);
252: PetscMalloc1(n, &plexNew->coneOrientations);
253: ISGetIndices(perm, &pperm);
254: PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
255: for (p = pStart; p < pEnd; ++p) {
256: PetscInt dof, off, offNew, d;
258: PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
259: PetscSectionGetOffset(plex->coneSection, p, &off);
260: PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
261: for (d = 0; d < dof; ++d) {
262: plexNew->cones[offNew+d] = pperm[plex->cones[off+d]];
263: plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
264: }
265: }
266: PetscSectionDestroy(&plexNew->supportSection);
267: PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
268: PetscSectionGetStorageSize(plexNew->supportSection, &n);
269: PetscMalloc1(n, &plexNew->supports);
270: PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
271: for (p = pStart; p < pEnd; ++p) {
272: PetscInt dof, off, offNew, d;
274: PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
275: PetscSectionGetOffset(plex->supportSection, p, &off);
276: PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
277: for (d = 0; d < dof; ++d) {
278: plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
279: }
280: }
281: ISRestoreIndices(perm, &pperm);
282: }
283: return(0);
284: }