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