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