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