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