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