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