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