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