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 /*@ 1726 DMLabelPermuteValues - Permute the values in a label 1727 1728 Not collective 1729 1730 Input Parameters: 1731 + label - the `DMLabel` 1732 - permutation - the value permutation, permutation[old value] = new value 1733 1734 Output Parameter: 1735 . label - the `DMLabel` now with permuted values 1736 1737 Note: 1738 The modification is done in-place 1739 1740 Level: intermediate 1741 1742 .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1743 @*/ 1744 PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation) 1745 { 1746 PetscInt Nv, Np; 1747 1748 PetscFunctionBegin; 1749 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1750 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2); 1751 PetscCall(DMLabelGetNumValues(label, &Nv)); 1752 PetscCall(ISGetLocalSize(permutation, &Np)); 1753 PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv); 1754 if (PetscDefined(USE_DEBUG)) { 1755 PetscBool flg; 1756 PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg)); 1757 PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation"); 1758 } 1759 PetscCall(DMLabelRewriteValues(label, permutation)); 1760 PetscFunctionReturn(PETSC_SUCCESS); 1761 } 1762 1763 /*@ 1764 DMLabelRewriteValues - Permute the values in a label, but some may be omitted 1765 1766 Not collective 1767 1768 Input Parameters: 1769 + label - the `DMLabel` 1770 - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted 1771 1772 Output Parameter: 1773 . label - the `DMLabel` now with permuted values 1774 1775 Note: 1776 The modification is done in-place 1777 1778 Level: intermediate 1779 1780 .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1781 @*/ 1782 PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation) 1783 { 1784 const PetscInt *perm; 1785 PetscInt Nv, Np; 1786 1787 PetscFunctionBegin; 1788 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1789 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2); 1790 PetscCall(DMLabelMakeAllValid_Private(label)); 1791 PetscCall(DMLabelGetNumValues(label, &Nv)); 1792 PetscCall(ISGetLocalSize(permutation, &Np)); 1793 PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv); 1794 PetscCall(ISGetIndices(permutation, &perm)); 1795 for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]]; 1796 PetscCall(ISRestoreIndices(permutation, &perm)); 1797 PetscFunctionReturn(PETSC_SUCCESS); 1798 } 1799 1800 static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata) 1801 { 1802 MPI_Comm comm; 1803 PetscInt s, l, nroots, nleaves, offset, size; 1804 PetscInt *remoteOffsets, *rootStrata, *rootIdx; 1805 PetscSection rootSection; 1806 PetscSF labelSF; 1807 1808 PetscFunctionBegin; 1809 if (label) PetscCall(DMLabelMakeAllValid_Private(label)); 1810 PetscCall(PetscObjectGetComm((PetscObject)sf, &comm)); 1811 /* Build a section of stratum values per point, generate the according SF 1812 and distribute point-wise stratum values to leaves. */ 1813 PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL)); 1814 PetscCall(PetscSectionCreate(comm, &rootSection)); 1815 PetscCall(PetscSectionSetChart(rootSection, 0, nroots)); 1816 if (label) { 1817 for (s = 0; s < label->numStrata; ++s) { 1818 const PetscInt *points; 1819 1820 PetscCall(ISGetIndices(label->points[s], &points)); 1821 for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1)); 1822 PetscCall(ISRestoreIndices(label->points[s], &points)); 1823 } 1824 } 1825 PetscCall(PetscSectionSetUp(rootSection)); 1826 /* Create a point-wise array of stratum values */ 1827 PetscCall(PetscSectionGetStorageSize(rootSection, &size)); 1828 PetscCall(PetscMalloc1(size, &rootStrata)); 1829 PetscCall(PetscCalloc1(nroots, &rootIdx)); 1830 if (label) { 1831 for (s = 0; s < label->numStrata; ++s) { 1832 const PetscInt *points; 1833 1834 PetscCall(ISGetIndices(label->points[s], &points)); 1835 for (l = 0; l < label->stratumSizes[s]; l++) { 1836 const PetscInt p = points[l]; 1837 PetscCall(PetscSectionGetOffset(rootSection, p, &offset)); 1838 rootStrata[offset + rootIdx[p]++] = label->stratumValues[s]; 1839 } 1840 PetscCall(ISRestoreIndices(label->points[s], &points)); 1841 } 1842 } 1843 /* Build SF that maps label points to remote processes */ 1844 PetscCall(PetscSectionCreate(comm, leafSection)); 1845 PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection)); 1846 PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF)); 1847 PetscCall(PetscFree(remoteOffsets)); 1848 /* Send the strata for each point over the derived SF */ 1849 PetscCall(PetscSectionGetStorageSize(*leafSection, &size)); 1850 PetscCall(PetscMalloc1(size, leafStrata)); 1851 PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE)); 1852 PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE)); 1853 /* Clean up */ 1854 PetscCall(PetscFree(rootStrata)); 1855 PetscCall(PetscFree(rootIdx)); 1856 PetscCall(PetscSectionDestroy(&rootSection)); 1857 PetscCall(PetscSFDestroy(&labelSF)); 1858 PetscFunctionReturn(PETSC_SUCCESS); 1859 } 1860 1861 /*@ 1862 DMLabelDistribute - Create a new label pushed forward over the `PetscSF` 1863 1864 Collective 1865 1866 Input Parameters: 1867 + label - the `DMLabel` 1868 - sf - the map from old to new distribution 1869 1870 Output Parameter: 1871 . labelNew - the new redistributed label 1872 1873 Level: intermediate 1874 1875 .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()` 1876 @*/ 1877 PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew) 1878 { 1879 MPI_Comm comm; 1880 PetscSection leafSection; 1881 PetscInt p, pStart, pEnd, s, size, dof, offset, stratum; 1882 PetscInt *leafStrata, *strataIdx; 1883 PetscInt **points; 1884 const char *lname = NULL; 1885 char *name; 1886 PetscMPIInt nameSize; 1887 PetscHSetI stratumHash; 1888 size_t len = 0; 1889 PetscMPIInt rank; 1890 1891 PetscFunctionBegin; 1892 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 1893 if (label) { 1894 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1895 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 1896 PetscCall(DMLabelMakeAllValid_Private(label)); 1897 } 1898 PetscCall(PetscObjectGetComm((PetscObject)sf, &comm)); 1899 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1900 /* Bcast name */ 1901 if (rank == 0) { 1902 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1903 PetscCall(PetscStrlen(lname, &len)); 1904 } 1905 PetscCall(PetscMPIIntCast(len, &nameSize)); 1906 PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm)); 1907 PetscCall(PetscMalloc1(nameSize + 1, &name)); 1908 if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1)); 1909 PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm)); 1910 PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew)); 1911 PetscCall(PetscFree(name)); 1912 /* Bcast defaultValue */ 1913 if (rank == 0) (*labelNew)->defaultValue = label->defaultValue; 1914 PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm)); 1915 /* Distribute stratum values over the SF and get the point mapping on the receiver */ 1916 PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata)); 1917 /* Determine received stratum values and initialise new label*/ 1918 PetscCall(PetscHSetICreate(&stratumHash)); 1919 PetscCall(PetscSectionGetStorageSize(leafSection, &size)); 1920 for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p])); 1921 PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata)); 1922 PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS)); 1923 for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE; 1924 PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues)); 1925 /* Turn leafStrata into indices rather than stratum values */ 1926 offset = 0; 1927 PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues)); 1928 PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues)); 1929 for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s)); 1930 for (p = 0; p < size; ++p) { 1931 for (s = 0; s < (*labelNew)->numStrata; ++s) { 1932 if (leafStrata[p] == (*labelNew)->stratumValues[s]) { 1933 leafStrata[p] = s; 1934 break; 1935 } 1936 } 1937 } 1938 /* Rebuild the point strata on the receiver */ 1939 PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes)); 1940 PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd)); 1941 for (p = pStart; p < pEnd; p++) { 1942 PetscCall(PetscSectionGetDof(leafSection, p, &dof)); 1943 PetscCall(PetscSectionGetOffset(leafSection, p, &offset)); 1944 for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++; 1945 } 1946 PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht)); 1947 PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points)); 1948 PetscCall(PetscCalloc1((*labelNew)->numStrata, &points)); 1949 for (s = 0; s < (*labelNew)->numStrata; ++s) { 1950 PetscCall(PetscHSetICreate(&(*labelNew)->ht[s])); 1951 PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s])); 1952 } 1953 /* Insert points into new strata */ 1954 PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx)); 1955 PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd)); 1956 for (p = pStart; p < pEnd; p++) { 1957 PetscCall(PetscSectionGetDof(leafSection, p, &dof)); 1958 PetscCall(PetscSectionGetOffset(leafSection, p, &offset)); 1959 for (s = 0; s < dof; s++) { 1960 stratum = leafStrata[offset + s]; 1961 points[stratum][strataIdx[stratum]++] = p; 1962 } 1963 } 1964 for (s = 0; s < (*labelNew)->numStrata; s++) { 1965 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s]))); 1966 PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices")); 1967 } 1968 PetscCall(PetscFree(points)); 1969 PetscCall(PetscHSetIDestroy(&stratumHash)); 1970 PetscCall(PetscFree(leafStrata)); 1971 PetscCall(PetscFree(strataIdx)); 1972 PetscCall(PetscSectionDestroy(&leafSection)); 1973 PetscFunctionReturn(PETSC_SUCCESS); 1974 } 1975 1976 /*@ 1977 DMLabelGather - Gather all label values from leafs into roots 1978 1979 Collective 1980 1981 Input Parameters: 1982 + label - the `DMLabel` 1983 - sf - the `PetscSF` communication map 1984 1985 Output Parameter: 1986 . labelNew - the new `DMLabel` with localised leaf values 1987 1988 Level: developer 1989 1990 Note: 1991 This is the inverse operation to `DMLabelDistribute()`. 1992 1993 .seealso: `DMLabel`, `DM`, `DMLabelDistribute()` 1994 @*/ 1995 PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew) 1996 { 1997 MPI_Comm comm; 1998 PetscSection rootSection; 1999 PetscSF sfLabel; 2000 PetscSFNode *rootPoints, *leafPoints; 2001 PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset; 2002 const PetscInt *rootDegree, *ilocal; 2003 PetscInt *rootStrata; 2004 const char *lname; 2005 char *name; 2006 PetscMPIInt nameSize; 2007 size_t len = 0; 2008 PetscMPIInt rank, size; 2009 2010 PetscFunctionBegin; 2011 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2012 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 2013 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2014 PetscCall(PetscObjectGetComm((PetscObject)sf, &comm)); 2015 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2016 PetscCallMPI(MPI_Comm_size(comm, &size)); 2017 /* Bcast name */ 2018 if (rank == 0) { 2019 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 2020 PetscCall(PetscStrlen(lname, &len)); 2021 } 2022 PetscCall(PetscMPIIntCast(len, &nameSize)); 2023 PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm)); 2024 PetscCall(PetscMalloc1(nameSize + 1, &name)); 2025 if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1)); 2026 PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm)); 2027 PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew)); 2028 PetscCall(PetscFree(name)); 2029 /* Gather rank/index pairs of leaves into local roots to build 2030 an inverse, multi-rooted SF. Note that this ignores local leaf 2031 indexing due to the use of the multiSF in PetscSFGather. */ 2032 PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL)); 2033 PetscCall(PetscMalloc1(nroots, &leafPoints)); 2034 for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1; 2035 for (p = 0; p < nleaves; p++) { 2036 PetscInt ilp = ilocal ? ilocal[p] : p; 2037 2038 leafPoints[ilp].index = ilp; 2039 leafPoints[ilp].rank = rank; 2040 } 2041 PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree)); 2042 PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree)); 2043 for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p]; 2044 PetscCall(PetscMalloc1(nmultiroots, &rootPoints)); 2045 PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints)); 2046 PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints)); 2047 PetscCall(PetscSFCreate(comm, &sfLabel)); 2048 PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER)); 2049 /* Migrate label over inverted SF to pull stratum values at leaves into roots. */ 2050 PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata)); 2051 /* Rebuild the point strata on the receiver */ 2052 for (p = 0, idx = 0; p < nroots; p++) { 2053 for (d = 0; d < rootDegree[p]; d++) { 2054 PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof)); 2055 PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset)); 2056 for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s])); 2057 } 2058 idx += rootDegree[p]; 2059 } 2060 PetscCall(PetscFree(leafPoints)); 2061 PetscCall(PetscFree(rootStrata)); 2062 PetscCall(PetscSectionDestroy(&rootSection)); 2063 PetscCall(PetscSFDestroy(&sfLabel)); 2064 PetscFunctionReturn(PETSC_SUCCESS); 2065 } 2066 2067 static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[]) 2068 { 2069 const PetscInt *degree; 2070 const PetscInt *points; 2071 PetscInt Nr, r, Nl, l, val, defVal; 2072 2073 PetscFunctionBegin; 2074 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 2075 /* Add in leaves */ 2076 PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL)); 2077 for (l = 0; l < Nl; ++l) { 2078 PetscCall(DMLabelGetValue(label, points[l], &val)); 2079 if (val != defVal) valArray[points[l]] = val; 2080 } 2081 /* Add in shared roots */ 2082 PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree)); 2083 PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree)); 2084 for (r = 0; r < Nr; ++r) { 2085 if (degree[r]) { 2086 PetscCall(DMLabelGetValue(label, r, &val)); 2087 if (val != defVal) valArray[r] = val; 2088 } 2089 } 2090 PetscFunctionReturn(PETSC_SUCCESS); 2091 } 2092 2093 static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx) 2094 { 2095 const PetscInt *degree; 2096 const PetscInt *points; 2097 PetscInt Nr, r, Nl, l, val, defVal; 2098 2099 PetscFunctionBegin; 2100 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 2101 /* Read out leaves */ 2102 PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL)); 2103 for (l = 0; l < Nl; ++l) { 2104 const PetscInt p = points[l]; 2105 const PetscInt cval = valArray[p]; 2106 2107 if (cval != defVal) { 2108 PetscCall(DMLabelGetValue(label, p, &val)); 2109 if (val == defVal) { 2110 PetscCall(DMLabelSetValue(label, p, cval)); 2111 if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx)); 2112 } 2113 } 2114 } 2115 /* Read out shared roots */ 2116 PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree)); 2117 PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree)); 2118 for (r = 0; r < Nr; ++r) { 2119 if (degree[r]) { 2120 const PetscInt cval = valArray[r]; 2121 2122 if (cval != defVal) { 2123 PetscCall(DMLabelGetValue(label, r, &val)); 2124 if (val == defVal) { 2125 PetscCall(DMLabelSetValue(label, r, cval)); 2126 if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx)); 2127 } 2128 } 2129 } 2130 } 2131 PetscFunctionReturn(PETSC_SUCCESS); 2132 } 2133 2134 /*@ 2135 DMLabelPropagateBegin - Setup a cycle of label propagation 2136 2137 Collective 2138 2139 Input Parameters: 2140 + label - The `DMLabel` to propagate across processes 2141 - sf - The `PetscSF` describing parallel layout of the label points 2142 2143 Level: intermediate 2144 2145 .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()` 2146 @*/ 2147 PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf) 2148 { 2149 PetscInt Nr, r, defVal; 2150 PetscMPIInt size; 2151 2152 PetscFunctionBegin; 2153 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2154 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size)); 2155 if (size > 1) { 2156 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 2157 PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL)); 2158 if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray)); 2159 for (r = 0; r < Nr; ++r) label->propArray[r] = defVal; 2160 } 2161 PetscFunctionReturn(PETSC_SUCCESS); 2162 } 2163 2164 /*@ 2165 DMLabelPropagateEnd - Tear down a cycle of label propagation 2166 2167 Collective 2168 2169 Input Parameters: 2170 + label - The `DMLabel` to propagate across processes 2171 - pointSF - The `PetscSF` describing parallel layout of the label points 2172 2173 Level: intermediate 2174 2175 .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()` 2176 @*/ 2177 PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF) 2178 { 2179 PetscFunctionBegin; 2180 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2181 PetscCall(PetscFree(label->propArray)); 2182 label->propArray = NULL; 2183 PetscFunctionReturn(PETSC_SUCCESS); 2184 } 2185 2186 /*@C 2187 DMLabelPropagatePush - Tear down a cycle of label propagation 2188 2189 Collective 2190 2191 Input Parameters: 2192 + label - The `DMLabel` to propagate across processes 2193 . pointSF - The `PetscSF` describing parallel layout of the label points 2194 . markPoint - An optional callback that is called when a point is marked, or `NULL` 2195 - ctx - An optional user context for the callback, or `NULL` 2196 2197 Calling sequence of `markPoint`: 2198 + label - The `DMLabel` 2199 . p - The point being marked 2200 . val - The label value for `p` 2201 - ctx - An optional user context 2202 2203 Level: intermediate 2204 2205 .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()` 2206 @*/ 2207 PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, void *ctx), void *ctx) 2208 { 2209 PetscInt *valArray = label->propArray, Nr; 2210 PetscMPIInt size; 2211 2212 PetscFunctionBegin; 2213 PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered"); 2214 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size)); 2215 PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL)); 2216 if (size > 1 && Nr >= 0) { 2217 /* Communicate marked edges 2218 The current implementation allocates an array the size of the number of root. We put the label values into the 2219 array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent. 2220 2221 TODO: We could use in-place communication with a different SF 2222 We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has 2223 already been marked. If not, it might have been handled on the process in this round, but we add it anyway. 2224 2225 In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new 2226 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 2227 edge to the queue. 2228 */ 2229 PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray)); 2230 PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX)); 2231 PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX)); 2232 PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE)); 2233 PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE)); 2234 PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx)); 2235 } 2236 PetscFunctionReturn(PETSC_SUCCESS); 2237 } 2238 2239 /*@ 2240 DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label 2241 2242 Not Collective 2243 2244 Input Parameter: 2245 . label - the `DMLabel` 2246 2247 Output Parameters: 2248 + section - the section giving offsets for each stratum 2249 - is - An `IS` containing all the label points 2250 2251 Level: developer 2252 2253 .seealso: `DMLabel`, `DM`, `DMLabelDistribute()` 2254 @*/ 2255 PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is) 2256 { 2257 IS vIS; 2258 const PetscInt *values; 2259 PetscInt *points; 2260 PetscInt nV, vS = 0, vE = 0, v, N; 2261 2262 PetscFunctionBegin; 2263 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2264 PetscCall(DMLabelGetNumValues(label, &nV)); 2265 PetscCall(DMLabelGetValueIS(label, &vIS)); 2266 PetscCall(ISGetIndices(vIS, &values)); 2267 if (nV) { 2268 vS = values[0]; 2269 vE = values[0] + 1; 2270 } 2271 for (v = 1; v < nV; ++v) { 2272 vS = PetscMin(vS, values[v]); 2273 vE = PetscMax(vE, values[v] + 1); 2274 } 2275 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section)); 2276 PetscCall(PetscSectionSetChart(*section, vS, vE)); 2277 for (v = 0; v < nV; ++v) { 2278 PetscInt n; 2279 2280 PetscCall(DMLabelGetStratumSize(label, values[v], &n)); 2281 PetscCall(PetscSectionSetDof(*section, values[v], n)); 2282 } 2283 PetscCall(PetscSectionSetUp(*section)); 2284 PetscCall(PetscSectionGetStorageSize(*section, &N)); 2285 PetscCall(PetscMalloc1(N, &points)); 2286 for (v = 0; v < nV; ++v) { 2287 IS is; 2288 const PetscInt *spoints; 2289 PetscInt dof, off, p; 2290 2291 PetscCall(PetscSectionGetDof(*section, values[v], &dof)); 2292 PetscCall(PetscSectionGetOffset(*section, values[v], &off)); 2293 PetscCall(DMLabelGetStratumIS(label, values[v], &is)); 2294 PetscCall(ISGetIndices(is, &spoints)); 2295 for (p = 0; p < dof; ++p) points[off + p] = spoints[p]; 2296 PetscCall(ISRestoreIndices(is, &spoints)); 2297 PetscCall(ISDestroy(&is)); 2298 } 2299 PetscCall(ISRestoreIndices(vIS, &values)); 2300 PetscCall(ISDestroy(&vIS)); 2301 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is)); 2302 PetscFunctionReturn(PETSC_SUCCESS); 2303 } 2304 2305 /*@C 2306 DMLabelRegister - Adds a new label component implementation 2307 2308 Not Collective 2309 2310 Input Parameters: 2311 + name - The name of a new user-defined creation routine 2312 - create_func - The creation routine itself 2313 2314 Notes: 2315 `DMLabelRegister()` may be called multiple times to add several user-defined labels 2316 2317 Example Usage: 2318 .vb 2319 DMLabelRegister("my_label", MyLabelCreate); 2320 .ve 2321 2322 Then, your label type can be chosen with the procedural interface via 2323 .vb 2324 DMLabelCreate(MPI_Comm, DMLabel *); 2325 DMLabelSetType(DMLabel, "my_label"); 2326 .ve 2327 or at runtime via the option 2328 .vb 2329 -dm_label_type my_label 2330 .ve 2331 2332 Level: advanced 2333 2334 .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()` 2335 @*/ 2336 PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel)) 2337 { 2338 PetscFunctionBegin; 2339 PetscCall(DMInitializePackage()); 2340 PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func)); 2341 PetscFunctionReturn(PETSC_SUCCESS); 2342 } 2343 2344 PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel); 2345 PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel); 2346 2347 /*@C 2348 DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package. 2349 2350 Not Collective 2351 2352 Level: advanced 2353 2354 .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()` 2355 @*/ 2356 PetscErrorCode DMLabelRegisterAll(void) 2357 { 2358 PetscFunctionBegin; 2359 if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 2360 DMLabelRegisterAllCalled = PETSC_TRUE; 2361 2362 PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete)); 2363 PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral)); 2364 PetscFunctionReturn(PETSC_SUCCESS); 2365 } 2366 2367 /*@C 2368 DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`. 2369 2370 Level: developer 2371 2372 .seealso: `DMLabel`, `DM`, `PetscInitialize()` 2373 @*/ 2374 PetscErrorCode DMLabelRegisterDestroy(void) 2375 { 2376 PetscFunctionBegin; 2377 PetscCall(PetscFunctionListDestroy(&DMLabelList)); 2378 DMLabelRegisterAllCalled = PETSC_FALSE; 2379 PetscFunctionReturn(PETSC_SUCCESS); 2380 } 2381 2382 /*@ 2383 DMLabelSetType - Sets the particular implementation for a label. 2384 2385 Collective 2386 2387 Input Parameters: 2388 + label - The label 2389 - method - The name of the label type 2390 2391 Options Database Key: 2392 . -dm_label_type <type> - Sets the label type; use -help for a list of available types or see `DMLabelType` 2393 2394 Level: intermediate 2395 2396 .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()` 2397 @*/ 2398 PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method) 2399 { 2400 PetscErrorCode (*r)(DMLabel); 2401 PetscBool match; 2402 2403 PetscFunctionBegin; 2404 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2405 PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match)); 2406 if (match) PetscFunctionReturn(PETSC_SUCCESS); 2407 2408 PetscCall(DMLabelRegisterAll()); 2409 PetscCall(PetscFunctionListFind(DMLabelList, method, &r)); 2410 PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method); 2411 2412 PetscTryTypeMethod(label, destroy); 2413 PetscCall(PetscMemzero(label->ops, sizeof(*label->ops))); 2414 PetscCall(PetscObjectChangeTypeName((PetscObject)label, method)); 2415 PetscCall((*r)(label)); 2416 PetscFunctionReturn(PETSC_SUCCESS); 2417 } 2418 2419 /*@ 2420 DMLabelGetType - Gets the type name (as a string) from the label. 2421 2422 Not Collective 2423 2424 Input Parameter: 2425 . label - The `DMLabel` 2426 2427 Output Parameter: 2428 . type - The `DMLabel` type name 2429 2430 Level: intermediate 2431 2432 .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()` 2433 @*/ 2434 PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type) 2435 { 2436 PetscFunctionBegin; 2437 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2438 PetscAssertPointer(type, 2); 2439 PetscCall(DMLabelRegisterAll()); 2440 *type = ((PetscObject)label)->type_name; 2441 PetscFunctionReturn(PETSC_SUCCESS); 2442 } 2443 2444 static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label) 2445 { 2446 PetscFunctionBegin; 2447 label->ops->view = DMLabelView_Concrete; 2448 label->ops->setup = NULL; 2449 label->ops->duplicate = DMLabelDuplicate_Concrete; 2450 label->ops->getstratumis = DMLabelGetStratumIS_Concrete; 2451 PetscFunctionReturn(PETSC_SUCCESS); 2452 } 2453 2454 PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label) 2455 { 2456 PetscFunctionBegin; 2457 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2458 PetscCall(DMLabelInitialize_Concrete(label)); 2459 PetscFunctionReturn(PETSC_SUCCESS); 2460 } 2461 2462 /*@ 2463 PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using 2464 the local section and an `PetscSF` describing the section point overlap. 2465 2466 Collective 2467 2468 Input Parameters: 2469 + s - The `PetscSection` for the local field layout 2470 . sf - The `PetscSF` describing parallel layout of the section points 2471 . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs 2472 . label - The label specifying the points 2473 - labelValue - The label stratum specifying the points 2474 2475 Output Parameter: 2476 . gsection - The `PetscSection` for the global field layout 2477 2478 Level: developer 2479 2480 Note: 2481 This gives negative sizes and offsets to points not owned by this process 2482 2483 .seealso: `DMLabel`, `DM`, `PetscSectionCreate()` 2484 @*/ 2485 PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection) 2486 { 2487 PetscInt *neg = NULL, *tmpOff = NULL; 2488 PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots; 2489 2490 PetscFunctionBegin; 2491 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2492 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 2493 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4); 2494 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection)); 2495 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 2496 PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd)); 2497 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 2498 if (nroots >= 0) { 2499 PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart); 2500 PetscCall(PetscCalloc1(nroots, &neg)); 2501 if (nroots > pEnd - pStart) { 2502 PetscCall(PetscCalloc1(nroots, &tmpOff)); 2503 } else { 2504 tmpOff = &(*gsection)->atlasDof[-pStart]; 2505 } 2506 } 2507 /* Mark ghost points with negative dof */ 2508 for (p = pStart; p < pEnd; ++p) { 2509 PetscInt value; 2510 2511 PetscCall(DMLabelGetValue(label, p, &value)); 2512 if (value != labelValue) continue; 2513 PetscCall(PetscSectionGetDof(s, p, &dof)); 2514 PetscCall(PetscSectionSetDof(*gsection, p, dof)); 2515 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 2516 if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof)); 2517 if (neg) neg[p] = -(dof + 1); 2518 } 2519 PetscCall(PetscSectionSetUpBC(*gsection)); 2520 if (nroots >= 0) { 2521 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2522 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2523 if (nroots > pEnd - pStart) { 2524 for (p = pStart; p < pEnd; ++p) { 2525 if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p]; 2526 } 2527 } 2528 } 2529 /* Calculate new sizes, get process offset, and calculate point offsets */ 2530 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 2531 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0; 2532 (*gsection)->atlasOff[p] = off; 2533 off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0; 2534 } 2535 PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s))); 2536 globalOff -= off; 2537 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 2538 (*gsection)->atlasOff[p] += globalOff; 2539 if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1); 2540 } 2541 /* Put in negative offsets for ghost points */ 2542 if (nroots >= 0) { 2543 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2544 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 2545 if (nroots > pEnd - pStart) { 2546 for (p = pStart; p < pEnd; ++p) { 2547 if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p]; 2548 } 2549 } 2550 } 2551 if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff)); 2552 PetscCall(PetscFree(neg)); 2553 PetscFunctionReturn(PETSC_SUCCESS); 2554 } 2555 2556 typedef struct _n_PetscSectionSym_Label { 2557 DMLabel label; 2558 PetscCopyMode *modes; 2559 PetscInt *sizes; 2560 const PetscInt ***perms; 2561 const PetscScalar ***rots; 2562 PetscInt (*minMaxOrients)[2]; 2563 PetscInt numStrata; /* numStrata is only increasing, functions as a state */ 2564 } PetscSectionSym_Label; 2565 2566 static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym) 2567 { 2568 PetscInt i, j; 2569 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2570 2571 PetscFunctionBegin; 2572 for (i = 0; i <= sl->numStrata; i++) { 2573 if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) { 2574 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2575 if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j])); 2576 if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j])); 2577 } 2578 if (sl->perms[i]) { 2579 const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]]; 2580 2581 PetscCall(PetscFree(perms)); 2582 } 2583 if (sl->rots[i]) { 2584 const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]]; 2585 2586 PetscCall(PetscFree(rots)); 2587 } 2588 } 2589 } 2590 PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients)); 2591 PetscCall(DMLabelDestroy(&sl->label)); 2592 sl->numStrata = 0; 2593 PetscFunctionReturn(PETSC_SUCCESS); 2594 } 2595 2596 static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym) 2597 { 2598 PetscFunctionBegin; 2599 PetscCall(PetscSectionSymLabelReset(sym)); 2600 PetscCall(PetscFree(sym->data)); 2601 PetscFunctionReturn(PETSC_SUCCESS); 2602 } 2603 2604 static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer) 2605 { 2606 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2607 PetscBool isAscii; 2608 DMLabel label = sl->label; 2609 const char *name; 2610 2611 PetscFunctionBegin; 2612 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii)); 2613 if (isAscii) { 2614 PetscInt i, j, k; 2615 PetscViewerFormat format; 2616 2617 PetscCall(PetscViewerGetFormat(viewer, &format)); 2618 if (label) { 2619 PetscCall(PetscViewerGetFormat(viewer, &format)); 2620 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2621 PetscCall(PetscViewerASCIIPushTab(viewer)); 2622 PetscCall(DMLabelView(label, viewer)); 2623 PetscCall(PetscViewerASCIIPopTab(viewer)); 2624 } else { 2625 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2626 PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name)); 2627 } 2628 } else { 2629 PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n")); 2630 } 2631 PetscCall(PetscViewerASCIIPushTab(viewer)); 2632 for (i = 0; i <= sl->numStrata; i++) { 2633 PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue; 2634 2635 if (!(sl->perms[i] || sl->rots[i])) { 2636 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i])); 2637 } else { 2638 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i])); 2639 PetscCall(PetscViewerASCIIPushTab(viewer)); 2640 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1])); 2641 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2642 PetscCall(PetscViewerASCIIPushTab(viewer)); 2643 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2644 if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) { 2645 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j)); 2646 } else { 2647 PetscInt tab; 2648 2649 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j)); 2650 PetscCall(PetscViewerASCIIPushTab(viewer)); 2651 PetscCall(PetscViewerASCIIGetTab(viewer, &tab)); 2652 if (sl->perms[i] && sl->perms[i][j]) { 2653 PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:")); 2654 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 2655 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k])); 2656 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 2657 PetscCall(PetscViewerASCIISetTab(viewer, tab)); 2658 } 2659 if (sl->rots[i] && sl->rots[i][j]) { 2660 PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: ")); 2661 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 2662 #if defined(PETSC_USE_COMPLEX) 2663 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]))); 2664 #else 2665 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k])); 2666 #endif 2667 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 2668 PetscCall(PetscViewerASCIISetTab(viewer, tab)); 2669 } 2670 PetscCall(PetscViewerASCIIPopTab(viewer)); 2671 } 2672 } 2673 PetscCall(PetscViewerASCIIPopTab(viewer)); 2674 } 2675 PetscCall(PetscViewerASCIIPopTab(viewer)); 2676 } 2677 } 2678 PetscCall(PetscViewerASCIIPopTab(viewer)); 2679 } 2680 PetscFunctionReturn(PETSC_SUCCESS); 2681 } 2682 2683 /*@ 2684 PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries 2685 2686 Logically 2687 2688 Input Parameters: 2689 + sym - the section symmetries 2690 - label - the `DMLabel` describing the types of points 2691 2692 Level: developer: 2693 2694 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()` 2695 @*/ 2696 PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label) 2697 { 2698 PetscSectionSym_Label *sl; 2699 2700 PetscFunctionBegin; 2701 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2702 sl = (PetscSectionSym_Label *)sym->data; 2703 if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym)); 2704 if (label) { 2705 sl->label = label; 2706 PetscCall(PetscObjectReference((PetscObject)label)); 2707 PetscCall(DMLabelGetNumValues(label, &sl->numStrata)); 2708 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)); 2709 PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode))); 2710 PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt))); 2711 PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **))); 2712 PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **))); 2713 PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2]))); 2714 } 2715 PetscFunctionReturn(PETSC_SUCCESS); 2716 } 2717 2718 /*@C 2719 PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum 2720 2721 Logically Collective 2722 2723 Input Parameters: 2724 + sym - the section symmetries 2725 - stratum - the stratum value in the label that we are assigning symmetries for 2726 2727 Output Parameters: 2728 + size - the number of dofs for points in the `stratum` of the label 2729 . minOrient - the smallest orientation for a point in this `stratum` 2730 . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`)) 2731 . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity 2732 - 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 2733 2734 Level: developer 2735 2736 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()` 2737 @*/ 2738 PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots) 2739 { 2740 PetscSectionSym_Label *sl; 2741 const char *name; 2742 PetscInt i; 2743 2744 PetscFunctionBegin; 2745 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2746 sl = (PetscSectionSym_Label *)sym->data; 2747 PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2748 for (i = 0; i <= sl->numStrata; i++) { 2749 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2750 2751 if (stratum == value) break; 2752 } 2753 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2754 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name); 2755 if (size) { 2756 PetscAssertPointer(size, 3); 2757 *size = sl->sizes[i]; 2758 } 2759 if (minOrient) { 2760 PetscAssertPointer(minOrient, 4); 2761 *minOrient = sl->minMaxOrients[i][0]; 2762 } 2763 if (maxOrient) { 2764 PetscAssertPointer(maxOrient, 5); 2765 *maxOrient = sl->minMaxOrients[i][1]; 2766 } 2767 if (perms) { 2768 PetscAssertPointer(perms, 6); 2769 *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]); 2770 } 2771 if (rots) { 2772 PetscAssertPointer(rots, 7); 2773 *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]); 2774 } 2775 PetscFunctionReturn(PETSC_SUCCESS); 2776 } 2777 2778 /*@C 2779 PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum 2780 2781 Logically 2782 2783 Input Parameters: 2784 + sym - the section symmetries 2785 . stratum - the stratum value in the label that we are assigning symmetries for 2786 . size - the number of dofs for points in the `stratum` of the label 2787 . minOrient - the smallest orientation for a point in this `stratum` 2788 . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`)) 2789 . mode - how `sym` should copy the `perms` and `rots` arrays 2790 . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity 2791 - 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 2792 2793 Level: developer 2794 2795 .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()` 2796 @*/ 2797 PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots) 2798 { 2799 PetscSectionSym_Label *sl; 2800 const char *name; 2801 PetscInt i, j, k; 2802 2803 PetscFunctionBegin; 2804 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2805 sl = (PetscSectionSym_Label *)sym->data; 2806 PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2807 for (i = 0; i <= sl->numStrata; i++) { 2808 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2809 2810 if (stratum == value) break; 2811 } 2812 PetscCall(PetscObjectGetName((PetscObject)sl->label, &name)); 2813 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name); 2814 sl->sizes[i] = size; 2815 sl->modes[i] = mode; 2816 sl->minMaxOrients[i][0] = minOrient; 2817 sl->minMaxOrients[i][1] = maxOrient; 2818 if (mode == PETSC_COPY_VALUES) { 2819 if (perms) { 2820 PetscInt **ownPerms; 2821 2822 PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms)); 2823 for (j = 0; j < maxOrient - minOrient; j++) { 2824 if (perms[j]) { 2825 PetscCall(PetscMalloc1(size, &ownPerms[j])); 2826 for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k]; 2827 } 2828 } 2829 sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient]; 2830 } 2831 if (rots) { 2832 PetscScalar **ownRots; 2833 2834 PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots)); 2835 for (j = 0; j < maxOrient - minOrient; j++) { 2836 if (rots[j]) { 2837 PetscCall(PetscMalloc1(size, &ownRots[j])); 2838 for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k]; 2839 } 2840 } 2841 sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient]; 2842 } 2843 } else { 2844 sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient); 2845 sl->rots[i] = PetscSafePointerPlusOffset(rots, -minOrient); 2846 } 2847 PetscFunctionReturn(PETSC_SUCCESS); 2848 } 2849 2850 static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots) 2851 { 2852 PetscInt i, j, numStrata; 2853 PetscSectionSym_Label *sl; 2854 DMLabel label; 2855 2856 PetscFunctionBegin; 2857 sl = (PetscSectionSym_Label *)sym->data; 2858 numStrata = sl->numStrata; 2859 label = sl->label; 2860 for (i = 0; i < numPoints; i++) { 2861 PetscInt point = points[2 * i]; 2862 PetscInt ornt = points[2 * i + 1]; 2863 2864 for (j = 0; j < numStrata; j++) { 2865 if (label->validIS[j]) { 2866 PetscInt k; 2867 2868 PetscCall(ISLocate(label->points[j], point, &k)); 2869 if (k >= 0) break; 2870 } else { 2871 PetscBool has; 2872 2873 PetscCall(PetscHSetIHas(label->ht[j], point, &has)); 2874 if (has) break; 2875 } 2876 } 2877 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], 2878 j < numStrata ? label->stratumValues[j] : label->defaultValue); 2879 if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL; 2880 if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL; 2881 } 2882 PetscFunctionReturn(PETSC_SUCCESS); 2883 } 2884 2885 static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym) 2886 { 2887 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data; 2888 IS valIS; 2889 const PetscInt *values; 2890 PetscInt Nv, v; 2891 2892 PetscFunctionBegin; 2893 PetscCall(DMLabelGetNumValues(sl->label, &Nv)); 2894 PetscCall(DMLabelGetValueIS(sl->label, &valIS)); 2895 PetscCall(ISGetIndices(valIS, &values)); 2896 for (v = 0; v < Nv; ++v) { 2897 const PetscInt val = values[v]; 2898 PetscInt size, minOrient, maxOrient; 2899 const PetscInt **perms; 2900 const PetscScalar **rots; 2901 2902 PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots)); 2903 PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots)); 2904 } 2905 PetscCall(ISDestroy(&valIS)); 2906 PetscFunctionReturn(PETSC_SUCCESS); 2907 } 2908 2909 static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym) 2910 { 2911 PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data; 2912 DMLabel dlabel; 2913 2914 PetscFunctionBegin; 2915 PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel)); 2916 PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym)); 2917 PetscCall(DMLabelDestroy(&dlabel)); 2918 PetscCall(PetscSectionSymCopy(sym, *dsym)); 2919 PetscFunctionReturn(PETSC_SUCCESS); 2920 } 2921 2922 PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym) 2923 { 2924 PetscSectionSym_Label *sl; 2925 2926 PetscFunctionBegin; 2927 PetscCall(PetscNew(&sl)); 2928 sym->ops->getpoints = PetscSectionSymGetPoints_Label; 2929 sym->ops->distribute = PetscSectionSymDistribute_Label; 2930 sym->ops->copy = PetscSectionSymCopy_Label; 2931 sym->ops->view = PetscSectionSymView_Label; 2932 sym->ops->destroy = PetscSectionSymDestroy_Label; 2933 sym->data = (void *)sl; 2934 PetscFunctionReturn(PETSC_SUCCESS); 2935 } 2936 2937 /*@ 2938 PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label 2939 2940 Collective 2941 2942 Input Parameters: 2943 + comm - the MPI communicator for the new symmetry 2944 - label - the label defining the strata 2945 2946 Output Parameter: 2947 . sym - the section symmetries 2948 2949 Level: developer 2950 2951 .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()` 2952 @*/ 2953 PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym) 2954 { 2955 PetscFunctionBegin; 2956 PetscCall(DMInitializePackage()); 2957 PetscCall(PetscSectionSymCreate(comm, sym)); 2958 PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL)); 2959 PetscCall(PetscSectionSymLabelSetLabel(*sym, label)); 2960 PetscFunctionReturn(PETSC_SUCCESS); 2961 } 2962