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 . sf - 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 $ PetscErrorCode markPoint(DMLabel label, PetscInt p, PetscInt val, void *ctx); 2090 + label - The `DMLabel` 2091 . p - The point being marked 2092 . pointSF - The label value for `p` 2093 - ctx - An optional user context 2094 2095 Level: intermediate 2096 2097 .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()` 2098 @*/ 2099 PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx) 2100 { 2101 PetscInt *valArray = label->propArray, Nr; 2102 PetscMPIInt size; 2103 2104 PetscFunctionBegin; 2105 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2106 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size)); 2107 PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL)); 2108 if (size > 1 && Nr >= 0) { 2109 /* Communicate marked edges 2110 The current implementation allocates an array the size of the number of root. We put the label values into the 2111 array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent. 2112 2113 TODO: We could use in-place communication with a different SF 2114 We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has 2115 already been marked. If not, it might have been handled on the process in this round, but we add it anyway. 2116 2117 In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new 2118 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 2119 edge to the queue. 2120 */ 2121 PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray)); 2122 PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX)); 2123 PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX)); 2124 PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE)); 2125 PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE)); 2126 PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx)); 2127 } 2128 PetscFunctionReturn(PETSC_SUCCESS); 2129 } 2130 2131 /*@ 2132 DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label 2133 2134 Not Collective 2135 2136 Input Parameter: 2137 . label - the `DMLabel` 2138 2139 Output Parameters: 2140 + section - the section giving offsets for each stratum 2141 - is - An `IS` containing all the label points 2142 2143 Level: developer 2144 2145 .seealso: `DMLabel`, `DM`, `DMLabelDistribute()` 2146 @*/ 2147 PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is) 2148 { 2149 IS vIS; 2150 const PetscInt *values; 2151 PetscInt *points; 2152 PetscInt nV, vS = 0, vE = 0, v, N; 2153 2154 PetscFunctionBegin; 2155 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2156 PetscCall(DMLabelGetNumValues(label, &nV)); 2157 PetscCall(DMLabelGetValueIS(label, &vIS)); 2158 PetscCall(ISGetIndices(vIS, &values)); 2159 if (nV) { 2160 vS = values[0]; 2161 vE = values[0] + 1; 2162 } 2163 for (v = 1; v < nV; ++v) { 2164 vS = PetscMin(vS, values[v]); 2165 vE = PetscMax(vE, values[v] + 1); 2166 } 2167 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section)); 2168 PetscCall(PetscSectionSetChart(*section, vS, vE)); 2169 for (v = 0; v < nV; ++v) { 2170 PetscInt n; 2171 2172 PetscCall(DMLabelGetStratumSize(label, values[v], &n)); 2173 PetscCall(PetscSectionSetDof(*section, values[v], n)); 2174 } 2175 PetscCall(PetscSectionSetUp(*section)); 2176 PetscCall(PetscSectionGetStorageSize(*section, &N)); 2177 PetscCall(PetscMalloc1(N, &points)); 2178 for (v = 0; v < nV; ++v) { 2179 IS is; 2180 const PetscInt *spoints; 2181 PetscInt dof, off, p; 2182 2183 PetscCall(PetscSectionGetDof(*section, values[v], &dof)); 2184 PetscCall(PetscSectionGetOffset(*section, values[v], &off)); 2185 PetscCall(DMLabelGetStratumIS(label, values[v], &is)); 2186 PetscCall(ISGetIndices(is, &spoints)); 2187 for (p = 0; p < dof; ++p) points[off + p] = spoints[p]; 2188 PetscCall(ISRestoreIndices(is, &spoints)); 2189 PetscCall(ISDestroy(&is)); 2190 } 2191 PetscCall(ISRestoreIndices(vIS, &values)); 2192 PetscCall(ISDestroy(&vIS)); 2193 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is)); 2194 PetscFunctionReturn(PETSC_SUCCESS); 2195 } 2196 2197 /*@C 2198 DMLabelRegister - Adds a new label component implementation 2199 2200 Not Collective 2201 2202 Input Parameters: 2203 + name - The name of a new user-defined creation routine 2204 - create_func - The creation routine itself 2205 2206 Notes: 2207 `DMLabelRegister()` may be called multiple times to add several user-defined labels 2208 2209 Example Usage: 2210 .vb 2211 DMLabelRegister("my_label", MyLabelCreate); 2212 .ve 2213 2214 Then, your label type can be chosen with the procedural interface via 2215 .vb 2216 DMLabelCreate(MPI_Comm, DMLabel *); 2217 DMLabelSetType(DMLabel, "my_label"); 2218 .ve 2219 or at runtime via the option 2220 .vb 2221 -dm_label_type my_label 2222 .ve 2223 2224 Level: advanced 2225 2226 .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()` 2227 @*/ 2228 PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel)) 2229 { 2230 PetscFunctionBegin; 2231 PetscCall(DMInitializePackage()); 2232 PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func)); 2233 PetscFunctionReturn(PETSC_SUCCESS); 2234 } 2235 2236 PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel); 2237 PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel); 2238 2239 /*@C 2240 DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package. 2241 2242 Not Collective 2243 2244 Level: advanced 2245 2246 .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()` 2247 @*/ 2248 PetscErrorCode DMLabelRegisterAll(void) 2249 { 2250 PetscFunctionBegin; 2251 if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 2252 DMLabelRegisterAllCalled = PETSC_TRUE; 2253 2254 PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete)); 2255 PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral)); 2256 PetscFunctionReturn(PETSC_SUCCESS); 2257 } 2258 2259 /*@C 2260 DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`. 2261 2262 Level: developer 2263 2264 .seealso: `DMLabel`, `DM`, `PetscInitialize()` 2265 @*/ 2266 PetscErrorCode DMLabelRegisterDestroy(void) 2267 { 2268 PetscFunctionBegin; 2269 PetscCall(PetscFunctionListDestroy(&DMLabelList)); 2270 DMLabelRegisterAllCalled = PETSC_FALSE; 2271 PetscFunctionReturn(PETSC_SUCCESS); 2272 } 2273 2274 /*@C 2275 DMLabelSetType - Sets the particular implementation for a label. 2276 2277 Collective 2278 2279 Input Parameters: 2280 + label - The label 2281 - method - The name of the label type 2282 2283 Options Database Key: 2284 . -dm_label_type <type> - Sets the label type; use -help for a list of available types or see `DMLabelType` 2285 2286 Level: intermediate 2287 2288 .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()` 2289 @*/ 2290 PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method) 2291 { 2292 PetscErrorCode (*r)(DMLabel); 2293 PetscBool match; 2294 2295 PetscFunctionBegin; 2296 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2297 PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match)); 2298 if (match) PetscFunctionReturn(PETSC_SUCCESS); 2299 2300 PetscCall(DMLabelRegisterAll()); 2301 PetscCall(PetscFunctionListFind(DMLabelList, method, &r)); 2302 PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method); 2303 2304 PetscTryTypeMethod(label, destroy); 2305 PetscCall(PetscMemzero(label->ops, sizeof(*label->ops))); 2306 PetscCall(PetscObjectChangeTypeName((PetscObject)label, method)); 2307 PetscCall((*r)(label)); 2308 PetscFunctionReturn(PETSC_SUCCESS); 2309 } 2310 2311 /*@C 2312 DMLabelGetType - Gets the type name (as a string) from the label. 2313 2314 Not Collective 2315 2316 Input Parameter: 2317 . label - The `DMLabel` 2318 2319 Output Parameter: 2320 . type - The `DMLabel` type name 2321 2322 Level: intermediate 2323 2324 .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()` 2325 @*/ 2326 PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type) 2327 { 2328 PetscFunctionBegin; 2329 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2330 PetscAssertPointer(type, 2); 2331 PetscCall(DMLabelRegisterAll()); 2332 *type = ((PetscObject)label)->type_name; 2333 PetscFunctionReturn(PETSC_SUCCESS); 2334 } 2335 2336 static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label) 2337 { 2338 PetscFunctionBegin; 2339 label->ops->view = DMLabelView_Concrete; 2340 label->ops->setup = NULL; 2341 label->ops->duplicate = DMLabelDuplicate_Concrete; 2342 label->ops->getstratumis = DMLabelGetStratumIS_Concrete; 2343 PetscFunctionReturn(PETSC_SUCCESS); 2344 } 2345 2346 PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label) 2347 { 2348 PetscFunctionBegin; 2349 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2350 PetscCall(DMLabelInitialize_Concrete(label)); 2351 PetscFunctionReturn(PETSC_SUCCESS); 2352 } 2353 2354 /*@ 2355 PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using 2356 the local section and an `PetscSF` describing the section point overlap. 2357 2358 Collective 2359 2360 Input Parameters: 2361 + s - The `PetscSection` for the local field layout 2362 . sf - The `PetscSF` describing parallel layout of the section points 2363 . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs 2364 . label - The label specifying the points 2365 - labelValue - The label stratum specifying the points 2366 2367 Output Parameter: 2368 . gsection - The `PetscSection` for the global field layout 2369 2370 Level: developer 2371 2372 Note: 2373 This gives negative sizes and offsets to points not owned by this process 2374 2375 .seealso: `DMLabel`, `DM`, `PetscSectionCreate()` 2376 @*/ 2377 PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection) 2378 { 2379 PetscInt *neg = NULL, *tmpOff = NULL; 2380 PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots; 2381 2382 PetscFunctionBegin; 2383 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2384 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 2385 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4); 2386 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection)); 2387 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 2388 PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd)); 2389 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 2390 if (nroots >= 0) { 2391 PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart); 2392 PetscCall(PetscCalloc1(nroots, &neg)); 2393 if (nroots > pEnd - pStart) { 2394 PetscCall(PetscCalloc1(nroots, &tmpOff)); 2395 } else { 2396 tmpOff = &(*gsection)->atlasDof[-pStart]; 2397 } 2398 } 2399 /* Mark ghost points with negative dof */ 2400 for (p = pStart; p < pEnd; ++p) { 2401 PetscInt value; 2402 2403 PetscCall(DMLabelGetValue(label, p, &value)); 2404 if (value != labelValue) continue; 2405 PetscCall(PetscSectionGetDof(s, p, &dof)); 2406 PetscCall(PetscSectionSetDof(*gsection, p, dof)); 2407 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 2408 if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof)); 2409 if (neg) neg[p] = -(dof + 1); 2410 } 2411 PetscCall(PetscSectionSetUpBC(*gsection)); 2412 if (nroots >= 0) { 2413 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2414 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2415 if (nroots > pEnd - pStart) { 2416 for (p = pStart; p < pEnd; ++p) { 2417 if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p]; 2418 } 2419 } 2420 } 2421 /* Calculate new sizes, get process offset, and calculate point offsets */ 2422 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 2423 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0; 2424 (*gsection)->atlasOff[p] = off; 2425 off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0; 2426 } 2427 PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s))); 2428 globalOff -= off; 2429 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 2430 (*gsection)->atlasOff[p] += globalOff; 2431 if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1); 2432 } 2433 /* Put in negative offsets for ghost points */ 2434 if (nroots >= 0) { 2435 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2436 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2437 if (nroots > pEnd - pStart) { 2438 for (p = pStart; p < pEnd; ++p) { 2439 if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p]; 2440 } 2441 } 2442 } 2443 if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff)); 2444 PetscCall(PetscFree(neg)); 2445 PetscFunctionReturn(PETSC_SUCCESS); 2446 } 2447 2448 typedef struct _n_PetscSectionSym_Label { 2449 DMLabel label; 2450 PetscCopyMode *modes; 2451 PetscInt *sizes; 2452 const PetscInt ***perms; 2453 const PetscScalar ***rots; 2454 PetscInt (*minMaxOrients)[2]; 2455 PetscInt numStrata; /* numStrata is only increasing, functions as a state */ 2456 } PetscSectionSym_Label; 2457 2458 static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym) 2459 { 2460 PetscInt i, j; 2461 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2462 2463 PetscFunctionBegin; 2464 for (i = 0; i <= sl->numStrata; i++) { 2465 if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) { 2466 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2467 if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j])); 2468 if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j])); 2469 } 2470 if (sl->perms[i]) { 2471 const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]]; 2472 2473 PetscCall(PetscFree(perms)); 2474 } 2475 if (sl->rots[i]) { 2476 const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]]; 2477 2478 PetscCall(PetscFree(rots)); 2479 } 2480 } 2481 } 2482 PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients)); 2483 PetscCall(DMLabelDestroy(&sl->label)); 2484 sl->numStrata = 0; 2485 PetscFunctionReturn(PETSC_SUCCESS); 2486 } 2487 2488 static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym) 2489 { 2490 PetscFunctionBegin; 2491 PetscCall(PetscSectionSymLabelReset(sym)); 2492 PetscCall(PetscFree(sym->data)); 2493 PetscFunctionReturn(PETSC_SUCCESS); 2494 } 2495 2496 static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer) 2497 { 2498 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2499 PetscBool isAscii; 2500 DMLabel label = sl->label; 2501 const char *name; 2502 2503 PetscFunctionBegin; 2504 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii)); 2505 if (isAscii) { 2506 PetscInt i, j, k; 2507 PetscViewerFormat format; 2508 2509 PetscCall(PetscViewerGetFormat(viewer, &format)); 2510 if (label) { 2511 PetscCall(PetscViewerGetFormat(viewer, &format)); 2512 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2513 PetscCall(PetscViewerASCIIPushTab(viewer)); 2514 PetscCall(DMLabelView(label, viewer)); 2515 PetscCall(PetscViewerASCIIPopTab(viewer)); 2516 } else { 2517 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2518 PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name)); 2519 } 2520 } else { 2521 PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n")); 2522 } 2523 PetscCall(PetscViewerASCIIPushTab(viewer)); 2524 for (i = 0; i <= sl->numStrata; i++) { 2525 PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue; 2526 2527 if (!(sl->perms[i] || sl->rots[i])) { 2528 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i])); 2529 } else { 2530 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i])); 2531 PetscCall(PetscViewerASCIIPushTab(viewer)); 2532 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1])); 2533 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2534 PetscCall(PetscViewerASCIIPushTab(viewer)); 2535 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2536 if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) { 2537 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j)); 2538 } else { 2539 PetscInt tab; 2540 2541 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j)); 2542 PetscCall(PetscViewerASCIIPushTab(viewer)); 2543 PetscCall(PetscViewerASCIIGetTab(viewer, &tab)); 2544 if (sl->perms[i] && sl->perms[i][j]) { 2545 PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:")); 2546 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 2547 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k])); 2548 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 2549 PetscCall(PetscViewerASCIISetTab(viewer, tab)); 2550 } 2551 if (sl->rots[i] && sl->rots[i][j]) { 2552 PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: ")); 2553 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 2554 #if defined(PETSC_USE_COMPLEX) 2555 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]))); 2556 #else 2557 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k])); 2558 #endif 2559 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 2560 PetscCall(PetscViewerASCIISetTab(viewer, tab)); 2561 } 2562 PetscCall(PetscViewerASCIIPopTab(viewer)); 2563 } 2564 } 2565 PetscCall(PetscViewerASCIIPopTab(viewer)); 2566 } 2567 PetscCall(PetscViewerASCIIPopTab(viewer)); 2568 } 2569 } 2570 PetscCall(PetscViewerASCIIPopTab(viewer)); 2571 } 2572 PetscFunctionReturn(PETSC_SUCCESS); 2573 } 2574 2575 /*@ 2576 PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries 2577 2578 Logically 2579 2580 Input Parameters: 2581 + sym - the section symmetries 2582 - label - the `DMLabel` describing the types of points 2583 2584 Level: developer: 2585 2586 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()` 2587 @*/ 2588 PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label) 2589 { 2590 PetscSectionSym_Label *sl; 2591 2592 PetscFunctionBegin; 2593 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2594 sl = (PetscSectionSym_Label *)sym->data; 2595 if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym)); 2596 if (label) { 2597 sl->label = label; 2598 PetscCall(PetscObjectReference((PetscObject)label)); 2599 PetscCall(DMLabelGetNumValues(label, &sl->numStrata)); 2600 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)); 2601 PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode))); 2602 PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt))); 2603 PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **))); 2604 PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **))); 2605 PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2]))); 2606 } 2607 PetscFunctionReturn(PETSC_SUCCESS); 2608 } 2609 2610 /*@C 2611 PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum 2612 2613 Logically Collective 2614 2615 Input Parameters: 2616 + sym - the section symmetries 2617 - stratum - the stratum value in the label that we are assigning symmetries for 2618 2619 Output Parameters: 2620 + size - the number of dofs for points in the `stratum` of the label 2621 . minOrient - the smallest orientation for a point in this `stratum` 2622 . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`)) 2623 . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity 2624 - 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 2625 2626 Level: developer 2627 2628 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()` 2629 @*/ 2630 PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots) 2631 { 2632 PetscSectionSym_Label *sl; 2633 const char *name; 2634 PetscInt i; 2635 2636 PetscFunctionBegin; 2637 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2638 sl = (PetscSectionSym_Label *)sym->data; 2639 PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2640 for (i = 0; i <= sl->numStrata; i++) { 2641 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2642 2643 if (stratum == value) break; 2644 } 2645 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2646 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name); 2647 if (size) { 2648 PetscAssertPointer(size, 3); 2649 *size = sl->sizes[i]; 2650 } 2651 if (minOrient) { 2652 PetscAssertPointer(minOrient, 4); 2653 *minOrient = sl->minMaxOrients[i][0]; 2654 } 2655 if (maxOrient) { 2656 PetscAssertPointer(maxOrient, 5); 2657 *maxOrient = sl->minMaxOrients[i][1]; 2658 } 2659 if (perms) { 2660 PetscAssertPointer(perms, 6); 2661 *perms = sl->perms[i] ? &sl->perms[i][sl->minMaxOrients[i][0]] : NULL; 2662 } 2663 if (rots) { 2664 PetscAssertPointer(rots, 7); 2665 *rots = sl->rots[i] ? &sl->rots[i][sl->minMaxOrients[i][0]] : NULL; 2666 } 2667 PetscFunctionReturn(PETSC_SUCCESS); 2668 } 2669 2670 /*@C 2671 PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum 2672 2673 Logically 2674 2675 Input Parameters: 2676 + sym - the section symmetries 2677 . stratum - the stratum value in the label that we are assigning symmetries for 2678 . size - the number of dofs for points in the `stratum` of the label 2679 . minOrient - the smallest orientation for a point in this `stratum` 2680 . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`)) 2681 . mode - how `sym` should copy the `perms` and `rots` arrays 2682 . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity 2683 - 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 2684 2685 Level: developer 2686 2687 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()` 2688 @*/ 2689 PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots) 2690 { 2691 PetscSectionSym_Label *sl; 2692 const char *name; 2693 PetscInt i, j, k; 2694 2695 PetscFunctionBegin; 2696 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2697 sl = (PetscSectionSym_Label *)sym->data; 2698 PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2699 for (i = 0; i <= sl->numStrata; i++) { 2700 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2701 2702 if (stratum == value) break; 2703 } 2704 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2705 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name); 2706 sl->sizes[i] = size; 2707 sl->modes[i] = mode; 2708 sl->minMaxOrients[i][0] = minOrient; 2709 sl->minMaxOrients[i][1] = maxOrient; 2710 if (mode == PETSC_COPY_VALUES) { 2711 if (perms) { 2712 PetscInt **ownPerms; 2713 2714 PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms)); 2715 for (j = 0; j < maxOrient - minOrient; j++) { 2716 if (perms[j]) { 2717 PetscCall(PetscMalloc1(size, &ownPerms[j])); 2718 for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k]; 2719 } 2720 } 2721 sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient]; 2722 } 2723 if (rots) { 2724 PetscScalar **ownRots; 2725 2726 PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots)); 2727 for (j = 0; j < maxOrient - minOrient; j++) { 2728 if (rots[j]) { 2729 PetscCall(PetscMalloc1(size, &ownRots[j])); 2730 for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k]; 2731 } 2732 } 2733 sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient]; 2734 } 2735 } else { 2736 sl->perms[i] = perms ? &perms[-minOrient] : NULL; 2737 sl->rots[i] = rots ? &rots[-minOrient] : NULL; 2738 } 2739 PetscFunctionReturn(PETSC_SUCCESS); 2740 } 2741 2742 static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots) 2743 { 2744 PetscInt i, j, numStrata; 2745 PetscSectionSym_Label *sl; 2746 DMLabel label; 2747 2748 PetscFunctionBegin; 2749 sl = (PetscSectionSym_Label *)sym->data; 2750 numStrata = sl->numStrata; 2751 label = sl->label; 2752 for (i = 0; i < numPoints; i++) { 2753 PetscInt point = points[2 * i]; 2754 PetscInt ornt = points[2 * i + 1]; 2755 2756 for (j = 0; j < numStrata; j++) { 2757 if (label->validIS[j]) { 2758 PetscInt k; 2759 2760 PetscCall(ISLocate(label->points[j], point, &k)); 2761 if (k >= 0) break; 2762 } else { 2763 PetscBool has; 2764 2765 PetscCall(PetscHSetIHas(label->ht[j], point, &has)); 2766 if (has) break; 2767 } 2768 } 2769 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], 2770 j < numStrata ? label->stratumValues[j] : label->defaultValue); 2771 if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL; 2772 if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL; 2773 } 2774 PetscFunctionReturn(PETSC_SUCCESS); 2775 } 2776 2777 static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym) 2778 { 2779 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data; 2780 IS valIS; 2781 const PetscInt *values; 2782 PetscInt Nv, v; 2783 2784 PetscFunctionBegin; 2785 PetscCall(DMLabelGetNumValues(sl->label, &Nv)); 2786 PetscCall(DMLabelGetValueIS(sl->label, &valIS)); 2787 PetscCall(ISGetIndices(valIS, &values)); 2788 for (v = 0; v < Nv; ++v) { 2789 const PetscInt val = values[v]; 2790 PetscInt size, minOrient, maxOrient; 2791 const PetscInt **perms; 2792 const PetscScalar **rots; 2793 2794 PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots)); 2795 PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots)); 2796 } 2797 PetscCall(ISDestroy(&valIS)); 2798 PetscFunctionReturn(PETSC_SUCCESS); 2799 } 2800 2801 static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym) 2802 { 2803 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2804 DMLabel dlabel; 2805 2806 PetscFunctionBegin; 2807 PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel)); 2808 PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym)); 2809 PetscCall(DMLabelDestroy(&dlabel)); 2810 PetscCall(PetscSectionSymCopy(sym, *dsym)); 2811 PetscFunctionReturn(PETSC_SUCCESS); 2812 } 2813 2814 PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym) 2815 { 2816 PetscSectionSym_Label *sl; 2817 2818 PetscFunctionBegin; 2819 PetscCall(PetscNew(&sl)); 2820 sym->ops->getpoints = PetscSectionSymGetPoints_Label; 2821 sym->ops->distribute = PetscSectionSymDistribute_Label; 2822 sym->ops->copy = PetscSectionSymCopy_Label; 2823 sym->ops->view = PetscSectionSymView_Label; 2824 sym->ops->destroy = PetscSectionSymDestroy_Label; 2825 sym->data = (void *)sl; 2826 PetscFunctionReturn(PETSC_SUCCESS); 2827 } 2828 2829 /*@ 2830 PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label 2831 2832 Collective 2833 2834 Input Parameters: 2835 + comm - the MPI communicator for the new symmetry 2836 - label - the label defining the strata 2837 2838 Output Parameter: 2839 . sym - the section symmetries 2840 2841 Level: developer 2842 2843 .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()` 2844 @*/ 2845 PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym) 2846 { 2847 PetscFunctionBegin; 2848 PetscCall(DMInitializePackage()); 2849 PetscCall(PetscSectionSymCreate(comm, sym)); 2850 PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL)); 2851 PetscCall(PetscSectionSymLabelSetLabel(*sym, label)); 2852 PetscFunctionReturn(PETSC_SUCCESS); 2853 } 2854