1 #include <petscdm.h> 2 #include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/ 3 #include <petsc/private/sectionimpl.h> /*I "petscsection.h" I*/ 4 #include <petscsf.h> 5 #include <petscsection.h> 6 7 PetscFunctionList DMLabelList = NULL; 8 PetscBool DMLabelRegisterAllCalled = PETSC_FALSE; 9 10 /*@ 11 DMLabelCreate - Create a `DMLabel` object, which is a multimap 12 13 Collective 14 15 Input Parameters: 16 + comm - The communicator, usually `PETSC_COMM_SELF` 17 - name - The label name 18 19 Output Parameter: 20 . label - The `DMLabel` 21 22 Level: beginner 23 24 Notes: 25 The label name is actually usually the `PetscObject` name. 26 One can get/set it with `PetscObjectGetName()`/`PetscObjectSetName()`. 27 28 .seealso: `DMLabel`, `DM`, `DMLabelDestroy()` 29 @*/ 30 PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label) 31 { 32 PetscFunctionBegin; 33 PetscAssertPointer(label, 3); 34 PetscCall(DMInitializePackage()); 35 36 PetscCall(PetscHeaderCreate(*label, DMLABEL_CLASSID, "DMLabel", "DMLabel", "DM", comm, DMLabelDestroy, DMLabelView)); 37 (*label)->numStrata = 0; 38 (*label)->defaultValue = -1; 39 (*label)->stratumValues = NULL; 40 (*label)->validIS = NULL; 41 (*label)->stratumSizes = NULL; 42 (*label)->points = NULL; 43 (*label)->ht = NULL; 44 (*label)->pStart = -1; 45 (*label)->pEnd = -1; 46 (*label)->bt = NULL; 47 PetscCall(PetscHMapICreate(&(*label)->hmap)); 48 PetscCall(PetscObjectSetName((PetscObject)*label, name)); 49 PetscCall(DMLabelSetType(*label, DMLABELCONCRETE)); 50 PetscFunctionReturn(PETSC_SUCCESS); 51 } 52 53 /*@ 54 DMLabelSetUp - SetUp a `DMLabel` object 55 56 Collective 57 58 Input Parameters: 59 . label - The `DMLabel` 60 61 Level: intermediate 62 63 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()` 64 @*/ 65 PetscErrorCode DMLabelSetUp(DMLabel label) 66 { 67 PetscFunctionBegin; 68 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 69 PetscTryTypeMethod(label, setup); 70 PetscFunctionReturn(PETSC_SUCCESS); 71 } 72 73 /* 74 DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format 75 76 Not collective 77 78 Input parameter: 79 + label - The `DMLabel` 80 - v - The stratum value 81 82 Output parameter: 83 . label - The `DMLabel` with stratum in sorted list format 84 85 Level: developer 86 87 .seealso: `DMLabel`, `DM`, `DMLabelCreate()` 88 */ 89 static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v) 90 { 91 IS is; 92 PetscInt off = 0, *pointArray, p; 93 94 PetscFunctionBegin; 95 if ((PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS); 96 PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeValid_Private", v); 97 PetscCall(PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v])); 98 PetscCall(PetscMalloc1(label->stratumSizes[v], &pointArray)); 99 PetscCall(PetscHSetIGetElems(label->ht[v], &off, pointArray)); 100 PetscCall(PetscHSetIClear(label->ht[v])); 101 PetscCall(PetscSortInt(label->stratumSizes[v], pointArray)); 102 if (label->bt) { 103 for (p = 0; p < label->stratumSizes[v]; ++p) { 104 const PetscInt point = pointArray[p]; 105 PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd); 106 PetscCall(PetscBTSet(label->bt, point - label->pStart)); 107 } 108 } 109 if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v] - 1] == pointArray[0] + label->stratumSizes[v] - 1) { 110 PetscCall(ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is)); 111 PetscCall(PetscFree(pointArray)); 112 } else { 113 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is)); 114 } 115 PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, PETSC_TRUE)); 116 PetscCall(PetscObjectSetName((PetscObject)is, "indices")); 117 label->points[v] = is; 118 label->validIS[v] = PETSC_TRUE; 119 PetscCall(PetscObjectStateIncrease((PetscObject)label)); 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 /* 124 DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format 125 126 Not Collective 127 128 Input parameter: 129 . label - The `DMLabel` 130 131 Output parameter: 132 . label - The `DMLabel` with all strata in sorted list format 133 134 Level: developer 135 136 .seealso: `DMLabel`, `DM`, `DMLabelCreate()` 137 */ 138 static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label) 139 { 140 PetscInt v; 141 142 PetscFunctionBegin; 143 for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeValid_Private(label, v)); 144 PetscFunctionReturn(PETSC_SUCCESS); 145 } 146 147 /* 148 DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format 149 150 Not Collective 151 152 Input parameter: 153 + label - The `DMLabel` 154 - v - The stratum value 155 156 Output parameter: 157 . label - The `DMLabel` with stratum in hash format 158 159 Level: developer 160 161 .seealso: `DMLabel`, `DM`, `DMLabelCreate()` 162 */ 163 static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v) 164 { 165 PetscInt p; 166 const PetscInt *points; 167 168 PetscFunctionBegin; 169 if ((PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS); 170 PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeInvalid_Private", v); 171 if (label->points[v]) { 172 PetscCall(ISGetIndices(label->points[v], &points)); 173 for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p])); 174 PetscCall(ISRestoreIndices(label->points[v], &points)); 175 PetscCall(ISDestroy(&label->points[v])); 176 } 177 label->validIS[v] = PETSC_FALSE; 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 PetscErrorCode DMLabelMakeAllInvalid_Internal(DMLabel label) 182 { 183 PetscInt v; 184 185 PetscFunctionBegin; 186 for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeInvalid_Private(label, v)); 187 PetscFunctionReturn(PETSC_SUCCESS); 188 } 189 190 #if !defined(DMLABEL_LOOKUP_THRESHOLD) 191 #define DMLABEL_LOOKUP_THRESHOLD 16 192 #endif 193 194 PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index) 195 { 196 PetscInt v; 197 198 PetscFunctionBegin; 199 *index = -1; 200 if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD || label->readonly) { 201 for (v = 0; v < label->numStrata; ++v) 202 if (label->stratumValues[v] == value) { 203 *index = v; 204 break; 205 } 206 } else { 207 PetscCall(PetscHMapIGet(label->hmap, value, index)); 208 } 209 if (PetscDefined(USE_DEBUG) && !label->readonly) { /* Check strata hash map consistency */ 210 PetscInt len, loc = -1; 211 PetscCall(PetscHMapIGetSize(label->hmap, &len)); 212 PetscCheck(len == label->numStrata, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size"); 213 if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) { 214 PetscCall(PetscHMapIGet(label->hmap, value, &loc)); 215 } else { 216 for (v = 0; v < label->numStrata; ++v) 217 if (label->stratumValues[v] == value) { 218 loc = v; 219 break; 220 } 221 } 222 PetscCheck(loc == *index, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup"); 223 } 224 PetscFunctionReturn(PETSC_SUCCESS); 225 } 226 227 static inline PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index) 228 { 229 PetscInt v; 230 PetscInt *tmpV; 231 PetscInt *tmpS; 232 PetscHSetI *tmpH, ht; 233 IS *tmpP, is; 234 PetscBool *tmpB; 235 PetscHMapI hmap = label->hmap; 236 237 PetscFunctionBegin; 238 v = label->numStrata; 239 tmpV = label->stratumValues; 240 tmpS = label->stratumSizes; 241 tmpH = label->ht; 242 tmpP = label->points; 243 tmpB = label->validIS; 244 { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */ 245 PetscInt *oldV = tmpV; 246 PetscInt *oldS = tmpS; 247 PetscHSetI *oldH = tmpH; 248 IS *oldP = tmpP; 249 PetscBool *oldB = tmpB; 250 PetscCall(PetscMalloc((v + 1) * sizeof(*tmpV), &tmpV)); 251 PetscCall(PetscMalloc((v + 1) * sizeof(*tmpS), &tmpS)); 252 PetscCall(PetscCalloc((v + 1) * sizeof(*tmpH), &tmpH)); 253 PetscCall(PetscCalloc((v + 1) * sizeof(*tmpP), &tmpP)); 254 PetscCall(PetscMalloc((v + 1) * sizeof(*tmpB), &tmpB)); 255 PetscCall(PetscArraycpy(tmpV, oldV, v)); 256 PetscCall(PetscArraycpy(tmpS, oldS, v)); 257 PetscCall(PetscArraycpy(tmpH, oldH, v)); 258 PetscCall(PetscArraycpy(tmpP, oldP, v)); 259 PetscCall(PetscArraycpy(tmpB, oldB, v)); 260 PetscCall(PetscFree(oldV)); 261 PetscCall(PetscFree(oldS)); 262 PetscCall(PetscFree(oldH)); 263 PetscCall(PetscFree(oldP)); 264 PetscCall(PetscFree(oldB)); 265 } 266 label->numStrata = v + 1; 267 label->stratumValues = tmpV; 268 label->stratumSizes = tmpS; 269 label->ht = tmpH; 270 label->points = tmpP; 271 label->validIS = tmpB; 272 PetscCall(PetscHSetICreate(&ht)); 273 PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is)); 274 PetscCall(PetscHMapISet(hmap, value, v)); 275 tmpV[v] = value; 276 tmpS[v] = 0; 277 tmpH[v] = ht; 278 tmpP[v] = is; 279 tmpB[v] = PETSC_TRUE; 280 PetscCall(PetscObjectStateIncrease((PetscObject)label)); 281 *index = v; 282 PetscFunctionReturn(PETSC_SUCCESS); 283 } 284 285 static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index) 286 { 287 PetscFunctionBegin; 288 PetscCall(DMLabelLookupStratum(label, value, index)); 289 if (*index < 0) PetscCall(DMLabelNewStratum(label, value, index)); 290 PetscFunctionReturn(PETSC_SUCCESS); 291 } 292 293 PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size) 294 { 295 PetscFunctionBegin; 296 *size = 0; 297 if (v < 0) PetscFunctionReturn(PETSC_SUCCESS); 298 if (label->readonly || label->validIS[v]) { 299 *size = label->stratumSizes[v]; 300 } else { 301 PetscCall(PetscHSetIGetSize(label->ht[v], size)); 302 } 303 PetscFunctionReturn(PETSC_SUCCESS); 304 } 305 306 /*@ 307 DMLabelAddStratum - Adds a new stratum value in a `DMLabel` 308 309 Input Parameters: 310 + label - The `DMLabel` 311 - value - The stratum value 312 313 Level: beginner 314 315 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()` 316 @*/ 317 PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value) 318 { 319 PetscInt v; 320 321 PetscFunctionBegin; 322 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 323 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 324 PetscCall(DMLabelLookupAddStratum(label, value, &v)); 325 PetscFunctionReturn(PETSC_SUCCESS); 326 } 327 328 /*@ 329 DMLabelAddStrata - Adds new stratum values in a `DMLabel` 330 331 Not Collective 332 333 Input Parameters: 334 + label - The `DMLabel` 335 . numStrata - The number of stratum values 336 - stratumValues - The stratum values 337 338 Level: beginner 339 340 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()` 341 @*/ 342 PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[]) 343 { 344 PetscInt *values, v; 345 346 PetscFunctionBegin; 347 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 348 if (numStrata) PetscAssertPointer(stratumValues, 3); 349 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 350 PetscCall(PetscMalloc1(numStrata, &values)); 351 PetscCall(PetscArraycpy(values, stratumValues, numStrata)); 352 PetscCall(PetscSortRemoveDupsInt(&numStrata, values)); 353 if (!label->numStrata) { /* Fast preallocation */ 354 PetscInt *tmpV; 355 PetscInt *tmpS; 356 PetscHSetI *tmpH, ht; 357 IS *tmpP, is; 358 PetscBool *tmpB; 359 PetscHMapI hmap = label->hmap; 360 361 PetscCall(PetscMalloc1(numStrata, &tmpV)); 362 PetscCall(PetscMalloc1(numStrata, &tmpS)); 363 PetscCall(PetscCalloc1(numStrata, &tmpH)); 364 PetscCall(PetscCalloc1(numStrata, &tmpP)); 365 PetscCall(PetscMalloc1(numStrata, &tmpB)); 366 label->numStrata = numStrata; 367 label->stratumValues = tmpV; 368 label->stratumSizes = tmpS; 369 label->ht = tmpH; 370 label->points = tmpP; 371 label->validIS = tmpB; 372 for (v = 0; v < numStrata; ++v) { 373 PetscCall(PetscHSetICreate(&ht)); 374 PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is)); 375 PetscCall(PetscHMapISet(hmap, values[v], v)); 376 tmpV[v] = values[v]; 377 tmpS[v] = 0; 378 tmpH[v] = ht; 379 tmpP[v] = is; 380 tmpB[v] = PETSC_TRUE; 381 } 382 PetscCall(PetscObjectStateIncrease((PetscObject)label)); 383 } else { 384 for (v = 0; v < numStrata; ++v) PetscCall(DMLabelAddStratum(label, values[v])); 385 } 386 PetscCall(PetscFree(values)); 387 PetscFunctionReturn(PETSC_SUCCESS); 388 } 389 390 /*@ 391 DMLabelAddStrataIS - Adds new stratum values in a `DMLabel` 392 393 Not Collective 394 395 Input Parameters: 396 + label - The `DMLabel` 397 - valueIS - Index set with stratum values 398 399 Level: beginner 400 401 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()` 402 @*/ 403 PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS) 404 { 405 PetscInt numStrata; 406 const PetscInt *stratumValues; 407 408 PetscFunctionBegin; 409 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 410 PetscValidHeaderSpecific(valueIS, IS_CLASSID, 2); 411 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 412 PetscCall(ISGetLocalSize(valueIS, &numStrata)); 413 PetscCall(ISGetIndices(valueIS, &stratumValues)); 414 PetscCall(DMLabelAddStrata(label, numStrata, stratumValues)); 415 PetscFunctionReturn(PETSC_SUCCESS); 416 } 417 418 static PetscErrorCode DMLabelView_Concrete_Ascii(DMLabel label, PetscViewer viewer) 419 { 420 PetscInt v; 421 PetscMPIInt rank; 422 423 PetscFunctionBegin; 424 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank)); 425 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 426 if (label) { 427 const char *name; 428 429 PetscCall(PetscObjectGetName((PetscObject)label, &name)); 430 PetscCall(PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name)); 431 if (label->bt) PetscCall(PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", label->pStart, label->pEnd)); 432 for (v = 0; v < label->numStrata; ++v) { 433 const PetscInt value = label->stratumValues[v]; 434 const PetscInt *points; 435 PetscInt p; 436 437 PetscCall(ISGetIndices(label->points[v], &points)); 438 for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value)); 439 PetscCall(ISRestoreIndices(label->points[v], &points)); 440 } 441 } 442 PetscCall(PetscViewerFlush(viewer)); 443 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 444 PetscFunctionReturn(PETSC_SUCCESS); 445 } 446 447 static PetscErrorCode DMLabelView_Concrete(DMLabel label, PetscViewer viewer) 448 { 449 PetscBool isascii; 450 451 PetscFunctionBegin; 452 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 453 if (isascii) PetscCall(DMLabelView_Concrete_Ascii(label, viewer)); 454 PetscFunctionReturn(PETSC_SUCCESS); 455 } 456 457 /*@ 458 DMLabelView - View the label 459 460 Collective 461 462 Input Parameters: 463 + label - The `DMLabel` 464 - viewer - The `PetscViewer` 465 466 Level: intermediate 467 468 .seealso: `DMLabel`, `PetscViewer`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()` 469 @*/ 470 PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer) 471 { 472 PetscFunctionBegin; 473 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 474 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer)); 475 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 476 PetscCall(DMLabelMakeAllValid_Private(label)); 477 PetscUseTypeMethod(label, view, viewer); 478 PetscFunctionReturn(PETSC_SUCCESS); 479 } 480 481 /*@ 482 DMLabelViewFromOptions - View a `DMLabel` in a particular way based on a request in the options database 483 484 Collective 485 486 Input Parameters: 487 + label - the `DMLabel` object 488 . obj - optional object that provides the prefix for the options database (if `NULL` then the prefix in `obj` is used) 489 - name - option string that is used to activate viewing 490 491 Level: intermediate 492 493 Note: 494 See `PetscObjectViewFromOptions()` for a list of values that can be provided in the options database to determine how the `DMLabel` is viewed 495 496 .seealso: [](ch_dmbase), `DMLabel`, `DMLabelView()`, `PetscObjectViewFromOptions()`, `DMLabelCreate()` 497 @*/ 498 PetscErrorCode DMLabelViewFromOptions(DMLabel label, PeOp PetscObject obj, const char name[]) 499 { 500 PetscFunctionBegin; 501 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 502 PetscCall(PetscObjectViewFromOptions((PetscObject)label, obj, name)); 503 PetscFunctionReturn(PETSC_SUCCESS); 504 } 505 506 /*@ 507 DMLabelReset - Destroys internal data structures in a `DMLabel` 508 509 Not Collective 510 511 Input Parameter: 512 . label - The `DMLabel` 513 514 Level: beginner 515 516 .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`, `DMLabelCreate()` 517 @*/ 518 PetscErrorCode DMLabelReset(DMLabel label) 519 { 520 PetscInt v; 521 522 PetscFunctionBegin; 523 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 524 for (v = 0; v < label->numStrata; ++v) { 525 if (label->ht[v]) PetscCall(PetscHSetIDestroy(&label->ht[v])); 526 PetscCall(ISDestroy(&label->points[v])); 527 } 528 label->numStrata = 0; 529 PetscCall(PetscFree(label->stratumValues)); 530 PetscCall(PetscFree(label->stratumSizes)); 531 PetscCall(PetscFree(label->ht)); 532 PetscCall(PetscFree(label->points)); 533 PetscCall(PetscFree(label->validIS)); 534 PetscCall(PetscHMapIReset(label->hmap)); 535 label->pStart = -1; 536 label->pEnd = -1; 537 PetscCall(PetscBTDestroy(&label->bt)); 538 PetscFunctionReturn(PETSC_SUCCESS); 539 } 540 541 /*@ 542 DMLabelDestroy - Destroys a `DMLabel` 543 544 Collective 545 546 Input Parameter: 547 . label - The `DMLabel` 548 549 Level: beginner 550 551 .seealso: `DMLabel`, `DM`, `DMLabelReset()`, `DMLabelCreate()` 552 @*/ 553 PetscErrorCode DMLabelDestroy(DMLabel *label) 554 { 555 PetscFunctionBegin; 556 if (!*label) PetscFunctionReturn(PETSC_SUCCESS); 557 PetscValidHeaderSpecific(*label, DMLABEL_CLASSID, 1); 558 if (--((PetscObject)*label)->refct > 0) { 559 *label = NULL; 560 PetscFunctionReturn(PETSC_SUCCESS); 561 } 562 PetscCall(DMLabelReset(*label)); 563 PetscCall(PetscHMapIDestroy(&(*label)->hmap)); 564 PetscCall(PetscHeaderDestroy(label)); 565 PetscFunctionReturn(PETSC_SUCCESS); 566 } 567 568 static PetscErrorCode DMLabelDuplicate_Concrete(DMLabel label, DMLabel *labelnew) 569 { 570 PetscFunctionBegin; 571 for (PetscInt v = 0; v < label->numStrata; ++v) { 572 PetscCall(PetscHSetICreate(&(*labelnew)->ht[v])); 573 PetscCall(PetscObjectReference((PetscObject)label->points[v])); 574 (*labelnew)->points[v] = label->points[v]; 575 } 576 PetscCall(PetscHMapIDestroy(&(*labelnew)->hmap)); 577 PetscCall(PetscHMapIDuplicate(label->hmap, &(*labelnew)->hmap)); 578 PetscFunctionReturn(PETSC_SUCCESS); 579 } 580 581 /*@ 582 DMLabelDuplicate - Duplicates a `DMLabel` 583 584 Collective 585 586 Input Parameter: 587 . label - The `DMLabel` 588 589 Output Parameter: 590 . labelnew - new label 591 592 Level: intermediate 593 594 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()` 595 @*/ 596 PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew) 597 { 598 const char *name; 599 600 PetscFunctionBegin; 601 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 602 PetscCall(DMLabelMakeAllValid_Private(label)); 603 PetscCall(PetscObjectGetName((PetscObject)label, &name)); 604 PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew)); 605 606 (*labelnew)->numStrata = label->numStrata; 607 (*labelnew)->defaultValue = label->defaultValue; 608 (*labelnew)->readonly = label->readonly; 609 PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues)); 610 PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes)); 611 PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->ht)); 612 PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->points)); 613 PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->validIS)); 614 for (PetscInt v = 0; v < label->numStrata; ++v) { 615 (*labelnew)->stratumValues[v] = label->stratumValues[v]; 616 (*labelnew)->stratumSizes[v] = label->stratumSizes[v]; 617 (*labelnew)->validIS[v] = PETSC_TRUE; 618 } 619 (*labelnew)->pStart = -1; 620 (*labelnew)->pEnd = -1; 621 (*labelnew)->bt = NULL; 622 PetscUseTypeMethod(label, duplicate, labelnew); 623 PetscFunctionReturn(PETSC_SUCCESS); 624 } 625 626 /*@C 627 DMLabelCompare - Compare two `DMLabel` objects 628 629 Collective; No Fortran Support 630 631 Input Parameters: 632 + comm - Comm over which to compare labels 633 . l0 - First `DMLabel` 634 - l1 - Second `DMLabel` 635 636 Output Parameters: 637 + equal - (Optional) Flag whether the two labels are equal 638 - message - (Optional) Message describing the difference 639 640 Level: intermediate 641 642 Notes: 643 The output flag equal is the same on all processes. 644 If it is passed as `NULL` and difference is found, an error is thrown on all processes. 645 Make sure to pass `NULL` on all processes. 646 647 The output message is set independently on each rank. 648 It is set to `NULL` if no difference was found on the current rank. It must be freed by user. 649 If message is passed as `NULL` and difference is found, the difference description is printed to stderr in synchronized manner. 650 Make sure to pass `NULL` on all processes. 651 652 For the comparison, we ignore the order of stratum values, and strata with no points. 653 654 The communicator needs to be specified because currently `DMLabel` can live on `PETSC_COMM_SELF` even if the underlying `DM` is parallel. 655 656 Developer Note: 657 Fortran stub cannot be generated automatically because `message` must be freed with `PetscFree()` 658 659 .seealso: `DMLabel`, `DM`, `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()` 660 @*/ 661 PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message) 662 { 663 const char *name0, *name1; 664 char msg[PETSC_MAX_PATH_LEN] = ""; 665 PetscBool eq; 666 PetscMPIInt rank; 667 668 PetscFunctionBegin; 669 PetscValidHeaderSpecific(l0, DMLABEL_CLASSID, 2); 670 PetscValidHeaderSpecific(l1, DMLABEL_CLASSID, 3); 671 if (equal) PetscAssertPointer(equal, 4); 672 if (message) PetscAssertPointer(message, 5); 673 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 674 PetscCall(PetscObjectGetName((PetscObject)l0, &name0)); 675 PetscCall(PetscObjectGetName((PetscObject)l1, &name1)); 676 { 677 PetscInt v0, v1; 678 679 PetscCall(DMLabelGetDefaultValue(l0, &v0)); 680 PetscCall(DMLabelGetDefaultValue(l1, &v1)); 681 eq = (PetscBool)(v0 == v1); 682 if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Default value of DMLabel l0 \"%s\" = %" PetscInt_FMT " != %" PetscInt_FMT " = Default value of DMLabel l1 \"%s\"", name0, v0, v1, name1)); 683 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm)); 684 if (!eq) goto finish; 685 } 686 { 687 IS is0, is1; 688 689 PetscCall(DMLabelGetNonEmptyStratumValuesIS(l0, &is0)); 690 PetscCall(DMLabelGetNonEmptyStratumValuesIS(l1, &is1)); 691 PetscCall(ISEqual(is0, is1, &eq)); 692 PetscCall(ISDestroy(&is0)); 693 PetscCall(ISDestroy(&is1)); 694 if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1)); 695 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm)); 696 if (!eq) goto finish; 697 } 698 { 699 PetscInt i, nValues; 700 701 PetscCall(DMLabelGetNumValues(l0, &nValues)); 702 for (i = 0; i < nValues; i++) { 703 const PetscInt v = l0->stratumValues[i]; 704 PetscInt n; 705 IS is0, is1; 706 707 PetscCall(DMLabelGetStratumSize_Private(l0, i, &n)); 708 if (!n) continue; 709 PetscCall(DMLabelGetStratumIS(l0, v, &is0)); 710 PetscCall(DMLabelGetStratumIS(l1, v, &is1)); 711 PetscCall(ISEqualUnsorted(is0, is1, &eq)); 712 PetscCall(ISDestroy(&is0)); 713 PetscCall(ISDestroy(&is1)); 714 if (!eq) { 715 PetscCall(PetscSNPrintf(msg, sizeof(msg), "Stratum #%" PetscInt_FMT " with value %" PetscInt_FMT " contains different points in DMLabel l0 \"%s\" and DMLabel l1 \"%s\"", i, v, name0, name1)); 716 break; 717 } 718 } 719 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm)); 720 } 721 finish: 722 /* If message output arg not set, print to stderr */ 723 if (message) { 724 *message = NULL; 725 if (msg[0]) PetscCall(PetscStrallocpy(msg, message)); 726 } else { 727 if (msg[0]) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg)); 728 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR)); 729 } 730 /* If same output arg not ser and labels are not equal, throw error */ 731 if (equal) *equal = eq; 732 else PetscCheck(eq, comm, PETSC_ERR_ARG_INCOMP, "DMLabels l0 \"%s\" and l1 \"%s\" are not equal", name0, name1); 733 PetscFunctionReturn(PETSC_SUCCESS); 734 } 735 736 /*@ 737 DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds 738 739 Not Collective 740 741 Input Parameter: 742 . label - The `DMLabel` 743 744 Level: intermediate 745 746 .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()` 747 @*/ 748 PetscErrorCode DMLabelComputeIndex(DMLabel label) 749 { 750 PetscInt pStart = PETSC_INT_MAX, pEnd = -1, v; 751 752 PetscFunctionBegin; 753 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 754 PetscCall(DMLabelMakeAllValid_Private(label)); 755 for (v = 0; v < label->numStrata; ++v) { 756 const PetscInt *points; 757 PetscInt i; 758 759 PetscCall(ISGetIndices(label->points[v], &points)); 760 for (i = 0; i < label->stratumSizes[v]; ++i) { 761 const PetscInt point = points[i]; 762 763 pStart = PetscMin(point, pStart); 764 pEnd = PetscMax(point + 1, pEnd); 765 } 766 PetscCall(ISRestoreIndices(label->points[v], &points)); 767 } 768 label->pStart = pStart == PETSC_INT_MAX ? -1 : pStart; 769 label->pEnd = pEnd; 770 PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd)); 771 PetscFunctionReturn(PETSC_SUCCESS); 772 } 773 774 /*@ 775 DMLabelCreateIndex - Create an index structure for membership determination 776 777 Not Collective 778 779 Input Parameters: 780 + label - The `DMLabel` 781 . pStart - The smallest point 782 - pEnd - The largest point + 1 783 784 Level: intermediate 785 786 .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()` 787 @*/ 788 PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd) 789 { 790 PetscInt v; 791 792 PetscFunctionBegin; 793 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 794 PetscCall(DMLabelDestroyIndex(label)); 795 PetscCall(DMLabelMakeAllValid_Private(label)); 796 label->pStart = pStart; 797 label->pEnd = pEnd; 798 /* This can be hooked into SetValue(), ClearValue(), etc. for updating */ 799 PetscCall(PetscBTCreate(pEnd - pStart, &label->bt)); 800 for (v = 0; v < label->numStrata; ++v) { 801 IS pointIS; 802 const PetscInt *points; 803 PetscInt i; 804 805 PetscUseTypeMethod(label, getstratumis, v, &pointIS); 806 PetscCall(ISGetIndices(pointIS, &points)); 807 for (i = 0; i < label->stratumSizes[v]; ++i) { 808 const PetscInt point = points[i]; 809 810 PetscCheck(!(point < pStart) && !(point >= pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " in stratum %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->stratumValues[v], pStart, pEnd); 811 PetscCall(PetscBTSet(label->bt, point - pStart)); 812 } 813 PetscCall(ISRestoreIndices(label->points[v], &points)); 814 PetscCall(ISDestroy(&pointIS)); 815 } 816 PetscFunctionReturn(PETSC_SUCCESS); 817 } 818 819 /*@ 820 DMLabelDestroyIndex - Destroy the index structure 821 822 Not Collective 823 824 Input Parameter: 825 . label - the `DMLabel` 826 827 Level: intermediate 828 829 .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()` 830 @*/ 831 PetscErrorCode DMLabelDestroyIndex(DMLabel label) 832 { 833 PetscFunctionBegin; 834 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 835 label->pStart = -1; 836 label->pEnd = -1; 837 PetscCall(PetscBTDestroy(&label->bt)); 838 PetscFunctionReturn(PETSC_SUCCESS); 839 } 840 841 /*@ 842 DMLabelGetBounds - Return the smallest and largest point in the label 843 844 Not Collective 845 846 Input Parameter: 847 . label - the `DMLabel` 848 849 Output Parameters: 850 + pStart - The smallest point 851 - pEnd - The largest point + 1 852 853 Level: intermediate 854 855 Note: 856 This will compute an index for the label if one does not exist. 857 858 .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()` 859 @*/ 860 PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd) 861 { 862 PetscFunctionBegin; 863 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 864 if ((label->pStart == -1) && (label->pEnd == -1)) PetscCall(DMLabelComputeIndex(label)); 865 if (pStart) { 866 PetscAssertPointer(pStart, 2); 867 *pStart = label->pStart; 868 } 869 if (pEnd) { 870 PetscAssertPointer(pEnd, 3); 871 *pEnd = label->pEnd; 872 } 873 PetscFunctionReturn(PETSC_SUCCESS); 874 } 875 876 /*@ 877 DMLabelHasValue - Determine whether a label assigns the value to any point 878 879 Not Collective 880 881 Input Parameters: 882 + label - the `DMLabel` 883 - value - the value 884 885 Output Parameter: 886 . contains - Flag indicating whether the label maps this value to any point 887 888 Level: developer 889 890 .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()` 891 @*/ 892 PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains) 893 { 894 PetscInt v; 895 896 PetscFunctionBegin; 897 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 898 PetscAssertPointer(contains, 3); 899 PetscCall(DMLabelLookupStratum(label, value, &v)); 900 *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE; 901 PetscFunctionReturn(PETSC_SUCCESS); 902 } 903 904 /*@ 905 DMLabelHasPoint - Determine whether a label assigns a value to a point 906 907 Not Collective 908 909 Input Parameters: 910 + label - the `DMLabel` 911 - point - the point 912 913 Output Parameter: 914 . contains - Flag indicating whether the label maps this point to a value 915 916 Level: developer 917 918 Note: 919 The user must call `DMLabelCreateIndex()` before this function. 920 921 .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()` 922 @*/ 923 PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains) 924 { 925 PetscInt pStart, pEnd; 926 927 PetscFunctionBeginHot; 928 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 929 PetscAssertPointer(contains, 3); 930 /* DMLabelGetBounds() calls DMLabelCreateIndex() only if needed */ 931 PetscCall(DMLabelGetBounds(label, &pStart, &pEnd)); 932 PetscCall(DMLabelMakeAllValid_Private(label)); 933 *contains = point >= pStart && point < pEnd && (PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE); 934 PetscFunctionReturn(PETSC_SUCCESS); 935 } 936 937 /*@ 938 DMLabelStratumHasPoint - Return true if the stratum contains a point 939 940 Not Collective 941 942 Input Parameters: 943 + label - the `DMLabel` 944 . value - the stratum value 945 - point - the point 946 947 Output Parameter: 948 . contains - true if the stratum contains the point 949 950 Level: intermediate 951 952 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()` 953 @*/ 954 PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains) 955 { 956 PetscFunctionBeginHot; 957 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 958 PetscAssertPointer(contains, 4); 959 if (value == label->defaultValue) { 960 PetscInt pointVal; 961 962 PetscCall(DMLabelGetValue(label, point, &pointVal)); 963 *contains = (PetscBool)(pointVal == value); 964 } else { 965 PetscInt v; 966 967 PetscCall(DMLabelLookupStratum(label, value, &v)); 968 if (v >= 0) { 969 if (label->validIS[v] || label->readonly) { 970 IS is; 971 PetscInt i; 972 973 PetscUseTypeMethod(label, getstratumis, v, &is); 974 PetscCall(ISLocate(is, point, &i)); 975 PetscCall(ISDestroy(&is)); 976 *contains = (PetscBool)(i >= 0); 977 } else { 978 PetscCall(PetscHSetIHas(label->ht[v], point, contains)); 979 } 980 } else { // value is not present 981 *contains = PETSC_FALSE; 982 } 983 } 984 PetscFunctionReturn(PETSC_SUCCESS); 985 } 986 987 /*@ 988 DMLabelGetDefaultValue - Get the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value. 989 When a label is created, it is initialized to -1. 990 991 Not Collective 992 993 Input Parameter: 994 . label - a `DMLabel` object 995 996 Output Parameter: 997 . defaultValue - the default value 998 999 Level: beginner 1000 1001 .seealso: `DMLabel`, `DM`, `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()` 1002 @*/ 1003 PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue) 1004 { 1005 PetscFunctionBegin; 1006 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1007 *defaultValue = label->defaultValue; 1008 PetscFunctionReturn(PETSC_SUCCESS); 1009 } 1010 1011 /*@ 1012 DMLabelSetDefaultValue - Set the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value. 1013 When a label is created, it is initialized to -1. 1014 1015 Not Collective 1016 1017 Input Parameter: 1018 . label - a `DMLabel` object 1019 1020 Output Parameter: 1021 . defaultValue - the default value 1022 1023 Level: beginner 1024 1025 .seealso: `DMLabel`, `DM`, `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()` 1026 @*/ 1027 PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue) 1028 { 1029 PetscFunctionBegin; 1030 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1031 label->defaultValue = defaultValue; 1032 PetscFunctionReturn(PETSC_SUCCESS); 1033 } 1034 1035 /*@ 1036 DMLabelGetValue - Return the value a label assigns to a point, or the label's default value (which is initially -1, and can be changed with 1037 `DMLabelSetDefaultValue()`) 1038 1039 Not Collective 1040 1041 Input Parameters: 1042 + label - the `DMLabel` 1043 - point - the point 1044 1045 Output Parameter: 1046 . value - The point value, or the default value (-1 by default) 1047 1048 Level: intermediate 1049 1050 Note: 1051 A label may assign multiple values to a point. No guarantees are made about which value is returned in that case. 1052 Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum. 1053 1054 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()` 1055 @*/ 1056 PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value) 1057 { 1058 PetscInt v; 1059 1060 PetscFunctionBeginHot; 1061 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1062 PetscAssertPointer(value, 3); 1063 *value = label->defaultValue; 1064 for (v = 0; v < label->numStrata; ++v) { 1065 if (label->validIS[v] || label->readonly) { 1066 IS is; 1067 PetscInt i; 1068 1069 PetscUseTypeMethod(label, getstratumis, v, &is); 1070 PetscCall(ISLocate(label->points[v], point, &i)); 1071 PetscCall(ISDestroy(&is)); 1072 if (i >= 0) { 1073 *value = label->stratumValues[v]; 1074 break; 1075 } 1076 } else { 1077 PetscBool has; 1078 1079 PetscCall(PetscHSetIHas(label->ht[v], point, &has)); 1080 if (has) { 1081 *value = label->stratumValues[v]; 1082 break; 1083 } 1084 } 1085 } 1086 PetscFunctionReturn(PETSC_SUCCESS); 1087 } 1088 1089 /*@ 1090 DMLabelSetValue - Set the value a label assigns to a point. If the value is the same as the label's default value (which is initially -1, and can 1091 be changed with `DMLabelSetDefaultValue()` to something different), then this function will do nothing. 1092 1093 Not Collective 1094 1095 Input Parameters: 1096 + label - the `DMLabel` 1097 . point - the point 1098 - value - The point value 1099 1100 Level: intermediate 1101 1102 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()` 1103 @*/ 1104 PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value) 1105 { 1106 PetscInt v; 1107 1108 PetscFunctionBegin; 1109 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1110 /* Find label value, add new entry if needed */ 1111 if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS); 1112 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1113 PetscCall(DMLabelLookupAddStratum(label, value, &v)); 1114 /* Set key */ 1115 PetscCall(DMLabelMakeInvalid_Private(label, v)); 1116 PetscCall(PetscHSetIAdd(label->ht[v], point)); 1117 PetscFunctionReturn(PETSC_SUCCESS); 1118 } 1119 1120 /*@ 1121 DMLabelClearValue - Clear the value a label assigns to a point 1122 1123 Not Collective 1124 1125 Input Parameters: 1126 + label - the `DMLabel` 1127 . point - the point 1128 - value - The point value 1129 1130 Level: intermediate 1131 1132 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()` 1133 @*/ 1134 PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value) 1135 { 1136 PetscInt v; 1137 1138 PetscFunctionBegin; 1139 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1140 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1141 /* Find label value */ 1142 PetscCall(DMLabelLookupStratum(label, value, &v)); 1143 if (v < 0) PetscFunctionReturn(PETSC_SUCCESS); 1144 1145 if (label->bt && point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart)); 1146 1147 /* Delete key */ 1148 PetscCall(DMLabelMakeInvalid_Private(label, v)); 1149 PetscCall(PetscHSetIDel(label->ht[v], point)); 1150 PetscFunctionReturn(PETSC_SUCCESS); 1151 } 1152 1153 /*@ 1154 DMLabelInsertIS - Set all points in the `IS` to a value 1155 1156 Not Collective 1157 1158 Input Parameters: 1159 + label - the `DMLabel` 1160 . is - the point `IS` 1161 - value - The point value 1162 1163 Level: intermediate 1164 1165 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1166 @*/ 1167 PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value) 1168 { 1169 PetscInt v, n, p; 1170 const PetscInt *points; 1171 1172 PetscFunctionBegin; 1173 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1174 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 1175 /* Find label value, add new entry if needed */ 1176 if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS); 1177 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1178 PetscCall(DMLabelLookupAddStratum(label, value, &v)); 1179 /* Set keys */ 1180 PetscCall(DMLabelMakeInvalid_Private(label, v)); 1181 PetscCall(ISGetLocalSize(is, &n)); 1182 PetscCall(ISGetIndices(is, &points)); 1183 for (p = 0; p < n; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p])); 1184 PetscCall(ISRestoreIndices(is, &points)); 1185 PetscFunctionReturn(PETSC_SUCCESS); 1186 } 1187 1188 /*@ 1189 DMLabelGetNumValues - Get the number of values that the `DMLabel` takes 1190 1191 Not Collective 1192 1193 Input Parameter: 1194 . label - the `DMLabel` 1195 1196 Output Parameter: 1197 . numValues - the number of values 1198 1199 Level: intermediate 1200 1201 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1202 @*/ 1203 PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues) 1204 { 1205 PetscFunctionBegin; 1206 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1207 PetscAssertPointer(numValues, 2); 1208 *numValues = label->numStrata; 1209 PetscFunctionReturn(PETSC_SUCCESS); 1210 } 1211 1212 /*@ 1213 DMLabelGetValueIS - Get an `IS` of all values that the `DMlabel` takes 1214 1215 Not Collective 1216 1217 Input Parameter: 1218 . label - the `DMLabel` 1219 1220 Output Parameter: 1221 . values - the value `IS` 1222 1223 Level: intermediate 1224 1225 Notes: 1226 The output `IS` should be destroyed when no longer needed. 1227 Strata which are allocated but empty [`DMLabelGetStratumSize()` yields 0] are counted. 1228 If you need to count only nonempty strata, use `DMLabelGetNonEmptyStratumValuesIS()`. 1229 1230 .seealso: `DMLabel`, `DM`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1231 @*/ 1232 PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values) 1233 { 1234 PetscFunctionBegin; 1235 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1236 PetscAssertPointer(values, 2); 1237 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values)); 1238 PetscFunctionReturn(PETSC_SUCCESS); 1239 } 1240 1241 /*@ 1242 DMLabelGetValueBounds - Return the smallest and largest value in the label 1243 1244 Not Collective 1245 1246 Input Parameter: 1247 . label - the `DMLabel` 1248 1249 Output Parameters: 1250 + minValue - The smallest value 1251 - maxValue - The largest value 1252 1253 Level: intermediate 1254 1255 .seealso: `DMLabel`, `DM`, `DMLabelGetBounds()`, `DMLabelGetValue()`, `DMLabelSetValue()` 1256 @*/ 1257 PetscErrorCode DMLabelGetValueBounds(DMLabel label, PetscInt *minValue, PetscInt *maxValue) 1258 { 1259 PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN; 1260 1261 PetscFunctionBegin; 1262 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1263 for (PetscInt v = 0; v < label->numStrata; ++v) { 1264 min = PetscMin(min, label->stratumValues[v]); 1265 max = PetscMax(max, label->stratumValues[v]); 1266 } 1267 if (minValue) { 1268 PetscAssertPointer(minValue, 2); 1269 *minValue = min; 1270 } 1271 if (maxValue) { 1272 PetscAssertPointer(maxValue, 3); 1273 *maxValue = max; 1274 } 1275 PetscFunctionReturn(PETSC_SUCCESS); 1276 } 1277 1278 /*@ 1279 DMLabelGetNonEmptyStratumValuesIS - Get an `IS` of all values that the `DMlabel` takes 1280 1281 Not Collective 1282 1283 Input Parameter: 1284 . label - the `DMLabel` 1285 1286 Output Parameter: 1287 . values - the value `IS` 1288 1289 Level: intermediate 1290 1291 Notes: 1292 The output `IS` should be destroyed when no longer needed. 1293 This is similar to `DMLabelGetValueIS()` but counts only nonempty strata. 1294 1295 .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1296 @*/ 1297 PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values) 1298 { 1299 PetscInt i, j; 1300 PetscInt *valuesArr; 1301 1302 PetscFunctionBegin; 1303 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1304 PetscAssertPointer(values, 2); 1305 PetscCall(PetscMalloc1(label->numStrata, &valuesArr)); 1306 for (i = 0, j = 0; i < label->numStrata; i++) { 1307 PetscInt n; 1308 1309 PetscCall(DMLabelGetStratumSize_Private(label, i, &n)); 1310 if (n) valuesArr[j++] = label->stratumValues[i]; 1311 } 1312 if (j == label->numStrata) { 1313 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values)); 1314 } else { 1315 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values)); 1316 } 1317 PetscCall(PetscFree(valuesArr)); 1318 PetscFunctionReturn(PETSC_SUCCESS); 1319 } 1320 1321 /*@ 1322 DMLabelGetValueIndex - Get the index of a given value in the list of values for the `DMlabel`, or -1 if it is not present 1323 1324 Not Collective 1325 1326 Input Parameters: 1327 + label - the `DMLabel` 1328 - value - the value 1329 1330 Output Parameter: 1331 . index - the index of value in the list of values 1332 1333 Level: intermediate 1334 1335 .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1336 @*/ 1337 PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index) 1338 { 1339 PetscInt v; 1340 1341 PetscFunctionBegin; 1342 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1343 PetscAssertPointer(index, 3); 1344 /* Do not assume they are sorted */ 1345 for (v = 0; v < label->numStrata; ++v) 1346 if (label->stratumValues[v] == value) break; 1347 if (v >= label->numStrata) *index = -1; 1348 else *index = v; 1349 PetscFunctionReturn(PETSC_SUCCESS); 1350 } 1351 1352 /*@ 1353 DMLabelHasStratum - Determine whether points exist with the given value 1354 1355 Not Collective 1356 1357 Input Parameters: 1358 + label - the `DMLabel` 1359 - value - the stratum value 1360 1361 Output Parameter: 1362 . exists - Flag saying whether points exist 1363 1364 Level: intermediate 1365 1366 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1367 @*/ 1368 PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists) 1369 { 1370 PetscInt v; 1371 1372 PetscFunctionBegin; 1373 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1374 PetscAssertPointer(exists, 3); 1375 PetscCall(DMLabelLookupStratum(label, value, &v)); 1376 *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE; 1377 PetscFunctionReturn(PETSC_SUCCESS); 1378 } 1379 1380 /*@ 1381 DMLabelGetStratumSize - Get the size of a stratum 1382 1383 Not Collective 1384 1385 Input Parameters: 1386 + label - the `DMLabel` 1387 - value - the stratum value 1388 1389 Output Parameter: 1390 . size - The number of points in the stratum 1391 1392 Level: intermediate 1393 1394 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1395 @*/ 1396 PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size) 1397 { 1398 PetscInt v; 1399 1400 PetscFunctionBegin; 1401 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1402 PetscAssertPointer(size, 3); 1403 PetscCall(DMLabelLookupStratum(label, value, &v)); 1404 PetscCall(DMLabelGetStratumSize_Private(label, v, size)); 1405 PetscFunctionReturn(PETSC_SUCCESS); 1406 } 1407 1408 /*@ 1409 DMLabelGetStratumBounds - Get the largest and smallest point of a stratum 1410 1411 Not Collective 1412 1413 Input Parameters: 1414 + label - the `DMLabel` 1415 - value - the stratum value 1416 1417 Output Parameters: 1418 + start - the smallest point in the stratum 1419 - end - the largest point in the stratum 1420 1421 Level: intermediate 1422 1423 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1424 @*/ 1425 PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end) 1426 { 1427 IS is; 1428 PetscInt v, min, max; 1429 1430 PetscFunctionBegin; 1431 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1432 if (start) { 1433 PetscAssertPointer(start, 3); 1434 *start = -1; 1435 } 1436 if (end) { 1437 PetscAssertPointer(end, 4); 1438 *end = -1; 1439 } 1440 PetscCall(DMLabelLookupStratum(label, value, &v)); 1441 if (v < 0) PetscFunctionReturn(PETSC_SUCCESS); 1442 PetscCall(DMLabelMakeValid_Private(label, v)); 1443 if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS); 1444 PetscUseTypeMethod(label, getstratumis, v, &is); 1445 PetscCall(ISGetMinMax(is, &min, &max)); 1446 PetscCall(ISDestroy(&is)); 1447 if (start) *start = min; 1448 if (end) *end = max + 1; 1449 PetscFunctionReturn(PETSC_SUCCESS); 1450 } 1451 1452 static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS) 1453 { 1454 PetscFunctionBegin; 1455 PetscCall(PetscObjectReference((PetscObject)label->points[v])); 1456 *pointIS = label->points[v]; 1457 PetscFunctionReturn(PETSC_SUCCESS); 1458 } 1459 1460 /*@ 1461 DMLabelGetStratumIS - Get an `IS` with the stratum points 1462 1463 Not Collective 1464 1465 Input Parameters: 1466 + label - the `DMLabel` 1467 - value - the stratum value 1468 1469 Output Parameter: 1470 . points - The stratum points 1471 1472 Level: intermediate 1473 1474 Notes: 1475 The output `IS` should be destroyed when no longer needed. 1476 Returns `NULL` if the stratum is empty. 1477 1478 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1479 @*/ 1480 PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points) 1481 { 1482 PetscInt v; 1483 1484 PetscFunctionBegin; 1485 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1486 PetscAssertPointer(points, 3); 1487 *points = NULL; 1488 PetscCall(DMLabelLookupStratum(label, value, &v)); 1489 if (v < 0) PetscFunctionReturn(PETSC_SUCCESS); 1490 PetscCall(DMLabelMakeValid_Private(label, v)); 1491 PetscUseTypeMethod(label, getstratumis, v, points); 1492 PetscFunctionReturn(PETSC_SUCCESS); 1493 } 1494 1495 /*@ 1496 DMLabelSetStratumIS - Set the stratum points using an `IS` 1497 1498 Not Collective 1499 1500 Input Parameters: 1501 + label - the `DMLabel` 1502 . value - the stratum value 1503 - is - The stratum points 1504 1505 Level: intermediate 1506 1507 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1508 @*/ 1509 PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is) 1510 { 1511 PetscInt v; 1512 1513 PetscFunctionBegin; 1514 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1515 PetscValidHeaderSpecific(is, IS_CLASSID, 3); 1516 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1517 PetscCall(DMLabelLookupAddStratum(label, value, &v)); 1518 if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS); 1519 PetscCall(DMLabelClearStratum(label, value)); 1520 PetscCall(ISGetLocalSize(is, &label->stratumSizes[v])); 1521 PetscCall(PetscObjectReference((PetscObject)is)); 1522 PetscCall(ISDestroy(&label->points[v])); 1523 label->points[v] = is; 1524 label->validIS[v] = PETSC_TRUE; 1525 PetscCall(PetscObjectStateIncrease((PetscObject)label)); 1526 if (label->bt) { 1527 const PetscInt *points; 1528 PetscInt p; 1529 1530 PetscCall(ISGetIndices(is, &points)); 1531 for (p = 0; p < label->stratumSizes[v]; ++p) { 1532 const PetscInt point = points[p]; 1533 1534 PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd); 1535 PetscCall(PetscBTSet(label->bt, point - label->pStart)); 1536 } 1537 } 1538 PetscFunctionReturn(PETSC_SUCCESS); 1539 } 1540 1541 /*@ 1542 DMLabelClearStratum - Remove a stratum 1543 1544 Not Collective 1545 1546 Input Parameters: 1547 + label - the `DMLabel` 1548 - value - the stratum value 1549 1550 Level: intermediate 1551 1552 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1553 @*/ 1554 PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value) 1555 { 1556 PetscInt v; 1557 1558 PetscFunctionBegin; 1559 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1560 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1561 PetscCall(DMLabelLookupStratum(label, value, &v)); 1562 if (v < 0) PetscFunctionReturn(PETSC_SUCCESS); 1563 if (label->validIS[v]) { 1564 if (label->bt) { 1565 PetscInt i; 1566 const PetscInt *points; 1567 1568 PetscCall(ISGetIndices(label->points[v], &points)); 1569 for (i = 0; i < label->stratumSizes[v]; ++i) { 1570 const PetscInt point = points[i]; 1571 1572 if (point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart)); 1573 } 1574 PetscCall(ISRestoreIndices(label->points[v], &points)); 1575 } 1576 label->stratumSizes[v] = 0; 1577 PetscCall(ISDestroy(&label->points[v])); 1578 PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v])); 1579 PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices")); 1580 PetscCall(PetscObjectStateIncrease((PetscObject)label)); 1581 } else { 1582 PetscCall(PetscHSetIClear(label->ht[v])); 1583 } 1584 PetscFunctionReturn(PETSC_SUCCESS); 1585 } 1586 1587 /*@ 1588 DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value 1589 1590 Not Collective 1591 1592 Input Parameters: 1593 + label - The `DMLabel` 1594 . value - The label value for all points 1595 . pStart - The first point 1596 - pEnd - A point beyond all marked points 1597 1598 Level: intermediate 1599 1600 Note: 1601 The marks points are [`pStart`, `pEnd`), and only the bounds are stored. 1602 1603 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()` 1604 @*/ 1605 PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd) 1606 { 1607 IS pIS; 1608 1609 PetscFunctionBegin; 1610 PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS)); 1611 PetscCall(DMLabelSetStratumIS(label, value, pIS)); 1612 PetscCall(ISDestroy(&pIS)); 1613 PetscFunctionReturn(PETSC_SUCCESS); 1614 } 1615 1616 /*@ 1617 DMLabelGetStratumPointIndex - Get the index of a point in a given stratum 1618 1619 Not Collective 1620 1621 Input Parameters: 1622 + label - The `DMLabel` 1623 . value - The label value 1624 - p - A point with this value 1625 1626 Output Parameter: 1627 . index - The index of this point in the stratum, or -1 if the point is not in the stratum or the stratum does not exist 1628 1629 Level: intermediate 1630 1631 .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()` 1632 @*/ 1633 PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index) 1634 { 1635 IS pointIS; 1636 PetscInt v; 1637 1638 PetscFunctionBegin; 1639 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1640 PetscAssertPointer(index, 4); 1641 *index = -1; 1642 PetscCall(DMLabelLookupStratum(label, value, &v)); 1643 if (v < 0) PetscFunctionReturn(PETSC_SUCCESS); 1644 PetscCall(DMLabelMakeValid_Private(label, v)); 1645 PetscUseTypeMethod(label, getstratumis, v, &pointIS); 1646 PetscCall(ISLocate(pointIS, p, index)); 1647 PetscCall(ISDestroy(&pointIS)); 1648 PetscFunctionReturn(PETSC_SUCCESS); 1649 } 1650 1651 /*@ 1652 DMLabelFilter - Remove all points outside of [`start`, `end`) 1653 1654 Not Collective 1655 1656 Input Parameters: 1657 + label - the `DMLabel` 1658 . start - the first point kept 1659 - end - one more than the last point kept 1660 1661 Level: intermediate 1662 1663 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1664 @*/ 1665 PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end) 1666 { 1667 PetscInt v; 1668 1669 PetscFunctionBegin; 1670 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1671 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1672 PetscCall(DMLabelDestroyIndex(label)); 1673 PetscCall(DMLabelMakeAllValid_Private(label)); 1674 for (v = 0; v < label->numStrata; ++v) { 1675 PetscCall(ISGeneralFilter(label->points[v], start, end)); 1676 PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v])); 1677 } 1678 PetscCall(DMLabelCreateIndex(label, start, end)); 1679 PetscFunctionReturn(PETSC_SUCCESS); 1680 } 1681 1682 /*@ 1683 DMLabelPermute - Create a new label with permuted points 1684 1685 Not Collective 1686 1687 Input Parameters: 1688 + label - the `DMLabel` 1689 - permutation - the point permutation 1690 1691 Output Parameter: 1692 . labelNew - the new label containing the permuted points 1693 1694 Level: intermediate 1695 1696 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1697 @*/ 1698 PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew) 1699 { 1700 const PetscInt *perm; 1701 PetscInt numValues, numPoints, v, q; 1702 1703 PetscFunctionBegin; 1704 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1705 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2); 1706 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1707 PetscCall(DMLabelMakeAllValid_Private(label)); 1708 PetscCall(DMLabelDuplicate(label, labelNew)); 1709 PetscCall(DMLabelGetNumValues(*labelNew, &numValues)); 1710 PetscCall(ISGetLocalSize(permutation, &numPoints)); 1711 PetscCall(ISGetIndices(permutation, &perm)); 1712 for (v = 0; v < numValues; ++v) { 1713 const PetscInt size = (*labelNew)->stratumSizes[v]; 1714 const PetscInt *points; 1715 PetscInt *pointsNew; 1716 1717 PetscCall(ISGetIndices((*labelNew)->points[v], &points)); 1718 PetscCall(PetscCalloc1(size, &pointsNew)); 1719 for (q = 0; q < size; ++q) { 1720 const PetscInt point = points[q]; 1721 1722 PetscCheck(!(point < 0) && !(point >= numPoints), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ") for the remapping", point, numPoints); 1723 pointsNew[q] = perm[point]; 1724 } 1725 PetscCall(ISRestoreIndices((*labelNew)->points[v], &points)); 1726 PetscCall(PetscSortInt(size, pointsNew)); 1727 PetscCall(ISDestroy(&(*labelNew)->points[v])); 1728 if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) { 1729 PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v]))); 1730 PetscCall(PetscFree(pointsNew)); 1731 } else { 1732 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v]))); 1733 } 1734 PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices")); 1735 } 1736 PetscCall(ISRestoreIndices(permutation, &perm)); 1737 if (label->bt) { 1738 PetscCall(PetscBTDestroy(&label->bt)); 1739 PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd)); 1740 } 1741 PetscFunctionReturn(PETSC_SUCCESS); 1742 } 1743 1744 /*@ 1745 DMLabelPermuteValues - Permute the values in a label 1746 1747 Not collective 1748 1749 Input Parameters: 1750 + label - the `DMLabel` 1751 - permutation - the value permutation, permutation[old value] = new value 1752 1753 Output Parameter: 1754 . label - the `DMLabel` now with permuted values 1755 1756 Note: 1757 The modification is done in-place 1758 1759 Level: intermediate 1760 1761 .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1762 @*/ 1763 PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation) 1764 { 1765 PetscInt Nv, Np; 1766 1767 PetscFunctionBegin; 1768 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1769 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2); 1770 PetscCall(DMLabelGetNumValues(label, &Nv)); 1771 PetscCall(ISGetLocalSize(permutation, &Np)); 1772 PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv); 1773 if (PetscDefined(USE_DEBUG)) { 1774 PetscBool flg; 1775 PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg)); 1776 PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation"); 1777 } 1778 PetscCall(DMLabelRewriteValues(label, permutation)); 1779 PetscFunctionReturn(PETSC_SUCCESS); 1780 } 1781 1782 /*@ 1783 DMLabelRewriteValues - Permute the values in a label, but some may be omitted 1784 1785 Not collective 1786 1787 Input Parameters: 1788 + label - the `DMLabel` 1789 - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted 1790 1791 Output Parameter: 1792 . label - the `DMLabel` now with permuted values 1793 1794 Note: 1795 The modification is done in-place 1796 1797 Level: intermediate 1798 1799 .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1800 @*/ 1801 PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation) 1802 { 1803 const PetscInt *perm; 1804 PetscInt Nv, Np; 1805 1806 PetscFunctionBegin; 1807 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1808 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2); 1809 PetscCall(DMLabelMakeAllValid_Private(label)); 1810 PetscCall(DMLabelGetNumValues(label, &Nv)); 1811 PetscCall(ISGetLocalSize(permutation, &Np)); 1812 PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv); 1813 PetscCall(ISGetIndices(permutation, &perm)); 1814 for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]]; 1815 PetscCall(ISRestoreIndices(permutation, &perm)); 1816 PetscFunctionReturn(PETSC_SUCCESS); 1817 } 1818 1819 static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata) 1820 { 1821 MPI_Comm comm; 1822 PetscInt s, l, nroots, nleaves, offset, size; 1823 PetscInt *remoteOffsets, *rootStrata, *rootIdx; 1824 PetscSection rootSection; 1825 PetscSF labelSF; 1826 1827 PetscFunctionBegin; 1828 if (label) PetscCall(DMLabelMakeAllValid_Private(label)); 1829 PetscCall(PetscObjectGetComm((PetscObject)sf, &comm)); 1830 /* Build a section of stratum values per point, generate the according SF 1831 and distribute point-wise stratum values to leaves. */ 1832 PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL)); 1833 PetscCall(PetscSectionCreate(comm, &rootSection)); 1834 PetscCall(PetscSectionSetChart(rootSection, 0, nroots)); 1835 if (label) { 1836 for (s = 0; s < label->numStrata; ++s) { 1837 const PetscInt *points; 1838 1839 PetscCall(ISGetIndices(label->points[s], &points)); 1840 for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1)); 1841 PetscCall(ISRestoreIndices(label->points[s], &points)); 1842 } 1843 } 1844 PetscCall(PetscSectionSetUp(rootSection)); 1845 /* Create a point-wise array of stratum values */ 1846 PetscCall(PetscSectionGetStorageSize(rootSection, &size)); 1847 PetscCall(PetscMalloc1(size, &rootStrata)); 1848 PetscCall(PetscCalloc1(nroots, &rootIdx)); 1849 if (label) { 1850 for (s = 0; s < label->numStrata; ++s) { 1851 const PetscInt *points; 1852 1853 PetscCall(ISGetIndices(label->points[s], &points)); 1854 for (l = 0; l < label->stratumSizes[s]; l++) { 1855 const PetscInt p = points[l]; 1856 PetscCall(PetscSectionGetOffset(rootSection, p, &offset)); 1857 rootStrata[offset + rootIdx[p]++] = label->stratumValues[s]; 1858 } 1859 PetscCall(ISRestoreIndices(label->points[s], &points)); 1860 } 1861 } 1862 /* Build SF that maps label points to remote processes */ 1863 PetscCall(PetscSectionCreate(comm, leafSection)); 1864 PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection)); 1865 PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF)); 1866 PetscCall(PetscFree(remoteOffsets)); 1867 /* Send the strata for each point over the derived SF */ 1868 PetscCall(PetscSectionGetStorageSize(*leafSection, &size)); 1869 PetscCall(PetscMalloc1(size, leafStrata)); 1870 PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE)); 1871 PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE)); 1872 /* Clean up */ 1873 PetscCall(PetscFree(rootStrata)); 1874 PetscCall(PetscFree(rootIdx)); 1875 PetscCall(PetscSectionDestroy(&rootSection)); 1876 PetscCall(PetscSFDestroy(&labelSF)); 1877 PetscFunctionReturn(PETSC_SUCCESS); 1878 } 1879 1880 /*@ 1881 DMLabelDistribute - Create a new label pushed forward over the `PetscSF` 1882 1883 Collective 1884 1885 Input Parameters: 1886 + label - the `DMLabel` 1887 - sf - the map from old to new distribution 1888 1889 Output Parameter: 1890 . labelNew - the new redistributed label 1891 1892 Level: intermediate 1893 1894 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1895 @*/ 1896 PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew) 1897 { 1898 MPI_Comm comm; 1899 PetscSection leafSection; 1900 PetscInt p, pStart, pEnd, s, size, dof, offset, stratum; 1901 PetscInt *leafStrata, *strataIdx; 1902 PetscInt **points; 1903 const char *lname = NULL; 1904 char *name; 1905 PetscMPIInt nameSize; 1906 PetscHSetI stratumHash; 1907 size_t len = 0; 1908 PetscMPIInt rank; 1909 1910 PetscFunctionBegin; 1911 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 1912 if (label) { 1913 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1914 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1915 PetscCall(DMLabelMakeAllValid_Private(label)); 1916 } 1917 PetscCall(PetscObjectGetComm((PetscObject)sf, &comm)); 1918 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1919 /* Bcast name */ 1920 if (rank == 0) { 1921 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1922 PetscCall(PetscStrlen(lname, &len)); 1923 } 1924 PetscCall(PetscMPIIntCast(len, &nameSize)); 1925 PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm)); 1926 PetscCall(PetscMalloc1(nameSize + 1, &name)); 1927 if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1)); 1928 PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm)); 1929 PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew)); 1930 PetscCall(PetscFree(name)); 1931 /* Bcast defaultValue */ 1932 if (rank == 0) (*labelNew)->defaultValue = label->defaultValue; 1933 PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm)); 1934 /* Distribute stratum values over the SF and get the point mapping on the receiver */ 1935 PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata)); 1936 /* Determine received stratum values and initialise new label*/ 1937 PetscCall(PetscHSetICreate(&stratumHash)); 1938 PetscCall(PetscSectionGetStorageSize(leafSection, &size)); 1939 for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p])); 1940 PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata)); 1941 PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS)); 1942 for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE; 1943 PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues)); 1944 /* Turn leafStrata into indices rather than stratum values */ 1945 offset = 0; 1946 PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues)); 1947 PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues)); 1948 for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s)); 1949 for (p = 0; p < size; ++p) { 1950 for (s = 0; s < (*labelNew)->numStrata; ++s) { 1951 if (leafStrata[p] == (*labelNew)->stratumValues[s]) { 1952 leafStrata[p] = s; 1953 break; 1954 } 1955 } 1956 } 1957 /* Rebuild the point strata on the receiver */ 1958 PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes)); 1959 PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd)); 1960 for (p = pStart; p < pEnd; p++) { 1961 PetscCall(PetscSectionGetDof(leafSection, p, &dof)); 1962 PetscCall(PetscSectionGetOffset(leafSection, p, &offset)); 1963 for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++; 1964 } 1965 PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht)); 1966 PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points)); 1967 PetscCall(PetscCalloc1((*labelNew)->numStrata, &points)); 1968 for (s = 0; s < (*labelNew)->numStrata; ++s) { 1969 PetscCall(PetscHSetICreate(&(*labelNew)->ht[s])); 1970 PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s])); 1971 } 1972 /* Insert points into new strata */ 1973 PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx)); 1974 PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd)); 1975 for (p = pStart; p < pEnd; p++) { 1976 PetscCall(PetscSectionGetDof(leafSection, p, &dof)); 1977 PetscCall(PetscSectionGetOffset(leafSection, p, &offset)); 1978 for (s = 0; s < dof; s++) { 1979 stratum = leafStrata[offset + s]; 1980 points[stratum][strataIdx[stratum]++] = p; 1981 } 1982 } 1983 for (s = 0; s < (*labelNew)->numStrata; s++) { 1984 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s]))); 1985 PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices")); 1986 } 1987 PetscCall(PetscFree(points)); 1988 PetscCall(PetscHSetIDestroy(&stratumHash)); 1989 PetscCall(PetscFree(leafStrata)); 1990 PetscCall(PetscFree(strataIdx)); 1991 PetscCall(PetscSectionDestroy(&leafSection)); 1992 PetscFunctionReturn(PETSC_SUCCESS); 1993 } 1994 1995 /*@ 1996 DMLabelGather - Gather all label values from leafs into roots 1997 1998 Collective 1999 2000 Input Parameters: 2001 + label - the `DMLabel` 2002 - sf - the `PetscSF` communication map 2003 2004 Output Parameter: 2005 . labelNew - the new `DMLabel` with localised leaf values 2006 2007 Level: developer 2008 2009 Note: 2010 This is the inverse operation to `DMLabelDistribute()`. 2011 2012 .seealso: `DMLabel`, `DM`, `DMLabelDistribute()` 2013 @*/ 2014 PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew) 2015 { 2016 MPI_Comm comm; 2017 PetscSection rootSection; 2018 PetscSF sfLabel; 2019 PetscSFNode *rootPoints, *leafPoints; 2020 PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset; 2021 const PetscInt *rootDegree, *ilocal; 2022 PetscInt *rootStrata; 2023 const char *lname; 2024 char *name; 2025 PetscMPIInt nameSize; 2026 size_t len = 0; 2027 PetscMPIInt rank, size; 2028 2029 PetscFunctionBegin; 2030 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2031 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 2032 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2033 PetscCall(PetscObjectGetComm((PetscObject)sf, &comm)); 2034 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2035 PetscCallMPI(MPI_Comm_size(comm, &size)); 2036 /* Bcast name */ 2037 if (rank == 0) { 2038 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 2039 PetscCall(PetscStrlen(lname, &len)); 2040 } 2041 PetscCall(PetscMPIIntCast(len, &nameSize)); 2042 PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm)); 2043 PetscCall(PetscMalloc1(nameSize + 1, &name)); 2044 if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1)); 2045 PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm)); 2046 PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew)); 2047 PetscCall(PetscFree(name)); 2048 /* Gather rank/index pairs of leaves into local roots to build 2049 an inverse, multi-rooted SF. Note that this ignores local leaf 2050 indexing due to the use of the multiSF in PetscSFGather. */ 2051 PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL)); 2052 PetscCall(PetscMalloc1(nroots, &leafPoints)); 2053 for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1; 2054 for (p = 0; p < nleaves; p++) { 2055 PetscInt ilp = ilocal ? ilocal[p] : p; 2056 2057 leafPoints[ilp].index = ilp; 2058 leafPoints[ilp].rank = rank; 2059 } 2060 PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree)); 2061 PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree)); 2062 for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p]; 2063 PetscCall(PetscMalloc1(nmultiroots, &rootPoints)); 2064 PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints)); 2065 PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints)); 2066 PetscCall(PetscSFCreate(comm, &sfLabel)); 2067 PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER)); 2068 /* Migrate label over inverted SF to pull stratum values at leaves into roots. */ 2069 PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata)); 2070 /* Rebuild the point strata on the receiver */ 2071 for (p = 0, idx = 0; p < nroots; p++) { 2072 for (d = 0; d < rootDegree[p]; d++) { 2073 PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof)); 2074 PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset)); 2075 for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s])); 2076 } 2077 idx += rootDegree[p]; 2078 } 2079 PetscCall(PetscFree(leafPoints)); 2080 PetscCall(PetscFree(rootStrata)); 2081 PetscCall(PetscSectionDestroy(&rootSection)); 2082 PetscCall(PetscSFDestroy(&sfLabel)); 2083 PetscFunctionReturn(PETSC_SUCCESS); 2084 } 2085 2086 static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[]) 2087 { 2088 const PetscInt *degree; 2089 const PetscInt *points; 2090 PetscInt Nr, r, Nl, l, val, defVal; 2091 2092 PetscFunctionBegin; 2093 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 2094 /* Add in leaves */ 2095 PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL)); 2096 for (l = 0; l < Nl; ++l) { 2097 PetscCall(DMLabelGetValue(label, points[l], &val)); 2098 if (val != defVal) valArray[points[l]] = val; 2099 } 2100 /* Add in shared roots */ 2101 PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree)); 2102 PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree)); 2103 for (r = 0; r < Nr; ++r) { 2104 if (degree[r]) { 2105 PetscCall(DMLabelGetValue(label, r, &val)); 2106 if (val != defVal) valArray[r] = val; 2107 } 2108 } 2109 PetscFunctionReturn(PETSC_SUCCESS); 2110 } 2111 2112 static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx) 2113 { 2114 const PetscInt *degree; 2115 const PetscInt *points; 2116 PetscInt Nr, r, Nl, l, val, defVal; 2117 2118 PetscFunctionBegin; 2119 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 2120 /* Read out leaves */ 2121 PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL)); 2122 for (l = 0; l < Nl; ++l) { 2123 const PetscInt p = points[l]; 2124 const PetscInt cval = valArray[p]; 2125 2126 if (cval != defVal) { 2127 PetscCall(DMLabelGetValue(label, p, &val)); 2128 if (val == defVal) { 2129 PetscCall(DMLabelSetValue(label, p, cval)); 2130 if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx)); 2131 } 2132 } 2133 } 2134 /* Read out shared roots */ 2135 PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree)); 2136 PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree)); 2137 for (r = 0; r < Nr; ++r) { 2138 if (degree[r]) { 2139 const PetscInt cval = valArray[r]; 2140 2141 if (cval != defVal) { 2142 PetscCall(DMLabelGetValue(label, r, &val)); 2143 if (val == defVal) { 2144 PetscCall(DMLabelSetValue(label, r, cval)); 2145 if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx)); 2146 } 2147 } 2148 } 2149 } 2150 PetscFunctionReturn(PETSC_SUCCESS); 2151 } 2152 2153 /*@ 2154 DMLabelPropagateBegin - Setup a cycle of label propagation 2155 2156 Collective 2157 2158 Input Parameters: 2159 + label - The `DMLabel` to propagate across processes 2160 - sf - The `PetscSF` describing parallel layout of the label points 2161 2162 Level: intermediate 2163 2164 .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()` 2165 @*/ 2166 PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf) 2167 { 2168 PetscInt Nr, r, defVal; 2169 PetscMPIInt size; 2170 2171 PetscFunctionBegin; 2172 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2173 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size)); 2174 if (size > 1) { 2175 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 2176 PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL)); 2177 if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray)); 2178 for (r = 0; r < Nr; ++r) label->propArray[r] = defVal; 2179 } 2180 PetscFunctionReturn(PETSC_SUCCESS); 2181 } 2182 2183 /*@ 2184 DMLabelPropagateEnd - Tear down a cycle of label propagation 2185 2186 Collective 2187 2188 Input Parameters: 2189 + label - The `DMLabel` to propagate across processes 2190 - pointSF - The `PetscSF` describing parallel layout of the label points 2191 2192 Level: intermediate 2193 2194 .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()` 2195 @*/ 2196 PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF) 2197 { 2198 PetscFunctionBegin; 2199 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2200 PetscCall(PetscFree(label->propArray)); 2201 label->propArray = NULL; 2202 PetscFunctionReturn(PETSC_SUCCESS); 2203 } 2204 2205 /*@C 2206 DMLabelPropagatePush - Tear down a cycle of label propagation 2207 2208 Collective 2209 2210 Input Parameters: 2211 + label - The `DMLabel` to propagate across processes 2212 . pointSF - The `PetscSF` describing parallel layout of the label points 2213 . markPoint - An optional callback that is called when a point is marked, or `NULL` 2214 - ctx - An optional user context for the callback, or `NULL` 2215 2216 Calling sequence of `markPoint`: 2217 + label - The `DMLabel` 2218 . p - The point being marked 2219 . val - The label value for `p` 2220 - ctx - An optional user context 2221 2222 Level: intermediate 2223 2224 .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()` 2225 @*/ 2226 PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, void *ctx), void *ctx) 2227 { 2228 PetscInt *valArray = label->propArray, Nr; 2229 PetscMPIInt size; 2230 2231 PetscFunctionBegin; 2232 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2233 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size)); 2234 PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL)); 2235 if (size > 1 && Nr >= 0) { 2236 /* Communicate marked edges 2237 The current implementation allocates an array the size of the number of root. We put the label values into the 2238 array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent. 2239 2240 TODO: We could use in-place communication with a different SF 2241 We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has 2242 already been marked. If not, it might have been handled on the process in this round, but we add it anyway. 2243 2244 In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new 2245 values will have 1+0=1 and old values will have 1+1=2. Loop over these, resetting the values to 1, and adding any new 2246 edge to the queue. 2247 */ 2248 PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray)); 2249 PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX)); 2250 PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX)); 2251 PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE)); 2252 PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE)); 2253 PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx)); 2254 } 2255 PetscFunctionReturn(PETSC_SUCCESS); 2256 } 2257 2258 /*@ 2259 DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label 2260 2261 Not Collective 2262 2263 Input Parameter: 2264 . label - the `DMLabel` 2265 2266 Output Parameters: 2267 + section - the section giving offsets for each stratum 2268 - is - An `IS` containing all the label points 2269 2270 Level: developer 2271 2272 .seealso: `DMLabel`, `DM`, `DMLabelDistribute()` 2273 @*/ 2274 PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is) 2275 { 2276 IS vIS; 2277 const PetscInt *values; 2278 PetscInt *points; 2279 PetscInt nV, vS = 0, vE = 0, v, N; 2280 2281 PetscFunctionBegin; 2282 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2283 PetscCall(DMLabelGetNumValues(label, &nV)); 2284 PetscCall(DMLabelGetValueIS(label, &vIS)); 2285 PetscCall(ISGetIndices(vIS, &values)); 2286 if (nV) { 2287 vS = values[0]; 2288 vE = values[0] + 1; 2289 } 2290 for (v = 1; v < nV; ++v) { 2291 vS = PetscMin(vS, values[v]); 2292 vE = PetscMax(vE, values[v] + 1); 2293 } 2294 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section)); 2295 PetscCall(PetscSectionSetChart(*section, vS, vE)); 2296 for (v = 0; v < nV; ++v) { 2297 PetscInt n; 2298 2299 PetscCall(DMLabelGetStratumSize(label, values[v], &n)); 2300 PetscCall(PetscSectionSetDof(*section, values[v], n)); 2301 } 2302 PetscCall(PetscSectionSetUp(*section)); 2303 PetscCall(PetscSectionGetStorageSize(*section, &N)); 2304 PetscCall(PetscMalloc1(N, &points)); 2305 for (v = 0; v < nV; ++v) { 2306 IS is; 2307 const PetscInt *spoints; 2308 PetscInt dof, off, p; 2309 2310 PetscCall(PetscSectionGetDof(*section, values[v], &dof)); 2311 PetscCall(PetscSectionGetOffset(*section, values[v], &off)); 2312 PetscCall(DMLabelGetStratumIS(label, values[v], &is)); 2313 PetscCall(ISGetIndices(is, &spoints)); 2314 for (p = 0; p < dof; ++p) points[off + p] = spoints[p]; 2315 PetscCall(ISRestoreIndices(is, &spoints)); 2316 PetscCall(ISDestroy(&is)); 2317 } 2318 PetscCall(ISRestoreIndices(vIS, &values)); 2319 PetscCall(ISDestroy(&vIS)); 2320 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is)); 2321 PetscFunctionReturn(PETSC_SUCCESS); 2322 } 2323 2324 /*@C 2325 DMLabelRegister - Adds a new label component implementation 2326 2327 Not Collective 2328 2329 Input Parameters: 2330 + name - The name of a new user-defined creation routine 2331 - create_func - The creation routine itself 2332 2333 Notes: 2334 `DMLabelRegister()` may be called multiple times to add several user-defined labels 2335 2336 Example Usage: 2337 .vb 2338 DMLabelRegister("my_label", MyLabelCreate); 2339 .ve 2340 2341 Then, your label type can be chosen with the procedural interface via 2342 .vb 2343 DMLabelCreate(MPI_Comm, DMLabel *); 2344 DMLabelSetType(DMLabel, "my_label"); 2345 .ve 2346 or at runtime via the option 2347 .vb 2348 -dm_label_type my_label 2349 .ve 2350 2351 Level: advanced 2352 2353 .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()` 2354 @*/ 2355 PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel)) 2356 { 2357 PetscFunctionBegin; 2358 PetscCall(DMInitializePackage()); 2359 PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func)); 2360 PetscFunctionReturn(PETSC_SUCCESS); 2361 } 2362 2363 PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel); 2364 PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel); 2365 2366 /*@C 2367 DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package. 2368 2369 Not Collective 2370 2371 Level: advanced 2372 2373 .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()` 2374 @*/ 2375 PetscErrorCode DMLabelRegisterAll(void) 2376 { 2377 PetscFunctionBegin; 2378 if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 2379 DMLabelRegisterAllCalled = PETSC_TRUE; 2380 2381 PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete)); 2382 PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral)); 2383 PetscFunctionReturn(PETSC_SUCCESS); 2384 } 2385 2386 /*@C 2387 DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`. 2388 2389 Level: developer 2390 2391 .seealso: `DMLabel`, `DM`, `PetscInitialize()` 2392 @*/ 2393 PetscErrorCode DMLabelRegisterDestroy(void) 2394 { 2395 PetscFunctionBegin; 2396 PetscCall(PetscFunctionListDestroy(&DMLabelList)); 2397 DMLabelRegisterAllCalled = PETSC_FALSE; 2398 PetscFunctionReturn(PETSC_SUCCESS); 2399 } 2400 2401 /*@ 2402 DMLabelSetType - Sets the particular implementation for a label. 2403 2404 Collective 2405 2406 Input Parameters: 2407 + label - The label 2408 - method - The name of the label type 2409 2410 Options Database Key: 2411 . -dm_label_type <type> - Sets the label type; use -help for a list of available types or see `DMLabelType` 2412 2413 Level: intermediate 2414 2415 .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()` 2416 @*/ 2417 PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method) 2418 { 2419 PetscErrorCode (*r)(DMLabel); 2420 PetscBool match; 2421 2422 PetscFunctionBegin; 2423 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2424 PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match)); 2425 if (match) PetscFunctionReturn(PETSC_SUCCESS); 2426 2427 PetscCall(DMLabelRegisterAll()); 2428 PetscCall(PetscFunctionListFind(DMLabelList, method, &r)); 2429 PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method); 2430 2431 PetscTryTypeMethod(label, destroy); 2432 PetscCall(PetscMemzero(label->ops, sizeof(*label->ops))); 2433 PetscCall(PetscObjectChangeTypeName((PetscObject)label, method)); 2434 PetscCall((*r)(label)); 2435 PetscFunctionReturn(PETSC_SUCCESS); 2436 } 2437 2438 /*@ 2439 DMLabelGetType - Gets the type name (as a string) from the label. 2440 2441 Not Collective 2442 2443 Input Parameter: 2444 . label - The `DMLabel` 2445 2446 Output Parameter: 2447 . type - The `DMLabel` type name 2448 2449 Level: intermediate 2450 2451 .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()` 2452 @*/ 2453 PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type) 2454 { 2455 PetscFunctionBegin; 2456 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2457 PetscAssertPointer(type, 2); 2458 PetscCall(DMLabelRegisterAll()); 2459 *type = ((PetscObject)label)->type_name; 2460 PetscFunctionReturn(PETSC_SUCCESS); 2461 } 2462 2463 static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label) 2464 { 2465 PetscFunctionBegin; 2466 label->ops->view = DMLabelView_Concrete; 2467 label->ops->setup = NULL; 2468 label->ops->duplicate = DMLabelDuplicate_Concrete; 2469 label->ops->getstratumis = DMLabelGetStratumIS_Concrete; 2470 PetscFunctionReturn(PETSC_SUCCESS); 2471 } 2472 2473 PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label) 2474 { 2475 PetscFunctionBegin; 2476 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2477 PetscCall(DMLabelInitialize_Concrete(label)); 2478 PetscFunctionReturn(PETSC_SUCCESS); 2479 } 2480 2481 /*@ 2482 PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using 2483 the local section and an `PetscSF` describing the section point overlap. 2484 2485 Collective 2486 2487 Input Parameters: 2488 + s - The `PetscSection` for the local field layout 2489 . sf - The `PetscSF` describing parallel layout of the section points 2490 . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs 2491 . label - The label specifying the points 2492 - labelValue - The label stratum specifying the points 2493 2494 Output Parameter: 2495 . gsection - The `PetscSection` for the global field layout 2496 2497 Level: developer 2498 2499 Note: 2500 This gives negative sizes and offsets to points not owned by this process 2501 2502 .seealso: `DMLabel`, `DM`, `PetscSectionCreate()` 2503 @*/ 2504 PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection) 2505 { 2506 PetscInt *neg = NULL, *tmpOff = NULL; 2507 PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots; 2508 2509 PetscFunctionBegin; 2510 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2511 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 2512 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4); 2513 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection)); 2514 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 2515 PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd)); 2516 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 2517 if (nroots >= 0) { 2518 PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart); 2519 PetscCall(PetscCalloc1(nroots, &neg)); 2520 if (nroots > pEnd - pStart) { 2521 PetscCall(PetscCalloc1(nroots, &tmpOff)); 2522 } else { 2523 tmpOff = &(*gsection)->atlasDof[-pStart]; 2524 } 2525 } 2526 /* Mark ghost points with negative dof */ 2527 for (p = pStart; p < pEnd; ++p) { 2528 PetscInt value; 2529 2530 PetscCall(DMLabelGetValue(label, p, &value)); 2531 if (value != labelValue) continue; 2532 PetscCall(PetscSectionGetDof(s, p, &dof)); 2533 PetscCall(PetscSectionSetDof(*gsection, p, dof)); 2534 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 2535 if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof)); 2536 if (neg) neg[p] = -(dof + 1); 2537 } 2538 PetscCall(PetscSectionSetUpBC(*gsection)); 2539 if (nroots >= 0) { 2540 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2541 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2542 if (nroots > pEnd - pStart) { 2543 for (p = pStart; p < pEnd; ++p) { 2544 if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p]; 2545 } 2546 } 2547 } 2548 /* Calculate new sizes, get process offset, and calculate point offsets */ 2549 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 2550 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0; 2551 (*gsection)->atlasOff[p] = off; 2552 off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0; 2553 } 2554 PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s))); 2555 globalOff -= off; 2556 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 2557 (*gsection)->atlasOff[p] += globalOff; 2558 if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1); 2559 } 2560 /* Put in negative offsets for ghost points */ 2561 if (nroots >= 0) { 2562 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2563 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2564 if (nroots > pEnd - pStart) { 2565 for (p = pStart; p < pEnd; ++p) { 2566 if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p]; 2567 } 2568 } 2569 } 2570 if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff)); 2571 PetscCall(PetscFree(neg)); 2572 PetscFunctionReturn(PETSC_SUCCESS); 2573 } 2574 2575 typedef struct _n_PetscSectionSym_Label { 2576 DMLabel label; 2577 PetscCopyMode *modes; 2578 PetscInt *sizes; 2579 const PetscInt ***perms; 2580 const PetscScalar ***rots; 2581 PetscInt (*minMaxOrients)[2]; 2582 PetscInt numStrata; /* numStrata is only increasing, functions as a state */ 2583 } PetscSectionSym_Label; 2584 2585 static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym) 2586 { 2587 PetscInt i, j; 2588 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2589 2590 PetscFunctionBegin; 2591 for (i = 0; i <= sl->numStrata; i++) { 2592 if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) { 2593 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2594 if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j])); 2595 if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j])); 2596 } 2597 if (sl->perms[i]) { 2598 const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]]; 2599 2600 PetscCall(PetscFree(perms)); 2601 } 2602 if (sl->rots[i]) { 2603 const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]]; 2604 2605 PetscCall(PetscFree(rots)); 2606 } 2607 } 2608 } 2609 PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients)); 2610 PetscCall(DMLabelDestroy(&sl->label)); 2611 sl->numStrata = 0; 2612 PetscFunctionReturn(PETSC_SUCCESS); 2613 } 2614 2615 static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym) 2616 { 2617 PetscFunctionBegin; 2618 PetscCall(PetscSectionSymLabelReset(sym)); 2619 PetscCall(PetscFree(sym->data)); 2620 PetscFunctionReturn(PETSC_SUCCESS); 2621 } 2622 2623 static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer) 2624 { 2625 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2626 PetscBool isAscii; 2627 DMLabel label = sl->label; 2628 const char *name; 2629 2630 PetscFunctionBegin; 2631 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii)); 2632 if (isAscii) { 2633 PetscInt i, j, k; 2634 PetscViewerFormat format; 2635 2636 PetscCall(PetscViewerGetFormat(viewer, &format)); 2637 if (label) { 2638 PetscCall(PetscViewerGetFormat(viewer, &format)); 2639 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2640 PetscCall(PetscViewerASCIIPushTab(viewer)); 2641 PetscCall(DMLabelView(label, viewer)); 2642 PetscCall(PetscViewerASCIIPopTab(viewer)); 2643 } else { 2644 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2645 PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name)); 2646 } 2647 } else { 2648 PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n")); 2649 } 2650 PetscCall(PetscViewerASCIIPushTab(viewer)); 2651 for (i = 0; i <= sl->numStrata; i++) { 2652 PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue; 2653 2654 if (!(sl->perms[i] || sl->rots[i])) { 2655 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i])); 2656 } else { 2657 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i])); 2658 PetscCall(PetscViewerASCIIPushTab(viewer)); 2659 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1])); 2660 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2661 PetscCall(PetscViewerASCIIPushTab(viewer)); 2662 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2663 if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) { 2664 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j)); 2665 } else { 2666 PetscInt tab; 2667 2668 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j)); 2669 PetscCall(PetscViewerASCIIPushTab(viewer)); 2670 PetscCall(PetscViewerASCIIGetTab(viewer, &tab)); 2671 if (sl->perms[i] && sl->perms[i][j]) { 2672 PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:")); 2673 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 2674 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k])); 2675 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 2676 PetscCall(PetscViewerASCIISetTab(viewer, tab)); 2677 } 2678 if (sl->rots[i] && sl->rots[i][j]) { 2679 PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: ")); 2680 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 2681 #if defined(PETSC_USE_COMPLEX) 2682 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g+i*%+g", (double)PetscRealPart(sl->rots[i][j][k]), (double)PetscImaginaryPart(sl->rots[i][j][k]))); 2683 #else 2684 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k])); 2685 #endif 2686 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 2687 PetscCall(PetscViewerASCIISetTab(viewer, tab)); 2688 } 2689 PetscCall(PetscViewerASCIIPopTab(viewer)); 2690 } 2691 } 2692 PetscCall(PetscViewerASCIIPopTab(viewer)); 2693 } 2694 PetscCall(PetscViewerASCIIPopTab(viewer)); 2695 } 2696 } 2697 PetscCall(PetscViewerASCIIPopTab(viewer)); 2698 } 2699 PetscFunctionReturn(PETSC_SUCCESS); 2700 } 2701 2702 /*@ 2703 PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries 2704 2705 Logically 2706 2707 Input Parameters: 2708 + sym - the section symmetries 2709 - label - the `DMLabel` describing the types of points 2710 2711 Level: developer: 2712 2713 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()` 2714 @*/ 2715 PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label) 2716 { 2717 PetscSectionSym_Label *sl; 2718 2719 PetscFunctionBegin; 2720 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2721 sl = (PetscSectionSym_Label *)sym->data; 2722 if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym)); 2723 if (label) { 2724 sl->label = label; 2725 PetscCall(PetscObjectReference((PetscObject)label)); 2726 PetscCall(DMLabelGetNumValues(label, &sl->numStrata)); 2727 PetscCall(PetscMalloc5(sl->numStrata + 1, &sl->modes, sl->numStrata + 1, &sl->sizes, sl->numStrata + 1, &sl->perms, sl->numStrata + 1, &sl->rots, sl->numStrata + 1, &sl->minMaxOrients)); 2728 PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode))); 2729 PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt))); 2730 PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **))); 2731 PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **))); 2732 PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2]))); 2733 } 2734 PetscFunctionReturn(PETSC_SUCCESS); 2735 } 2736 2737 /*@C 2738 PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum 2739 2740 Logically Collective 2741 2742 Input Parameters: 2743 + sym - the section symmetries 2744 - stratum - the stratum value in the label that we are assigning symmetries for 2745 2746 Output Parameters: 2747 + size - the number of dofs for points in the `stratum` of the label 2748 . minOrient - the smallest orientation for a point in this `stratum` 2749 . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`)) 2750 . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity 2751 - rots - `NULL` if there are no rotations, or (`maxOrient` - `minOrient`) sets of rotations, one for each orientation. A `NULL` set of orientations is the identity 2752 2753 Level: developer 2754 2755 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()` 2756 @*/ 2757 PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots) 2758 { 2759 PetscSectionSym_Label *sl; 2760 const char *name; 2761 PetscInt i; 2762 2763 PetscFunctionBegin; 2764 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2765 sl = (PetscSectionSym_Label *)sym->data; 2766 PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2767 for (i = 0; i <= sl->numStrata; i++) { 2768 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2769 2770 if (stratum == value) break; 2771 } 2772 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2773 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name); 2774 if (size) { 2775 PetscAssertPointer(size, 3); 2776 *size = sl->sizes[i]; 2777 } 2778 if (minOrient) { 2779 PetscAssertPointer(minOrient, 4); 2780 *minOrient = sl->minMaxOrients[i][0]; 2781 } 2782 if (maxOrient) { 2783 PetscAssertPointer(maxOrient, 5); 2784 *maxOrient = sl->minMaxOrients[i][1]; 2785 } 2786 if (perms) { 2787 PetscAssertPointer(perms, 6); 2788 *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]); 2789 } 2790 if (rots) { 2791 PetscAssertPointer(rots, 7); 2792 *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]); 2793 } 2794 PetscFunctionReturn(PETSC_SUCCESS); 2795 } 2796 2797 /*@C 2798 PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum 2799 2800 Logically 2801 2802 Input Parameters: 2803 + sym - the section symmetries 2804 . stratum - the stratum value in the label that we are assigning symmetries for 2805 . size - the number of dofs for points in the `stratum` of the label 2806 . minOrient - the smallest orientation for a point in this `stratum` 2807 . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`)) 2808 . mode - how `sym` should copy the `perms` and `rots` arrays 2809 . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity 2810 - rots - `NULL` if there are no rotations, or (`maxOrient` - `minOrient`) sets of rotations, one for each orientation. A `NULL` set of orientations is the identity 2811 2812 Level: developer 2813 2814 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()` 2815 @*/ 2816 PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots) 2817 { 2818 PetscSectionSym_Label *sl; 2819 const char *name; 2820 PetscInt i, j, k; 2821 2822 PetscFunctionBegin; 2823 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2824 sl = (PetscSectionSym_Label *)sym->data; 2825 PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2826 for (i = 0; i <= sl->numStrata; i++) { 2827 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2828 2829 if (stratum == value) break; 2830 } 2831 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2832 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name); 2833 sl->sizes[i] = size; 2834 sl->modes[i] = mode; 2835 sl->minMaxOrients[i][0] = minOrient; 2836 sl->minMaxOrients[i][1] = maxOrient; 2837 if (mode == PETSC_COPY_VALUES) { 2838 if (perms) { 2839 PetscInt **ownPerms; 2840 2841 PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms)); 2842 for (j = 0; j < maxOrient - minOrient; j++) { 2843 if (perms[j]) { 2844 PetscCall(PetscMalloc1(size, &ownPerms[j])); 2845 for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k]; 2846 } 2847 } 2848 sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient]; 2849 } 2850 if (rots) { 2851 PetscScalar **ownRots; 2852 2853 PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots)); 2854 for (j = 0; j < maxOrient - minOrient; j++) { 2855 if (rots[j]) { 2856 PetscCall(PetscMalloc1(size, &ownRots[j])); 2857 for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k]; 2858 } 2859 } 2860 sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient]; 2861 } 2862 } else { 2863 sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient); 2864 sl->rots[i] = PetscSafePointerPlusOffset(rots, -minOrient); 2865 } 2866 PetscFunctionReturn(PETSC_SUCCESS); 2867 } 2868 2869 static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots) 2870 { 2871 PetscInt i, j, numStrata; 2872 PetscSectionSym_Label *sl; 2873 DMLabel label; 2874 2875 PetscFunctionBegin; 2876 sl = (PetscSectionSym_Label *)sym->data; 2877 numStrata = sl->numStrata; 2878 label = sl->label; 2879 for (i = 0; i < numPoints; i++) { 2880 PetscInt point = points[2 * i]; 2881 PetscInt ornt = points[2 * i + 1]; 2882 2883 for (j = 0; j < numStrata; j++) { 2884 if (label->validIS[j]) { 2885 PetscInt k; 2886 2887 PetscCall(ISLocate(label->points[j], point, &k)); 2888 if (k >= 0) break; 2889 } else { 2890 PetscBool has; 2891 2892 PetscCall(PetscHSetIHas(label->ht[j], point, &has)); 2893 if (has) break; 2894 } 2895 } 2896 PetscCheck(!(sl->minMaxOrients[j][1] > sl->minMaxOrients[j][0]) || !(ornt < sl->minMaxOrients[j][0] || ornt >= sl->minMaxOrients[j][1]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " orientation %" PetscInt_FMT " not in range [%" PetscInt_FMT ", %" PetscInt_FMT ") for stratum %" PetscInt_FMT, point, ornt, sl->minMaxOrients[j][0], sl->minMaxOrients[j][1], 2897 j < numStrata ? label->stratumValues[j] : label->defaultValue); 2898 if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL; 2899 if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL; 2900 } 2901 PetscFunctionReturn(PETSC_SUCCESS); 2902 } 2903 2904 static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym) 2905 { 2906 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data; 2907 IS valIS; 2908 const PetscInt *values; 2909 PetscInt Nv, v; 2910 2911 PetscFunctionBegin; 2912 PetscCall(DMLabelGetNumValues(sl->label, &Nv)); 2913 PetscCall(DMLabelGetValueIS(sl->label, &valIS)); 2914 PetscCall(ISGetIndices(valIS, &values)); 2915 for (v = 0; v < Nv; ++v) { 2916 const PetscInt val = values[v]; 2917 PetscInt size, minOrient, maxOrient; 2918 const PetscInt **perms; 2919 const PetscScalar **rots; 2920 2921 PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots)); 2922 PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots)); 2923 } 2924 PetscCall(ISDestroy(&valIS)); 2925 PetscFunctionReturn(PETSC_SUCCESS); 2926 } 2927 2928 static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym) 2929 { 2930 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2931 DMLabel dlabel; 2932 2933 PetscFunctionBegin; 2934 PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel)); 2935 PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym)); 2936 PetscCall(DMLabelDestroy(&dlabel)); 2937 PetscCall(PetscSectionSymCopy(sym, *dsym)); 2938 PetscFunctionReturn(PETSC_SUCCESS); 2939 } 2940 2941 PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym) 2942 { 2943 PetscSectionSym_Label *sl; 2944 2945 PetscFunctionBegin; 2946 PetscCall(PetscNew(&sl)); 2947 sym->ops->getpoints = PetscSectionSymGetPoints_Label; 2948 sym->ops->distribute = PetscSectionSymDistribute_Label; 2949 sym->ops->copy = PetscSectionSymCopy_Label; 2950 sym->ops->view = PetscSectionSymView_Label; 2951 sym->ops->destroy = PetscSectionSymDestroy_Label; 2952 sym->data = (void *)sl; 2953 PetscFunctionReturn(PETSC_SUCCESS); 2954 } 2955 2956 /*@ 2957 PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label 2958 2959 Collective 2960 2961 Input Parameters: 2962 + comm - the MPI communicator for the new symmetry 2963 - label - the label defining the strata 2964 2965 Output Parameter: 2966 . sym - the section symmetries 2967 2968 Level: developer 2969 2970 .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()` 2971 @*/ 2972 PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym) 2973 { 2974 PetscFunctionBegin; 2975 PetscCall(DMInitializePackage()); 2976 PetscCall(PetscSectionSymCreate(comm, sym)); 2977 PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL)); 2978 PetscCall(PetscSectionSymLabelSetLabel(*sym, label)); 2979 PetscFunctionReturn(PETSC_SUCCESS); 2980 } 2981