Actual source code: plexreorder.c

petsc-3.10.5 2019-03-28
Report Typos and Errors
  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;

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

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

 35:         if (newc < 0) {
 36:           perm[fMax]  = oldc;
 37:           iperm[oldc] = fMax++;
 38:         }
 39:       }
 40:     }
 41:     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);
 42:   }
 43:   *clperm    = perm;
 44:   *invclperm = iperm;
 45:   return(0);
 46: }

 48: /*@
 49:   DMPlexGetOrdering - Calculate a reordering of the mesh

 51:   Collective on DM

 53:   Input Parameter:
 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


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

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

 70:   Level: intermediate

 72: .keywords: mesh
 73: .seealso: MatGetOrdering()
 74: @*/
 75: PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
 76: {
 77:   PetscInt       numCells = 0;
 78:   PetscInt      *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *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:   /* Segregate */
 97:   if (label) {
 98:     IS              valueIS;
 99:     const PetscInt *values;
100:     PetscInt        numValues, numPoints = 0;
101:     PetscInt       *sperm, *vsize, *voff, v;

103:     DMLabelGetValueIS(label, &valueIS);
104:     ISSort(valueIS);
105:     ISGetLocalSize(valueIS, &numValues);
106:     ISGetIndices(valueIS, &values);
107:     PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff);
108:     for (v = 0; v < numValues; ++v) {
109:       DMLabelGetStratumSize(label, values[v], &vsize[v]);
110:       if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1];
111:       numPoints += vsize[v];
112:     }
113:     if (numPoints != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %D cells < %D total", numPoints, numCells);
114:     for (c = 0; c < numCells; ++c) {
115:       const PetscInt oldc = cperm[c];
116:       PetscInt       val, vloc;

118:       DMLabelGetValue(label, oldc, &val);
119:       if (val == -1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %D not present in label", oldc);
120:       PetscFindInt(val, numValues, values, &vloc);
121:       if (vloc < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %D not present label", val);
122:       sperm[voff[vloc+1]++] = oldc;
123:     }
124:     for (v = 0; v < numValues; ++v) {
125:       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]);
126:     }
127:     ISRestoreIndices(valueIS, &values);
128:     ISDestroy(&valueIS);
129:     PetscMemcpy(cperm, sperm, numCells * sizeof(PetscInt));
130:     PetscFree3(sperm, vsize, voff);
131:   }
132:   /* Construct closure */
133:   DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
134:   PetscFree3(cperm,mask,xls);
135:   PetscFree(clperm);
136:   /* Invert permutation */
137:   DMPlexGetChart(dm, &pStart, &pEnd);
138:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);
139:   return(0);
140: }

