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