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