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, Nr; 2010 PetscMPIInt size; 2011 2012 PetscFunctionBegin; 2013 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) pointSF), &size)); 2014 PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL)); 2015 if (size > 1 && Nr >= 0) { 2016 /* Communicate marked edges 2017 The current implementation allocates an array the size of the number of root. We put the label values into the 2018 array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent. 2019 2020 TODO: We could use in-place communication with a different SF 2021 We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has 2022 already been marked. If not, it might have been handled on the process in this round, but we add it anyway. 2023 2024 In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new 2025 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 2026 edge to the queue. 2027 */ 2028 PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray)); 2029 PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX)); 2030 PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX)); 2031 PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray,MPI_REPLACE)); 2032 PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray,MPI_REPLACE)); 2033 PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx)); 2034 } 2035 PetscFunctionReturn(0); 2036 } 2037 2038 /*@ 2039 DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label 2040 2041 Not collective 2042 2043 Input Parameter: 2044 . label - the DMLabel 2045 2046 Output Parameters: 2047 + section - the section giving offsets for each stratum 2048 - is - An IS containing all the label points 2049 2050 Level: developer 2051 2052 .seealso: DMLabelDistribute() 2053 @*/ 2054 PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is) 2055 { 2056 IS vIS; 2057 const PetscInt *values; 2058 PetscInt *points; 2059 PetscInt nV, vS = 0, vE = 0, v, N; 2060 2061 PetscFunctionBegin; 2062 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 2063 PetscCall(DMLabelGetNumValues(label, &nV)); 2064 PetscCall(DMLabelGetValueIS(label, &vIS)); 2065 PetscCall(ISGetIndices(vIS, &values)); 2066 if (nV) {vS = values[0]; vE = values[0]+1;} 2067 for (v = 1; v < nV; ++v) { 2068 vS = PetscMin(vS, values[v]); 2069 vE = PetscMax(vE, values[v]+1); 2070 } 2071 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section)); 2072 PetscCall(PetscSectionSetChart(*section, vS, vE)); 2073 for (v = 0; v < nV; ++v) { 2074 PetscInt n; 2075 2076 PetscCall(DMLabelGetStratumSize(label, values[v], &n)); 2077 PetscCall(PetscSectionSetDof(*section, values[v], n)); 2078 } 2079 PetscCall(PetscSectionSetUp(*section)); 2080 PetscCall(PetscSectionGetStorageSize(*section, &N)); 2081 PetscCall(PetscMalloc1(N, &points)); 2082 for (v = 0; v < nV; ++v) { 2083 IS is; 2084 const PetscInt *spoints; 2085 PetscInt dof, off, p; 2086 2087 PetscCall(PetscSectionGetDof(*section, values[v], &dof)); 2088 PetscCall(PetscSectionGetOffset(*section, values[v], &off)); 2089 PetscCall(DMLabelGetStratumIS(label, values[v], &is)); 2090 PetscCall(ISGetIndices(is, &spoints)); 2091 for (p = 0; p < dof; ++p) points[off+p] = spoints[p]; 2092 PetscCall(ISRestoreIndices(is, &spoints)); 2093 PetscCall(ISDestroy(&is)); 2094 } 2095 PetscCall(ISRestoreIndices(vIS, &values)); 2096 PetscCall(ISDestroy(&vIS)); 2097 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is)); 2098 PetscFunctionReturn(0); 2099 } 2100 2101 /*@ 2102 PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using 2103 the local section and an SF describing the section point overlap. 2104 2105 Collective on sf 2106 2107 Input Parameters: 2108 + s - The PetscSection for the local field layout 2109 . sf - The SF describing parallel layout of the section points 2110 . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs 2111 . label - The label specifying the points 2112 - labelValue - The label stratum specifying the points 2113 2114 Output Parameter: 2115 . gsection - The PetscSection for the global field layout 2116 2117 Note: This gives negative sizes and offsets to points not owned by this process 2118 2119 Level: developer 2120 2121 .seealso: PetscSectionCreate() 2122 @*/ 2123 PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection) 2124 { 2125 PetscInt *neg = NULL, *tmpOff = NULL; 2126 PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots; 2127 2128 PetscFunctionBegin; 2129 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2130 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 2131 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4); 2132 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection)); 2133 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 2134 PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd)); 2135 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 2136 if (nroots >= 0) { 2137 PetscCheck(nroots >= pEnd-pStart,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart); 2138 PetscCall(PetscCalloc1(nroots, &neg)); 2139 if (nroots > pEnd-pStart) { 2140 PetscCall(PetscCalloc1(nroots, &tmpOff)); 2141 } else { 2142 tmpOff = &(*gsection)->atlasDof[-pStart]; 2143 } 2144 } 2145 /* Mark ghost points with negative dof */ 2146 for (p = pStart; p < pEnd; ++p) { 2147 PetscInt value; 2148 2149 PetscCall(DMLabelGetValue(label, p, &value)); 2150 if (value != labelValue) continue; 2151 PetscCall(PetscSectionGetDof(s, p, &dof)); 2152 PetscCall(PetscSectionSetDof(*gsection, p, dof)); 2153 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 2154 if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof)); 2155 if (neg) neg[p] = -(dof+1); 2156 } 2157 PetscCall(PetscSectionSetUpBC(*gsection)); 2158 if (nroots >= 0) { 2159 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE)); 2160 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE)); 2161 if (nroots > pEnd-pStart) { 2162 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];} 2163 } 2164 } 2165 /* Calculate new sizes, get proccess offset, and calculate point offsets */ 2166 for (p = 0, off = 0; p < pEnd-pStart; ++p) { 2167 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0; 2168 (*gsection)->atlasOff[p] = off; 2169 off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0; 2170 } 2171 PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s))); 2172 globalOff -= off; 2173 for (p = 0, off = 0; p < pEnd-pStart; ++p) { 2174 (*gsection)->atlasOff[p] += globalOff; 2175 if (neg) neg[p] = -((*gsection)->atlasOff[p]+1); 2176 } 2177 /* Put in negative offsets for ghost points */ 2178 if (nroots >= 0) { 2179 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE)); 2180 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE)); 2181 if (nroots > pEnd-pStart) { 2182 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];} 2183 } 2184 } 2185 if (nroots >= 0 && nroots > pEnd-pStart) PetscCall(PetscFree(tmpOff)); 2186 PetscCall(PetscFree(neg)); 2187 PetscFunctionReturn(0); 2188 } 2189 2190 typedef struct _n_PetscSectionSym_Label 2191 { 2192 DMLabel label; 2193 PetscCopyMode *modes; 2194 PetscInt *sizes; 2195 const PetscInt ***perms; 2196 const PetscScalar ***rots; 2197 PetscInt (*minMaxOrients)[2]; 2198 PetscInt numStrata; /* numStrata is only increasing, functions as a state */ 2199 } PetscSectionSym_Label; 2200 2201 static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym) 2202 { 2203 PetscInt i, j; 2204 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data; 2205 2206 PetscFunctionBegin; 2207 for (i = 0; i <= sl->numStrata; i++) { 2208 if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) { 2209 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2210 if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j])); 2211 if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j])); 2212 } 2213 if (sl->perms[i]) { 2214 const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]]; 2215 2216 PetscCall(PetscFree(perms)); 2217 } 2218 if (sl->rots[i]) { 2219 const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]]; 2220 2221 PetscCall(PetscFree(rots)); 2222 } 2223 } 2224 } 2225 PetscCall(PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients)); 2226 PetscCall(DMLabelDestroy(&sl->label)); 2227 sl->numStrata = 0; 2228 PetscFunctionReturn(0); 2229 } 2230 2231 static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym) 2232 { 2233 PetscFunctionBegin; 2234 PetscCall(PetscSectionSymLabelReset(sym)); 2235 PetscCall(PetscFree(sym->data)); 2236 PetscFunctionReturn(0); 2237 } 2238 2239 static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer) 2240 { 2241 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data; 2242 PetscBool isAscii; 2243 DMLabel label = sl->label; 2244 const char *name; 2245 2246 PetscFunctionBegin; 2247 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii)); 2248 if (isAscii) { 2249 PetscInt i, j, k; 2250 PetscViewerFormat format; 2251 2252 PetscCall(PetscViewerGetFormat(viewer,&format)); 2253 if (label) { 2254 PetscCall(PetscViewerGetFormat(viewer,&format)); 2255 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2256 PetscCall(PetscViewerASCIIPushTab(viewer)); 2257 PetscCall(DMLabelView(label, viewer)); 2258 PetscCall(PetscViewerASCIIPopTab(viewer)); 2259 } else { 2260 PetscCall(PetscObjectGetName((PetscObject) sl->label, &name)); 2261 PetscCall(PetscViewerASCIIPrintf(viewer," Label '%s'\n",name)); 2262 } 2263 } else { 2264 PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n")); 2265 } 2266 PetscCall(PetscViewerASCIIPushTab(viewer)); 2267 for (i = 0; i <= sl->numStrata; i++) { 2268 PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue; 2269 2270 if (!(sl->perms[i] || sl->rots[i])) { 2271 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i])); 2272 } else { 2273 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i])); 2274 PetscCall(PetscViewerASCIIPushTab(viewer)); 2275 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1])); 2276 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2277 PetscCall(PetscViewerASCIIPushTab(viewer)); 2278 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2279 if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) { 2280 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j)); 2281 } else { 2282 PetscInt tab; 2283 2284 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j)); 2285 PetscCall(PetscViewerASCIIPushTab(viewer)); 2286 PetscCall(PetscViewerASCIIGetTab(viewer,&tab)); 2287 if (sl->perms[i] && sl->perms[i][j]) { 2288 PetscCall(PetscViewerASCIIPrintf(viewer,"Permutation:")); 2289 PetscCall(PetscViewerASCIISetTab(viewer,0)); 2290 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k])); 2291 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 2292 PetscCall(PetscViewerASCIISetTab(viewer,tab)); 2293 } 2294 if (sl->rots[i] && sl->rots[i][j]) { 2295 PetscCall(PetscViewerASCIIPrintf(viewer,"Rotations: ")); 2296 PetscCall(PetscViewerASCIISetTab(viewer,0)); 2297 #if defined(PETSC_USE_COMPLEX) 2298 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]))); 2299 #else 2300 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k])); 2301 #endif 2302 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 2303 PetscCall(PetscViewerASCIISetTab(viewer,tab)); 2304 } 2305 PetscCall(PetscViewerASCIIPopTab(viewer)); 2306 } 2307 } 2308 PetscCall(PetscViewerASCIIPopTab(viewer)); 2309 } 2310 PetscCall(PetscViewerASCIIPopTab(viewer)); 2311 } 2312 } 2313 PetscCall(PetscViewerASCIIPopTab(viewer)); 2314 } 2315 PetscFunctionReturn(0); 2316 } 2317 2318 /*@ 2319 PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries 2320 2321 Logically collective on sym 2322 2323 Input parameters: 2324 + sym - the section symmetries 2325 - label - the DMLabel describing the types of points 2326 2327 Level: developer: 2328 2329 .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms() 2330 @*/ 2331 PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label) 2332 { 2333 PetscSectionSym_Label *sl; 2334 2335 PetscFunctionBegin; 2336 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1); 2337 sl = (PetscSectionSym_Label *) sym->data; 2338 if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym)); 2339 if (label) { 2340 sl->label = label; 2341 PetscCall(PetscObjectReference((PetscObject) label)); 2342 PetscCall(DMLabelGetNumValues(label,&sl->numStrata)); 2343 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)); 2344 PetscCall(PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode))); 2345 PetscCall(PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt))); 2346 PetscCall(PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **))); 2347 PetscCall(PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **))); 2348 PetscCall(PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]))); 2349 } 2350 PetscFunctionReturn(0); 2351 } 2352 2353 /*@C 2354 PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum 2355 2356 Logically collective on sym 2357 2358 Input Parameters: 2359 + sym - the section symmetries 2360 - stratum - the stratum value in the label that we are assigning symmetries for 2361 2362 Output Parameters: 2363 + size - the number of dofs for points in the stratum of the label 2364 . minOrient - the smallest orientation for a point in this stratum 2365 . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient)) 2366 . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity 2367 - 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 2368 2369 Level: developer 2370 2371 .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel() 2372 @*/ 2373 PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots) 2374 { 2375 PetscSectionSym_Label *sl; 2376 const char *name; 2377 PetscInt i; 2378 2379 PetscFunctionBegin; 2380 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1); 2381 sl = (PetscSectionSym_Label *) sym->data; 2382 PetscCheck(sl->label, PetscObjectComm((PetscObject) sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2383 for (i = 0; i <= sl->numStrata; i++) { 2384 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2385 2386 if (stratum == value) break; 2387 } 2388 PetscCall(PetscObjectGetName((PetscObject) sl->label, &name)); 2389 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject) sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name); 2390 if (size) {PetscValidIntPointer(size, 3); *size = sl->sizes[i];} 2391 if (minOrient) {PetscValidIntPointer(minOrient, 4); *minOrient = sl->minMaxOrients[i][0];} 2392 if (maxOrient) {PetscValidIntPointer(maxOrient, 5); *maxOrient = sl->minMaxOrients[i][1];} 2393 if (perms) {PetscValidPointer(perms, 6); *perms = sl->perms[i] ? &sl->perms[i][sl->minMaxOrients[i][0]] : NULL;} 2394 if (rots) {PetscValidPointer(rots, 7); *rots = sl->rots[i] ? &sl->rots[i][sl->minMaxOrients[i][0]] : NULL;} 2395 PetscFunctionReturn(0); 2396 } 2397 2398 /*@C 2399 PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum 2400 2401 Logically collective on sym 2402 2403 InputParameters: 2404 + sym - the section symmetries 2405 . stratum - the stratum value in the label that we are assigning symmetries for 2406 . size - the number of dofs for points in the stratum of the label 2407 . minOrient - the smallest orientation for a point in this stratum 2408 . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient)) 2409 . mode - how sym should copy the perms and rots arrays 2410 . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity 2411 - 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 2412 2413 Level: developer 2414 2415 .seealso: PetscSectionSymLabelGetStratum(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel() 2416 @*/ 2417 PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots) 2418 { 2419 PetscSectionSym_Label *sl; 2420 const char *name; 2421 PetscInt i, j, k; 2422 2423 PetscFunctionBegin; 2424 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1); 2425 sl = (PetscSectionSym_Label *) sym->data; 2426 PetscCheck(sl->label, PetscObjectComm((PetscObject) sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2427 for (i = 0; i <= sl->numStrata; i++) { 2428 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2429 2430 if (stratum == value) break; 2431 } 2432 PetscCall(PetscObjectGetName((PetscObject) sl->label, &name)); 2433 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject) sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %D not found in label %s", stratum, name); 2434 sl->sizes[i] = size; 2435 sl->modes[i] = mode; 2436 sl->minMaxOrients[i][0] = minOrient; 2437 sl->minMaxOrients[i][1] = maxOrient; 2438 if (mode == PETSC_COPY_VALUES) { 2439 if (perms) { 2440 PetscInt **ownPerms; 2441 2442 PetscCall(PetscCalloc1(maxOrient - minOrient,&ownPerms)); 2443 for (j = 0; j < maxOrient-minOrient; j++) { 2444 if (perms[j]) { 2445 PetscCall(PetscMalloc1(size,&ownPerms[j])); 2446 for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];} 2447 } 2448 } 2449 sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient]; 2450 } 2451 if (rots) { 2452 PetscScalar **ownRots; 2453 2454 PetscCall(PetscCalloc1(maxOrient - minOrient,&ownRots)); 2455 for (j = 0; j < maxOrient-minOrient; j++) { 2456 if (rots[j]) { 2457 PetscCall(PetscMalloc1(size,&ownRots[j])); 2458 for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];} 2459 } 2460 } 2461 sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient]; 2462 } 2463 } else { 2464 sl->perms[i] = perms ? &perms[-minOrient] : NULL; 2465 sl->rots[i] = rots ? &rots[-minOrient] : NULL; 2466 } 2467 PetscFunctionReturn(0); 2468 } 2469 2470 static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots) 2471 { 2472 PetscInt i, j, numStrata; 2473 PetscSectionSym_Label *sl; 2474 DMLabel label; 2475 2476 PetscFunctionBegin; 2477 sl = (PetscSectionSym_Label *) sym->data; 2478 numStrata = sl->numStrata; 2479 label = sl->label; 2480 for (i = 0; i < numPoints; i++) { 2481 PetscInt point = points[2*i]; 2482 PetscInt ornt = points[2*i+1]; 2483 2484 for (j = 0; j < numStrata; j++) { 2485 if (label->validIS[j]) { 2486 PetscInt k; 2487 2488 PetscCall(ISLocate(label->points[j],point,&k)); 2489 if (k >= 0) break; 2490 } else { 2491 PetscBool has; 2492 2493 PetscCall(PetscHSetIHas(label->ht[j], point, &has)); 2494 if (has) break; 2495 } 2496 } 2497 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); 2498 if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;} 2499 if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;} 2500 } 2501 PetscFunctionReturn(0); 2502 } 2503 2504 static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym) 2505 { 2506 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) nsym->data; 2507 IS valIS; 2508 const PetscInt *values; 2509 PetscInt Nv, v; 2510 2511 PetscFunctionBegin; 2512 PetscCall(DMLabelGetNumValues(sl->label, &Nv)); 2513 PetscCall(DMLabelGetValueIS(sl->label, &valIS)); 2514 PetscCall(ISGetIndices(valIS, &values)); 2515 for (v = 0; v < Nv; ++v) { 2516 const PetscInt val = values[v]; 2517 PetscInt size, minOrient, maxOrient; 2518 const PetscInt **perms; 2519 const PetscScalar **rots; 2520 2521 PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots)); 2522 PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots)); 2523 } 2524 PetscCall(ISDestroy(&valIS)); 2525 PetscFunctionReturn(0); 2526 } 2527 2528 static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym) 2529 { 2530 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data; 2531 DMLabel dlabel; 2532 2533 PetscFunctionBegin; 2534 PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel)); 2535 PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject) sym), dlabel, dsym)); 2536 PetscCall(DMLabelDestroy(&dlabel)); 2537 PetscCall(PetscSectionSymCopy(sym, *dsym)); 2538 PetscFunctionReturn(0); 2539 } 2540 2541 PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym) 2542 { 2543 PetscSectionSym_Label *sl; 2544 2545 PetscFunctionBegin; 2546 PetscCall(PetscNewLog(sym,&sl)); 2547 sym->ops->getpoints = PetscSectionSymGetPoints_Label; 2548 sym->ops->distribute = PetscSectionSymDistribute_Label; 2549 sym->ops->copy = PetscSectionSymCopy_Label; 2550 sym->ops->view = PetscSectionSymView_Label; 2551 sym->ops->destroy = PetscSectionSymDestroy_Label; 2552 sym->data = (void *) sl; 2553 PetscFunctionReturn(0); 2554 } 2555 2556 /*@ 2557 PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label 2558 2559 Collective 2560 2561 Input Parameters: 2562 + comm - the MPI communicator for the new symmetry 2563 - label - the label defining the strata 2564 2565 Output Parameters: 2566 . sym - the section symmetries 2567 2568 Level: developer 2569 2570 .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms() 2571 @*/ 2572 PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym) 2573 { 2574 PetscFunctionBegin; 2575 PetscCall(DMInitializePackage()); 2576 PetscCall(PetscSectionSymCreate(comm,sym)); 2577 PetscCall(PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL)); 2578 PetscCall(PetscSectionSymLabelSetLabel(*sym,label)); 2579 PetscFunctionReturn(0); 2580 } 2581