1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/matorderimpl.h> /*I "petscmat.h" I*/ 3 #include <petsc/private/dmlabelimpl.h> 4 5 static PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm) 6 { 7 PetscInt *perm, *iperm; 8 PetscInt depth, d, pStart, pEnd, fStart, fMax, fEnd, p; 9 10 PetscFunctionBegin; 11 PetscCall(DMPlexGetDepth(dm, &depth)); 12 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 13 PetscCall(PetscMalloc1(pEnd-pStart,&perm)); 14 PetscCall(PetscMalloc1(pEnd-pStart,&iperm)); 15 for (p = pStart; p < pEnd; ++p) iperm[p] = -1; 16 for (d = depth; d > 0; --d) { 17 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 18 PetscCall(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; 23 24 if (d == depth) { 25 perm[p] = pperm[p]; 26 iperm[pperm[p]] = p; 27 } 28 point = perm[p]; 29 PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 30 PetscCall(DMPlexGetCone(dm, point, &cone)); 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 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); 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 78 PetscFunctionBegin; 79 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 80 PetscValidPointer(perm, 4); 81 PetscCall(DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency)); 82 PetscCall(PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls)); 83 if (numCells) { 84 /* Shift for Fortran numbering */ 85 for (i = 0; i < start[numCells]; ++i) ++adjacency[i]; 86 for (i = 0; i <= numCells; ++i) ++start[i]; 87 PetscCall(SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls)); 88 } 89 PetscCall(PetscFree(start)); 90 PetscCall(PetscFree(adjacency)); 91 /* Shift for Fortran numbering */ 92 for (c = 0; c < numCells; ++c) --cperm[c]; 93 /* Segregate */ 94 if (label) { 95 IS valueIS; 96 const PetscInt *valuesTmp; 97 PetscInt *values; 98 PetscInt numValues, numPoints = 0; 99 PetscInt *sperm, *vsize, *voff, v; 100 101 // Can't directly sort the valueIS, since it is a view into the DMLabel 102 PetscCall(DMLabelGetValueIS(label, &valueIS)); 103 PetscCall(ISGetLocalSize(valueIS, &numValues)); 104 PetscCall(ISGetIndices(valueIS, &valuesTmp)); 105 PetscCall(PetscCalloc4(numCells, &sperm, numValues, &values, numValues, &vsize, numValues+1,&voff)); 106 PetscCall(PetscArraycpy(values, valuesTmp, numValues)); 107 PetscCall(PetscSortInt(numValues, values)); 108 PetscCall(ISRestoreIndices(valueIS, &valuesTmp)); 109 PetscCall(ISDestroy(&valueIS)); 110 for (v = 0; v < numValues; ++v) { 111 PetscCall(DMLabelGetStratumSize(label, values[v], &vsize[v])); 112 if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1]; 113 numPoints += vsize[v]; 114 } 115 PetscCheck(numPoints == numCells, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %" PetscInt_FMT " cells != %" PetscInt_FMT " total cells", numPoints, numCells); 116 for (c = 0; c < numCells; ++c) { 117 const PetscInt oldc = cperm[c]; 118 PetscInt val, vloc; 119 120 PetscCall(DMLabelGetValue(label, oldc, &val)); 121 PetscCheck(val != -1,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " not present in label", oldc); 122 PetscCall(PetscFindInt(val, numValues, values, &vloc)); 123 PetscCheck(vloc >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " not present label", val); 124 sperm[voff[vloc+1]++] = oldc; 125 } 126 for (v = 0; v < numValues; ++v) { 127 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]); 128 } 129 PetscCall(PetscArraycpy(cperm, sperm, numCells)); 130 PetscCall(PetscFree4(sperm, values, vsize, voff)); 131 } 132 /* Construct closure */ 133 PetscCall(DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm)); 134 PetscCall(PetscFree3(cperm,mask,xls)); 135 PetscCall(PetscFree(clperm)); 136 /* Invert permutation */ 137 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 138 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm)); 139 PetscFunctionReturn(0); 140 } 141 142 /*@ 143 DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line 144 145 Collective on dm 146 147 Input Parameter: 148 . dm - The DMPlex object 149 150 Output Parameter: 151 . perm - The point permutation as an IS, perm[old point number] = new point number 152 153 Level: intermediate 154 155 .seealso: `DMPlexGetOrdering()`, `DMPlexPermute()`, `MatGetOrdering()` 156 @*/ 157 PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm) 158 { 159 PetscInt *points; 160 const PetscInt *support, *cone; 161 PetscInt dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0; 162 163 PetscFunctionBegin; 164 PetscCall(DMGetDimension(dm, &dim)); 165 PetscCheck(dim == 1, PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Input mesh must be one dimensional, not %" PetscInt_FMT, dim); 166 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 167 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 168 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 169 PetscCall(PetscMalloc1(pEnd-pStart, &points)); 170 for (c = cStart; c < cEnd; ++c) points[c] = c; 171 for (v = vStart; v < vEnd; ++v) points[v] = v; 172 for (v = vStart; v < vEnd; ++v) { 173 PetscCall(DMPlexGetSupportSize(dm, v, &suppSize)); 174 PetscCall(DMPlexGetSupport(dm, v, &support)); 175 if (suppSize == 1) {lastCell = support[0]; break;} 176 } 177 if (v < vEnd) { 178 PetscInt pos = cEnd; 179 180 points[v] = pos++; 181 while (lastCell >= cStart) { 182 PetscCall(DMPlexGetCone(dm, lastCell, &cone)); 183 if (cone[0] == v) v = cone[1]; 184 else v = cone[0]; 185 PetscCall(DMPlexGetSupport(dm, v, &support)); 186 PetscCall(DMPlexGetSupportSize(dm, v, &suppSize)); 187 if (suppSize == 1) {lastCell = -1;} 188 else { 189 if (support[0] == lastCell) lastCell = support[1]; 190 else lastCell = support[0]; 191 } 192 points[v] = pos++; 193 } 194 PetscCheck(pos == pEnd, PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Last vertex was %" PetscInt_FMT ", not %" PetscInt_FMT, pos, pEnd); 195 } 196 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, points, PETSC_OWN_POINTER, perm)); 197 PetscFunctionReturn(0); 198 } 199 200 static PetscErrorCode DMPlexRemapCoordinates_Private(IS perm, PetscSection cs, Vec coordinates, PetscSection *csNew, Vec *coordinatesNew) 201 { 202 PetscScalar *coords, *coordsNew; 203 const PetscInt *pperm; 204 PetscInt pStart, pEnd, p; 205 const char *name; 206 207 PetscFunctionBegin; 208 PetscCall(PetscSectionPermute(cs, perm, csNew)); 209 PetscCall(VecDuplicate(coordinates, coordinatesNew)); 210 PetscCall(PetscObjectGetName((PetscObject) coordinates, &name)); 211 PetscCall(PetscObjectSetName((PetscObject) *coordinatesNew, name)); 212 PetscCall(VecGetArray(coordinates, &coords)); 213 PetscCall(VecGetArray(*coordinatesNew, &coordsNew)); 214 PetscCall(PetscSectionGetChart(*csNew, &pStart, &pEnd)); 215 PetscCall(ISGetIndices(perm, &pperm)); 216 for (p = pStart; p < pEnd; ++p) { 217 PetscInt dof, off, offNew, d; 218 219 PetscCall(PetscSectionGetDof(*csNew, p, &dof)); 220 PetscCall(PetscSectionGetOffset(cs, p, &off)); 221 PetscCall(PetscSectionGetOffset(*csNew, pperm[p], &offNew)); 222 for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d]; 223 } 224 PetscCall(ISRestoreIndices(perm, &pperm)); 225 PetscCall(VecRestoreArray(coordinates, &coords)); 226 PetscCall(VecRestoreArray(*coordinatesNew, &coordsNew)); 227 PetscFunctionReturn(0); 228 } 229 230 /*@ 231 DMPlexPermute - Reorder the mesh according to the input permutation 232 233 Collective on dm 234 235 Input Parameters: 236 + dm - The DMPlex object 237 - perm - The point permutation, perm[old point number] = new point number 238 239 Output Parameter: 240 . pdm - The permuted DM 241 242 Level: intermediate 243 244 .seealso: `MatPermute()` 245 @*/ 246 PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm) 247 { 248 DM_Plex *plex = (DM_Plex *) dm->data, *plexNew; 249 PetscInt dim, cdim; 250 const char *name; 251 252 PetscFunctionBegin; 253 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 254 PetscValidHeaderSpecific(perm, IS_CLASSID, 2); 255 PetscValidPointer(pdm, 3); 256 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), pdm)); 257 PetscCall(DMSetType(*pdm, DMPLEX)); 258 PetscCall(PetscObjectGetName((PetscObject) dm, &name)); 259 PetscCall(PetscObjectSetName((PetscObject) *pdm, name)); 260 PetscCall(DMGetDimension(dm, &dim)); 261 PetscCall(DMSetDimension(*pdm, dim)); 262 PetscCall(DMGetCoordinateDim(dm, &cdim)); 263 PetscCall(DMSetCoordinateDim(*pdm, cdim)); 264 PetscCall(DMCopyDisc(dm, *pdm)); 265 if (dm->localSection) { 266 PetscSection section, sectionNew; 267 268 PetscCall(DMGetLocalSection(dm, §ion)); 269 PetscCall(PetscSectionPermute(section, perm, §ionNew)); 270 PetscCall(DMSetLocalSection(*pdm, sectionNew)); 271 PetscCall(PetscSectionDestroy(§ionNew)); 272 } 273 plexNew = (DM_Plex *) (*pdm)->data; 274 /* Ignore ltogmap, ltogmapb */ 275 /* Ignore sf, sectionSF */ 276 /* Ignore globalVertexNumbers, globalCellNumbers */ 277 /* Reorder labels */ 278 { 279 PetscInt numLabels, l; 280 DMLabel label, labelNew; 281 282 PetscCall(DMGetNumLabels(dm, &numLabels)); 283 for (l = 0; l < numLabels; ++l) { 284 PetscCall(DMGetLabelByNum(dm, l, &label)); 285 PetscCall(DMLabelPermute(label, perm, &labelNew)); 286 PetscCall(DMAddLabel(*pdm, labelNew)); 287 PetscCall(DMLabelDestroy(&labelNew)); 288 } 289 PetscCall(DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel)); 290 if (plex->subpointMap) PetscCall(DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap)); 291 } 292 if ((*pdm)->celltypeLabel) { 293 DMLabel ctLabel; 294 295 // Reset label for fast lookup 296 PetscCall(DMPlexGetCellTypeLabel(*pdm, &ctLabel)); 297 PetscCall(DMLabelMakeAllInvalid_Internal(ctLabel)); 298 } 299 /* Reorder topology */ 300 { 301 const PetscInt *pperm; 302 PetscInt n, pStart, pEnd, p; 303 304 PetscCall(PetscSectionDestroy(&plexNew->coneSection)); 305 PetscCall(PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection)); 306 PetscCall(PetscSectionGetStorageSize(plexNew->coneSection, &n)); 307 PetscCall(PetscMalloc1(n, &plexNew->cones)); 308 PetscCall(PetscMalloc1(n, &plexNew->coneOrientations)); 309 PetscCall(ISGetIndices(perm, &pperm)); 310 PetscCall(PetscSectionGetChart(plex->coneSection, &pStart, &pEnd)); 311 for (p = pStart; p < pEnd; ++p) { 312 PetscInt dof, off, offNew, d; 313 314 PetscCall(PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof)); 315 PetscCall(PetscSectionGetOffset(plex->coneSection, p, &off)); 316 PetscCall(PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew)); 317 for (d = 0; d < dof; ++d) { 318 plexNew->cones[offNew+d] = pperm[plex->cones[off+d]]; 319 plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d]; 320 } 321 } 322 PetscCall(PetscSectionDestroy(&plexNew->supportSection)); 323 PetscCall(PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection)); 324 PetscCall(PetscSectionGetStorageSize(plexNew->supportSection, &n)); 325 PetscCall(PetscMalloc1(n, &plexNew->supports)); 326 PetscCall(PetscSectionGetChart(plex->supportSection, &pStart, &pEnd)); 327 for (p = pStart; p < pEnd; ++p) { 328 PetscInt dof, off, offNew, d; 329 330 PetscCall(PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof)); 331 PetscCall(PetscSectionGetOffset(plex->supportSection, p, &off)); 332 PetscCall(PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew)); 333 for (d = 0; d < dof; ++d) { 334 plexNew->supports[offNew+d] = pperm[plex->supports[off+d]]; 335 } 336 } 337 PetscCall(ISRestoreIndices(perm, &pperm)); 338 } 339 /* Remap coordinates */ 340 { 341 DM cdm, cdmNew; 342 PetscSection cs, csNew; 343 Vec coordinates, coordinatesNew; 344 345 PetscCall(DMGetCoordinateSection(dm, &cs)); 346 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 347 PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew)); 348 PetscCall(DMSetCoordinateSection(*pdm, PETSC_DETERMINE, csNew)); 349 PetscCall(DMSetCoordinatesLocal(*pdm, coordinatesNew)); 350 PetscCall(PetscSectionDestroy(&csNew)); 351 PetscCall(VecDestroy(&coordinatesNew)); 352 353 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 354 if (cdm) { 355 PetscCall(DMGetCoordinateDM(*pdm, &cdm)); 356 PetscCall(DMClone(cdm, &cdmNew)); 357 PetscCall(DMSetCellCoordinateDM(*pdm, cdmNew)); 358 PetscCall(DMDestroy(&cdmNew)); 359 PetscCall(DMGetCellCoordinateSection(dm, &cs)); 360 PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates)); 361 PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew)); 362 PetscCall(DMSetCellCoordinateSection(*pdm, PETSC_DETERMINE, csNew)); 363 PetscCall(DMSetCellCoordinatesLocal(*pdm, coordinatesNew)); 364 PetscCall(PetscSectionDestroy(&csNew)); 365 PetscCall(VecDestroy(&coordinatesNew)); 366 } 367 } 368 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *pdm)); 369 (*pdm)->setupcalled = PETSC_TRUE; 370 PetscFunctionReturn(0); 371 } 372 373 PetscErrorCode DMPlexReorderSetDefault_Plex(DM dm, DMPlexReorderDefaultFlag reorder) 374 { 375 DM_Plex *mesh = (DM_Plex *) dm->data; 376 377 PetscFunctionBegin; 378 mesh->reorderDefault = reorder; 379 PetscFunctionReturn(0); 380 } 381 382 /*@ 383 DMPlexReorderSetDefault - Set flag indicating whether the DM should be reordered by default 384 385 Logically collective 386 387 Input Parameters: 388 + dm - The DM 389 - reorder - Flag for reordering 390 391 Level: intermediate 392 393 .seealso: `DMPlexReorderGetDefault()` 394 @*/ 395 PetscErrorCode DMPlexReorderSetDefault(DM dm, DMPlexReorderDefaultFlag reorder) 396 { 397 PetscFunctionBegin; 398 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 399 PetscTryMethod(dm,"DMPlexReorderSetDefault_C",(DM,DMPlexReorderDefaultFlag),(dm,reorder)); 400 PetscFunctionReturn(0); 401 } 402 403 PetscErrorCode DMPlexReorderGetDefault_Plex(DM dm, DMPlexReorderDefaultFlag *reorder) 404 { 405 DM_Plex *mesh = (DM_Plex *) dm->data; 406 407 PetscFunctionBegin; 408 *reorder = mesh->reorderDefault; 409 PetscFunctionReturn(0); 410 } 411 412 /*@ 413 DMPlexReorderGetDefault - Get flag indicating whether the DM should be reordered by default 414 415 Not collective 416 417 Input Parameter: 418 . dm - The DM 419 420 Output Parameter: 421 . reorder - Flag for reordering 422 423 Level: intermediate 424 425 .seealso: `DMPlexReorderSetDefault()` 426 @*/ 427 PetscErrorCode DMPlexReorderGetDefault(DM dm, DMPlexReorderDefaultFlag *reorder) 428 { 429 PetscFunctionBegin; 430 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 431 PetscValidPointer(reorder, 2); 432 PetscUseMethod(dm,"DMPlexReorderGetDefault_C",(DM,DMPlexReorderDefaultFlag*),(dm,reorder)); 433 PetscFunctionReturn(0); 434 } 435