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;

  9:   DMPlexGetDepth(dm, &depth);
 10:   DMPlexGetChart(dm, &pStart, &pEnd);
 11:   PetscMalloc1(pEnd-pStart,&perm);
 12:   PetscMalloc1(pEnd-pStart,&iperm);
 13:   for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
 14:   for (d = depth; d > 0; --d) {
 15:     DMPlexGetDepthStratum(dm, d,   &pStart, &pEnd);
 16:     DMPlexGetDepthStratum(dm, d-1, &fStart, &fEnd);
 17:     fMax = fStart;
 18:     for (p = pStart; p < pEnd; ++p) {
 19:       const PetscInt *cone;
 20:       PetscInt        point, coneSize, c;

 22:       if (d == depth) {
 23:         perm[p]         = pperm[p];
 24:         iperm[pperm[p]] = p;
 25:       }
 26:       point = perm[p];
 27:       DMPlexGetConeSize(dm, point, &coneSize);
 28:       DMPlexGetCone(dm, point, &cone);
 29:       for (c = 0; c < coneSize; ++c) {
 30:         const PetscInt oldc = cone[c];
 31:         const PetscInt newc = iperm[oldc];

 33:         if (newc < 0) {
 34:           perm[fMax]  = oldc;
 35:           iperm[oldc] = fMax++;
 36:         }
 37:       }
 38:     }
 40:   }
 41:   *clperm    = perm;
 42:   *invclperm = iperm;
 43:   return 0;
 44: }

 46: /*@
 47:   DMPlexGetOrdering - Calculate a reordering of the mesh

 49:   Collective on dm

 51:   Input Parameters:
 52: + dm - The DMPlex object
 53: . otype - type of reordering, one of the following:
 54: $     MATORDERINGNATURAL - Natural
 55: $     MATORDERINGND - Nested Dissection
 56: $     MATORDERING1WD - One-way Dissection
 57: $     MATORDERINGRCM - Reverse Cuthill-McKee
 58: $     MATORDERINGQMD - Quotient Minimum Degree
 59: - label - [Optional] Label used to segregate ordering into sets, or NULL

 61:   Output Parameter:
 62: . perm - The point permutation as an IS, perm[old point number] = new point number

 64:   Note: The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which
 65:   has different types of cells, and then loop over each set of reordered cells for assembly.

 67:   Level: intermediate

 69: .seealso: DMPlexPermute(), MatGetOrdering()
 70: @*/
 71: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
 72: {
 73:   PetscInt       numCells = 0;
 74:   PetscInt      *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i;

 78:   DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);
 79:   PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);
 80:   if (numCells) {
 81:     /* Shift for Fortran numbering */
 82:     for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
 83:     for (i = 0; i <= numCells; ++i)       ++start[i];
 84:     SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);
 85:   }
 86:   PetscFree(start);
 87:   PetscFree(adjacency);
 88:   /* Shift for Fortran numbering */
 89:   for (c = 0; c < numCells; ++c) --cperm[c];
 90:   /* Segregate */
 91:   if (label) {
 92:     IS              valueIS;
 93:     const PetscInt *values;
 94:     PetscInt        numValues, numPoints = 0;
 95:     PetscInt       *sperm, *vsize, *voff, v;

 97:     DMLabelGetValueIS(label, &valueIS);
 98:     ISSort(valueIS);
 99:     ISGetLocalSize(valueIS, &numValues);
100:     ISGetIndices(valueIS, &values);
101:     PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff);
102:     for (v = 0; v < numValues; ++v) {
103:       DMLabelGetStratumSize(label, values[v], &vsize[v]);
104:       if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1];
105:       numPoints += vsize[v];
106:     }
108:     for (c = 0; c < numCells; ++c) {
109:       const PetscInt oldc = cperm[c];
110:       PetscInt       val, vloc;

112:       DMLabelGetValue(label, oldc, &val);
114:       PetscFindInt(val, numValues, values, &vloc);
116:       sperm[voff[vloc+1]++] = oldc;
117:     }
118:     for (v = 0; v < numValues; ++v) {
120:     }
121:     ISRestoreIndices(valueIS, &values);
122:     ISDestroy(&valueIS);
123:     PetscArraycpy(cperm, sperm, numCells);
124:     PetscFree3(sperm, vsize, voff);
125:   }
126:   /* Construct closure */
127:   DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
128:   PetscFree3(cperm,mask,xls);
129:   PetscFree(clperm);
130:   /* Invert permutation */
131:   DMPlexGetChart(dm, &pStart, &pEnd);
132:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);
133:   return 0;
134: }

