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