Actual source code: plexreorder.c
petsc-3.5.4 2015-05-23
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: PetscMalloc((pEnd-pStart) * sizeof(PetscInt) ,&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: DMPlexGetDimension(dm, &dim);
139: DMPlexSetDimension(*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: DMLabel label = plex->labels, labelNew = NULL;
189: while (label) {
190: if (!plexNew->labels) {
191: DMLabelPermute(label, perm, &plexNew->labels);
192: labelNew = plexNew->labels;
193: } else {
194: DMLabelPermute(label, perm, &labelNew->next);
195: labelNew = labelNew->next;
196: }
197: label = label->next;
198: }
199: if (plex->subpointMap) {DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);}
200: }
201: /* Reorder topology */
202: {
203: const PetscInt *pperm;
204: PetscInt maxConeSize, maxSupportSize, n, pStart, pEnd, p;
206: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
207: plexNew->maxConeSize = maxConeSize;
208: plexNew->maxSupportSize = maxSupportSize;
209: PetscSectionDestroy(&plexNew->coneSection);
210: PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
211: PetscSectionGetStorageSize(plexNew->coneSection, &n);
212: PetscMalloc1(n, &plexNew->cones);
213: PetscMalloc1(n, &plexNew->coneOrientations);
214: ISGetIndices(perm, &pperm);
215: PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
216: for (p = pStart; p < pEnd; ++p) {
217: PetscInt dof, off, offNew, d;
219: PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
220: PetscSectionGetOffset(plex->coneSection, p, &off);
221: PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
222: for (d = 0; d < dof; ++d) {
223: plexNew->cones[offNew+d] = pperm[plex->cones[off+d]];
224: plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
225: }
226: }
227: PetscSectionDestroy(&plexNew->supportSection);
228: PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
229: PetscSectionGetStorageSize(plexNew->supportSection, &n);
230: PetscMalloc1(n, &plexNew->supports);
231: PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
232: for (p = pStart; p < pEnd; ++p) {
233: PetscInt dof, off, offNew, d;
235: PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
236: PetscSectionGetOffset(plex->supportSection, p, &off);
237: PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
238: for (d = 0; d < dof; ++d) {
239: plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
240: }
241: }
242: ISRestoreIndices(perm, &pperm);
243: }
244: return(0);
245: }