136: /*@
137:   DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line

139:   Collective on dm

141:   Input Parameter:
142: . dm - The DMPlex object

144:   Output Parameter:
145: . perm - The point permutation as an IS, perm[old point number] = new point number

147:   Level: intermediate

149: .seealso: DMPlexGetOrdering(), DMPlexPermute(), MatGetOrdering()
150: @*/
151: PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm)
152: {
153:   PetscInt       *points;
154:   const PetscInt *support, *cone;
155:   PetscInt        dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0;

157:   DMGetDimension(dm, &dim);
159:   DMPlexGetChart(dm, &pStart, &pEnd);
160:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
161:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
162:   PetscMalloc1(pEnd-pStart, &points);
163:   for (c = cStart; c < cEnd; ++c) points[c] = c;
164:   for (v = vStart; v < vEnd; ++v) points[v] = v;
165:   for (v = vStart; v < vEnd; ++v) {
166:     DMPlexGetSupportSize(dm, v, &suppSize);
167:     DMPlexGetSupport(dm, v, &support);
168:     if (suppSize == 1) {lastCell = support[0]; break;}
169:   }
170:   if (v < vEnd) {
171:     PetscInt pos = cEnd;

173:     points[v] = pos++;
174:     while (lastCell >= cStart) {
175:       DMPlexGetCone(dm, lastCell, &cone);
176:       if (cone[0] == v) v = cone[1];
177:       else              v = cone[0];
178:       DMPlexGetSupport(dm, v, &support);
179:       DMPlexGetSupportSize(dm, v, &suppSize);
180:       if (suppSize == 1) {lastCell = -1;}
181:       else {
182:         if (support[0] == lastCell) lastCell = support[1];
183:         else                        lastCell = support[0];
184:       }
185:       points[v] = pos++;
186:     }
188:   }
189:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, points, PETSC_OWN_POINTER, perm);
190:   return 0;
191: }