142: /*@
143:   DMPlexPermute - Reorder the mesh according to the input permutation

145:   Collective on DM

147:   Input Parameter:
148: + dm - The DMPlex object
149: - perm - The point permutation, perm[old point number] = new point number

151:   Output Parameter:
152: . pdm - The permuted DM

154:   Level: intermediate

156: .keywords: mesh
157: .seealso: MatPermute()
158: @*/
159: PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
160: {
161:   DM_Plex       *plex = (DM_Plex *) dm->data, *plexNew;
162:   PetscSection   section, sectionNew;
163:   PetscInt       dim;

170:   DMCreate(PetscObjectComm((PetscObject) dm), pdm);
171:   DMSetType(*pdm, DMPLEX);
172:   DMGetDimension(dm, &dim);
173:   DMSetDimension(*pdm, dim);
174:   DMGetSection(dm, &section);
175:   if (section) {
176:     PetscSectionPermute(section, perm, &sectionNew);
177:     DMSetSection(*pdm, sectionNew);
178:     PetscSectionDestroy(&sectionNew);
179:   }
180:   plexNew = (DM_Plex *) (*pdm)->data;
181:   /* Ignore ltogmap, ltogmapb */
182:   /* Ignore sf, defaultSF */
183:   /* Ignore globalVertexNumbers, globalCellNumbers */
184:   /* Remap coordinates */
185:   {
186:     DM              cdm, cdmNew;
187:     PetscSection    csection, csectionNew;
188:     Vec             coordinates, coordinatesNew;
189:     PetscScalar    *coords, *coordsNew;
190:     const PetscInt *pperm;
191:     PetscInt        pStart, pEnd, p;
192:     const char     *name;

194:     DMGetCoordinateDM(dm, &cdm);
195:     DMGetSection(cdm, &csection);
196:     PetscSectionPermute(csection, perm, &csectionNew);
197:     DMGetCoordinatesLocal(dm, &coordinates);
198:     VecDuplicate(coordinates, &coordinatesNew);
199:     PetscObjectGetName((PetscObject)coordinates,&name);
200:     PetscObjectSetName((PetscObject)coordinatesNew,name);
201:     VecGetArray(coordinates, &coords);
202:     VecGetArray(coordinatesNew, &coordsNew);
203:     PetscSectionGetChart(csectionNew, &pStart, &pEnd);
204:     ISGetIndices(perm, &pperm);
205:     for (p = pStart; p < pEnd; ++p) {
206:       PetscInt dof, off, offNew, d;

208:       PetscSectionGetDof(csectionNew, p, &dof);
209:       PetscSectionGetOffset(csection, p, &off);
210:       PetscSectionGetOffset(csectionNew, pperm[p], &offNew);
211:       for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
212:     }
213:     ISRestoreIndices(perm, &pperm);
214:     VecRestoreArray(coordinates, &coords);
215:     VecRestoreArray(coordinatesNew, &coordsNew);
216:     DMGetCoordinateDM(*pdm, &cdmNew);
217:     DMSetSection(cdmNew, csectionNew);
218:     DMSetCoordinatesLocal(*pdm, coordinatesNew);
219:     PetscSectionDestroy(&csectionNew);
220:     VecDestroy(&coordinatesNew);
221:   }
222:   /* Reorder labels */
223:   {
224:     PetscInt numLabels, l;
225:     DMLabel  label, labelNew;

227:     DMGetNumLabels(dm, &numLabels);
228:     for (l = numLabels-1; l >= 0; --l) {
229:       DMGetLabelByNum(dm, l, &label);
230:       DMLabelPermute(label, perm, &labelNew);
231:       DMAddLabel(*pdm, labelNew);
232:     }
233:     if (plex->subpointMap) {DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);}
234:   }
235:   /* Reorder topology */
236:   {
237:     const PetscInt *pperm;
238:     PetscInt        maxConeSize, maxSupportSize, n, pStart, pEnd, p;

240:     DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
241:     plexNew->maxConeSize    = maxConeSize;
242:     plexNew->maxSupportSize = maxSupportSize;
243:     PetscSectionDestroy(&plexNew->coneSection);
244:     PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);
245:     PetscSectionGetStorageSize(plexNew->coneSection, &n);
246:     PetscMalloc1(n, &plexNew->cones);
247:     PetscMalloc1(n, &plexNew->coneOrientations);
248:     ISGetIndices(perm, &pperm);
249:     PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);
250:     for (p = pStart; p < pEnd; ++p) {
251:       PetscInt dof, off, offNew, d;

253:       PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);
254:       PetscSectionGetOffset(plex->coneSection, p, &off);
255:       PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);
256:       for (d = 0; d < dof; ++d) {
257:         plexNew->cones[offNew+d]            = pperm[plex->cones[off+d]];
258:         plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
259:       }
260:     }
261:     PetscSectionDestroy(&plexNew->supportSection);
262:     PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);
263:     PetscSectionGetStorageSize(plexNew->supportSection, &n);
264:     PetscMalloc1(n, &plexNew->supports);
265:     PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);
266:     for (p = pStart; p < pEnd; ++p) {
267:       PetscInt dof, off, offNew, d;

269:       PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);
270:       PetscSectionGetOffset(plex->supportSection, p, &off);
271:       PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);
272:       for (d = 0; d < dof; ++d) {
273:         plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
274:       }
275:     }
276:     ISRestoreIndices(perm, &pperm);
277:   }
278:   return(0);
279: }