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