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