193: /*@
194:   DMPlexPermute - Reorder the mesh according to the input permutation

196:   Collective on dm

198:   Input Parameters:
199: + dm - The DMPlex object
200: - perm - The point permutation, perm[old point number] = new point number

202:   Output Parameter:
203: . pdm - The permuted DM

205:   Level: intermediate

207: .seealso: MatPermute()
208: @*/
209: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
210: {
211:   DM_Plex       *plex = (DM_Plex *) dm->data, *plexNew;
212:   PetscInt       dim, cdim;
213:   const char    *name;

218:   DMCreate(PetscObjectComm((PetscObject) dm), pdm);
219:   DMSetType(*pdm, DMPLEX);
220:   PetscObjectGetName((PetscObject) dm, &name);
221:   PetscObjectSetName((PetscObject) *pdm, name);
222:   DMGetDimension(dm, &dim);
223:   DMSetDimension(*pdm, dim);
224:   DMGetCoordinateDim(dm, &cdim);
225:   DMSetCoordinateDim(*pdm, cdim);
226:   DMCopyDisc(dm, *pdm);
227:   if (dm->localSection) {
228:     PetscSection section, sectionNew;

230:     DMGetLocalSection(dm, &section);
231:     PetscSectionPermute(section, perm, &sectionNew);
232:     DMSetLocalSection(*pdm, sectionNew);
233:     PetscSectionDestroy(&sectionNew);
234:   }
235:   plexNew = (DM_Plex *) (*pdm)->data;
236:   /* Ignore ltogmap, ltogmapb */
237:   /* Ignore sf, sectionSF */
238:   /* Ignore globalVertexNumbers, globalCellNumbers */
239:   /* Reorder labels */
240:   {
241:     PetscInt numLabels, l;
242:     DMLabel  label, labelNew;

244:     DMGetNumLabels(dm, &numLabels);
245:     for (l = 0; l < numLabels; ++l) {
246:       DMGetLabelByNum(dm, l, &label);
247:       DMLabelPermute(label, perm, &labelNew);
248:       DMAddLabel(*pdm, labelNew);
249:       DMLabelDestroy(&labelNew);
250:     }
251:     DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel);
252:     if (plex->subpointMap) DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);
253:   }
254:   /* Reorder topology */
255:   {
256:     const PetscInt *pperm;
257:     PetscInt        n, pStart, pEnd, p;

259:     PetscSectionDestroy(&plexNew->coneSection);
260:     PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
261:     PetscSectionGetStorageSize(plexNew->coneSection, &n);
262:     PetscMalloc1(n, &plexNew->cones);
263:     PetscMalloc1(n, &plexNew->coneOrientations);
264:     ISGetIndices(perm, &pperm);
265:     PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
266:     for (p = pStart; p < pEnd; ++p) {
267:       PetscInt dof, off, offNew, d;

269:       PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
270:       PetscSectionGetOffset(plex->coneSection, p, &off);
271:       PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
272:       for (d = 0; d < dof; ++d) {
273:         plexNew->cones[offNew+d]            = pperm[plex->cones[off+d]];
274:         plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
275:       }
276:     }
277:     PetscSectionDestroy(&plexNew->supportSection);
278:     PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
279:     PetscSectionGetStorageSize(plexNew->supportSection, &n);
280:     PetscMalloc1(n, &plexNew->supports);
281:     PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
282:     for (p = pStart; p < pEnd; ++p) {
283:       PetscInt dof, off, offNew, d;

285:       PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
286:       PetscSectionGetOffset(plex->supportSection, p, &off);
287:       PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
288:       for (d = 0; d < dof; ++d) {
289:         plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
290:       }
291:     }
292:     ISRestoreIndices(perm, &pperm);
293:   }
294:   /* Remap coordinates */
295:   {
296:     DM              cdm, cdmNew;
297:     PetscSection    csection, csectionNew;
298:     Vec             coordinates, coordinatesNew;
299:     PetscScalar    *coords, *coordsNew;
300:     const PetscInt *pperm;
301:     PetscInt        pStart, pEnd, p;
302:     const char     *name;

304:     DMGetCoordinateDM(dm, &cdm);
305:     DMGetLocalSection(cdm, &csection);
306:     PetscSectionPermute(csection, perm, &csectionNew);
307:     DMGetCoordinatesLocal(dm, &coordinates);
308:     VecDuplicate(coordinates, &coordinatesNew);
309:     PetscObjectGetName((PetscObject)coordinates,&name);
310:     PetscObjectSetName((PetscObject)coordinatesNew,name);
311:     VecGetArray(coordinates, &coords);
312:     VecGetArray(coordinatesNew, &coordsNew);
313:     PetscSectionGetChart(csectionNew, &pStart, &pEnd);
314:     ISGetIndices(perm, &pperm);
315:     for (p = pStart; p < pEnd; ++p) {
316:       PetscInt dof, off, offNew, d;

318:       PetscSectionGetDof(csectionNew, p, &dof);
319:       PetscSectionGetOffset(csection, p, &off);
320:       PetscSectionGetOffset(csectionNew, pperm[p], &offNew);
321:       for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
322:     }
323:     ISRestoreIndices(perm, &pperm);
324:     VecRestoreArray(coordinates, &coords);
325:     VecRestoreArray(coordinatesNew, &coordsNew);
326:     DMGetCoordinateDM(*pdm, &cdmNew);
327:     DMSetLocalSection(cdmNew, csectionNew);
328:     DMSetCoordinatesLocal(*pdm, coordinatesNew);
329:     PetscSectionDestroy(&csectionNew);
330:     VecDestroy(&coordinatesNew);
331:   }
332:   DMPlexCopy_Internal(dm, PETSC_TRUE, *pdm);
333:   (*pdm)->setupcalled = PETSC_TRUE;
334:   return 0;
335: }