xref: /petsc/src/dm/impls/plex/plexreorder.c (revision 6913077d92d45dc85e2cd558d8cf01c8c45125fc) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 #include <petsc/private/matorderimpl.h> /*I      "petscmat.h"      I*/
3 
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;
8   PetscErrorCode ierr;
9 
10   PetscFunctionBegin;
11   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
12   ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
13   ierr = PetscMalloc1(pEnd-pStart,&perm);CHKERRQ(ierr);
14   ierr = PetscMalloc1(pEnd-pStart,&iperm);CHKERRQ(ierr);
15   for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
16   for (d = depth; d > 0; --d) {
17     ierr = DMPlexGetDepthStratum(dm, d,   &pStart, &pEnd);CHKERRQ(ierr);
18     ierr = DMPlexGetDepthStratum(dm, d-1, &fStart, &fEnd);CHKERRQ(ierr);
19     fMax = fStart;
20     for (p = pStart; p < pEnd; ++p) {
21       const PetscInt *cone;
22       PetscInt        point, coneSize, c;
23 
24       if (d == depth) {
25         perm[p]         = pperm[p];
26         iperm[pperm[p]] = p;
27       }
28       point = perm[p];
29       ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr);
30       ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr);
31       for (c = 0; c < coneSize; ++c) {
32         const PetscInt oldc = cone[c];
33         const PetscInt newc = iperm[oldc];
34 
35         if (newc < 0) {
36           perm[fMax]  = oldc;
37           iperm[oldc] = fMax++;
38         }
39       }
40     }
41     PetscCheckFalse(fMax != fEnd,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   PetscFunctionReturn(0);
46 }
47 
48 /*@
49   DMPlexGetOrdering - Calculate a reordering of the mesh
50 
51   Collective on dm
52 
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
62 
63   Output Parameter:
64 . perm - The point permutation as an IS, perm[old point number] = new point number
65 
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.
68 
69   Level: intermediate
70 
71 .seealso: DMPlexPermute(), 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;
77   PetscErrorCode ierr;
78 
79   PetscFunctionBegin;
80   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
81   PetscValidPointer(perm, 4);
82   ierr = DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);CHKERRQ(ierr);
83   ierr = PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);CHKERRQ(ierr);
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     ierr = SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);CHKERRQ(ierr);
89   }
90   ierr = PetscFree(start);CHKERRQ(ierr);
91   ierr = PetscFree(adjacency);CHKERRQ(ierr);
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;
100 
101     ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr);
102     ierr = ISSort(valueIS);CHKERRQ(ierr);
103     ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr);
104     ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr);
105     ierr = PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff);CHKERRQ(ierr);
106     for (v = 0; v < numValues; ++v) {
107       ierr = DMLabelGetStratumSize(label, values[v], &vsize[v]);CHKERRQ(ierr);
108       if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1];
109       numPoints += vsize[v];
110     }
111     PetscCheckFalse(numPoints != numCells,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;
115 
116       ierr = DMLabelGetValue(label, oldc, &val);CHKERRQ(ierr);
117       PetscCheckFalse(val == -1,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %D not present in label", oldc);
118       ierr = PetscFindInt(val, numValues, values, &vloc);CHKERRQ(ierr);
119       PetscCheckFalse(vloc < 0,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       PetscCheckFalse(voff[v+1] - voff[v] != vsize[v],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     ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr);
126     ierr = ISDestroy(&valueIS);CHKERRQ(ierr);
127     ierr = PetscArraycpy(cperm, sperm, numCells);CHKERRQ(ierr);
128     ierr = PetscFree3(sperm, vsize, voff);CHKERRQ(ierr);
129   }
130   /* Construct closure */
131   ierr = DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);CHKERRQ(ierr);
132   ierr = PetscFree3(cperm,mask,xls);CHKERRQ(ierr);
133   ierr = PetscFree(clperm);CHKERRQ(ierr);
134   /* Invert permutation */
135   ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
136   ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);CHKERRQ(ierr);
137   PetscFunctionReturn(0);
138 }
139 
140 /*@
141   DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line
142 
143   Collective on dm
144 
145   Input Parameter:
146 . dm - The DMPlex object
147 
148   Output Parameter:
149 . perm - The point permutation as an IS, perm[old point number] = new point number
150 
151   Level: intermediate
152 
153 .seealso: DMPlexGetOrdering(), DMPlexPermute(), MatGetOrdering()
154 @*/
155 PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm)
156 {
157   PetscInt       *points;
158   const PetscInt *support, *cone;
159   PetscInt        dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0;
160   PetscErrorCode  ierr;
161 
162   PetscFunctionBegin;
163   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
164   PetscCheck(dim == 1, PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Input mesh must be one dimensional, not %" PetscInt_FMT, dim);
165   ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
166   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
167   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
168   ierr = PetscMalloc1(pEnd-pStart, &points);CHKERRQ(ierr);
169   for (c = cStart; c < cEnd; ++c) points[c] = c;
170   for (v = vStart; v < vEnd; ++v) points[v] = v;
171   for (v = vStart; v < vEnd; ++v) {
172     ierr = DMPlexGetSupportSize(dm, v, &suppSize);CHKERRQ(ierr);
173     ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr);
174     if (suppSize == 1) {lastCell = support[0]; break;}
175   }
176   if (v < vEnd) {
177     PetscInt pos = cEnd;
178 
179     points[v] = pos++;
180     while (lastCell >= cStart) {
181       ierr = DMPlexGetCone(dm, lastCell, &cone);CHKERRQ(ierr);
182       if (cone[0] == v) v = cone[1];
183       else              v = cone[0];
184       ierr = DMPlexGetSupport(dm, v, &support);CHKERRQ(ierr);
185       ierr = DMPlexGetSupportSize(dm, v, &suppSize);CHKERRQ(ierr);
186       if (suppSize == 1) {lastCell = -1;}
187       else {
188         if (support[0] == lastCell) lastCell = support[1];
189         else                        lastCell = support[0];
190       }
191       points[v] = pos++;
192     }
193     PetscCheck(pos == pEnd, PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Last vertex was %" PetscInt_FMT ", not %" PetscInt_FMT, pos, pEnd);
194   }
195   ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, points, PETSC_OWN_POINTER, perm);CHKERRQ(ierr);
196   PetscFunctionReturn(0);
197 }
198 
199 /*@
200   DMPlexPermute - Reorder the mesh according to the input permutation
201 
202   Collective on dm
203 
204   Input Parameters:
205 + dm - The DMPlex object
206 - perm - The point permutation, perm[old point number] = new point number
207 
208   Output Parameter:
209 . pdm - The permuted DM
210 
211   Level: intermediate
212 
213 .seealso: MatPermute()
214 @*/
215 PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
216 {
217   DM_Plex       *plex = (DM_Plex *) dm->data, *plexNew;
218   PetscInt       dim, cdim;
219   const char    *name;
220   PetscErrorCode ierr;
221 
222   PetscFunctionBegin;
223   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
224   PetscValidHeaderSpecific(perm, IS_CLASSID, 2);
225   PetscValidPointer(pdm, 3);
226   ierr = DMCreate(PetscObjectComm((PetscObject) dm), pdm);CHKERRQ(ierr);
227   ierr = DMSetType(*pdm, DMPLEX);CHKERRQ(ierr);
228   ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr);
229   ierr = PetscObjectSetName((PetscObject) *pdm, name);CHKERRQ(ierr);
230   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
231   ierr = DMSetDimension(*pdm, dim);CHKERRQ(ierr);
232   ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr);
233   ierr = DMSetCoordinateDim(*pdm, cdim);CHKERRQ(ierr);
234   ierr = DMCopyDisc(dm, *pdm);CHKERRQ(ierr);
235   if (dm->localSection) {
236     PetscSection section, sectionNew;
237 
238     ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
239     ierr = PetscSectionPermute(section, perm, &sectionNew);CHKERRQ(ierr);
240     ierr = DMSetLocalSection(*pdm, sectionNew);CHKERRQ(ierr);
241     ierr = PetscSectionDestroy(&sectionNew);CHKERRQ(ierr);
242   }
243   plexNew = (DM_Plex *) (*pdm)->data;
244   /* Ignore ltogmap, ltogmapb */
245   /* Ignore sf, sectionSF */
246   /* Ignore globalVertexNumbers, globalCellNumbers */
247   /* Reorder labels */
248   {
249     PetscInt numLabels, l;
250     DMLabel  label, labelNew;
251 
252     ierr = DMGetNumLabels(dm, &numLabels);CHKERRQ(ierr);
253     for (l = 0; l < numLabels; ++l) {
254       ierr = DMGetLabelByNum(dm, l, &label);CHKERRQ(ierr);
255       ierr = DMLabelPermute(label, perm, &labelNew);CHKERRQ(ierr);
256       ierr = DMAddLabel(*pdm, labelNew);CHKERRQ(ierr);
257       ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr);
258     }
259     ierr = DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel);CHKERRQ(ierr);
260     if (plex->subpointMap) {ierr = DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);CHKERRQ(ierr);}
261   }
262   /* Reorder topology */
263   {
264     const PetscInt *pperm;
265     PetscInt        maxConeSize, maxSupportSize, n, pStart, pEnd, p;
266 
267     ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
268     plexNew->maxConeSize    = maxConeSize;
269     plexNew->maxSupportSize = maxSupportSize;
270     ierr = PetscSectionDestroy(&plexNew->coneSection);CHKERRQ(ierr);
271     ierr = PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);CHKERRQ(ierr);
272     ierr = PetscSectionGetStorageSize(plexNew->coneSection, &n);CHKERRQ(ierr);
273     ierr = PetscMalloc1(n, &plexNew->cones);CHKERRQ(ierr);
274     ierr = PetscMalloc1(n, &plexNew->coneOrientations);CHKERRQ(ierr);
275     ierr = ISGetIndices(perm, &pperm);CHKERRQ(ierr);
276     ierr = PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);CHKERRQ(ierr);
277     for (p = pStart; p < pEnd; ++p) {
278       PetscInt dof, off, offNew, d;
279 
280       ierr = PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);CHKERRQ(ierr);
281       ierr = PetscSectionGetOffset(plex->coneSection, p, &off);CHKERRQ(ierr);
282       ierr = PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);CHKERRQ(ierr);
283       for (d = 0; d < dof; ++d) {
284         plexNew->cones[offNew+d]            = pperm[plex->cones[off+d]];
285         plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
286       }
287     }
288     ierr = PetscSectionDestroy(&plexNew->supportSection);CHKERRQ(ierr);
289     ierr = PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);CHKERRQ(ierr);
290     ierr = PetscSectionGetStorageSize(plexNew->supportSection, &n);CHKERRQ(ierr);
291     ierr = PetscMalloc1(n, &plexNew->supports);CHKERRQ(ierr);
292     ierr = PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);CHKERRQ(ierr);
293     for (p = pStart; p < pEnd; ++p) {
294       PetscInt dof, off, offNew, d;
295 
296       ierr = PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);CHKERRQ(ierr);
297       ierr = PetscSectionGetOffset(plex->supportSection, p, &off);CHKERRQ(ierr);
298       ierr = PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);CHKERRQ(ierr);
299       for (d = 0; d < dof; ++d) {
300         plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
301       }
302     }
303     ierr = ISRestoreIndices(perm, &pperm);CHKERRQ(ierr);
304   }
305   /* Remap coordinates */
306   {
307     DM              cdm, cdmNew;
308     PetscSection    csection, csectionNew;
309     Vec             coordinates, coordinatesNew;
310     PetscScalar    *coords, *coordsNew;
311     const PetscInt *pperm;
312     PetscInt        pStart, pEnd, p;
313     const char     *name;
314 
315     ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
316     ierr = DMGetLocalSection(cdm, &csection);CHKERRQ(ierr);
317     ierr = PetscSectionPermute(csection, perm, &csectionNew);CHKERRQ(ierr);
318     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
319     ierr = VecDuplicate(coordinates, &coordinatesNew);CHKERRQ(ierr);
320     ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr);
321     ierr = PetscObjectSetName((PetscObject)coordinatesNew,name);CHKERRQ(ierr);
322     ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
323     ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr);
324     ierr = PetscSectionGetChart(csectionNew, &pStart, &pEnd);CHKERRQ(ierr);
325     ierr = ISGetIndices(perm, &pperm);CHKERRQ(ierr);
326     for (p = pStart; p < pEnd; ++p) {
327       PetscInt dof, off, offNew, d;
328 
329       ierr = PetscSectionGetDof(csectionNew, p, &dof);CHKERRQ(ierr);
330       ierr = PetscSectionGetOffset(csection, p, &off);CHKERRQ(ierr);
331       ierr = PetscSectionGetOffset(csectionNew, pperm[p], &offNew);CHKERRQ(ierr);
332       for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
333     }
334     ierr = ISRestoreIndices(perm, &pperm);CHKERRQ(ierr);
335     ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
336     ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr);
337     ierr = DMGetCoordinateDM(*pdm, &cdmNew);CHKERRQ(ierr);
338     ierr = DMSetLocalSection(cdmNew, csectionNew);CHKERRQ(ierr);
339     ierr = DMSetCoordinatesLocal(*pdm, coordinatesNew);CHKERRQ(ierr);
340     ierr = PetscSectionDestroy(&csectionNew);CHKERRQ(ierr);
341     ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr);
342   }
343   ierr = DMPlexCopy_Internal(dm, PETSC_TRUE, *pdm);CHKERRQ(ierr);
344   (*pdm)->setupcalled = PETSC_TRUE;
345   PetscFunctionReturn(0);
346 }
347