Actual source code: plexreorder.c
petsc-3.6.1 2015-08-06
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
67: Output Parameter:
68: . perm - The point permutation as an IS
70: Level: intermediate
72: .keywords: mesh
73: .seealso: MatGetOrdering()
74: @*/
75: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, IS *perm)
76: {
77: PetscInt numCells = 0;
78: PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm, *invclperm, *mask, *xls, pStart, pEnd, c, i;
84: DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);
85: PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);
86: if (numCells) {
87: /* Shift for Fortran numbering */
88: for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
89: for (i = 0; i <= numCells; ++i) ++start[i];
90: SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);
91: }
92: PetscFree(start);
93: PetscFree(adjacency);
94: /* Shift for Fortran numbering */
95: for (c = 0; c < numCells; ++c) --cperm[c];
96: /* Construct closure */
97: DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
98: PetscFree3(cperm,mask,xls);
99: PetscFree(clperm);
100: /* Invert permutation */
101: DMPlexGetChart(dm, &pStart, &pEnd);
102: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);
103: return(0);
104: }
108: /*@
109: DMPlexPermute - Reorder the mesh according to the input permutation
111: Collective on DM
113: Input Parameter:
114: + dm - The DMPlex object
115: - perm - The point permutation
117: Output Parameter:
118: . pdm - The permuted DM
120: Level: intermediate
122: .keywords: mesh
123: .seealso: MatPermute()
124: @*/
125: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
126: {
127: DM_Plex *plex = (DM_Plex *) dm->data, *plexNew;
128: PetscSection section, sectionNew;
129: PetscInt dim;
136: DMCreate(PetscObjectComm((PetscObject) dm), pdm);
137: DMSetType(*pdm, DMPLEX);
138: DMGetDimension(dm, &dim);
139: DMSetDimension(*pdm, dim);
140: DMGetDefaultSection(dm, §ion);
141: if (section) {
142: PetscSectionPermute(section, perm, §ionNew);
143: DMSetDefaultSection(*pdm, sectionNew);
144: PetscSectionDestroy(§ionNew);
145: }
146: plexNew = (DM_Plex *) (*pdm)->data;
147: /* Ignore ltogmap, ltogmapb */
148: /* Ignore sf, defaultSF */
149: /* Ignore globalVertexNumbers, globalCellNumbers */
150: /* Remap coordinates */
151: {
152: DM cdm, cdmNew;
153: PetscSection csection, csectionNew;
154: Vec coordinates, coordinatesNew;
155: PetscScalar *coords, *coordsNew;
156: const PetscInt *pperm;
157: PetscInt pStart, pEnd, p;
159: DMGetCoordinateDM(dm, &cdm);
160: DMGetDefaultSection(cdm, &csection);
161: PetscSectionPermute(csection, perm, &csectionNew);
162: DMGetCoordinatesLocal(dm, &coordinates);
163: VecDuplicate(coordinates, &coordinatesNew);
164: VecGetArray(coordinates, &coords);
165: VecGetArray(coordinatesNew, &coordsNew);
166: PetscSectionGetChart(csectionNew, &pStart, &pEnd);
167: ISGetIndices(perm, &pperm);
168: for (p = pStart; p < pEnd; ++p) {
169: PetscInt dof, off, offNew, d;
171: PetscSectionGetDof(csectionNew, p, &dof);
172: PetscSectionGetOffset(csection, p, &off);
173: PetscSectionGetOffset(csectionNew, pperm[p], &offNew);
174: for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
175: }
176: ISRestoreIndices(perm, &pperm);
177: VecRestoreArray(coordinates, &coords);
178: VecRestoreArray(coordinatesNew, &coordsNew);
179: DMGetCoordinateDM(*pdm, &cdmNew);
180: DMSetDefaultSection(cdmNew, csectionNew);
181: DMSetCoordinatesLocal(*pdm, coordinatesNew);
182: PetscSectionDestroy(&csectionNew);
183: VecDestroy(&coordinatesNew);
184: }
185: /* Reorder labels */
186: {
187: PetscInt numLabels, l;
188: DMLabel label, labelNew;
190: DMPlexGetNumLabels(dm, &numLabels);
191: for (l = numLabels-1; l >= 0; --l) {
192: DMPlexGetLabelByNum(dm, l, &label);
193: DMLabelPermute(label, perm, &labelNew);
194: DMPlexAddLabel(*pdm, labelNew);
195: }
196: if (plex->subpointMap) {DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);}
197: }
198: /* Reorder topology */
199: {
200: const PetscInt *pperm;
201: PetscInt maxConeSize, maxSupportSize, n, pStart, pEnd, p;
203: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
204: plexNew->maxConeSize = maxConeSize;
205: plexNew->maxSupportSize = maxSupportSize;
206: PetscSectionDestroy(&plexNew->coneSection);
207: PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
208: PetscSectionGetStorageSize(plexNew->coneSection, &n);
209: PetscMalloc1(n, &plexNew->cones);
210: PetscMalloc1(n, &plexNew->coneOrientations);
211: ISGetIndices(perm, &pperm);
212: PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
213: for (p = pStart; p < pEnd; ++p) {
214: PetscInt dof, off, offNew, d;
216: PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
217: PetscSectionGetOffset(plex->coneSection, p, &off);
218: PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
219: for (d = 0; d < dof; ++d) {
220: plexNew->cones[offNew+d] = pperm[plex->cones[off+d]];
221: plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
222: }
223: }
224: PetscSectionDestroy(&plexNew->supportSection);
225: PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
226: PetscSectionGetStorageSize(plexNew->supportSection, &n);
227: PetscMalloc1(n, &plexNew->supports);
228: PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
229: for (p = pStart; p < pEnd; ++p) {
230: PetscInt dof, off, offNew, d;
232: PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
233: PetscSectionGetOffset(plex->supportSection, p, &off);
234: PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
235: for (d = 0; d < dof; ++d) {
236: plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
237: }
238: }
239: ISRestoreIndices(perm, &pperm);
240: }
241: return(0);
242: }