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 9 PetscFunctionBegin; 10 PetscCall(DMPlexGetDepth(dm, &depth)); 11 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 12 PetscCall(PetscMalloc1(pEnd-pStart,&perm)); 13 PetscCall(PetscMalloc1(pEnd-pStart,&iperm)); 14 for (p = pStart; p < pEnd; ++p) iperm[p] = -1; 15 for (d = depth; d > 0; --d) { 16 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 17 PetscCall(DMPlexGetDepthStratum(dm, d-1, &fStart, &fEnd)); 18 fMax = fStart; 19 for (p = pStart; p < pEnd; ++p) { 20 const PetscInt *cone; 21 PetscInt point, coneSize, c; 22 23 if (d == depth) { 24 perm[p] = pperm[p]; 25 iperm[pperm[p]] = p; 26 } 27 point = perm[p]; 28 PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 29 PetscCall(DMPlexGetCone(dm, point, &cone)); 30 for (c = 0; c < coneSize; ++c) { 31 const PetscInt oldc = cone[c]; 32 const PetscInt newc = iperm[oldc]; 33 34 if (newc < 0) { 35 perm[fMax] = oldc; 36 iperm[oldc] = fMax++; 37 } 38 } 39 } 40 PetscCheck(fMax == fEnd,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %" PetscInt_FMT " faces %" PetscInt_FMT " does not match permuted number %" PetscInt_FMT, d, fEnd-fStart, fMax-fStart); 41 } 42 *clperm = perm; 43 *invclperm = iperm; 44 PetscFunctionReturn(0); 45 } 46 47 /*@ 48 DMPlexGetOrdering - Calculate a reordering of the mesh 49 50 Collective on dm 51 52 Input Parameters: 53 + dm - The DMPlex object 54 . otype - type of reordering, one of the following: 55 $ MATORDERINGNATURAL - Natural 56 $ MATORDERINGND - Nested Dissection 57 $ MATORDERING1WD - One-way Dissection 58 $ MATORDERINGRCM - Reverse Cuthill-McKee 59 $ MATORDERINGQMD - Quotient Minimum Degree 60 - label - [Optional] Label used to segregate ordering into sets, or NULL 61 62 Output Parameter: 63 . perm - The point permutation as an IS, perm[old point number] = new point number 64 65 Note: The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which 66 has different types of cells, and then loop over each set of reordered cells for assembly. 67 68 Level: intermediate 69 70 .seealso: `DMPlexPermute()`, `MatGetOrdering()` 71 @*/ 72 PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm) 73 { 74 PetscInt numCells = 0; 75 PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i; 76 77 PetscFunctionBegin; 78 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 79 PetscValidPointer(perm, 4); 80 PetscCall(DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency)); 81 PetscCall(PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls)); 82 if (numCells) { 83 /* Shift for Fortran numbering */ 84 for (i = 0; i < start[numCells]; ++i) ++adjacency[i]; 85 for (i = 0; i <= numCells; ++i) ++start[i]; 86 PetscCall(SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls)); 87 } 88 PetscCall(PetscFree(start)); 89 PetscCall(PetscFree(adjacency)); 90 /* Shift for Fortran numbering */ 91 for (c = 0; c < numCells; ++c) --cperm[c]; 92 /* Segregate */ 93 if (label) { 94 IS valueIS; 95 const PetscInt *values; 96 PetscInt numValues, numPoints = 0; 97 PetscInt *sperm, *vsize, *voff, v; 98 99 PetscCall(DMLabelGetValueIS(label, &valueIS)); 100 PetscCall(ISSort(valueIS)); 101 PetscCall(ISGetLocalSize(valueIS, &numValues)); 102 PetscCall(ISGetIndices(valueIS, &values)); 103 PetscCall(PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff)); 104 for (v = 0; v < numValues; ++v) { 105 PetscCall(DMLabelGetStratumSize(label, values[v], &vsize[v])); 106 if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1]; 107 numPoints += vsize[v]; 108 } 109 PetscCheck(numPoints == numCells, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %" PetscInt_FMT " cells != %" PetscInt_FMT " total cells", numPoints, numCells); 110 for (c = 0; c < numCells; ++c) { 111 const PetscInt oldc = cperm[c]; 112 PetscInt val, vloc; 113 114 PetscCall(DMLabelGetValue(label, oldc, &val)); 115 PetscCheck(val != -1,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " not present in label", oldc); 116 PetscCall(PetscFindInt(val, numValues, values, &vloc)); 117 PetscCheck(vloc >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " not present label", val); 118 sperm[voff[vloc+1]++] = oldc; 119 } 120 for (v = 0; v < numValues; ++v) { 121 PetscCheck(voff[v+1] - voff[v] == vsize[v],PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of %" PetscInt_FMT " values found is %" PetscInt_FMT " != %" PetscInt_FMT, values[v], voff[v+1] - voff[v], vsize[v]); 122 } 123 PetscCall(ISRestoreIndices(valueIS, &values)); 124 PetscCall(ISDestroy(&valueIS)); 125 PetscCall(PetscArraycpy(cperm, sperm, numCells)); 126 PetscCall(PetscFree3(sperm, vsize, voff)); 127 } 128 /* Construct closure */ 129 PetscCall(DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm)); 130 PetscCall(PetscFree3(cperm,mask,xls)); 131 PetscCall(PetscFree(clperm)); 132 /* Invert permutation */ 133 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 134 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm)); 135 PetscFunctionReturn(0); 136 } 137 138 /*@ 139 DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line 140 141 Collective on dm 142 143 Input Parameter: 144 . dm - The DMPlex object 145 146 Output Parameter: 147 . perm - The point permutation as an IS, perm[old point number] = new point number 148 149 Level: intermediate 150 151 .seealso: `DMPlexGetOrdering()`, `DMPlexPermute()`, `MatGetOrdering()` 152 @*/ 153 PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm) 154 { 155 PetscInt *points; 156 const PetscInt *support, *cone; 157 PetscInt dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0; 158 159 PetscFunctionBegin; 160 PetscCall(DMGetDimension(dm, &dim)); 161 PetscCheck(dim == 1, PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Input mesh must be one dimensional, not %" PetscInt_FMT, dim); 162 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 163 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 164 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 165 PetscCall(PetscMalloc1(pEnd-pStart, &points)); 166 for (c = cStart; c < cEnd; ++c) points[c] = c; 167 for (v = vStart; v < vEnd; ++v) points[v] = v; 168 for (v = vStart; v < vEnd; ++v) { 169 PetscCall(DMPlexGetSupportSize(dm, v, &suppSize)); 170 PetscCall(DMPlexGetSupport(dm, v, &support)); 171 if (suppSize == 1) {lastCell = support[0]; break;} 172 } 173 if (v < vEnd) { 174 PetscInt pos = cEnd; 175 176 points[v] = pos++; 177 while (lastCell >= cStart) { 178 PetscCall(DMPlexGetCone(dm, lastCell, &cone)); 179 if (cone[0] == v) v = cone[1]; 180 else v = cone[0]; 181 PetscCall(DMPlexGetSupport(dm, v, &support)); 182 PetscCall(DMPlexGetSupportSize(dm, v, &suppSize)); 183 if (suppSize == 1) {lastCell = -1;} 184 else { 185 if (support[0] == lastCell) lastCell = support[1]; 186 else lastCell = support[0]; 187 } 188 points[v] = pos++; 189 } 190 PetscCheck(pos == pEnd, PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Last vertex was %" PetscInt_FMT ", not %" PetscInt_FMT, pos, pEnd); 191 } 192 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, points, PETSC_OWN_POINTER, perm)); 193 PetscFunctionReturn(0); 194 } 195 196 /*@ 197 DMPlexPermute - Reorder the mesh according to the input permutation 198 199 Collective on dm 200 201 Input Parameters: 202 + dm - The DMPlex object 203 - perm - The point permutation, perm[old point number] = new point number 204 205 Output Parameter: 206 . pdm - The permuted DM 207 208 Level: intermediate 209 210 .seealso: `MatPermute()` 211 @*/ 212 PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm) 213 { 214 DM_Plex *plex = (DM_Plex *) dm->data, *plexNew; 215 PetscInt dim, cdim; 216 const char *name; 217 218 PetscFunctionBegin; 219 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 220 PetscValidHeaderSpecific(perm, IS_CLASSID, 2); 221 PetscValidPointer(pdm, 3); 222 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), pdm)); 223 PetscCall(DMSetType(*pdm, DMPLEX)); 224 PetscCall(PetscObjectGetName((PetscObject) dm, &name)); 225 PetscCall(PetscObjectSetName((PetscObject) *pdm, name)); 226 PetscCall(DMGetDimension(dm, &dim)); 227 PetscCall(DMSetDimension(*pdm, dim)); 228 PetscCall(DMGetCoordinateDim(dm, &cdim)); 229 PetscCall(DMSetCoordinateDim(*pdm, cdim)); 230 PetscCall(DMCopyDisc(dm, *pdm)); 231 if (dm->localSection) { 232 PetscSection section, sectionNew; 233 234 PetscCall(DMGetLocalSection(dm, §ion)); 235 PetscCall(PetscSectionPermute(section, perm, §ionNew)); 236 PetscCall(DMSetLocalSection(*pdm, sectionNew)); 237 PetscCall(PetscSectionDestroy(§ionNew)); 238 } 239 plexNew = (DM_Plex *) (*pdm)->data; 240 /* Ignore ltogmap, ltogmapb */ 241 /* Ignore sf, sectionSF */ 242 /* Ignore globalVertexNumbers, globalCellNumbers */ 243 /* Reorder labels */ 244 { 245 PetscInt numLabels, l; 246 DMLabel label, labelNew; 247 248 PetscCall(DMGetNumLabels(dm, &numLabels)); 249 for (l = 0; l < numLabels; ++l) { 250 PetscCall(DMGetLabelByNum(dm, l, &label)); 251 PetscCall(DMLabelPermute(label, perm, &labelNew)); 252 PetscCall(DMAddLabel(*pdm, labelNew)); 253 PetscCall(DMLabelDestroy(&labelNew)); 254 } 255 PetscCall(DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel)); 256 if (plex->subpointMap) PetscCall(DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap)); 257 } 258 /* Reorder topology */ 259 { 260 const PetscInt *pperm; 261 PetscInt n, pStart, pEnd, p; 262 263 PetscCall(PetscSectionDestroy(&plexNew->coneSection)); 264 PetscCall(PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection)); 265 PetscCall(PetscSectionGetStorageSize(plexNew->coneSection, &n)); 266 PetscCall(PetscMalloc1(n, &plexNew->cones)); 267 PetscCall(PetscMalloc1(n, &plexNew->coneOrientations)); 268 PetscCall(ISGetIndices(perm, &pperm)); 269 PetscCall(PetscSectionGetChart(plex->coneSection, &pStart, &pEnd)); 270 for (p = pStart; p < pEnd; ++p) { 271 PetscInt dof, off, offNew, d; 272 273 PetscCall(PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof)); 274 PetscCall(PetscSectionGetOffset(plex->coneSection, p, &off)); 275 PetscCall(PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew)); 276 for (d = 0; d < dof; ++d) { 277 plexNew->cones[offNew+d] = pperm[plex->cones[off+d]]; 278 plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d]; 279 } 280 } 281 PetscCall(PetscSectionDestroy(&plexNew->supportSection)); 282 PetscCall(PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection)); 283 PetscCall(PetscSectionGetStorageSize(plexNew->supportSection, &n)); 284 PetscCall(PetscMalloc1(n, &plexNew->supports)); 285 PetscCall(PetscSectionGetChart(plex->supportSection, &pStart, &pEnd)); 286 for (p = pStart; p < pEnd; ++p) { 287 PetscInt dof, off, offNew, d; 288 289 PetscCall(PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof)); 290 PetscCall(PetscSectionGetOffset(plex->supportSection, p, &off)); 291 PetscCall(PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew)); 292 for (d = 0; d < dof; ++d) { 293 plexNew->supports[offNew+d] = pperm[plex->supports[off+d]]; 294 } 295 } 296 PetscCall(ISRestoreIndices(perm, &pperm)); 297 } 298 /* Remap coordinates */ 299 { 300 DM cdm, cdmNew; 301 PetscSection csection, csectionNew; 302 Vec coordinates, coordinatesNew; 303 PetscScalar *coords, *coordsNew; 304 const PetscInt *pperm; 305 PetscInt pStart, pEnd, p; 306 const char *name; 307 308 PetscCall(DMGetCoordinateDM(dm, &cdm)); 309 PetscCall(DMGetLocalSection(cdm, &csection)); 310 PetscCall(PetscSectionPermute(csection, perm, &csectionNew)); 311 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 312 PetscCall(VecDuplicate(coordinates, &coordinatesNew)); 313 PetscCall(PetscObjectGetName((PetscObject)coordinates,&name)); 314 PetscCall(PetscObjectSetName((PetscObject)coordinatesNew,name)); 315 PetscCall(VecGetArray(coordinates, &coords)); 316 PetscCall(VecGetArray(coordinatesNew, &coordsNew)); 317 PetscCall(PetscSectionGetChart(csectionNew, &pStart, &pEnd)); 318 PetscCall(ISGetIndices(perm, &pperm)); 319 for (p = pStart; p < pEnd; ++p) { 320 PetscInt dof, off, offNew, d; 321 322 PetscCall(PetscSectionGetDof(csectionNew, p, &dof)); 323 PetscCall(PetscSectionGetOffset(csection, p, &off)); 324 PetscCall(PetscSectionGetOffset(csectionNew, pperm[p], &offNew)); 325 for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d]; 326 } 327 PetscCall(ISRestoreIndices(perm, &pperm)); 328 PetscCall(VecRestoreArray(coordinates, &coords)); 329 PetscCall(VecRestoreArray(coordinatesNew, &coordsNew)); 330 PetscCall(DMGetCoordinateDM(*pdm, &cdmNew)); 331 PetscCall(DMSetLocalSection(cdmNew, csectionNew)); 332 PetscCall(DMSetCoordinatesLocal(*pdm, coordinatesNew)); 333 PetscCall(PetscSectionDestroy(&csectionNew)); 334 PetscCall(VecDestroy(&coordinatesNew)); 335 } 336 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *pdm)); 337 (*pdm)->setupcalled = PETSC_TRUE; 338 PetscFunctionReturn(0); 339 } 340