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