1 /* 2 This file contains routines for basic section object implementation. 3 */ 4 5 #include <petsc/private/sectionimpl.h> /*I "petscsection.h" I*/ 6 #include <petscsf.h> 7 8 PetscClassId PETSC_SECTION_CLASSID; 9 10 /*@ 11 PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default. 12 13 Collective 14 15 Input Parameters: 16 + comm - the MPI communicator 17 - s - pointer to the section 18 19 Level: beginner 20 21 Notes: 22 Typical calling sequence 23 $ PetscSectionCreate(MPI_Comm,PetscSection *); 24 $ PetscSectionSetNumFields(PetscSection, numFields); 25 $ PetscSectionSetChart(PetscSection,low,high); 26 $ PetscSectionSetDof(PetscSection,point,numdof); 27 $ PetscSectionSetUp(PetscSection); 28 $ PetscSectionGetOffset(PetscSection,point,PetscInt *); 29 $ PetscSectionDestroy(PetscSection); 30 31 The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementations; it is 32 recommended they not be used in user codes unless you really gain something in their use. 33 34 .seealso: `PetscSection`, `PetscSectionDestroy()` 35 @*/ 36 PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s) 37 { 38 PetscFunctionBegin; 39 PetscValidPointer(s, 2); 40 PetscCall(ISInitializePackage()); 41 42 PetscCall(PetscHeaderCreate(*s, PETSC_SECTION_CLASSID, "PetscSection", "Section", "IS", comm, PetscSectionDestroy, PetscSectionView)); 43 44 (*s)->pStart = -1; 45 (*s)->pEnd = -1; 46 (*s)->perm = NULL; 47 (*s)->pointMajor = PETSC_TRUE; 48 (*s)->includesConstraints = PETSC_TRUE; 49 (*s)->atlasDof = NULL; 50 (*s)->atlasOff = NULL; 51 (*s)->bc = NULL; 52 (*s)->bcIndices = NULL; 53 (*s)->setup = PETSC_FALSE; 54 (*s)->numFields = 0; 55 (*s)->fieldNames = NULL; 56 (*s)->field = NULL; 57 (*s)->useFieldOff = PETSC_FALSE; 58 (*s)->compNames = NULL; 59 (*s)->clObj = NULL; 60 (*s)->clHash = NULL; 61 (*s)->clSection = NULL; 62 (*s)->clPoints = NULL; 63 PetscCall(PetscSectionInvalidateMaxDof_Internal(*s)); 64 PetscFunctionReturn(0); 65 } 66 67 /*@ 68 PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection 69 70 Collective 71 72 Input Parameter: 73 . section - the PetscSection 74 75 Output Parameter: 76 . newSection - the copy 77 78 Level: intermediate 79 80 .seealso: `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()` 81 @*/ 82 PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection) 83 { 84 PetscSectionSym sym; 85 IS perm; 86 PetscInt numFields, f, c, pStart, pEnd, p; 87 88 PetscFunctionBegin; 89 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 90 PetscValidHeaderSpecific(newSection, PETSC_SECTION_CLASSID, 2); 91 PetscCall(PetscSectionReset(newSection)); 92 PetscCall(PetscSectionGetNumFields(section, &numFields)); 93 if (numFields) PetscCall(PetscSectionSetNumFields(newSection, numFields)); 94 for (f = 0; f < numFields; ++f) { 95 const char *fieldName = NULL, *compName = NULL; 96 PetscInt numComp = 0; 97 98 PetscCall(PetscSectionGetFieldName(section, f, &fieldName)); 99 PetscCall(PetscSectionSetFieldName(newSection, f, fieldName)); 100 PetscCall(PetscSectionGetFieldComponents(section, f, &numComp)); 101 PetscCall(PetscSectionSetFieldComponents(newSection, f, numComp)); 102 for (c = 0; c < numComp; ++c) { 103 PetscCall(PetscSectionGetComponentName(section, f, c, &compName)); 104 PetscCall(PetscSectionSetComponentName(newSection, f, c, compName)); 105 } 106 PetscCall(PetscSectionGetFieldSym(section, f, &sym)); 107 PetscCall(PetscSectionSetFieldSym(newSection, f, sym)); 108 } 109 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 110 PetscCall(PetscSectionSetChart(newSection, pStart, pEnd)); 111 PetscCall(PetscSectionGetPermutation(section, &perm)); 112 PetscCall(PetscSectionSetPermutation(newSection, perm)); 113 PetscCall(PetscSectionGetSym(section, &sym)); 114 PetscCall(PetscSectionSetSym(newSection, sym)); 115 for (p = pStart; p < pEnd; ++p) { 116 PetscInt dof, cdof, fcdof = 0; 117 118 PetscCall(PetscSectionGetDof(section, p, &dof)); 119 PetscCall(PetscSectionSetDof(newSection, p, dof)); 120 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 121 if (cdof) PetscCall(PetscSectionSetConstraintDof(newSection, p, cdof)); 122 for (f = 0; f < numFields; ++f) { 123 PetscCall(PetscSectionGetFieldDof(section, p, f, &dof)); 124 PetscCall(PetscSectionSetFieldDof(newSection, p, f, dof)); 125 if (cdof) { 126 PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof)); 127 if (fcdof) PetscCall(PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof)); 128 } 129 } 130 } 131 PetscCall(PetscSectionSetUp(newSection)); 132 for (p = pStart; p < pEnd; ++p) { 133 PetscInt off, cdof, fcdof = 0; 134 const PetscInt *cInd; 135 136 /* Must set offsets in case they do not agree with the prefix sums */ 137 PetscCall(PetscSectionGetOffset(section, p, &off)); 138 PetscCall(PetscSectionSetOffset(newSection, p, off)); 139 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 140 if (cdof) { 141 PetscCall(PetscSectionGetConstraintIndices(section, p, &cInd)); 142 PetscCall(PetscSectionSetConstraintIndices(newSection, p, cInd)); 143 for (f = 0; f < numFields; ++f) { 144 PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof)); 145 if (fcdof) { 146 PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &cInd)); 147 PetscCall(PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd)); 148 } 149 } 150 } 151 } 152 PetscFunctionReturn(0); 153 } 154 155 /*@ 156 PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection 157 158 Collective 159 160 Input Parameter: 161 . section - the PetscSection 162 163 Output Parameter: 164 . newSection - the copy 165 166 Level: beginner 167 168 .seealso: `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()` 169 @*/ 170 PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection) 171 { 172 PetscFunctionBegin; 173 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 174 PetscValidPointer(newSection, 2); 175 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection)); 176 PetscCall(PetscSectionCopy(section, *newSection)); 177 PetscFunctionReturn(0); 178 } 179 180 /*@ 181 PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database 182 183 Collective 184 185 Input Parameter: 186 . section - the PetscSection 187 188 Options Database: 189 . -petscsection_point_major - PETSC_TRUE for point-major order 190 191 Level: intermediate 192 193 .seealso: `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()` 194 @*/ 195 PetscErrorCode PetscSectionSetFromOptions(PetscSection s) 196 { 197 PetscFunctionBegin; 198 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 199 PetscObjectOptionsBegin((PetscObject)s); 200 PetscCall(PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL)); 201 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 202 PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject)); 203 PetscOptionsEnd(); 204 PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view")); 205 PetscFunctionReturn(0); 206 } 207 208 /*@ 209 PetscSectionCompare - Compares two sections 210 211 Collective 212 213 Input Parameters: 214 + s1 - the first PetscSection 215 - s2 - the second PetscSection 216 217 Output Parameter: 218 . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise 219 220 Level: intermediate 221 222 Notes: 223 Field names are disregarded. 224 225 .seealso: `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()` 226 @*/ 227 PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent) 228 { 229 PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2; 230 const PetscInt *idx1, *idx2; 231 IS perm1, perm2; 232 PetscBool flg; 233 PetscMPIInt mflg; 234 235 PetscFunctionBegin; 236 PetscValidHeaderSpecific(s1, PETSC_SECTION_CLASSID, 1); 237 PetscValidHeaderSpecific(s2, PETSC_SECTION_CLASSID, 2); 238 PetscValidBoolPointer(congruent, 3); 239 flg = PETSC_FALSE; 240 241 PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg)); 242 if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) { 243 *congruent = PETSC_FALSE; 244 PetscFunctionReturn(0); 245 } 246 247 PetscCall(PetscSectionGetChart(s1, &pStart, &pEnd)); 248 PetscCall(PetscSectionGetChart(s2, &n1, &n2)); 249 if (pStart != n1 || pEnd != n2) goto not_congruent; 250 251 PetscCall(PetscSectionGetPermutation(s1, &perm1)); 252 PetscCall(PetscSectionGetPermutation(s2, &perm2)); 253 if (perm1 && perm2) { 254 PetscCall(ISEqual(perm1, perm2, congruent)); 255 if (!(*congruent)) goto not_congruent; 256 } else if (perm1 != perm2) goto not_congruent; 257 258 for (p = pStart; p < pEnd; ++p) { 259 PetscCall(PetscSectionGetOffset(s1, p, &n1)); 260 PetscCall(PetscSectionGetOffset(s2, p, &n2)); 261 if (n1 != n2) goto not_congruent; 262 263 PetscCall(PetscSectionGetDof(s1, p, &n1)); 264 PetscCall(PetscSectionGetDof(s2, p, &n2)); 265 if (n1 != n2) goto not_congruent; 266 267 PetscCall(PetscSectionGetConstraintDof(s1, p, &ncdof)); 268 PetscCall(PetscSectionGetConstraintDof(s2, p, &n2)); 269 if (ncdof != n2) goto not_congruent; 270 271 PetscCall(PetscSectionGetConstraintIndices(s1, p, &idx1)); 272 PetscCall(PetscSectionGetConstraintIndices(s2, p, &idx2)); 273 PetscCall(PetscArraycmp(idx1, idx2, ncdof, congruent)); 274 if (!(*congruent)) goto not_congruent; 275 } 276 277 PetscCall(PetscSectionGetNumFields(s1, &nfields)); 278 PetscCall(PetscSectionGetNumFields(s2, &n2)); 279 if (nfields != n2) goto not_congruent; 280 281 for (f = 0; f < nfields; ++f) { 282 PetscCall(PetscSectionGetFieldComponents(s1, f, &n1)); 283 PetscCall(PetscSectionGetFieldComponents(s2, f, &n2)); 284 if (n1 != n2) goto not_congruent; 285 286 for (p = pStart; p < pEnd; ++p) { 287 PetscCall(PetscSectionGetFieldOffset(s1, p, f, &n1)); 288 PetscCall(PetscSectionGetFieldOffset(s2, p, f, &n2)); 289 if (n1 != n2) goto not_congruent; 290 291 PetscCall(PetscSectionGetFieldDof(s1, p, f, &n1)); 292 PetscCall(PetscSectionGetFieldDof(s2, p, f, &n2)); 293 if (n1 != n2) goto not_congruent; 294 295 PetscCall(PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof)); 296 PetscCall(PetscSectionGetFieldConstraintDof(s2, p, f, &n2)); 297 if (nfcdof != n2) goto not_congruent; 298 299 PetscCall(PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1)); 300 PetscCall(PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2)); 301 PetscCall(PetscArraycmp(idx1, idx2, nfcdof, congruent)); 302 if (!(*congruent)) goto not_congruent; 303 } 304 } 305 306 flg = PETSC_TRUE; 307 not_congruent: 308 PetscCall(MPIU_Allreduce(&flg, congruent, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1))); 309 PetscFunctionReturn(0); 310 } 311 312 /*@ 313 PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined. 314 315 Not Collective 316 317 Input Parameter: 318 . s - the PetscSection 319 320 Output Parameter: 321 . numFields - the number of fields defined, or 0 if none were defined 322 323 Level: intermediate 324 325 .seealso: `PetscSectionSetNumFields()` 326 @*/ 327 PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields) 328 { 329 PetscFunctionBegin; 330 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 331 PetscValidIntPointer(numFields, 2); 332 *numFields = s->numFields; 333 PetscFunctionReturn(0); 334 } 335 336 /*@ 337 PetscSectionSetNumFields - Sets the number of fields. 338 339 Not Collective 340 341 Input Parameters: 342 + s - the PetscSection 343 - numFields - the number of fields 344 345 Level: intermediate 346 347 .seealso: `PetscSectionGetNumFields()` 348 @*/ 349 PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields) 350 { 351 PetscInt f; 352 353 PetscFunctionBegin; 354 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 355 PetscCheck(numFields > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %" PetscInt_FMT " must be positive", numFields); 356 PetscCall(PetscSectionReset(s)); 357 358 s->numFields = numFields; 359 PetscCall(PetscMalloc1(s->numFields, &s->numFieldComponents)); 360 PetscCall(PetscMalloc1(s->numFields, &s->fieldNames)); 361 PetscCall(PetscMalloc1(s->numFields, &s->compNames)); 362 PetscCall(PetscMalloc1(s->numFields, &s->field)); 363 for (f = 0; f < s->numFields; ++f) { 364 char name[64]; 365 366 s->numFieldComponents[f] = 1; 367 368 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f])); 369 PetscCall(PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f)); 370 PetscCall(PetscStrallocpy(name, (char **)&s->fieldNames[f])); 371 PetscCall(PetscSNPrintf(name, 64, "Component_0")); 372 PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f])); 373 PetscCall(PetscStrallocpy(name, (char **)&s->compNames[f][0])); 374 } 375 PetscFunctionReturn(0); 376 } 377 378 /*@C 379 PetscSectionGetFieldName - Returns the name of a field in the PetscSection 380 381 Not Collective 382 383 Input Parameters: 384 + s - the PetscSection 385 - field - the field number 386 387 Output Parameter: 388 . fieldName - the field name 389 390 Level: intermediate 391 392 .seealso: `PetscSectionSetFieldName()` 393 @*/ 394 PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[]) 395 { 396 PetscFunctionBegin; 397 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 398 PetscValidPointer(fieldName, 3); 399 PetscSectionCheckValidField(field, s->numFields); 400 *fieldName = s->fieldNames[field]; 401 PetscFunctionReturn(0); 402 } 403 404 /*@C 405 PetscSectionSetFieldName - Sets the name of a field in the PetscSection 406 407 Not Collective 408 409 Input Parameters: 410 + s - the PetscSection 411 . field - the field number 412 - fieldName - the field name 413 414 Level: intermediate 415 416 .seealso: `PetscSectionGetFieldName()` 417 @*/ 418 PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[]) 419 { 420 PetscFunctionBegin; 421 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 422 if (fieldName) PetscValidCharPointer(fieldName, 3); 423 PetscSectionCheckValidField(field, s->numFields); 424 PetscCall(PetscFree(s->fieldNames[field])); 425 PetscCall(PetscStrallocpy(fieldName, (char **)&s->fieldNames[field])); 426 PetscFunctionReturn(0); 427 } 428 429 /*@C 430 PetscSectionGetComponentName - Gets the name of a field component in the PetscSection 431 432 Not Collective 433 434 Input Parameters: 435 + s - the PetscSection 436 . field - the field number 437 . comp - the component number 438 - compName - the component name 439 440 Level: intermediate 441 442 .seealso: `PetscSectionSetComponentName()` 443 @*/ 444 PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[]) 445 { 446 PetscFunctionBegin; 447 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 448 PetscValidPointer(compName, 4); 449 PetscSectionCheckValidField(field, s->numFields); 450 PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]); 451 *compName = s->compNames[field][comp]; 452 PetscFunctionReturn(0); 453 } 454 455 /*@C 456 PetscSectionSetComponentName - Sets the name of a field component in the PetscSection 457 458 Not Collective 459 460 Input Parameters: 461 + s - the PetscSection 462 . field - the field number 463 . comp - the component number 464 - compName - the component name 465 466 Level: intermediate 467 468 .seealso: `PetscSectionGetComponentName()` 469 @*/ 470 PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[]) 471 { 472 PetscFunctionBegin; 473 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 474 if (compName) PetscValidCharPointer(compName, 4); 475 PetscSectionCheckValidField(field, s->numFields); 476 PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]); 477 PetscCall(PetscFree(s->compNames[field][comp])); 478 PetscCall(PetscStrallocpy(compName, (char **)&s->compNames[field][comp])); 479 PetscFunctionReturn(0); 480 } 481 482 /*@ 483 PetscSectionGetFieldComponents - Returns the number of field components for the given field. 484 485 Not Collective 486 487 Input Parameters: 488 + s - the PetscSection 489 - field - the field number 490 491 Output Parameter: 492 . numComp - the number of field components 493 494 Level: intermediate 495 496 .seealso: `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()` 497 @*/ 498 PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp) 499 { 500 PetscFunctionBegin; 501 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 502 PetscValidIntPointer(numComp, 3); 503 PetscSectionCheckValidField(field, s->numFields); 504 *numComp = s->numFieldComponents[field]; 505 PetscFunctionReturn(0); 506 } 507 508 /*@ 509 PetscSectionSetFieldComponents - Sets the number of field components for the given field. 510 511 Not Collective 512 513 Input Parameters: 514 + s - the PetscSection 515 . field - the field number 516 - numComp - the number of field components 517 518 Level: intermediate 519 520 .seealso: `PetscSectionGetFieldComponents()`, `PetscSectionGetNumFields()` 521 @*/ 522 PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp) 523 { 524 PetscInt c; 525 526 PetscFunctionBegin; 527 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 528 PetscSectionCheckValidField(field, s->numFields); 529 if (s->compNames) { 530 for (c = 0; c < s->numFieldComponents[field]; ++c) PetscCall(PetscFree(s->compNames[field][c])); 531 PetscCall(PetscFree(s->compNames[field])); 532 } 533 534 s->numFieldComponents[field] = numComp; 535 if (numComp) { 536 PetscCall(PetscMalloc1(numComp, (char ***)&s->compNames[field])); 537 for (c = 0; c < numComp; ++c) { 538 char name[64]; 539 540 PetscCall(PetscSNPrintf(name, 64, "%" PetscInt_FMT, c)); 541 PetscCall(PetscStrallocpy(name, (char **)&s->compNames[field][c])); 542 } 543 } 544 PetscFunctionReturn(0); 545 } 546 547 /*@ 548 PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie. 549 550 Not Collective 551 552 Input Parameter: 553 . s - the PetscSection 554 555 Output Parameters: 556 + pStart - the first point 557 - pEnd - one past the last point 558 559 Level: intermediate 560 561 .seealso: `PetscSectionSetChart()`, `PetscSectionCreate()` 562 @*/ 563 PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd) 564 { 565 PetscFunctionBegin; 566 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 567 if (pStart) *pStart = s->pStart; 568 if (pEnd) *pEnd = s->pEnd; 569 PetscFunctionReturn(0); 570 } 571 572 /*@ 573 PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie. 574 575 Not Collective 576 577 Input Parameters: 578 + s - the PetscSection 579 . pStart - the first point 580 - pEnd - one past the last point 581 582 Level: intermediate 583 584 .seealso: `PetscSectionGetChart()`, `PetscSectionCreate()` 585 @*/ 586 PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd) 587 { 588 PetscInt f; 589 590 PetscFunctionBegin; 591 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 592 if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(0); 593 /* Cannot Reset() because it destroys field information */ 594 s->setup = PETSC_FALSE; 595 PetscCall(PetscSectionDestroy(&s->bc)); 596 PetscCall(PetscFree(s->bcIndices)); 597 PetscCall(PetscFree2(s->atlasDof, s->atlasOff)); 598 599 s->pStart = pStart; 600 s->pEnd = pEnd; 601 PetscCall(PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff)); 602 PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart)); 603 for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd)); 604 PetscFunctionReturn(0); 605 } 606 607 /*@ 608 PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL 609 610 Not Collective 611 612 Input Parameter: 613 . s - the PetscSection 614 615 Output Parameters: 616 . perm - The permutation as an IS 617 618 Level: intermediate 619 620 .seealso: `PetscSectionSetPermutation()`, `PetscSectionCreate()` 621 @*/ 622 PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm) 623 { 624 PetscFunctionBegin; 625 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 626 if (perm) { 627 PetscValidPointer(perm, 2); 628 *perm = s->perm; 629 } 630 PetscFunctionReturn(0); 631 } 632 633 /*@ 634 PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart) 635 636 Not Collective 637 638 Input Parameters: 639 + s - the PetscSection 640 - perm - the permutation of points 641 642 Level: intermediate 643 644 .seealso: `PetscSectionGetPermutation()`, `PetscSectionCreate()` 645 @*/ 646 PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm) 647 { 648 PetscFunctionBegin; 649 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 650 if (perm) PetscValidHeaderSpecific(perm, IS_CLASSID, 2); 651 PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup"); 652 if (s->perm != perm) { 653 PetscCall(ISDestroy(&s->perm)); 654 if (perm) { 655 s->perm = perm; 656 PetscCall(PetscObjectReference((PetscObject)s->perm)); 657 } 658 } 659 PetscFunctionReturn(0); 660 } 661 662 /*@ 663 PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major 664 665 Not Collective 666 667 Input Parameter: 668 . s - the PetscSection 669 670 Output Parameter: 671 . pm - the flag for point major ordering 672 673 Level: intermediate 674 675 .seealso: `PetscSectionSetPointMajor()` 676 @*/ 677 PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm) 678 { 679 PetscFunctionBegin; 680 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 681 PetscValidBoolPointer(pm, 2); 682 *pm = s->pointMajor; 683 PetscFunctionReturn(0); 684 } 685 686 /*@ 687 PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major 688 689 Not Collective 690 691 Input Parameters: 692 + s - the PetscSection 693 - pm - the flag for point major ordering 694 695 Level: intermediate 696 697 .seealso: `PetscSectionGetPointMajor()` 698 @*/ 699 PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm) 700 { 701 PetscFunctionBegin; 702 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 703 PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup"); 704 s->pointMajor = pm; 705 PetscFunctionReturn(0); 706 } 707 708 /*@ 709 PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets 710 711 Not Collective 712 713 Input Parameter: 714 . s - the PetscSection 715 716 Output Parameter: 717 . includesConstraints - the flag indicating if constrained dofs were included when computing offsets 718 719 Level: intermediate 720 721 .seealso: `PetscSectionSetIncludesConstraints()` 722 @*/ 723 PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints) 724 { 725 PetscFunctionBegin; 726 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 727 PetscValidBoolPointer(includesConstraints, 2); 728 *includesConstraints = s->includesConstraints; 729 PetscFunctionReturn(0); 730 } 731 732 /*@ 733 PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets 734 735 Not Collective 736 737 Input Parameters: 738 + s - the PetscSection 739 - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets 740 741 Level: intermediate 742 743 .seealso: `PetscSectionGetIncludesConstraints()` 744 @*/ 745 PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints) 746 { 747 PetscFunctionBegin; 748 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 749 PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up"); 750 s->includesConstraints = includesConstraints; 751 PetscFunctionReturn(0); 752 } 753 754 /*@ 755 PetscSectionGetDof - Return the number of degrees of freedom associated with a given point. 756 757 Not Collective 758 759 Input Parameters: 760 + s - the PetscSection 761 - point - the point 762 763 Output Parameter: 764 . numDof - the number of dof 765 766 Level: intermediate 767 768 .seealso: `PetscSectionSetDof()`, `PetscSectionCreate()` 769 @*/ 770 PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof) 771 { 772 PetscFunctionBeginHot; 773 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 774 PetscValidIntPointer(numDof, 3); 775 PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd); 776 *numDof = s->atlasDof[point - s->pStart]; 777 PetscFunctionReturn(0); 778 } 779 780 /*@ 781 PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point. 782 783 Not Collective 784 785 Input Parameters: 786 + s - the PetscSection 787 . point - the point 788 - numDof - the number of dof 789 790 Level: intermediate 791 792 .seealso: `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()` 793 @*/ 794 PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof) 795 { 796 PetscFunctionBegin; 797 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 798 PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd); 799 s->atlasDof[point - s->pStart] = numDof; 800 PetscCall(PetscSectionInvalidateMaxDof_Internal(s)); 801 PetscFunctionReturn(0); 802 } 803 804 /*@ 805 PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point. 806 807 Not Collective 808 809 Input Parameters: 810 + s - the PetscSection 811 . point - the point 812 - numDof - the number of additional dof 813 814 Level: intermediate 815 816 .seealso: `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()` 817 @*/ 818 PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof) 819 { 820 PetscFunctionBeginHot; 821 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 822 PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd); 823 s->atlasDof[point - s->pStart] += numDof; 824 PetscCall(PetscSectionInvalidateMaxDof_Internal(s)); 825 PetscFunctionReturn(0); 826 } 827 828 /*@ 829 PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point. 830 831 Not Collective 832 833 Input Parameters: 834 + s - the PetscSection 835 . point - the point 836 - field - the field 837 838 Output Parameter: 839 . numDof - the number of dof 840 841 Level: intermediate 842 843 .seealso: `PetscSectionSetFieldDof()`, `PetscSectionCreate()` 844 @*/ 845 PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof) 846 { 847 PetscFunctionBegin; 848 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 849 PetscValidIntPointer(numDof, 4); 850 PetscSectionCheckValidField(field, s->numFields); 851 PetscCall(PetscSectionGetDof(s->field[field], point, numDof)); 852 PetscFunctionReturn(0); 853 } 854 855 /*@ 856 PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point. 857 858 Not Collective 859 860 Input Parameters: 861 + s - the PetscSection 862 . point - the point 863 . field - the field 864 - numDof - the number of dof 865 866 Level: intermediate 867 868 .seealso: `PetscSectionGetFieldDof()`, `PetscSectionCreate()` 869 @*/ 870 PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof) 871 { 872 PetscFunctionBegin; 873 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 874 PetscSectionCheckValidField(field, s->numFields); 875 PetscCall(PetscSectionSetDof(s->field[field], point, numDof)); 876 PetscFunctionReturn(0); 877 } 878 879 /*@ 880 PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point. 881 882 Not Collective 883 884 Input Parameters: 885 + s - the PetscSection 886 . point - the point 887 . field - the field 888 - numDof - the number of dof 889 890 Level: intermediate 891 892 .seealso: `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()` 893 @*/ 894 PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof) 895 { 896 PetscFunctionBegin; 897 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 898 PetscSectionCheckValidField(field, s->numFields); 899 PetscCall(PetscSectionAddDof(s->field[field], point, numDof)); 900 PetscFunctionReturn(0); 901 } 902 903 /*@ 904 PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point. 905 906 Not Collective 907 908 Input Parameters: 909 + s - the PetscSection 910 - point - the point 911 912 Output Parameter: 913 . numDof - the number of dof which are fixed by constraints 914 915 Level: intermediate 916 917 .seealso: `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()` 918 @*/ 919 PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof) 920 { 921 PetscFunctionBegin; 922 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 923 PetscValidIntPointer(numDof, 3); 924 if (s->bc) { 925 PetscCall(PetscSectionGetDof(s->bc, point, numDof)); 926 } else *numDof = 0; 927 PetscFunctionReturn(0); 928 } 929 930 /*@ 931 PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point. 932 933 Not Collective 934 935 Input Parameters: 936 + s - the PetscSection 937 . point - the point 938 - numDof - the number of dof which are fixed by constraints 939 940 Level: intermediate 941 942 .seealso: `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()` 943 @*/ 944 PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof) 945 { 946 PetscFunctionBegin; 947 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 948 if (numDof) { 949 PetscCall(PetscSectionCheckConstraints_Private(s)); 950 PetscCall(PetscSectionSetDof(s->bc, point, numDof)); 951 } 952 PetscFunctionReturn(0); 953 } 954 955 /*@ 956 PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point. 957 958 Not Collective 959 960 Input Parameters: 961 + s - the PetscSection 962 . point - the point 963 - numDof - the number of additional dof which are fixed by constraints 964 965 Level: intermediate 966 967 .seealso: `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()` 968 @*/ 969 PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof) 970 { 971 PetscFunctionBegin; 972 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 973 if (numDof) { 974 PetscCall(PetscSectionCheckConstraints_Private(s)); 975 PetscCall(PetscSectionAddDof(s->bc, point, numDof)); 976 } 977 PetscFunctionReturn(0); 978 } 979 980 /*@ 981 PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point. 982 983 Not Collective 984 985 Input Parameters: 986 + s - the PetscSection 987 . point - the point 988 - field - the field 989 990 Output Parameter: 991 . numDof - the number of dof which are fixed by constraints 992 993 Level: intermediate 994 995 .seealso: `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()` 996 @*/ 997 PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof) 998 { 999 PetscFunctionBegin; 1000 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1001 PetscValidIntPointer(numDof, 4); 1002 PetscSectionCheckValidField(field, s->numFields); 1003 PetscCall(PetscSectionGetConstraintDof(s->field[field], point, numDof)); 1004 PetscFunctionReturn(0); 1005 } 1006 1007 /*@ 1008 PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point. 1009 1010 Not Collective 1011 1012 Input Parameters: 1013 + s - the PetscSection 1014 . point - the point 1015 . field - the field 1016 - numDof - the number of dof which are fixed by constraints 1017 1018 Level: intermediate 1019 1020 .seealso: `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()` 1021 @*/ 1022 PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof) 1023 { 1024 PetscFunctionBegin; 1025 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1026 PetscSectionCheckValidField(field, s->numFields); 1027 PetscCall(PetscSectionSetConstraintDof(s->field[field], point, numDof)); 1028 PetscFunctionReturn(0); 1029 } 1030 1031 /*@ 1032 PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point. 1033 1034 Not Collective 1035 1036 Input Parameters: 1037 + s - the PetscSection 1038 . point - the point 1039 . field - the field 1040 - numDof - the number of additional dof which are fixed by constraints 1041 1042 Level: intermediate 1043 1044 .seealso: `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()` 1045 @*/ 1046 PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof) 1047 { 1048 PetscFunctionBegin; 1049 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1050 PetscSectionCheckValidField(field, s->numFields); 1051 PetscCall(PetscSectionAddConstraintDof(s->field[field], point, numDof)); 1052 PetscFunctionReturn(0); 1053 } 1054 1055 /*@ 1056 PetscSectionSetUpBC - Setup the subsections describing boundary conditions. 1057 1058 Not Collective 1059 1060 Input Parameter: 1061 . s - the PetscSection 1062 1063 Level: advanced 1064 1065 .seealso: `PetscSectionSetUp()`, `PetscSectionCreate()` 1066 @*/ 1067 PetscErrorCode PetscSectionSetUpBC(PetscSection s) 1068 { 1069 PetscFunctionBegin; 1070 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1071 if (s->bc) { 1072 const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1; 1073 1074 PetscCall(PetscSectionSetUp(s->bc)); 1075 PetscCall(PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices)); 1076 } 1077 PetscFunctionReturn(0); 1078 } 1079 1080 /*@ 1081 PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point. 1082 1083 Not Collective 1084 1085 Input Parameter: 1086 . s - the PetscSection 1087 1088 Level: intermediate 1089 1090 .seealso: `PetscSectionCreate()` 1091 @*/ 1092 PetscErrorCode PetscSectionSetUp(PetscSection s) 1093 { 1094 const PetscInt *pind = NULL; 1095 PetscInt offset = 0, foff, p, f; 1096 1097 PetscFunctionBegin; 1098 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1099 if (s->setup) PetscFunctionReturn(0); 1100 s->setup = PETSC_TRUE; 1101 /* Set offsets and field offsets for all points */ 1102 /* Assume that all fields have the same chart */ 1103 PetscCheck(s->includesConstraints, PETSC_COMM_SELF, PETSC_ERR_SUP, "PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE"); 1104 if (s->perm) PetscCall(ISGetIndices(s->perm, &pind)); 1105 if (s->pointMajor) { 1106 for (p = 0; p < s->pEnd - s->pStart; ++p) { 1107 const PetscInt q = pind ? pind[p] : p; 1108 1109 /* Set point offset */ 1110 s->atlasOff[q] = offset; 1111 offset += s->atlasDof[q]; 1112 /* Set field offset */ 1113 for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) { 1114 PetscSection sf = s->field[f]; 1115 1116 sf->atlasOff[q] = foff; 1117 foff += sf->atlasDof[q]; 1118 } 1119 } 1120 } else { 1121 /* Set field offsets for all points */ 1122 for (f = 0; f < s->numFields; ++f) { 1123 PetscSection sf = s->field[f]; 1124 1125 for (p = 0; p < s->pEnd - s->pStart; ++p) { 1126 const PetscInt q = pind ? pind[p] : p; 1127 1128 sf->atlasOff[q] = offset; 1129 offset += sf->atlasDof[q]; 1130 } 1131 } 1132 /* Disable point offsets since these are unused */ 1133 for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1; 1134 } 1135 if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind)); 1136 /* Setup BC sections */ 1137 PetscCall(PetscSectionSetUpBC(s)); 1138 for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetUpBC(s->field[f])); 1139 PetscFunctionReturn(0); 1140 } 1141 1142 /*@ 1143 PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart 1144 1145 Not Collective 1146 1147 Input Parameters: 1148 . s - the PetscSection 1149 1150 Output Parameter: 1151 . maxDof - the maximum dof 1152 1153 Level: intermediate 1154 1155 Notes: 1156 The returned number is up-to-date without need for PetscSectionSetUp(). 1157 1158 Developer Notes: 1159 The returned number is calculated lazily and stashed. 1160 A call to PetscSectionInvalidateMaxDof_Internal() invalidates the stashed value. 1161 PetscSectionInvalidateMaxDof_Internal() is called in PetscSectionSetDof(), PetscSectionAddDof() and PetscSectionReset(). 1162 It should also be called every time atlasDof is modified directly. 1163 1164 .seealso: `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()` 1165 @*/ 1166 PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof) 1167 { 1168 PetscInt p; 1169 1170 PetscFunctionBegin; 1171 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1172 PetscValidIntPointer(maxDof, 2); 1173 if (s->maxDof == PETSC_MIN_INT) { 1174 s->maxDof = 0; 1175 for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]); 1176 } 1177 *maxDof = s->maxDof; 1178 PetscFunctionReturn(0); 1179 } 1180 1181 /*@ 1182 PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom. 1183 1184 Not Collective 1185 1186 Input Parameter: 1187 . s - the PetscSection 1188 1189 Output Parameter: 1190 . size - the size of an array which can hold all the dofs 1191 1192 Level: intermediate 1193 1194 .seealso: `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()` 1195 @*/ 1196 PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size) 1197 { 1198 PetscInt p, n = 0; 1199 1200 PetscFunctionBegin; 1201 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1202 PetscValidIntPointer(size, 2); 1203 for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0; 1204 *size = n; 1205 PetscFunctionReturn(0); 1206 } 1207 1208 /*@ 1209 PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom. 1210 1211 Not Collective 1212 1213 Input Parameter: 1214 . s - the PetscSection 1215 1216 Output Parameter: 1217 . size - the size of an array which can hold all unconstrained dofs 1218 1219 Level: intermediate 1220 1221 .seealso: `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()` 1222 @*/ 1223 PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size) 1224 { 1225 PetscInt p, n = 0; 1226 1227 PetscFunctionBegin; 1228 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1229 PetscValidIntPointer(size, 2); 1230 for (p = 0; p < s->pEnd - s->pStart; ++p) { 1231 const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0; 1232 n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0; 1233 } 1234 *size = n; 1235 PetscFunctionReturn(0); 1236 } 1237 1238 /*@ 1239 PetscSectionCreateGlobalSection - Create a section describing the global field layout using 1240 the local section and an SF describing the section point overlap. 1241 1242 Input Parameters: 1243 + s - The PetscSection for the local field layout 1244 . sf - The SF describing parallel layout of the section points (leaves are unowned local points) 1245 . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs 1246 - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points 1247 1248 Output Parameter: 1249 . gsection - The PetscSection for the global field layout 1250 1251 Note: This gives negative sizes and offsets to points not owned by this process 1252 1253 Level: intermediate 1254 1255 .seealso: `PetscSectionCreate()` 1256 @*/ 1257 PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection) 1258 { 1259 PetscSection gs; 1260 const PetscInt *pind = NULL; 1261 PetscInt *recv = NULL, *neg = NULL; 1262 PetscInt pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf; 1263 PetscInt numFields, f, numComponents; 1264 1265 PetscFunctionBegin; 1266 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1267 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 1268 PetscValidLogicalCollectiveBool(s, includeConstraints, 3); 1269 PetscValidLogicalCollectiveBool(s, localOffsets, 4); 1270 PetscValidPointer(gsection, 5); 1271 PetscCheck(s->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering"); 1272 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs)); 1273 PetscCall(PetscSectionGetNumFields(s, &numFields)); 1274 if (numFields > 0) PetscCall(PetscSectionSetNumFields(gs, numFields)); 1275 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1276 PetscCall(PetscSectionSetChart(gs, pStart, pEnd)); 1277 gs->includesConstraints = includeConstraints; 1278 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 1279 nlocal = nroots; /* The local/leaf space matches global/root space */ 1280 /* Must allocate for all points visible to SF, which may be more than this section */ 1281 if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */ 1282 PetscCall(PetscSFGetLeafRange(sf, NULL, &maxleaf)); 1283 PetscCheck(nroots >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %" PetscInt_FMT " < pEnd %" PetscInt_FMT, nroots, pEnd); 1284 PetscCheck(maxleaf < nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %" PetscInt_FMT " >= nroots %" PetscInt_FMT, maxleaf, nroots); 1285 PetscCall(PetscMalloc2(nroots, &neg, nlocal, &recv)); 1286 PetscCall(PetscArrayzero(neg, nroots)); 1287 } 1288 /* Mark all local points with negative dof */ 1289 for (p = pStart; p < pEnd; ++p) { 1290 PetscCall(PetscSectionGetDof(s, p, &dof)); 1291 PetscCall(PetscSectionSetDof(gs, p, dof)); 1292 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 1293 if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(gs, p, cdof)); 1294 if (neg) neg[p] = -(dof + 1); 1295 } 1296 PetscCall(PetscSectionSetUpBC(gs)); 1297 if (gs->bcIndices) PetscCall(PetscArraycpy(gs->bcIndices, s->bcIndices, gs->bc->atlasOff[gs->bc->pEnd - gs->bc->pStart - 1] + gs->bc->atlasDof[gs->bc->pEnd - gs->bc->pStart - 1])); 1298 if (nroots >= 0) { 1299 PetscCall(PetscArrayzero(recv, nlocal)); 1300 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE)); 1301 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE)); 1302 for (p = pStart; p < pEnd; ++p) { 1303 if (recv[p] < 0) { 1304 gs->atlasDof[p - pStart] = recv[p]; 1305 PetscCall(PetscSectionGetDof(s, p, &dof)); 1306 PetscCheck(-(recv[p] + 1) == dof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %" PetscInt_FMT " for point %" PetscInt_FMT " is not the unconstrained %" PetscInt_FMT, -(recv[p] + 1), p, dof); 1307 } 1308 } 1309 } 1310 /* Calculate new sizes, get process offset, and calculate point offsets */ 1311 if (s->perm) PetscCall(ISGetIndices(s->perm, &pind)); 1312 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 1313 const PetscInt q = pind ? pind[p] : p; 1314 1315 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0; 1316 gs->atlasOff[q] = off; 1317 off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0; 1318 } 1319 if (!localOffsets) { 1320 PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s))); 1321 globalOff -= off; 1322 } 1323 for (p = pStart, off = 0; p < pEnd; ++p) { 1324 gs->atlasOff[p - pStart] += globalOff; 1325 if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1); 1326 } 1327 if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind)); 1328 /* Put in negative offsets for ghost points */ 1329 if (nroots >= 0) { 1330 PetscCall(PetscArrayzero(recv, nlocal)); 1331 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE)); 1332 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE)); 1333 for (p = pStart; p < pEnd; ++p) { 1334 if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p]; 1335 } 1336 } 1337 PetscCall(PetscFree2(neg, recv)); 1338 /* Set field dofs/offsets/constraints */ 1339 for (f = 0; f < numFields; ++f) { 1340 gs->field[f]->includesConstraints = includeConstraints; 1341 PetscCall(PetscSectionGetFieldComponents(s, f, &numComponents)); 1342 PetscCall(PetscSectionSetFieldComponents(gs, f, numComponents)); 1343 } 1344 for (p = pStart; p < pEnd; ++p) { 1345 PetscCall(PetscSectionGetOffset(gs, p, &off)); 1346 for (f = 0, foff = off; f < numFields; ++f) { 1347 PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof)); 1348 if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetFieldConstraintDof(gs, p, f, cdof)); 1349 PetscCall(PetscSectionGetFieldDof(s, p, f, &dof)); 1350 PetscCall(PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof)); 1351 PetscCall(PetscSectionSetFieldOffset(gs, p, f, foff)); 1352 PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof)); 1353 foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof); 1354 } 1355 } 1356 for (f = 0; f < numFields; ++f) { 1357 PetscSection gfs = gs->field[f]; 1358 1359 PetscCall(PetscSectionSetUpBC(gfs)); 1360 if (gfs->bcIndices) PetscCall(PetscArraycpy(gfs->bcIndices, s->field[f]->bcIndices, gfs->bc->atlasOff[gfs->bc->pEnd - gfs->bc->pStart - 1] + gfs->bc->atlasDof[gfs->bc->pEnd - gfs->bc->pStart - 1])); 1361 } 1362 gs->setup = PETSC_TRUE; 1363 PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view")); 1364 *gsection = gs; 1365 PetscFunctionReturn(0); 1366 } 1367 1368 /*@ 1369 PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using 1370 the local section and an SF describing the section point overlap. 1371 1372 Input Parameters: 1373 + s - The PetscSection for the local field layout 1374 . sf - The SF describing parallel layout of the section points 1375 . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs 1376 . numExcludes - The number of exclusion ranges 1377 - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs 1378 1379 Output Parameter: 1380 . gsection - The PetscSection for the global field layout 1381 1382 Note: This gives negative sizes and offsets to points not owned by this process 1383 1384 Level: advanced 1385 1386 .seealso: `PetscSectionCreate()` 1387 @*/ 1388 PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection) 1389 { 1390 const PetscInt *pind = NULL; 1391 PetscInt *neg = NULL, *tmpOff = NULL; 1392 PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots; 1393 1394 PetscFunctionBegin; 1395 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1396 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 1397 PetscValidPointer(gsection, 6); 1398 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection)); 1399 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1400 PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd)); 1401 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 1402 if (nroots >= 0) { 1403 PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart); 1404 PetscCall(PetscCalloc1(nroots, &neg)); 1405 if (nroots > pEnd - pStart) { 1406 PetscCall(PetscCalloc1(nroots, &tmpOff)); 1407 } else { 1408 tmpOff = &(*gsection)->atlasDof[-pStart]; 1409 } 1410 } 1411 /* Mark ghost points with negative dof */ 1412 for (p = pStart; p < pEnd; ++p) { 1413 for (e = 0; e < numExcludes; ++e) { 1414 if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) { 1415 PetscCall(PetscSectionSetDof(*gsection, p, 0)); 1416 break; 1417 } 1418 } 1419 if (e < numExcludes) continue; 1420 PetscCall(PetscSectionGetDof(s, p, &dof)); 1421 PetscCall(PetscSectionSetDof(*gsection, p, dof)); 1422 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 1423 if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof)); 1424 if (neg) neg[p] = -(dof + 1); 1425 } 1426 PetscCall(PetscSectionSetUpBC(*gsection)); 1427 if (nroots >= 0) { 1428 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 1429 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 1430 if (nroots > pEnd - pStart) { 1431 for (p = pStart; p < pEnd; ++p) { 1432 if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p]; 1433 } 1434 } 1435 } 1436 /* Calculate new sizes, get proccess offset, and calculate point offsets */ 1437 if (s->perm) PetscCall(ISGetIndices(s->perm, &pind)); 1438 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 1439 const PetscInt q = pind ? pind[p] : p; 1440 1441 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0; 1442 (*gsection)->atlasOff[q] = off; 1443 off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0; 1444 } 1445 PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s))); 1446 globalOff -= off; 1447 for (p = 0, off = 0; p < pEnd - pStart; ++p) { 1448 (*gsection)->atlasOff[p] += globalOff; 1449 if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1); 1450 } 1451 if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind)); 1452 /* Put in negative offsets for ghost points */ 1453 if (nroots >= 0) { 1454 if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart]; 1455 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 1456 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE)); 1457 if (nroots > pEnd - pStart) { 1458 for (p = pStart; p < pEnd; ++p) { 1459 if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p]; 1460 } 1461 } 1462 } 1463 if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff)); 1464 PetscCall(PetscFree(neg)); 1465 PetscFunctionReturn(0); 1466 } 1467 1468 /*@ 1469 PetscSectionGetPointLayout - Get the PetscLayout associated with the section points 1470 1471 Collective on comm 1472 1473 Input Parameters: 1474 + comm - The MPI_Comm 1475 - s - The PetscSection 1476 1477 Output Parameter: 1478 . layout - The point layout for the section 1479 1480 Note: This is usually called for the default global section. 1481 1482 Level: advanced 1483 1484 .seealso: `PetscSectionGetValueLayout()`, `PetscSectionCreate()` 1485 @*/ 1486 PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout) 1487 { 1488 PetscInt pStart, pEnd, p, localSize = 0; 1489 1490 PetscFunctionBegin; 1491 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1492 for (p = pStart; p < pEnd; ++p) { 1493 PetscInt dof; 1494 1495 PetscCall(PetscSectionGetDof(s, p, &dof)); 1496 if (dof >= 0) ++localSize; 1497 } 1498 PetscCall(PetscLayoutCreate(comm, layout)); 1499 PetscCall(PetscLayoutSetLocalSize(*layout, localSize)); 1500 PetscCall(PetscLayoutSetBlockSize(*layout, 1)); 1501 PetscCall(PetscLayoutSetUp(*layout)); 1502 PetscFunctionReturn(0); 1503 } 1504 1505 /*@ 1506 PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs. 1507 1508 Collective on comm 1509 1510 Input Parameters: 1511 + comm - The MPI_Comm 1512 - s - The PetscSection 1513 1514 Output Parameter: 1515 . layout - The dof layout for the section 1516 1517 Note: This is usually called for the default global section. 1518 1519 Level: advanced 1520 1521 .seealso: `PetscSectionGetPointLayout()`, `PetscSectionCreate()` 1522 @*/ 1523 PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout) 1524 { 1525 PetscInt pStart, pEnd, p, localSize = 0; 1526 1527 PetscFunctionBegin; 1528 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2); 1529 PetscValidPointer(layout, 3); 1530 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1531 for (p = pStart; p < pEnd; ++p) { 1532 PetscInt dof, cdof; 1533 1534 PetscCall(PetscSectionGetDof(s, p, &dof)); 1535 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 1536 if (dof - cdof > 0) localSize += dof - cdof; 1537 } 1538 PetscCall(PetscLayoutCreate(comm, layout)); 1539 PetscCall(PetscLayoutSetLocalSize(*layout, localSize)); 1540 PetscCall(PetscLayoutSetBlockSize(*layout, 1)); 1541 PetscCall(PetscLayoutSetUp(*layout)); 1542 PetscFunctionReturn(0); 1543 } 1544 1545 /*@ 1546 PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point. 1547 1548 Not Collective 1549 1550 Input Parameters: 1551 + s - the PetscSection 1552 - point - the point 1553 1554 Output Parameter: 1555 . offset - the offset 1556 1557 Level: intermediate 1558 1559 .seealso: `PetscSectionGetFieldOffset()`, `PetscSectionCreate()` 1560 @*/ 1561 PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset) 1562 { 1563 PetscFunctionBegin; 1564 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1565 PetscValidIntPointer(offset, 3); 1566 if (PetscDefined(USE_DEBUG)) PetscCheck(!(point < s->pStart) && !(point >= s->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd); 1567 *offset = s->atlasOff[point - s->pStart]; 1568 PetscFunctionReturn(0); 1569 } 1570 1571 /*@ 1572 PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point. 1573 1574 Not Collective 1575 1576 Input Parameters: 1577 + s - the PetscSection 1578 . point - the point 1579 - offset - the offset 1580 1581 Note: The user usually does not call this function, but uses PetscSectionSetUp() 1582 1583 Level: intermediate 1584 1585 .seealso: `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()` 1586 @*/ 1587 PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset) 1588 { 1589 PetscFunctionBegin; 1590 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1591 PetscCheck(!(point < s->pStart) && !(point >= s->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd); 1592 s->atlasOff[point - s->pStart] = offset; 1593 PetscFunctionReturn(0); 1594 } 1595 1596 /*@ 1597 PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point. 1598 1599 Not Collective 1600 1601 Input Parameters: 1602 + s - the PetscSection 1603 . point - the point 1604 - field - the field 1605 1606 Output Parameter: 1607 . offset - the offset 1608 1609 Level: intermediate 1610 1611 .seealso: `PetscSectionGetOffset()`, `PetscSectionCreate()` 1612 @*/ 1613 PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset) 1614 { 1615 PetscFunctionBegin; 1616 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1617 PetscValidIntPointer(offset, 4); 1618 PetscSectionCheckValidField(field, s->numFields); 1619 PetscCall(PetscSectionGetOffset(s->field[field], point, offset)); 1620 PetscFunctionReturn(0); 1621 } 1622 1623 /*@ 1624 PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given field at a point. 1625 1626 Not Collective 1627 1628 Input Parameters: 1629 + s - the PetscSection 1630 . point - the point 1631 . field - the field 1632 - offset - the offset 1633 1634 Note: The user usually does not call this function, but uses PetscSectionSetUp() 1635 1636 Level: intermediate 1637 1638 .seealso: `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()` 1639 @*/ 1640 PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset) 1641 { 1642 PetscFunctionBegin; 1643 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1644 PetscSectionCheckValidField(field, s->numFields); 1645 PetscCall(PetscSectionSetOffset(s->field[field], point, offset)); 1646 PetscFunctionReturn(0); 1647 } 1648 1649 /*@ 1650 PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point. 1651 1652 Not Collective 1653 1654 Input Parameters: 1655 + s - the PetscSection 1656 . point - the point 1657 - field - the field 1658 1659 Output Parameter: 1660 . offset - the offset 1661 1662 Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for 1663 this point, what number is the first dof with this field. 1664 1665 Level: advanced 1666 1667 .seealso: `PetscSectionGetOffset()`, `PetscSectionCreate()` 1668 @*/ 1669 PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset) 1670 { 1671 PetscInt off, foff; 1672 1673 PetscFunctionBegin; 1674 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1675 PetscValidIntPointer(offset, 4); 1676 PetscSectionCheckValidField(field, s->numFields); 1677 PetscCall(PetscSectionGetOffset(s, point, &off)); 1678 PetscCall(PetscSectionGetOffset(s->field[field], point, &foff)); 1679 *offset = foff - off; 1680 PetscFunctionReturn(0); 1681 } 1682 1683 /*@ 1684 PetscSectionGetOffsetRange - Return the full range of offsets [start, end) 1685 1686 Not Collective 1687 1688 Input Parameter: 1689 . s - the PetscSection 1690 1691 Output Parameters: 1692 + start - the minimum offset 1693 - end - one more than the maximum offset 1694 1695 Level: intermediate 1696 1697 .seealso: `PetscSectionGetOffset()`, `PetscSectionCreate()` 1698 @*/ 1699 PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end) 1700 { 1701 PetscInt os = 0, oe = 0, pStart, pEnd, p; 1702 1703 PetscFunctionBegin; 1704 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1705 if (s->atlasOff) { 1706 os = s->atlasOff[0]; 1707 oe = s->atlasOff[0]; 1708 } 1709 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1710 for (p = 0; p < pEnd - pStart; ++p) { 1711 PetscInt dof = s->atlasDof[p], off = s->atlasOff[p]; 1712 1713 if (off >= 0) { 1714 os = PetscMin(os, off); 1715 oe = PetscMax(oe, off + dof); 1716 } 1717 } 1718 if (start) *start = os; 1719 if (end) *end = oe; 1720 PetscFunctionReturn(0); 1721 } 1722 1723 /*@ 1724 PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields 1725 1726 Collective 1727 1728 Input Parameters: 1729 + s - the PetscSection 1730 . len - the number of subfields 1731 - fields - the subfield numbers 1732 1733 Output Parameter: 1734 . subs - the subsection 1735 1736 Note: The section offsets now refer to a new, smaller vector. 1737 1738 Level: advanced 1739 1740 .seealso: `PetscSectionCreateSupersection()`, `PetscSectionCreate()` 1741 @*/ 1742 PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs) 1743 { 1744 PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0; 1745 1746 PetscFunctionBegin; 1747 if (!len) PetscFunctionReturn(0); 1748 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1749 PetscValidIntPointer(fields, 3); 1750 PetscValidPointer(subs, 4); 1751 PetscCall(PetscSectionGetNumFields(s, &nF)); 1752 PetscCheck(len <= nF, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "Number of requested fields %" PetscInt_FMT " greater than number of fields %" PetscInt_FMT, len, nF); 1753 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs)); 1754 PetscCall(PetscSectionSetNumFields(*subs, len)); 1755 for (f = 0; f < len; ++f) { 1756 const char *name = NULL; 1757 PetscInt numComp = 0; 1758 1759 PetscCall(PetscSectionGetFieldName(s, fields[f], &name)); 1760 PetscCall(PetscSectionSetFieldName(*subs, f, name)); 1761 PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp)); 1762 PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp)); 1763 for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) { 1764 PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name)); 1765 PetscCall(PetscSectionSetComponentName(*subs, f, c, name)); 1766 } 1767 } 1768 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1769 PetscCall(PetscSectionSetChart(*subs, pStart, pEnd)); 1770 for (p = pStart; p < pEnd; ++p) { 1771 PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0; 1772 1773 for (f = 0; f < len; ++f) { 1774 PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof)); 1775 PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof)); 1776 PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof)); 1777 if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof)); 1778 dof += fdof; 1779 cdof += cfdof; 1780 } 1781 PetscCall(PetscSectionSetDof(*subs, p, dof)); 1782 if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof)); 1783 maxCdof = PetscMax(cdof, maxCdof); 1784 } 1785 PetscCall(PetscSectionSetUp(*subs)); 1786 if (maxCdof) { 1787 PetscInt *indices; 1788 1789 PetscCall(PetscMalloc1(maxCdof, &indices)); 1790 for (p = pStart; p < pEnd; ++p) { 1791 PetscInt cdof; 1792 1793 PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof)); 1794 if (cdof) { 1795 const PetscInt *oldIndices = NULL; 1796 PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0; 1797 1798 for (f = 0; f < len; ++f) { 1799 PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof)); 1800 PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof)); 1801 PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices)); 1802 PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices)); 1803 for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff; 1804 numConst += cfdof; 1805 fOff += fdof; 1806 } 1807 PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices)); 1808 } 1809 } 1810 PetscCall(PetscFree(indices)); 1811 } 1812 PetscFunctionReturn(0); 1813 } 1814 1815 /*@ 1816 PetscSectionCreateSupersection - Create a new, larger section composed of the input sections 1817 1818 Collective 1819 1820 Input Parameters: 1821 + s - the input sections 1822 - len - the number of input sections 1823 1824 Output Parameter: 1825 . supers - the supersection 1826 1827 Note: The section offsets now refer to a new, larger vector. 1828 1829 Level: advanced 1830 1831 .seealso: `PetscSectionCreateSubsection()`, `PetscSectionCreate()` 1832 @*/ 1833 PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers) 1834 { 1835 PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i; 1836 1837 PetscFunctionBegin; 1838 if (!len) PetscFunctionReturn(0); 1839 for (i = 0; i < len; ++i) { 1840 PetscInt nf, pStarti, pEndi; 1841 1842 PetscCall(PetscSectionGetNumFields(s[i], &nf)); 1843 PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi)); 1844 pStart = PetscMin(pStart, pStarti); 1845 pEnd = PetscMax(pEnd, pEndi); 1846 Nf += nf; 1847 } 1848 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers)); 1849 PetscCall(PetscSectionSetNumFields(*supers, Nf)); 1850 for (i = 0, f = 0; i < len; ++i) { 1851 PetscInt nf, fi, ci; 1852 1853 PetscCall(PetscSectionGetNumFields(s[i], &nf)); 1854 for (fi = 0; fi < nf; ++fi, ++f) { 1855 const char *name = NULL; 1856 PetscInt numComp = 0; 1857 1858 PetscCall(PetscSectionGetFieldName(s[i], fi, &name)); 1859 PetscCall(PetscSectionSetFieldName(*supers, f, name)); 1860 PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp)); 1861 PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp)); 1862 for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) { 1863 PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name)); 1864 PetscCall(PetscSectionSetComponentName(*supers, f, ci, name)); 1865 } 1866 } 1867 } 1868 PetscCall(PetscSectionSetChart(*supers, pStart, pEnd)); 1869 for (p = pStart; p < pEnd; ++p) { 1870 PetscInt dof = 0, cdof = 0; 1871 1872 for (i = 0, f = 0; i < len; ++i) { 1873 PetscInt nf, fi, pStarti, pEndi; 1874 PetscInt fdof = 0, cfdof = 0; 1875 1876 PetscCall(PetscSectionGetNumFields(s[i], &nf)); 1877 PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi)); 1878 if ((p < pStarti) || (p >= pEndi)) continue; 1879 for (fi = 0; fi < nf; ++fi, ++f) { 1880 PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof)); 1881 PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof)); 1882 PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof)); 1883 if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof)); 1884 dof += fdof; 1885 cdof += cfdof; 1886 } 1887 } 1888 PetscCall(PetscSectionSetDof(*supers, p, dof)); 1889 if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof)); 1890 maxCdof = PetscMax(cdof, maxCdof); 1891 } 1892 PetscCall(PetscSectionSetUp(*supers)); 1893 if (maxCdof) { 1894 PetscInt *indices; 1895 1896 PetscCall(PetscMalloc1(maxCdof, &indices)); 1897 for (p = pStart; p < pEnd; ++p) { 1898 PetscInt cdof; 1899 1900 PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof)); 1901 if (cdof) { 1902 PetscInt dof, numConst = 0, fOff = 0; 1903 1904 for (i = 0, f = 0; i < len; ++i) { 1905 const PetscInt *oldIndices = NULL; 1906 PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc; 1907 1908 PetscCall(PetscSectionGetNumFields(s[i], &nf)); 1909 PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi)); 1910 if ((p < pStarti) || (p >= pEndi)) continue; 1911 for (fi = 0; fi < nf; ++fi, ++f) { 1912 PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof)); 1913 PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof)); 1914 PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices)); 1915 for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc]; 1916 PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst])); 1917 for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff; 1918 numConst += cfdof; 1919 } 1920 PetscCall(PetscSectionGetDof(s[i], p, &dof)); 1921 fOff += dof; 1922 } 1923 PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices)); 1924 } 1925 } 1926 PetscCall(PetscFree(indices)); 1927 } 1928 PetscFunctionReturn(0); 1929 } 1930 1931 PetscErrorCode PetscSectionCreateSubplexSection_Internal(PetscSection s, IS subpointMap, PetscBool renumberPoints, PetscSection *subs) 1932 { 1933 const PetscInt *points = NULL, *indices = NULL; 1934 PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp; 1935 1936 PetscFunctionBegin; 1937 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1938 PetscValidHeaderSpecific(subpointMap, IS_CLASSID, 2); 1939 PetscValidPointer(subs, 4); 1940 PetscCall(PetscSectionGetNumFields(s, &numFields)); 1941 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs)); 1942 if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields)); 1943 for (f = 0; f < numFields; ++f) { 1944 const char *name = NULL; 1945 PetscInt numComp = 0; 1946 1947 PetscCall(PetscSectionGetFieldName(s, f, &name)); 1948 PetscCall(PetscSectionSetFieldName(*subs, f, name)); 1949 PetscCall(PetscSectionGetFieldComponents(s, f, &numComp)); 1950 PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp)); 1951 for (c = 0; c < s->numFieldComponents[f]; ++c) { 1952 PetscCall(PetscSectionGetComponentName(s, f, c, &name)); 1953 PetscCall(PetscSectionSetComponentName(*subs, f, c, name)); 1954 } 1955 } 1956 /* For right now, we do not try to squeeze the subchart */ 1957 if (subpointMap) { 1958 PetscCall(ISGetSize(subpointMap, &numSubpoints)); 1959 PetscCall(ISGetIndices(subpointMap, &points)); 1960 } 1961 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 1962 if (renumberPoints) { 1963 spStart = 0; 1964 spEnd = numSubpoints; 1965 } else { 1966 PetscCall(ISGetMinMax(subpointMap, &spStart, &spEnd)); 1967 ++spEnd; 1968 } 1969 PetscCall(PetscSectionSetChart(*subs, spStart, spEnd)); 1970 for (p = pStart; p < pEnd; ++p) { 1971 PetscInt dof, cdof, fdof = 0, cfdof = 0; 1972 1973 PetscCall(PetscFindInt(p, numSubpoints, points, &subp)); 1974 if (subp < 0) continue; 1975 if (!renumberPoints) subp = p; 1976 for (f = 0; f < numFields; ++f) { 1977 PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof)); 1978 PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof)); 1979 PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof)); 1980 if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof)); 1981 } 1982 PetscCall(PetscSectionGetDof(s, p, &dof)); 1983 PetscCall(PetscSectionSetDof(*subs, subp, dof)); 1984 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 1985 if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof)); 1986 } 1987 PetscCall(PetscSectionSetUp(*subs)); 1988 /* Change offsets to original offsets */ 1989 for (p = pStart; p < pEnd; ++p) { 1990 PetscInt off, foff = 0; 1991 1992 PetscCall(PetscFindInt(p, numSubpoints, points, &subp)); 1993 if (subp < 0) continue; 1994 if (!renumberPoints) subp = p; 1995 for (f = 0; f < numFields; ++f) { 1996 PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff)); 1997 PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff)); 1998 } 1999 PetscCall(PetscSectionGetOffset(s, p, &off)); 2000 PetscCall(PetscSectionSetOffset(*subs, subp, off)); 2001 } 2002 /* Copy constraint indices */ 2003 for (subp = spStart; subp < spEnd; ++subp) { 2004 PetscInt cdof; 2005 2006 PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof)); 2007 if (cdof) { 2008 for (f = 0; f < numFields; ++f) { 2009 PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices)); 2010 PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices)); 2011 } 2012 PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices)); 2013 PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices)); 2014 } 2015 } 2016 if (subpointMap) PetscCall(ISRestoreIndices(subpointMap, &points)); 2017 PetscFunctionReturn(0); 2018 } 2019 2020 /*@ 2021 PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh 2022 2023 Collective on s 2024 2025 Input Parameters: 2026 + s - the PetscSection 2027 - subpointMap - a sorted list of points in the original mesh which are in the submesh 2028 2029 Output Parameter: 2030 . subs - the subsection 2031 2032 Note: 2033 The points are renumbered from 0, and the section offsets now refer to a new, smaller vector. 2034 2035 Level: advanced 2036 2037 .seealso: `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()` 2038 @*/ 2039 PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs) 2040 { 2041 PetscFunctionBegin; 2042 PetscCall(PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_TRUE, subs)); 2043 PetscFunctionReturn(0); 2044 } 2045 2046 /*@ 2047 PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh 2048 2049 Collective on s 2050 2051 Input Parameters: 2052 + s - the PetscSection 2053 - subpointMap - a sorted list of points in the original mesh which are in the subdomain 2054 2055 Output Parameter: 2056 . subs - the subsection 2057 2058 Note: 2059 The point numbers remain the same, but the section offsets now refer to a new, smaller vector. 2060 2061 Level: advanced 2062 2063 .seealso: `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()` 2064 @*/ 2065 PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs) 2066 { 2067 PetscFunctionBegin; 2068 PetscCall(PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_FALSE, subs)); 2069 PetscFunctionReturn(0); 2070 } 2071 2072 static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer) 2073 { 2074 PetscInt p; 2075 PetscMPIInt rank; 2076 2077 PetscFunctionBegin; 2078 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank)); 2079 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 2080 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank)); 2081 for (p = 0; p < s->pEnd - s->pStart; ++p) { 2082 if ((s->bc) && (s->bc->atlasDof[p] > 0)) { 2083 PetscInt b; 2084 2085 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p])); 2086 if (s->bcIndices) { 2087 for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b])); 2088 } 2089 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n")); 2090 } else { 2091 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p])); 2092 } 2093 } 2094 PetscCall(PetscViewerFlush(viewer)); 2095 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 2096 if (s->sym) { 2097 PetscCall(PetscViewerASCIIPushTab(viewer)); 2098 PetscCall(PetscSectionSymView(s->sym, viewer)); 2099 PetscCall(PetscViewerASCIIPopTab(viewer)); 2100 } 2101 PetscFunctionReturn(0); 2102 } 2103 2104 /*@C 2105 PetscSectionViewFromOptions - View from Options 2106 2107 Collective 2108 2109 Input Parameters: 2110 + A - the PetscSection object to view 2111 . obj - Optional object 2112 - name - command line option 2113 2114 Level: intermediate 2115 .seealso: `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()` 2116 @*/ 2117 PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[]) 2118 { 2119 PetscFunctionBegin; 2120 PetscValidHeaderSpecific(A, PETSC_SECTION_CLASSID, 1); 2121 PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name)); 2122 PetscFunctionReturn(0); 2123 } 2124 2125 /*@C 2126 PetscSectionView - Views a PetscSection 2127 2128 Collective 2129 2130 Input Parameters: 2131 + s - the PetscSection object to view 2132 - v - the viewer 2133 2134 Note: 2135 PetscSectionView(), when viewer is of type PETSCVIEWERHDF5, only saves 2136 distribution independent data, such as dofs, offsets, constraint dofs, 2137 and constraint indices. Points that have negative dofs, for instance, 2138 are not saved as they represent points owned by other processes. 2139 Point numbering and rank assignment is currently not stored. 2140 The saved section can be loaded with PetscSectionLoad(). 2141 2142 Level: beginner 2143 2144 .seealso `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()` 2145 @*/ 2146 PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer) 2147 { 2148 PetscBool isascii, ishdf5; 2149 PetscInt f; 2150 2151 PetscFunctionBegin; 2152 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2153 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer)); 2154 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 2155 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 2156 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5)); 2157 if (isascii) { 2158 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer)); 2159 if (s->numFields) { 2160 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields)); 2161 for (f = 0; f < s->numFields; ++f) { 2162 PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f])); 2163 PetscCall(PetscSectionView_ASCII(s->field[f], viewer)); 2164 } 2165 } else { 2166 PetscCall(PetscSectionView_ASCII(s, viewer)); 2167 } 2168 } else if (ishdf5) { 2169 #if PetscDefined(HAVE_HDF5) 2170 PetscCall(PetscSectionView_HDF5_Internal(s, viewer)); 2171 #else 2172 SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 2173 #endif 2174 } 2175 PetscFunctionReturn(0); 2176 } 2177 2178 /*@C 2179 PetscSectionLoad - Loads a PetscSection 2180 2181 Collective 2182 2183 Input Parameters: 2184 + s - the PetscSection object to load 2185 - v - the viewer 2186 2187 Note: 2188 PetscSectionLoad(), when viewer is of type PETSCVIEWERHDF5, loads 2189 a section saved with PetscSectionView(). The number of processes 2190 used here (N) does not need to be the same as that used when saving. 2191 After calling this function, the chart of s on rank i will be set 2192 to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of 2193 saved section points. 2194 2195 Level: beginner 2196 2197 .seealso `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()` 2198 @*/ 2199 PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer) 2200 { 2201 PetscBool ishdf5; 2202 2203 PetscFunctionBegin; 2204 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2205 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 2206 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5)); 2207 if (ishdf5) { 2208 #if PetscDefined(HAVE_HDF5) 2209 PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer)); 2210 PetscFunctionReturn(0); 2211 #else 2212 SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 2213 #endif 2214 } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name); 2215 } 2216 2217 static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section) 2218 { 2219 PetscSectionClosurePermVal clVal; 2220 2221 PetscFunctionBegin; 2222 if (!section->clHash) PetscFunctionReturn(0); 2223 kh_foreach_value(section->clHash, clVal, { 2224 PetscCall(PetscFree(clVal.perm)); 2225 PetscCall(PetscFree(clVal.invPerm)); 2226 }); 2227 kh_destroy(ClPerm, section->clHash); 2228 section->clHash = NULL; 2229 PetscFunctionReturn(0); 2230 } 2231 2232 /*@ 2233 PetscSectionReset - Frees all section data. 2234 2235 Not Collective 2236 2237 Input Parameters: 2238 . s - the PetscSection 2239 2240 Level: beginner 2241 2242 .seealso: `PetscSection`, `PetscSectionCreate()` 2243 @*/ 2244 PetscErrorCode PetscSectionReset(PetscSection s) 2245 { 2246 PetscInt f, c; 2247 2248 PetscFunctionBegin; 2249 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2250 for (f = 0; f < s->numFields; ++f) { 2251 PetscCall(PetscSectionDestroy(&s->field[f])); 2252 PetscCall(PetscFree(s->fieldNames[f])); 2253 for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c])); 2254 PetscCall(PetscFree(s->compNames[f])); 2255 } 2256 PetscCall(PetscFree(s->numFieldComponents)); 2257 PetscCall(PetscFree(s->fieldNames)); 2258 PetscCall(PetscFree(s->compNames)); 2259 PetscCall(PetscFree(s->field)); 2260 PetscCall(PetscSectionDestroy(&s->bc)); 2261 PetscCall(PetscFree(s->bcIndices)); 2262 PetscCall(PetscFree2(s->atlasDof, s->atlasOff)); 2263 PetscCall(PetscSectionDestroy(&s->clSection)); 2264 PetscCall(ISDestroy(&s->clPoints)); 2265 PetscCall(ISDestroy(&s->perm)); 2266 PetscCall(PetscSectionResetClosurePermutation(s)); 2267 PetscCall(PetscSectionSymDestroy(&s->sym)); 2268 PetscCall(PetscSectionDestroy(&s->clSection)); 2269 PetscCall(ISDestroy(&s->clPoints)); 2270 PetscCall(PetscSectionInvalidateMaxDof_Internal(s)); 2271 s->pStart = -1; 2272 s->pEnd = -1; 2273 s->maxDof = 0; 2274 s->setup = PETSC_FALSE; 2275 s->numFields = 0; 2276 s->clObj = NULL; 2277 PetscFunctionReturn(0); 2278 } 2279 2280 /*@ 2281 PetscSectionDestroy - Frees a section object and frees its range if that exists. 2282 2283 Not Collective 2284 2285 Input Parameters: 2286 . s - the PetscSection 2287 2288 Level: beginner 2289 2290 .seealso: `PetscSection`, `PetscSectionCreate()` 2291 @*/ 2292 PetscErrorCode PetscSectionDestroy(PetscSection *s) 2293 { 2294 PetscFunctionBegin; 2295 if (!*s) PetscFunctionReturn(0); 2296 PetscValidHeaderSpecific(*s, PETSC_SECTION_CLASSID, 1); 2297 if (--((PetscObject)(*s))->refct > 0) { 2298 *s = NULL; 2299 PetscFunctionReturn(0); 2300 } 2301 PetscCall(PetscSectionReset(*s)); 2302 PetscCall(PetscHeaderDestroy(s)); 2303 PetscFunctionReturn(0); 2304 } 2305 2306 PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values) 2307 { 2308 const PetscInt p = point - s->pStart; 2309 2310 PetscFunctionBegin; 2311 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2); 2312 *values = &baseArray[s->atlasOff[p]]; 2313 PetscFunctionReturn(0); 2314 } 2315 2316 PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode) 2317 { 2318 PetscInt *array; 2319 const PetscInt p = point - s->pStart; 2320 const PetscInt orientation = 0; /* Needs to be included for use in closure operations */ 2321 PetscInt cDim = 0; 2322 2323 PetscFunctionBegin; 2324 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2); 2325 PetscCall(PetscSectionGetConstraintDof(s, p, &cDim)); 2326 array = &baseArray[s->atlasOff[p]]; 2327 if (!cDim) { 2328 if (orientation >= 0) { 2329 const PetscInt dim = s->atlasDof[p]; 2330 PetscInt i; 2331 2332 if (mode == INSERT_VALUES) { 2333 for (i = 0; i < dim; ++i) array[i] = values[i]; 2334 } else { 2335 for (i = 0; i < dim; ++i) array[i] += values[i]; 2336 } 2337 } else { 2338 PetscInt offset = 0; 2339 PetscInt j = -1, field, i; 2340 2341 for (field = 0; field < s->numFields; ++field) { 2342 const PetscInt dim = s->field[field]->atlasDof[p]; 2343 2344 for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset]; 2345 offset += dim; 2346 } 2347 } 2348 } else { 2349 if (orientation >= 0) { 2350 const PetscInt dim = s->atlasDof[p]; 2351 PetscInt cInd = 0, i; 2352 const PetscInt *cDof; 2353 2354 PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof)); 2355 if (mode == INSERT_VALUES) { 2356 for (i = 0; i < dim; ++i) { 2357 if ((cInd < cDim) && (i == cDof[cInd])) { 2358 ++cInd; 2359 continue; 2360 } 2361 array[i] = values[i]; 2362 } 2363 } else { 2364 for (i = 0; i < dim; ++i) { 2365 if ((cInd < cDim) && (i == cDof[cInd])) { 2366 ++cInd; 2367 continue; 2368 } 2369 array[i] += values[i]; 2370 } 2371 } 2372 } else { 2373 const PetscInt *cDof; 2374 PetscInt offset = 0; 2375 PetscInt cOffset = 0; 2376 PetscInt j = 0, field; 2377 2378 PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof)); 2379 for (field = 0; field < s->numFields; ++field) { 2380 const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */ 2381 const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */ 2382 const PetscInt sDim = dim - tDim; 2383 PetscInt cInd = 0, i, k; 2384 2385 for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) { 2386 if ((cInd < sDim) && (j == cDof[cInd + cOffset])) { 2387 ++cInd; 2388 continue; 2389 } 2390 array[j] = values[k]; 2391 } 2392 offset += dim; 2393 cOffset += dim - tDim; 2394 } 2395 } 2396 } 2397 PetscFunctionReturn(0); 2398 } 2399 2400 /*@C 2401 PetscSectionHasConstraints - Determine whether a section has constrained dofs 2402 2403 Not Collective 2404 2405 Input Parameter: 2406 . s - The PetscSection 2407 2408 Output Parameter: 2409 . hasConstraints - flag indicating that the section has constrained dofs 2410 2411 Level: intermediate 2412 2413 .seealso: `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection` 2414 @*/ 2415 PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints) 2416 { 2417 PetscFunctionBegin; 2418 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2419 PetscValidBoolPointer(hasConstraints, 2); 2420 *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE; 2421 PetscFunctionReturn(0); 2422 } 2423 2424 /*@C 2425 PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained 2426 2427 Not Collective 2428 2429 Input Parameters: 2430 + s - The PetscSection 2431 - point - The point 2432 2433 Output Parameter: 2434 . indices - The constrained dofs 2435 2436 Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90() 2437 2438 Level: intermediate 2439 2440 .seealso: `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection` 2441 @*/ 2442 PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices) 2443 { 2444 PetscFunctionBegin; 2445 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2446 if (s->bc) { 2447 PetscCall(VecIntGetValuesSection(s->bcIndices, s->bc, point, indices)); 2448 } else *indices = NULL; 2449 PetscFunctionReturn(0); 2450 } 2451 2452 /*@C 2453 PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained 2454 2455 Not Collective 2456 2457 Input Parameters: 2458 + s - The PetscSection 2459 . point - The point 2460 - indices - The constrained dofs 2461 2462 Note: The Fortran is PetscSectionSetConstraintIndicesF90() 2463 2464 Level: intermediate 2465 2466 .seealso: `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection` 2467 @*/ 2468 PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[]) 2469 { 2470 PetscFunctionBegin; 2471 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2472 if (s->bc) { 2473 const PetscInt dof = s->atlasDof[point]; 2474 const PetscInt cdof = s->bc->atlasDof[point]; 2475 PetscInt d; 2476 2477 for (d = 0; d < cdof; ++d) PetscCheck(indices[d] < dof, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " dof %" PetscInt_FMT ", invalid constraint index[%" PetscInt_FMT "]: %" PetscInt_FMT, point, dof, d, indices[d]); 2478 PetscCall(VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES)); 2479 } 2480 PetscFunctionReturn(0); 2481 } 2482 2483 /*@C 2484 PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained 2485 2486 Not Collective 2487 2488 Input Parameters: 2489 + s - The PetscSection 2490 . field - The field number 2491 - point - The point 2492 2493 Output Parameter: 2494 . indices - The constrained dofs sorted in ascending order 2495 2496 Notes: 2497 The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by PetscSectionGetConstraintDof(). 2498 2499 Fortran Note: 2500 In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90() 2501 2502 Level: intermediate 2503 2504 .seealso: `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection` 2505 @*/ 2506 PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices) 2507 { 2508 PetscFunctionBegin; 2509 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2510 PetscValidPointer(indices, 4); 2511 PetscSectionCheckValidField(field, s->numFields); 2512 PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices)); 2513 PetscFunctionReturn(0); 2514 } 2515 2516 /*@C 2517 PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained 2518 2519 Not Collective 2520 2521 Input Parameters: 2522 + s - The PetscSection 2523 . point - The point 2524 . field - The field number 2525 - indices - The constrained dofs 2526 2527 Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90() 2528 2529 Level: intermediate 2530 2531 .seealso: `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection` 2532 @*/ 2533 PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[]) 2534 { 2535 PetscFunctionBegin; 2536 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2537 if (PetscDefined(USE_DEBUG)) { 2538 PetscInt nfdof; 2539 2540 PetscCall(PetscSectionGetFieldConstraintDof(s, point, field, &nfdof)); 2541 if (nfdof) PetscValidIntPointer(indices, 4); 2542 } 2543 PetscSectionCheckValidField(field, s->numFields); 2544 PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices)); 2545 PetscFunctionReturn(0); 2546 } 2547 2548 /*@ 2549 PetscSectionPermute - Reorder the section according to the input point permutation 2550 2551 Collective 2552 2553 Input Parameters: 2554 + section - The PetscSection object 2555 - perm - The point permutation, old point p becomes new point perm[p] 2556 2557 Output Parameter: 2558 . sectionNew - The permuted PetscSection 2559 2560 Level: intermediate 2561 2562 .seealso: `MatPermute()` 2563 @*/ 2564 PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew) 2565 { 2566 PetscSection s = section, sNew; 2567 const PetscInt *perm; 2568 PetscInt numFields, f, c, numPoints, pStart, pEnd, p; 2569 2570 PetscFunctionBegin; 2571 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 2572 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2); 2573 PetscValidPointer(sectionNew, 3); 2574 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew)); 2575 PetscCall(PetscSectionGetNumFields(s, &numFields)); 2576 if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields)); 2577 for (f = 0; f < numFields; ++f) { 2578 const char *name; 2579 PetscInt numComp; 2580 2581 PetscCall(PetscSectionGetFieldName(s, f, &name)); 2582 PetscCall(PetscSectionSetFieldName(sNew, f, name)); 2583 PetscCall(PetscSectionGetFieldComponents(s, f, &numComp)); 2584 PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp)); 2585 for (c = 0; c < s->numFieldComponents[f]; ++c) { 2586 PetscCall(PetscSectionGetComponentName(s, f, c, &name)); 2587 PetscCall(PetscSectionSetComponentName(sNew, f, c, name)); 2588 } 2589 } 2590 PetscCall(ISGetLocalSize(permutation, &numPoints)); 2591 PetscCall(ISGetIndices(permutation, &perm)); 2592 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 2593 PetscCall(PetscSectionSetChart(sNew, pStart, pEnd)); 2594 PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd); 2595 for (p = pStart; p < pEnd; ++p) { 2596 PetscInt dof, cdof; 2597 2598 PetscCall(PetscSectionGetDof(s, p, &dof)); 2599 PetscCall(PetscSectionSetDof(sNew, perm[p], dof)); 2600 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 2601 if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof)); 2602 for (f = 0; f < numFields; ++f) { 2603 PetscCall(PetscSectionGetFieldDof(s, p, f, &dof)); 2604 PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof)); 2605 PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof)); 2606 if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof)); 2607 } 2608 } 2609 PetscCall(PetscSectionSetUp(sNew)); 2610 for (p = pStart; p < pEnd; ++p) { 2611 const PetscInt *cind; 2612 PetscInt cdof; 2613 2614 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 2615 if (cdof) { 2616 PetscCall(PetscSectionGetConstraintIndices(s, p, &cind)); 2617 PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind)); 2618 } 2619 for (f = 0; f < numFields; ++f) { 2620 PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof)); 2621 if (cdof) { 2622 PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind)); 2623 PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind)); 2624 } 2625 } 2626 } 2627 PetscCall(ISRestoreIndices(permutation, &perm)); 2628 *sectionNew = sNew; 2629 PetscFunctionReturn(0); 2630 } 2631 2632 /*@ 2633 PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section 2634 2635 Collective 2636 2637 Input Parameters: 2638 + section - The PetscSection 2639 . obj - A PetscObject which serves as the key for this index 2640 . clSection - Section giving the size of the closure of each point 2641 - clPoints - IS giving the points in each closure 2642 2643 Note: We compress out closure points with no dofs in this section 2644 2645 Level: advanced 2646 2647 .seealso: `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()` 2648 @*/ 2649 PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints) 2650 { 2651 PetscFunctionBegin; 2652 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 2653 PetscValidHeaderSpecific(clSection, PETSC_SECTION_CLASSID, 3); 2654 PetscValidHeaderSpecific(clPoints, IS_CLASSID, 4); 2655 if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section)); 2656 section->clObj = obj; 2657 PetscCall(PetscObjectReference((PetscObject)clSection)); 2658 PetscCall(PetscObjectReference((PetscObject)clPoints)); 2659 PetscCall(PetscSectionDestroy(§ion->clSection)); 2660 PetscCall(ISDestroy(§ion->clPoints)); 2661 section->clSection = clSection; 2662 section->clPoints = clPoints; 2663 PetscFunctionReturn(0); 2664 } 2665 2666 /*@ 2667 PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section 2668 2669 Collective 2670 2671 Input Parameters: 2672 + section - The PetscSection 2673 - obj - A PetscObject which serves as the key for this index 2674 2675 Output Parameters: 2676 + clSection - Section giving the size of the closure of each point 2677 - clPoints - IS giving the points in each closure 2678 2679 Note: We compress out closure points with no dofs in this section 2680 2681 Level: advanced 2682 2683 .seealso: `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()` 2684 @*/ 2685 PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints) 2686 { 2687 PetscFunctionBegin; 2688 if (section->clObj == obj) { 2689 if (clSection) *clSection = section->clSection; 2690 if (clPoints) *clPoints = section->clPoints; 2691 } else { 2692 if (clSection) *clSection = NULL; 2693 if (clPoints) *clPoints = NULL; 2694 } 2695 PetscFunctionReturn(0); 2696 } 2697 2698 PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm) 2699 { 2700 PetscInt i; 2701 khiter_t iter; 2702 int new_entry; 2703 PetscSectionClosurePermKey key = {depth, clSize}; 2704 PetscSectionClosurePermVal *val; 2705 2706 PetscFunctionBegin; 2707 if (section->clObj != obj) { 2708 PetscCall(PetscSectionDestroy(§ion->clSection)); 2709 PetscCall(ISDestroy(§ion->clPoints)); 2710 } 2711 section->clObj = obj; 2712 if (!section->clHash) PetscCall(PetscClPermCreate(§ion->clHash)); 2713 iter = kh_put(ClPerm, section->clHash, key, &new_entry); 2714 val = &kh_val(section->clHash, iter); 2715 if (!new_entry) { 2716 PetscCall(PetscFree(val->perm)); 2717 PetscCall(PetscFree(val->invPerm)); 2718 } 2719 if (mode == PETSC_COPY_VALUES) { 2720 PetscCall(PetscMalloc1(clSize, &val->perm)); 2721 PetscCall(PetscArraycpy(val->perm, clPerm, clSize)); 2722 } else if (mode == PETSC_OWN_POINTER) { 2723 val->perm = clPerm; 2724 } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays"); 2725 PetscCall(PetscMalloc1(clSize, &val->invPerm)); 2726 for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i; 2727 PetscFunctionReturn(0); 2728 } 2729 2730 /*@ 2731 PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex. 2732 2733 Not Collective 2734 2735 Input Parameters: 2736 + section - The PetscSection 2737 . obj - A PetscObject which serves as the key for this index (usually a DM) 2738 . depth - Depth of points on which to apply the given permutation 2739 - perm - Permutation of the cell dof closure 2740 2741 Note: 2742 The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a 2743 mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for 2744 each topology and degree. 2745 2746 This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for 2747 exotic/enriched spaces on mixed topology meshes. 2748 2749 Level: intermediate 2750 2751 .seealso: `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode` 2752 @*/ 2753 PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm) 2754 { 2755 const PetscInt *clPerm = NULL; 2756 PetscInt clSize = 0; 2757 2758 PetscFunctionBegin; 2759 if (perm) { 2760 PetscCall(ISGetLocalSize(perm, &clSize)); 2761 PetscCall(ISGetIndices(perm, &clPerm)); 2762 } 2763 PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm)); 2764 if (perm) PetscCall(ISRestoreIndices(perm, &clPerm)); 2765 PetscFunctionReturn(0); 2766 } 2767 2768 PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[]) 2769 { 2770 PetscFunctionBegin; 2771 if (section->clObj == obj) { 2772 PetscSectionClosurePermKey k = {depth, size}; 2773 PetscSectionClosurePermVal v; 2774 PetscCall(PetscClPermGet(section->clHash, k, &v)); 2775 if (perm) *perm = v.perm; 2776 } else { 2777 if (perm) *perm = NULL; 2778 } 2779 PetscFunctionReturn(0); 2780 } 2781 2782 /*@ 2783 PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex. 2784 2785 Not Collective 2786 2787 Input Parameters: 2788 + section - The PetscSection 2789 . obj - A PetscObject which serves as the key for this index (usually a DM) 2790 . depth - Depth stratum on which to obtain closure permutation 2791 - clSize - Closure size to be permuted (e.g., may vary with element topology and degree) 2792 2793 Output Parameter: 2794 . perm - The dof closure permutation 2795 2796 Note: 2797 The user must destroy the IS that is returned. 2798 2799 Level: intermediate 2800 2801 .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()` 2802 @*/ 2803 PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm) 2804 { 2805 const PetscInt *clPerm; 2806 2807 PetscFunctionBegin; 2808 PetscCall(PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm)); 2809 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm)); 2810 PetscFunctionReturn(0); 2811 } 2812 2813 PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[]) 2814 { 2815 PetscFunctionBegin; 2816 if (section->clObj == obj && section->clHash) { 2817 PetscSectionClosurePermKey k = {depth, size}; 2818 PetscSectionClosurePermVal v; 2819 PetscCall(PetscClPermGet(section->clHash, k, &v)); 2820 if (perm) *perm = v.invPerm; 2821 } else { 2822 if (perm) *perm = NULL; 2823 } 2824 PetscFunctionReturn(0); 2825 } 2826 2827 /*@ 2828 PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex. 2829 2830 Not Collective 2831 2832 Input Parameters: 2833 + section - The PetscSection 2834 . obj - A PetscObject which serves as the key for this index (usually a DM) 2835 . depth - Depth stratum on which to obtain closure permutation 2836 - clSize - Closure size to be permuted (e.g., may vary with element topology and degree) 2837 2838 Output Parameters: 2839 . perm - The dof closure permutation 2840 2841 Note: 2842 The user must destroy the IS that is returned. 2843 2844 Level: intermediate 2845 2846 .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()` 2847 @*/ 2848 PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm) 2849 { 2850 const PetscInt *clPerm; 2851 2852 PetscFunctionBegin; 2853 PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm)); 2854 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm)); 2855 PetscFunctionReturn(0); 2856 } 2857 2858 /*@ 2859 PetscSectionGetField - Get the subsection associated with a single field 2860 2861 Input Parameters: 2862 + s - The PetscSection 2863 - field - The field number 2864 2865 Output Parameter: 2866 . subs - The subsection for the given field 2867 2868 Level: intermediate 2869 2870 .seealso: `PetscSectionSetNumFields()` 2871 @*/ 2872 PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs) 2873 { 2874 PetscFunctionBegin; 2875 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2876 PetscValidPointer(subs, 3); 2877 PetscSectionCheckValidField(field, s->numFields); 2878 *subs = s->field[field]; 2879 PetscFunctionReturn(0); 2880 } 2881 2882 PetscClassId PETSC_SECTION_SYM_CLASSID; 2883 PetscFunctionList PetscSectionSymList = NULL; 2884 2885 /*@ 2886 PetscSectionSymCreate - Creates an empty PetscSectionSym object. 2887 2888 Collective 2889 2890 Input Parameter: 2891 . comm - the MPI communicator 2892 2893 Output Parameter: 2894 . sym - pointer to the new set of symmetries 2895 2896 Level: developer 2897 2898 .seealso: `PetscSectionSym`, `PetscSectionSymDestroy()` 2899 @*/ 2900 PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym) 2901 { 2902 PetscFunctionBegin; 2903 PetscValidPointer(sym, 2); 2904 PetscCall(ISInitializePackage()); 2905 PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView)); 2906 PetscFunctionReturn(0); 2907 } 2908 2909 /*@C 2910 PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation. 2911 2912 Collective 2913 2914 Input Parameters: 2915 + sym - The section symmetry object 2916 - method - The name of the section symmetry type 2917 2918 Level: developer 2919 2920 .seealso: `PetscSectionSymGetType()`, `PetscSectionSymCreate()` 2921 @*/ 2922 PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method) 2923 { 2924 PetscErrorCode (*r)(PetscSectionSym); 2925 PetscBool match; 2926 2927 PetscFunctionBegin; 2928 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2929 PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match)); 2930 if (match) PetscFunctionReturn(0); 2931 2932 PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r)); 2933 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method); 2934 PetscTryTypeMethod(sym, destroy); 2935 sym->ops->destroy = NULL; 2936 2937 PetscCall((*r)(sym)); 2938 PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method)); 2939 PetscFunctionReturn(0); 2940 } 2941 2942 /*@C 2943 PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym. 2944 2945 Not Collective 2946 2947 Input Parameter: 2948 . sym - The section symmetry 2949 2950 Output Parameter: 2951 . type - The index set type name 2952 2953 Level: developer 2954 2955 .seealso: `PetscSectionSymSetType()`, `PetscSectionSymCreate()` 2956 @*/ 2957 PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type) 2958 { 2959 PetscFunctionBegin; 2960 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 2961 PetscValidPointer(type, 2); 2962 *type = ((PetscObject)sym)->type_name; 2963 PetscFunctionReturn(0); 2964 } 2965 2966 /*@C 2967 PetscSectionSymRegister - Adds a new section symmetry implementation 2968 2969 Not Collective 2970 2971 Input Parameters: 2972 + name - The name of a new user-defined creation routine 2973 - create_func - The creation routine itself 2974 2975 Notes: 2976 PetscSectionSymRegister() may be called multiple times to add several user-defined vectors 2977 2978 Level: developer 2979 2980 .seealso: `PetscSectionSymCreate()`, `PetscSectionSymSetType()` 2981 @*/ 2982 PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym)) 2983 { 2984 PetscFunctionBegin; 2985 PetscCall(ISInitializePackage()); 2986 PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function)); 2987 PetscFunctionReturn(0); 2988 } 2989 2990 /*@ 2991 PetscSectionSymDestroy - Destroys a section symmetry. 2992 2993 Collective 2994 2995 Input Parameters: 2996 . sym - the section symmetry 2997 2998 Level: developer 2999 3000 .seealso: `PetscSectionSymCreate()`, `PetscSectionSymDestroy()` 3001 @*/ 3002 PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym) 3003 { 3004 SymWorkLink link, next; 3005 3006 PetscFunctionBegin; 3007 if (!*sym) PetscFunctionReturn(0); 3008 PetscValidHeaderSpecific((*sym), PETSC_SECTION_SYM_CLASSID, 1); 3009 if (--((PetscObject)(*sym))->refct > 0) { 3010 *sym = NULL; 3011 PetscFunctionReturn(0); 3012 } 3013 if ((*sym)->ops->destroy) PetscCall((*(*sym)->ops->destroy)(*sym)); 3014 PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out"); 3015 for (link = (*sym)->workin; link; link = next) { 3016 PetscInt **perms = (PetscInt **)link->perms; 3017 PetscScalar **rots = (PetscScalar **)link->rots; 3018 PetscCall(PetscFree2(perms, rots)); 3019 next = link->next; 3020 PetscCall(PetscFree(link)); 3021 } 3022 (*sym)->workin = NULL; 3023 PetscCall(PetscHeaderDestroy(sym)); 3024 PetscFunctionReturn(0); 3025 } 3026 3027 /*@C 3028 PetscSectionSymView - Displays a section symmetry 3029 3030 Collective 3031 3032 Input Parameters: 3033 + sym - the index set 3034 - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF. 3035 3036 Level: developer 3037 3038 .seealso: `PetscViewerASCIIOpen()` 3039 @*/ 3040 PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer) 3041 { 3042 PetscFunctionBegin; 3043 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 3044 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer)); 3045 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 3046 PetscCheckSameComm(sym, 1, viewer, 2); 3047 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer)); 3048 PetscTryTypeMethod(sym, view, viewer); 3049 PetscFunctionReturn(0); 3050 } 3051 3052 /*@ 3053 PetscSectionSetSym - Set the symmetries for the data referred to by the section 3054 3055 Collective 3056 3057 Input Parameters: 3058 + section - the section describing data layout 3059 - sym - the symmetry describing the affect of orientation on the access of the data 3060 3061 Level: developer 3062 3063 .seealso: `PetscSectionGetSym()`, `PetscSectionSymCreate()` 3064 @*/ 3065 PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym) 3066 { 3067 PetscFunctionBegin; 3068 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 3069 PetscCall(PetscSectionSymDestroy(&(section->sym))); 3070 if (sym) { 3071 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 2); 3072 PetscCheckSameComm(section, 1, sym, 2); 3073 PetscCall(PetscObjectReference((PetscObject)sym)); 3074 } 3075 section->sym = sym; 3076 PetscFunctionReturn(0); 3077 } 3078 3079 /*@ 3080 PetscSectionGetSym - Get the symmetries for the data referred to by the section 3081 3082 Not Collective 3083 3084 Input Parameters: 3085 . section - the section describing data layout 3086 3087 Output Parameters: 3088 . sym - the symmetry describing the affect of orientation on the access of the data 3089 3090 Level: developer 3091 3092 .seealso: `PetscSectionSetSym()`, `PetscSectionSymCreate()` 3093 @*/ 3094 PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym) 3095 { 3096 PetscFunctionBegin; 3097 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 3098 *sym = section->sym; 3099 PetscFunctionReturn(0); 3100 } 3101 3102 /*@ 3103 PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section 3104 3105 Collective 3106 3107 Input Parameters: 3108 + section - the section describing data layout 3109 . field - the field number 3110 - sym - the symmetry describing the affect of orientation on the access of the data 3111 3112 Level: developer 3113 3114 .seealso: `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()` 3115 @*/ 3116 PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym) 3117 { 3118 PetscFunctionBegin; 3119 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 3120 PetscSectionCheckValidField(field, section->numFields); 3121 PetscCall(PetscSectionSetSym(section->field[field], sym)); 3122 PetscFunctionReturn(0); 3123 } 3124 3125 /*@ 3126 PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section 3127 3128 Collective 3129 3130 Input Parameters: 3131 + section - the section describing data layout 3132 - field - the field number 3133 3134 Output Parameters: 3135 . sym - the symmetry describing the affect of orientation on the access of the data 3136 3137 Level: developer 3138 3139 .seealso: `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()` 3140 @*/ 3141 PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym) 3142 { 3143 PetscFunctionBegin; 3144 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 3145 PetscSectionCheckValidField(field, section->numFields); 3146 *sym = section->field[field]->sym; 3147 PetscFunctionReturn(0); 3148 } 3149 3150 /*@C 3151 PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations. 3152 3153 Not Collective 3154 3155 Input Parameters: 3156 + section - the section 3157 . numPoints - the number of points 3158 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3159 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3160 context, see DMPlexGetConeOrientation()). 3161 3162 Output Parameters: 3163 + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity). 3164 - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all 3165 identity). 3166 3167 Example of usage, gathering dofs into a local array (lArray) from a section array (sArray): 3168 .vb 3169 const PetscInt **perms; 3170 const PetscScalar **rots; 3171 PetscInt lOffset; 3172 3173 PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots); 3174 for (i = 0, lOffset = 0; i < numPoints; i++) { 3175 PetscInt point = points[2*i], dof, sOffset; 3176 const PetscInt *perm = perms ? perms[i] : NULL; 3177 const PetscScalar *rot = rots ? rots[i] : NULL; 3178 3179 PetscSectionGetDof(section,point,&dof); 3180 PetscSectionGetOffset(section,point,&sOffset); 3181 3182 if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}} 3183 else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}} 3184 if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }} 3185 lOffset += dof; 3186 } 3187 PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots); 3188 .ve 3189 3190 Example of usage, adding dofs into a section array (sArray) from a local array (lArray): 3191 .vb 3192 const PetscInt **perms; 3193 const PetscScalar **rots; 3194 PetscInt lOffset; 3195 3196 PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots); 3197 for (i = 0, lOffset = 0; i < numPoints; i++) { 3198 PetscInt point = points[2*i], dof, sOffset; 3199 const PetscInt *perm = perms ? perms[i] : NULL; 3200 const PetscScalar *rot = rots ? rots[i] : NULL; 3201 3202 PetscSectionGetDof(section,point,&dof); 3203 PetscSectionGetOffset(section,point,&sOff); 3204 3205 if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}} 3206 else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}} 3207 offset += dof; 3208 } 3209 PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots); 3210 .ve 3211 3212 Level: developer 3213 3214 .seealso: `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()` 3215 @*/ 3216 PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3217 { 3218 PetscSectionSym sym; 3219 3220 PetscFunctionBegin; 3221 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 3222 if (numPoints) PetscValidIntPointer(points, 3); 3223 if (perms) *perms = NULL; 3224 if (rots) *rots = NULL; 3225 sym = section->sym; 3226 if (sym && (perms || rots)) { 3227 SymWorkLink link; 3228 3229 if (sym->workin) { 3230 link = sym->workin; 3231 sym->workin = sym->workin->next; 3232 } else { 3233 PetscCall(PetscNew(&link)); 3234 } 3235 if (numPoints > link->numPoints) { 3236 PetscInt **perms = (PetscInt **)link->perms; 3237 PetscScalar **rots = (PetscScalar **)link->rots; 3238 PetscCall(PetscFree2(perms, rots)); 3239 PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots)); 3240 link->numPoints = numPoints; 3241 } 3242 link->next = sym->workout; 3243 sym->workout = link; 3244 PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints)); 3245 PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints)); 3246 PetscCall((*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots)); 3247 if (perms) *perms = link->perms; 3248 if (rots) *rots = link->rots; 3249 } 3250 PetscFunctionReturn(0); 3251 } 3252 3253 /*@C 3254 PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms() 3255 3256 Not Collective 3257 3258 Input Parameters: 3259 + section - the section 3260 . numPoints - the number of points 3261 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3262 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3263 context, see DMPlexGetConeOrientation()). 3264 3265 Output Parameters: 3266 + perms - The permutations for the given orientations: set to NULL at conclusion 3267 - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion 3268 3269 Level: developer 3270 3271 .seealso: `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()` 3272 @*/ 3273 PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3274 { 3275 PetscSectionSym sym; 3276 3277 PetscFunctionBegin; 3278 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 3279 sym = section->sym; 3280 if (sym && (perms || rots)) { 3281 SymWorkLink *p, link; 3282 3283 for (p = &sym->workout; (link = *p); p = &link->next) { 3284 if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) { 3285 *p = link->next; 3286 link->next = sym->workin; 3287 sym->workin = link; 3288 if (perms) *perms = NULL; 3289 if (rots) *rots = NULL; 3290 PetscFunctionReturn(0); 3291 } 3292 } 3293 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out"); 3294 } 3295 PetscFunctionReturn(0); 3296 } 3297 3298 /*@C 3299 PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations. 3300 3301 Not Collective 3302 3303 Input Parameters: 3304 + section - the section 3305 . field - the field of the section 3306 . numPoints - the number of points 3307 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3308 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3309 context, see DMPlexGetConeOrientation()). 3310 3311 Output Parameters: 3312 + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity). 3313 - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all 3314 identity). 3315 3316 Level: developer 3317 3318 .seealso: `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()` 3319 @*/ 3320 PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3321 { 3322 PetscFunctionBegin; 3323 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 3324 PetscCheck(field <= section->numFields, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "field %" PetscInt_FMT " greater than number of fields (%" PetscInt_FMT ") in section", field, section->numFields); 3325 PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots)); 3326 PetscFunctionReturn(0); 3327 } 3328 3329 /*@C 3330 PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms() 3331 3332 Not Collective 3333 3334 Input Parameters: 3335 + section - the section 3336 . field - the field number 3337 . numPoints - the number of points 3338 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3339 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3340 context, see DMPlexGetConeOrientation()). 3341 3342 Output Parameters: 3343 + perms - The permutations for the given orientations: set to NULL at conclusion 3344 - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion 3345 3346 Level: developer 3347 3348 .seealso: `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()` 3349 @*/ 3350 PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3351 { 3352 PetscFunctionBegin; 3353 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 3354 PetscCheck(field <= section->numFields, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "field %" PetscInt_FMT " greater than number of fields (%" PetscInt_FMT ") in section", field, section->numFields); 3355 PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots)); 3356 PetscFunctionReturn(0); 3357 } 3358 3359 /*@ 3360 PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible 3361 3362 Not Collective 3363 3364 Input Parameter: 3365 . sym - the PetscSectionSym 3366 3367 Output Parameter: 3368 . nsym - the equivalent symmetries 3369 3370 Level: developer 3371 3372 .seealso: `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()` 3373 @*/ 3374 PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym) 3375 { 3376 PetscFunctionBegin; 3377 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 3378 PetscValidHeaderSpecific(nsym, PETSC_SECTION_SYM_CLASSID, 2); 3379 PetscTryTypeMethod(sym, copy, nsym); 3380 PetscFunctionReturn(0); 3381 } 3382 3383 /*@ 3384 PetscSectionSymDistribute - Distribute the symmetries in accordance with the input SF 3385 3386 Collective 3387 3388 Input Parameters: 3389 + sym - the PetscSectionSym 3390 - migrationSF - the distribution map from roots to leaves 3391 3392 Output Parameters: 3393 . dsym - the redistributed symmetries 3394 3395 Level: developer 3396 3397 .seealso: `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()` 3398 @*/ 3399 PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym) 3400 { 3401 PetscFunctionBegin; 3402 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID, 1); 3403 PetscValidHeaderSpecific(migrationSF, PETSCSF_CLASSID, 2); 3404 PetscValidPointer(dsym, 3); 3405 PetscTryTypeMethod(sym, distribute, migrationSF, dsym); 3406 PetscFunctionReturn(0); 3407 } 3408 3409 /*@ 3410 PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset 3411 3412 Not Collective 3413 3414 Input Parameter: 3415 . s - the global PetscSection 3416 3417 Output Parameters: 3418 . flg - the flag 3419 3420 Level: developer 3421 3422 .seealso: `PetscSectionSetChart()`, `PetscSectionCreate()` 3423 @*/ 3424 PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg) 3425 { 3426 PetscFunctionBegin; 3427 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 3428 *flg = s->useFieldOff; 3429 PetscFunctionReturn(0); 3430 } 3431 3432 /*@ 3433 PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset 3434 3435 Not Collective 3436 3437 Input Parameters: 3438 + s - the global PetscSection 3439 - flg - the flag 3440 3441 Level: developer 3442 3443 .seealso: `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()` 3444 @*/ 3445 PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg) 3446 { 3447 PetscFunctionBegin; 3448 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 3449 s->useFieldOff = flg; 3450 PetscFunctionReturn(0); 3451 } 3452 3453 #define PetscSectionExpandPoints_Loop(TYPE) \ 3454 { \ 3455 PetscInt i, n, o0, o1, size; \ 3456 TYPE *a0 = (TYPE *)origArray, *a1; \ 3457 PetscCall(PetscSectionGetStorageSize(s, &size)); \ 3458 PetscCall(PetscMalloc1(size, &a1)); \ 3459 for (i = 0; i < npoints; i++) { \ 3460 PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \ 3461 PetscCall(PetscSectionGetOffset(s, i, &o1)); \ 3462 PetscCall(PetscSectionGetDof(s, i, &n)); \ 3463 PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n *unitsize)); \ 3464 } \ 3465 *newArray = (void *)a1; \ 3466 } 3467 3468 /*@ 3469 PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points. 3470 3471 Not Collective 3472 3473 Input Parameters: 3474 + origSection - the PetscSection describing the layout of the array 3475 . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL) 3476 . origArray - the array; its size must be equal to the storage size of origSection 3477 - points - IS with points to extract; its indices must lie in the chart of origSection 3478 3479 Output Parameters: 3480 + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs) 3481 - newArray - the array of the extracted DOFs; its size is the storage size of newSection 3482 3483 Level: developer 3484 3485 .seealso: `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()` 3486 @*/ 3487 PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[]) 3488 { 3489 PetscSection s; 3490 const PetscInt *points_; 3491 PetscInt i, n, npoints, pStart, pEnd; 3492 PetscMPIInt unitsize; 3493 3494 PetscFunctionBegin; 3495 PetscValidHeaderSpecific(origSection, PETSC_SECTION_CLASSID, 1); 3496 PetscValidPointer(origArray, 3); 3497 PetscValidHeaderSpecific(points, IS_CLASSID, 4); 3498 if (newSection) PetscValidPointer(newSection, 5); 3499 if (newArray) PetscValidPointer(newArray, 6); 3500 PetscCallMPI(MPI_Type_size(dataType, &unitsize)); 3501 PetscCall(ISGetLocalSize(points, &npoints)); 3502 PetscCall(ISGetIndices(points, &points_)); 3503 PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd)); 3504 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s)); 3505 PetscCall(PetscSectionSetChart(s, 0, npoints)); 3506 for (i = 0; i < npoints; i++) { 3507 PetscCheck(points_[i] >= pStart && points_[i] < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " (index %" PetscInt_FMT ") in input IS out of input section's chart", points_[i], i); 3508 PetscCall(PetscSectionGetDof(origSection, points_[i], &n)); 3509 PetscCall(PetscSectionSetDof(s, i, n)); 3510 } 3511 PetscCall(PetscSectionSetUp(s)); 3512 if (newArray) { 3513 if (dataType == MPIU_INT) { 3514 PetscSectionExpandPoints_Loop(PetscInt); 3515 } else if (dataType == MPIU_SCALAR) { 3516 PetscSectionExpandPoints_Loop(PetscScalar); 3517 } else if (dataType == MPIU_REAL) { 3518 PetscSectionExpandPoints_Loop(PetscReal); 3519 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype"); 3520 } 3521 if (newSection) { 3522 *newSection = s; 3523 } else { 3524 PetscCall(PetscSectionDestroy(&s)); 3525 } 3526 PetscCall(ISRestoreIndices(points, &points_)); 3527 PetscFunctionReturn(0); 3528 } 3529