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