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