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