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 PetscCheckFalse((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 PetscCheckFalse(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 PetscCheckFalse(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 PetscCheckFalse((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 PetscCheckFalse((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 PetscCheckFalse((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 PetscCheckFalse((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 PetscCheckFalse((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 PetscCheckFalse((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 /*@ 1868 DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label 1869 1870 Not collective 1871 1872 Input Parameter: 1873 . label - the DMLabel 1874 1875 Output Parameters: 1876 + section - the section giving offsets for each stratum 1877 - is - An IS containing all the label points 1878 1879 Level: developer 1880 1881 .seealso: DMLabelDistribute() 1882 @*/ 1883 PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is) 1884 { 1885 IS vIS; 1886 const PetscInt *values; 1887 PetscInt *points; 1888 PetscInt nV, vS = 0, vE = 0, v, N; 1889 1890 PetscFunctionBegin; 1891 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1892 PetscCall(DMLabelGetNumValues(label, &nV)); 1893 PetscCall(DMLabelGetValueIS(label, &vIS)); 1894 PetscCall(ISGetIndices(vIS, &values)); 1895 if (nV) {vS = values[0]; vE = values[0]+1;} 1896 for (v = 1; v < nV; ++v) { 1897 vS = PetscMin(vS, values[v]); 1898 vE = PetscMax(vE, values[v]+1); 1899 } 1900 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section)); 1901 PetscCall(PetscSectionSetChart(*section, vS, vE)); 1902 for (v = 0; v < nV; ++v) { 1903 PetscInt n; 1904 1905 PetscCall(DMLabelGetStratumSize(label, values[v], &n)); 1906 PetscCall(PetscSectionSetDof(*section, values[v], n)); 1907 } 1908 PetscCall(PetscSectionSetUp(*section)); 1909 PetscCall(PetscSectionGetStorageSize(*section, &N)); 1910 PetscCall(PetscMalloc1(N, &points)); 1911 for (v = 0; v < nV; ++v) { 1912 IS is; 1913 const PetscInt *spoints; 1914 PetscInt dof, off, p; 1915 1916 PetscCall(PetscSectionGetDof(*section, values[v], &dof)); 1917 PetscCall(PetscSectionGetOffset(*section, values[v], &off)); 1918 PetscCall(DMLabelGetStratumIS(label, values[v], &is)); 1919 PetscCall(ISGetIndices(is, &spoints)); 1920 for (p = 0; p < dof; ++p) points[off+p] = spoints[p]; 1921 PetscCall(ISRestoreIndices(is, &spoints)); 1922 PetscCall(ISDestroy(&is)); 1923 } 1924 PetscCall(ISRestoreIndices(vIS, &values)); 1925 PetscCall(ISDestroy(&vIS)); 1926 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is)); 1927 PetscFunctionReturn(0); 1928 } 1929 1930 /*@ 1931 PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using 1932 the local section and an SF describing the section point overlap. 1933 1934 Collective on sf 1935 1936 Input Parameters: 1937 + s - The PetscSection for the local field layout 1938 . sf - The SF describing parallel layout of the section points 1939 . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs 1940 . label - The label specifying the points 1941 - labelValue - The label stratum specifying the points 1942 1943 Output Parameter: 1944 . gsection - The PetscSection for the global field layout 1945 1946 Note: This gives negative sizes and offsets to points not owned by this process 1947 1948 Level: developer 1949 1950 .seealso: PetscSectionCreate() 1951 @*/ 1952 PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection) 1953 { 1954 PetscInt *neg = NULL, *tmpOff = NULL; 1955 PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots; 1956 1957 PetscFunctionBegin; 1958 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1959 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 1960 PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4); 1961 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection)); 1962 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1963 PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd)); 1964 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 1965 if (nroots >= 0) { 1966 PetscCheckFalse(nroots < pEnd-pStart,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart); 1967 PetscCall(PetscCalloc1(nroots, &neg)); 1968 if (nroots > pEnd-pStart) { 1969 PetscCall(PetscCalloc1(nroots, &tmpOff)); 1970 } else { 1971 tmpOff = &(*gsection)->atlasDof[-pStart]; 1972 } 1973 } 1974 /* Mark ghost points with negative dof */ 1975 for (p = pStart; p < pEnd; ++p) { 1976 PetscInt value; 1977 1978 PetscCall(DMLabelGetValue(label, p, &value)); 1979 if (value != labelValue) continue; 1980 PetscCall(PetscSectionGetDof(s, p, &dof)); 1981 PetscCall(PetscSectionSetDof(*gsection, p, dof)); 1982 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 1983 if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof)); 1984 if (neg) neg[p] = -(dof+1); 1985 } 1986 PetscCall(PetscSectionSetUpBC(*gsection)); 1987 if (nroots >= 0) { 1988 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE)); 1989 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE)); 1990 if (nroots > pEnd-pStart) { 1991 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];} 1992 } 1993 } 1994 /* Calculate new sizes, get proccess offset, and calculate point offsets */ 1995 for (p = 0, off = 0; p < pEnd-pStart; ++p) { 1996 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0; 1997 (*gsection)->atlasOff[p] = off; 1998 off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0; 1999 } 2000 PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s))); 2001 globalOff -= off; 2002 for (p = 0, off = 0; p < pEnd-pStart; ++p) { 2003 (*gsection)->atlasOff[p] += globalOff; 2004 if (neg) neg[p] = -((*gsection)->atlasOff[p]+1); 2005 } 2006 /* Put in negative offsets for ghost points */ 2007 if (nroots >= 0) { 2008 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE)); 2009 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE)); 2010 if (nroots > pEnd-pStart) { 2011 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];} 2012 } 2013 } 2014 if (nroots >= 0 && nroots > pEnd-pStart) PetscCall(PetscFree(tmpOff)); 2015 PetscCall(PetscFree(neg)); 2016 PetscFunctionReturn(0); 2017 } 2018 2019 typedef struct _n_PetscSectionSym_Label 2020 { 2021 DMLabel label; 2022 PetscCopyMode *modes; 2023 PetscInt *sizes; 2024 const PetscInt ***perms; 2025 const PetscScalar ***rots; 2026 PetscInt (*minMaxOrients)[2]; 2027 PetscInt numStrata; /* numStrata is only increasing, functions as a state */ 2028 } PetscSectionSym_Label; 2029 2030 static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym) 2031 { 2032 PetscInt i, j; 2033 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data; 2034 2035 PetscFunctionBegin; 2036 for (i = 0; i <= sl->numStrata; i++) { 2037 if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) { 2038 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2039 if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j])); 2040 if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j])); 2041 } 2042 if (sl->perms[i]) { 2043 const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]]; 2044 2045 PetscCall(PetscFree(perms)); 2046 } 2047 if (sl->rots[i]) { 2048 const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]]; 2049 2050 PetscCall(PetscFree(rots)); 2051 } 2052 } 2053 } 2054 PetscCall(PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients)); 2055 PetscCall(DMLabelDestroy(&sl->label)); 2056 sl->numStrata = 0; 2057 PetscFunctionReturn(0); 2058 } 2059 2060 static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym) 2061 { 2062 PetscFunctionBegin; 2063 PetscCall(PetscSectionSymLabelReset(sym)); 2064 PetscCall(PetscFree(sym->data)); 2065 PetscFunctionReturn(0); 2066 } 2067 2068 static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer) 2069 { 2070 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data; 2071 PetscBool isAscii; 2072 DMLabel label = sl->label; 2073 const char *name; 2074 2075 PetscFunctionBegin; 2076 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii)); 2077 if (isAscii) { 2078 PetscInt i, j, k; 2079 PetscViewerFormat format; 2080 2081 PetscCall(PetscViewerGetFormat(viewer,&format)); 2082 if (label) { 2083 PetscCall(PetscViewerGetFormat(viewer,&format)); 2084 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2085 PetscCall(PetscViewerASCIIPushTab(viewer)); 2086 PetscCall(DMLabelView(label, viewer)); 2087 PetscCall(PetscViewerASCIIPopTab(viewer)); 2088 } else { 2089 PetscCall(PetscObjectGetName((PetscObject) sl->label, &name)); 2090 PetscCall(PetscViewerASCIIPrintf(viewer," Label '%s'\n",name)); 2091 } 2092 } else { 2093 PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n")); 2094 } 2095 PetscCall(PetscViewerASCIIPushTab(viewer)); 2096 for (i = 0; i <= sl->numStrata; i++) { 2097 PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue; 2098 2099 if (!(sl->perms[i] || sl->rots[i])) { 2100 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i])); 2101 } else { 2102 PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i])); 2103 PetscCall(PetscViewerASCIIPushTab(viewer)); 2104 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1])); 2105 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2106 PetscCall(PetscViewerASCIIPushTab(viewer)); 2107 for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) { 2108 if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) { 2109 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j)); 2110 } else { 2111 PetscInt tab; 2112 2113 PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j)); 2114 PetscCall(PetscViewerASCIIPushTab(viewer)); 2115 PetscCall(PetscViewerASCIIGetTab(viewer,&tab)); 2116 if (sl->perms[i] && sl->perms[i][j]) { 2117 PetscCall(PetscViewerASCIIPrintf(viewer,"Permutation:")); 2118 PetscCall(PetscViewerASCIISetTab(viewer,0)); 2119 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k])); 2120 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 2121 PetscCall(PetscViewerASCIISetTab(viewer,tab)); 2122 } 2123 if (sl->rots[i] && sl->rots[i][j]) { 2124 PetscCall(PetscViewerASCIIPrintf(viewer,"Rotations: ")); 2125 PetscCall(PetscViewerASCIISetTab(viewer,0)); 2126 #if defined(PETSC_USE_COMPLEX) 2127 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]))); 2128 #else 2129 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k])); 2130 #endif 2131 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 2132 PetscCall(PetscViewerASCIISetTab(viewer,tab)); 2133 } 2134 PetscCall(PetscViewerASCIIPopTab(viewer)); 2135 } 2136 } 2137 PetscCall(PetscViewerASCIIPopTab(viewer)); 2138 } 2139 PetscCall(PetscViewerASCIIPopTab(viewer)); 2140 } 2141 } 2142 PetscCall(PetscViewerASCIIPopTab(viewer)); 2143 } 2144 PetscFunctionReturn(0); 2145 } 2146 2147 /*@ 2148 PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries 2149 2150 Logically collective on sym 2151 2152 Input parameters: 2153 + sym - the section symmetries 2154 - label - the DMLabel describing the types of points 2155 2156 Level: developer: 2157 2158 .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms() 2159 @*/ 2160 PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label) 2161 { 2162 PetscSectionSym_Label *sl; 2163 2164 PetscFunctionBegin; 2165 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1); 2166 sl = (PetscSectionSym_Label *) sym->data; 2167 if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym)); 2168 if (label) { 2169 sl->label = label; 2170 PetscCall(PetscObjectReference((PetscObject) label)); 2171 PetscCall(DMLabelGetNumValues(label,&sl->numStrata)); 2172 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)); 2173 PetscCall(PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode))); 2174 PetscCall(PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt))); 2175 PetscCall(PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **))); 2176 PetscCall(PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **))); 2177 PetscCall(PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]))); 2178 } 2179 PetscFunctionReturn(0); 2180 } 2181 2182 /*@C 2183 PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum 2184 2185 Logically collective on sym 2186 2187 Input Parameters: 2188 + sym - the section symmetries 2189 - stratum - the stratum value in the label that we are assigning symmetries for 2190 2191 Output Parameters: 2192 + size - the number of dofs for points in the stratum of the label 2193 . minOrient - the smallest orientation for a point in this stratum 2194 . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient)) 2195 . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity 2196 - 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 2197 2198 Level: developer 2199 2200 .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel() 2201 @*/ 2202 PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots) 2203 { 2204 PetscSectionSym_Label *sl; 2205 const char *name; 2206 PetscInt i; 2207 2208 PetscFunctionBegin; 2209 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1); 2210 sl = (PetscSectionSym_Label *) sym->data; 2211 PetscCheck(sl->label, PetscObjectComm((PetscObject) sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2212 for (i = 0; i <= sl->numStrata; i++) { 2213 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2214 2215 if (stratum == value) break; 2216 } 2217 PetscCall(PetscObjectGetName((PetscObject) sl->label, &name)); 2218 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject) sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name); 2219 if (size) {PetscValidIntPointer(size, 3); *size = sl->sizes[i];} 2220 if (minOrient) {PetscValidIntPointer(minOrient, 4); *minOrient = sl->minMaxOrients[i][0];} 2221 if (maxOrient) {PetscValidIntPointer(maxOrient, 5); *maxOrient = sl->minMaxOrients[i][1];} 2222 if (perms) {PetscValidPointer(perms, 6); *perms = sl->perms[i] ? &sl->perms[i][sl->minMaxOrients[i][0]] : NULL;} 2223 if (rots) {PetscValidPointer(rots, 7); *rots = sl->rots[i] ? &sl->rots[i][sl->minMaxOrients[i][0]] : NULL;} 2224 PetscFunctionReturn(0); 2225 } 2226 2227 /*@C 2228 PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum 2229 2230 Logically collective on sym 2231 2232 InputParameters: 2233 + sym - the section symmetries 2234 . stratum - the stratum value in the label that we are assigning symmetries for 2235 . size - the number of dofs for points in the stratum of the label 2236 . minOrient - the smallest orientation for a point in this stratum 2237 . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient)) 2238 . mode - how sym should copy the perms and rots arrays 2239 . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity 2240 - 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 2241 2242 Level: developer 2243 2244 .seealso: PetscSectionSymLabelGetStratum(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel() 2245 @*/ 2246 PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots) 2247 { 2248 PetscSectionSym_Label *sl; 2249 const char *name; 2250 PetscInt i, j, k; 2251 2252 PetscFunctionBegin; 2253 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1); 2254 sl = (PetscSectionSym_Label *) sym->data; 2255 PetscCheck(sl->label, PetscObjectComm((PetscObject) sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet"); 2256 for (i = 0; i <= sl->numStrata; i++) { 2257 PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue; 2258 2259 if (stratum == value) break; 2260 } 2261 PetscCall(PetscObjectGetName((PetscObject) sl->label, &name)); 2262 PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject) sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %D not found in label %s", stratum, name); 2263 sl->sizes[i] = size; 2264 sl->modes[i] = mode; 2265 sl->minMaxOrients[i][0] = minOrient; 2266 sl->minMaxOrients[i][1] = maxOrient; 2267 if (mode == PETSC_COPY_VALUES) { 2268 if (perms) { 2269 PetscInt **ownPerms; 2270 2271 PetscCall(PetscCalloc1(maxOrient - minOrient,&ownPerms)); 2272 for (j = 0; j < maxOrient-minOrient; j++) { 2273 if (perms[j]) { 2274 PetscCall(PetscMalloc1(size,&ownPerms[j])); 2275 for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];} 2276 } 2277 } 2278 sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient]; 2279 } 2280 if (rots) { 2281 PetscScalar **ownRots; 2282 2283 PetscCall(PetscCalloc1(maxOrient - minOrient,&ownRots)); 2284 for (j = 0; j < maxOrient-minOrient; j++) { 2285 if (rots[j]) { 2286 PetscCall(PetscMalloc1(size,&ownRots[j])); 2287 for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];} 2288 } 2289 } 2290 sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient]; 2291 } 2292 } else { 2293 sl->perms[i] = perms ? &perms[-minOrient] : NULL; 2294 sl->rots[i] = rots ? &rots[-minOrient] : NULL; 2295 } 2296 PetscFunctionReturn(0); 2297 } 2298 2299 static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots) 2300 { 2301 PetscInt i, j, numStrata; 2302 PetscSectionSym_Label *sl; 2303 DMLabel label; 2304 2305 PetscFunctionBegin; 2306 sl = (PetscSectionSym_Label *) sym->data; 2307 numStrata = sl->numStrata; 2308 label = sl->label; 2309 for (i = 0; i < numPoints; i++) { 2310 PetscInt point = points[2*i]; 2311 PetscInt ornt = points[2*i+1]; 2312 2313 for (j = 0; j < numStrata; j++) { 2314 if (label->validIS[j]) { 2315 PetscInt k; 2316 2317 PetscCall(ISLocate(label->points[j],point,&k)); 2318 if (k >= 0) break; 2319 } else { 2320 PetscBool has; 2321 2322 PetscCall(PetscHSetIHas(label->ht[j], point, &has)); 2323 if (has) break; 2324 } 2325 } 2326 PetscCheckFalse((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); 2327 if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;} 2328 if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;} 2329 } 2330 PetscFunctionReturn(0); 2331 } 2332 2333 static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym) 2334 { 2335 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) nsym->data; 2336 IS valIS; 2337 const PetscInt *values; 2338 PetscInt Nv, v; 2339 2340 PetscFunctionBegin; 2341 PetscCall(DMLabelGetNumValues(sl->label, &Nv)); 2342 PetscCall(DMLabelGetValueIS(sl->label, &valIS)); 2343 PetscCall(ISGetIndices(valIS, &values)); 2344 for (v = 0; v < Nv; ++v) { 2345 const PetscInt val = values[v]; 2346 PetscInt size, minOrient, maxOrient; 2347 const PetscInt **perms; 2348 const PetscScalar **rots; 2349 2350 PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots)); 2351 PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots)); 2352 } 2353 PetscCall(ISDestroy(&valIS)); 2354 PetscFunctionReturn(0); 2355 } 2356 2357 static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym) 2358 { 2359 PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data; 2360 DMLabel dlabel; 2361 2362 PetscFunctionBegin; 2363 PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel)); 2364 PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject) sym), dlabel, dsym)); 2365 PetscCall(DMLabelDestroy(&dlabel)); 2366 PetscCall(PetscSectionSymCopy(sym, *dsym)); 2367 PetscFunctionReturn(0); 2368 } 2369 2370 PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym) 2371 { 2372 PetscSectionSym_Label *sl; 2373 2374 PetscFunctionBegin; 2375 PetscCall(PetscNewLog(sym,&sl)); 2376 sym->ops->getpoints = PetscSectionSymGetPoints_Label; 2377 sym->ops->distribute = PetscSectionSymDistribute_Label; 2378 sym->ops->copy = PetscSectionSymCopy_Label; 2379 sym->ops->view = PetscSectionSymView_Label; 2380 sym->ops->destroy = PetscSectionSymDestroy_Label; 2381 sym->data = (void *) sl; 2382 PetscFunctionReturn(0); 2383 } 2384 2385 /*@ 2386 PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label 2387 2388 Collective 2389 2390 Input Parameters: 2391 + comm - the MPI communicator for the new symmetry 2392 - label - the label defining the strata 2393 2394 Output Parameters: 2395 . sym - the section symmetries 2396 2397 Level: developer 2398 2399 .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms() 2400 @*/ 2401 PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym) 2402 { 2403 PetscFunctionBegin; 2404 PetscCall(DMInitializePackage()); 2405 PetscCall(PetscSectionSymCreate(comm,sym)); 2406 PetscCall(PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL)); 2407 PetscCall(PetscSectionSymLabelSetLabel(*sym,label)); 2408 PetscFunctionReturn(0); 2409 } 2410