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