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