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(s->bc->atlasOff[last] + s->bc->atlasDof[last], &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 PetscSectionView - Views a PetscSection 1952 1953 Collective on PetscSection 1954 1955 Input Parameters: 1956 + s - the PetscSection object to view 1957 - v - the viewer 1958 1959 Level: beginner 1960 1961 .seealso PetscSectionCreate(), PetscSectionDestroy() 1962 @*/ 1963 PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer) 1964 { 1965 PetscBool isascii; 1966 PetscInt f; 1967 PetscErrorCode ierr; 1968 1969 PetscFunctionBegin; 1970 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 1971 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);CHKERRQ(ierr);} 1972 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1973 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 1974 if (isascii) { 1975 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);CHKERRQ(ierr); 1976 if (s->numFields) { 1977 ierr = PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);CHKERRQ(ierr); 1978 for (f = 0; f < s->numFields; ++f) { 1979 ierr = PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);CHKERRQ(ierr); 1980 ierr = PetscSectionView_ASCII(s->field[f], viewer);CHKERRQ(ierr); 1981 } 1982 } else { 1983 ierr = PetscSectionView_ASCII(s, viewer);CHKERRQ(ierr); 1984 } 1985 } 1986 PetscFunctionReturn(0); 1987 } 1988 1989 /*@ 1990 PetscSectionReset - Frees all section data. 1991 1992 Not collective 1993 1994 Input Parameters: 1995 . s - the PetscSection 1996 1997 Level: beginner 1998 1999 .seealso: PetscSection, PetscSectionCreate() 2000 @*/ 2001 PetscErrorCode PetscSectionReset(PetscSection s) 2002 { 2003 PetscInt f; 2004 PetscErrorCode ierr; 2005 2006 PetscFunctionBegin; 2007 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2008 ierr = PetscFree(s->numFieldComponents);CHKERRQ(ierr); 2009 for (f = 0; f < s->numFields; ++f) { 2010 ierr = PetscSectionDestroy(&s->field[f]);CHKERRQ(ierr); 2011 ierr = PetscFree(s->fieldNames[f]);CHKERRQ(ierr); 2012 } 2013 ierr = PetscFree(s->fieldNames);CHKERRQ(ierr); 2014 ierr = PetscFree(s->field);CHKERRQ(ierr); 2015 ierr = PetscSectionDestroy(&s->bc);CHKERRQ(ierr); 2016 ierr = PetscFree(s->bcIndices);CHKERRQ(ierr); 2017 ierr = PetscFree2(s->atlasDof, s->atlasOff);CHKERRQ(ierr); 2018 ierr = PetscSectionDestroy(&s->clSection);CHKERRQ(ierr); 2019 ierr = ISDestroy(&s->clPoints);CHKERRQ(ierr); 2020 ierr = ISDestroy(&s->perm);CHKERRQ(ierr); 2021 ierr = PetscFree(s->clPerm);CHKERRQ(ierr); 2022 ierr = PetscFree(s->clInvPerm);CHKERRQ(ierr); 2023 ierr = PetscSectionSymDestroy(&s->sym);CHKERRQ(ierr); 2024 2025 s->pStart = -1; 2026 s->pEnd = -1; 2027 s->maxDof = 0; 2028 s->setup = PETSC_FALSE; 2029 s->numFields = 0; 2030 s->clObj = NULL; 2031 PetscFunctionReturn(0); 2032 } 2033 2034 /*@ 2035 PetscSectionDestroy - Frees a section object and frees its range if that exists. 2036 2037 Not collective 2038 2039 Input Parameters: 2040 . s - the PetscSection 2041 2042 Level: beginner 2043 2044 .seealso: PetscSection, PetscSectionCreate() 2045 @*/ 2046 PetscErrorCode PetscSectionDestroy(PetscSection *s) 2047 { 2048 PetscErrorCode ierr; 2049 2050 PetscFunctionBegin; 2051 if (!*s) PetscFunctionReturn(0); 2052 PetscValidHeaderSpecific((*s),PETSC_SECTION_CLASSID,1); 2053 if (--((PetscObject)(*s))->refct > 0) { 2054 *s = NULL; 2055 PetscFunctionReturn(0); 2056 } 2057 ierr = PetscSectionReset(*s);CHKERRQ(ierr); 2058 ierr = PetscHeaderDestroy(s);CHKERRQ(ierr); 2059 PetscFunctionReturn(0); 2060 } 2061 2062 PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values) 2063 { 2064 const PetscInt p = point - s->pStart; 2065 2066 PetscFunctionBegin; 2067 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2); 2068 *values = &baseArray[s->atlasOff[p]]; 2069 PetscFunctionReturn(0); 2070 } 2071 2072 PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode) 2073 { 2074 PetscInt *array; 2075 const PetscInt p = point - s->pStart; 2076 const PetscInt orientation = 0; /* Needs to be included for use in closure operations */ 2077 PetscInt cDim = 0; 2078 PetscErrorCode ierr; 2079 2080 PetscFunctionBegin; 2081 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2); 2082 ierr = PetscSectionGetConstraintDof(s, p, &cDim);CHKERRQ(ierr); 2083 array = &baseArray[s->atlasOff[p]]; 2084 if (!cDim) { 2085 if (orientation >= 0) { 2086 const PetscInt dim = s->atlasDof[p]; 2087 PetscInt i; 2088 2089 if (mode == INSERT_VALUES) { 2090 for (i = 0; i < dim; ++i) array[i] = values[i]; 2091 } else { 2092 for (i = 0; i < dim; ++i) array[i] += values[i]; 2093 } 2094 } else { 2095 PetscInt offset = 0; 2096 PetscInt j = -1, field, i; 2097 2098 for (field = 0; field < s->numFields; ++field) { 2099 const PetscInt dim = s->field[field]->atlasDof[p]; 2100 2101 for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset]; 2102 offset += dim; 2103 } 2104 } 2105 } else { 2106 if (orientation >= 0) { 2107 const PetscInt dim = s->atlasDof[p]; 2108 PetscInt cInd = 0, i; 2109 const PetscInt *cDof; 2110 2111 ierr = PetscSectionGetConstraintIndices(s, point, &cDof);CHKERRQ(ierr); 2112 if (mode == INSERT_VALUES) { 2113 for (i = 0; i < dim; ++i) { 2114 if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;} 2115 array[i] = values[i]; 2116 } 2117 } else { 2118 for (i = 0; i < dim; ++i) { 2119 if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;} 2120 array[i] += values[i]; 2121 } 2122 } 2123 } else { 2124 const PetscInt *cDof; 2125 PetscInt offset = 0; 2126 PetscInt cOffset = 0; 2127 PetscInt j = 0, field; 2128 2129 ierr = PetscSectionGetConstraintIndices(s, point, &cDof);CHKERRQ(ierr); 2130 for (field = 0; field < s->numFields; ++field) { 2131 const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */ 2132 const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */ 2133 const PetscInt sDim = dim - tDim; 2134 PetscInt cInd = 0, i,k; 2135 2136 for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) { 2137 if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;} 2138 array[j] = values[k]; 2139 } 2140 offset += dim; 2141 cOffset += dim - tDim; 2142 } 2143 } 2144 } 2145 PetscFunctionReturn(0); 2146 } 2147 2148 /*@C 2149 PetscSectionHasConstraints - Determine whether a section has constrained dofs 2150 2151 Not collective 2152 2153 Input Parameter: 2154 . s - The PetscSection 2155 2156 Output Parameter: 2157 . hasConstraints - flag indicating that the section has constrained dofs 2158 2159 Level: intermediate 2160 2161 .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection 2162 @*/ 2163 PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints) 2164 { 2165 PetscFunctionBegin; 2166 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2167 PetscValidPointer(hasConstraints, 2); 2168 *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE; 2169 PetscFunctionReturn(0); 2170 } 2171 2172 /*@C 2173 PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained 2174 2175 Not collective 2176 2177 Input Parameters: 2178 + s - The PetscSection 2179 - point - The point 2180 2181 Output Parameter: 2182 . indices - The constrained dofs 2183 2184 Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90() 2185 2186 Level: intermediate 2187 2188 .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection 2189 @*/ 2190 PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices) 2191 { 2192 PetscErrorCode ierr; 2193 2194 PetscFunctionBegin; 2195 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2196 if (s->bc) { 2197 ierr = VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);CHKERRQ(ierr); 2198 } else *indices = NULL; 2199 PetscFunctionReturn(0); 2200 } 2201 2202 /*@C 2203 PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained 2204 2205 Not collective 2206 2207 Input Parameters: 2208 + s - The PetscSection 2209 . point - The point 2210 - indices - The constrained dofs 2211 2212 Note: The Fortran is PetscSectionSetConstraintIndicesF90() 2213 2214 Level: intermediate 2215 2216 .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection 2217 @*/ 2218 PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[]) 2219 { 2220 PetscErrorCode ierr; 2221 2222 PetscFunctionBegin; 2223 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2224 if (s->bc) { 2225 ierr = VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);CHKERRQ(ierr); 2226 } 2227 PetscFunctionReturn(0); 2228 } 2229 2230 /*@C 2231 PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained 2232 2233 Not collective 2234 2235 Input Parameters: 2236 + s - The PetscSection 2237 .field - The field number 2238 - point - The point 2239 2240 Output Parameter: 2241 . indices - The constrained dofs 2242 2243 Note: In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90() 2244 2245 Level: intermediate 2246 2247 .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection 2248 @*/ 2249 PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices) 2250 { 2251 PetscErrorCode ierr; 2252 2253 PetscFunctionBegin; 2254 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2255 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); 2256 ierr = PetscSectionGetConstraintIndices(s->field[field], point, indices);CHKERRQ(ierr); 2257 PetscFunctionReturn(0); 2258 } 2259 2260 /*@C 2261 PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained 2262 2263 Not collective 2264 2265 Input Parameters: 2266 + s - The PetscSection 2267 . point - The point 2268 . field - The field number 2269 - indices - The constrained dofs 2270 2271 Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90() 2272 2273 Level: intermediate 2274 2275 .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection 2276 @*/ 2277 PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[]) 2278 { 2279 PetscErrorCode ierr; 2280 2281 PetscFunctionBegin; 2282 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 2283 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); 2284 ierr = PetscSectionSetConstraintIndices(s->field[field], point, indices);CHKERRQ(ierr); 2285 PetscFunctionReturn(0); 2286 } 2287 2288 /*@ 2289 PetscSectionPermute - Reorder the section according to the input point permutation 2290 2291 Collective on PetscSection 2292 2293 Input Parameter: 2294 + section - The PetscSection object 2295 - perm - The point permutation, old point p becomes new point perm[p] 2296 2297 Output Parameter: 2298 . sectionNew - The permuted PetscSection 2299 2300 Level: intermediate 2301 2302 .seealso: MatPermute() 2303 @*/ 2304 PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew) 2305 { 2306 PetscSection s = section, sNew; 2307 const PetscInt *perm; 2308 PetscInt numFields, f, numPoints, pStart, pEnd, p; 2309 PetscErrorCode ierr; 2310 2311 PetscFunctionBegin; 2312 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1); 2313 PetscValidHeaderSpecific(permutation, IS_CLASSID, 2); 2314 PetscValidPointer(sectionNew, 3); 2315 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);CHKERRQ(ierr); 2316 ierr = PetscSectionGetNumFields(s, &numFields);CHKERRQ(ierr); 2317 if (numFields) {ierr = PetscSectionSetNumFields(sNew, numFields);CHKERRQ(ierr);} 2318 for (f = 0; f < numFields; ++f) { 2319 const char *name; 2320 PetscInt numComp; 2321 2322 ierr = PetscSectionGetFieldName(s, f, &name);CHKERRQ(ierr); 2323 ierr = PetscSectionSetFieldName(sNew, f, name);CHKERRQ(ierr); 2324 ierr = PetscSectionGetFieldComponents(s, f, &numComp);CHKERRQ(ierr); 2325 ierr = PetscSectionSetFieldComponents(sNew, f, numComp);CHKERRQ(ierr); 2326 } 2327 ierr = ISGetLocalSize(permutation, &numPoints);CHKERRQ(ierr); 2328 ierr = ISGetIndices(permutation, &perm);CHKERRQ(ierr); 2329 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 2330 ierr = PetscSectionSetChart(sNew, pStart, pEnd);CHKERRQ(ierr); 2331 if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd); 2332 for (p = pStart; p < pEnd; ++p) { 2333 PetscInt dof, cdof; 2334 2335 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr); 2336 ierr = PetscSectionSetDof(sNew, perm[p], dof);CHKERRQ(ierr); 2337 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 2338 if (cdof) {ierr = PetscSectionSetConstraintDof(sNew, perm[p], cdof);CHKERRQ(ierr);} 2339 for (f = 0; f < numFields; ++f) { 2340 ierr = PetscSectionGetFieldDof(s, p, f, &dof);CHKERRQ(ierr); 2341 ierr = PetscSectionSetFieldDof(sNew, perm[p], f, dof);CHKERRQ(ierr); 2342 ierr = PetscSectionGetFieldConstraintDof(s, p, f, &cdof);CHKERRQ(ierr); 2343 if (cdof) {ierr = PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);CHKERRQ(ierr);} 2344 } 2345 } 2346 ierr = PetscSectionSetUp(sNew);CHKERRQ(ierr); 2347 for (p = pStart; p < pEnd; ++p) { 2348 const PetscInt *cind; 2349 PetscInt cdof; 2350 2351 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 2352 if (cdof) { 2353 ierr = PetscSectionGetConstraintIndices(s, p, &cind);CHKERRQ(ierr); 2354 ierr = PetscSectionSetConstraintIndices(sNew, perm[p], cind);CHKERRQ(ierr); 2355 } 2356 for (f = 0; f < numFields; ++f) { 2357 ierr = PetscSectionGetFieldConstraintDof(s, p, f, &cdof);CHKERRQ(ierr); 2358 if (cdof) { 2359 ierr = PetscSectionGetFieldConstraintIndices(s, p, f, &cind);CHKERRQ(ierr); 2360 ierr = PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);CHKERRQ(ierr); 2361 } 2362 } 2363 } 2364 ierr = ISRestoreIndices(permutation, &perm);CHKERRQ(ierr); 2365 *sectionNew = sNew; 2366 PetscFunctionReturn(0); 2367 } 2368 2369 /* TODO: the next three functions should be moved to sf/utils */ 2370 #include <petsc/private/sfimpl.h> 2371 2372 /*@C 2373 PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF 2374 2375 Collective on sf 2376 2377 Input Parameters: 2378 + sf - The SF 2379 - rootSection - Section defined on root space 2380 2381 Output Parameters: 2382 + remoteOffsets - root offsets in leaf storage, or NULL 2383 - leafSection - Section defined on the leaf space 2384 2385 Level: advanced 2386 2387 .seealso: PetscSFCreate() 2388 @*/ 2389 PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection) 2390 { 2391 PetscSF embedSF; 2392 const PetscInt *indices; 2393 IS selected; 2394 PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f; 2395 PetscBool sub,lsub; 2396 PetscErrorCode ierr; 2397 2398 PetscFunctionBegin; 2399 ierr = PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);CHKERRQ(ierr); 2400 ierr = PetscSectionGetNumFields(rootSection, &numFields);CHKERRQ(ierr); 2401 if (numFields) {ierr = PetscSectionSetNumFields(leafSection, numFields);CHKERRQ(ierr);} 2402 for (f = 0; f < numFields; ++f) { 2403 PetscInt numComp = 0; 2404 ierr = PetscSectionGetFieldComponents(rootSection, f, &numComp);CHKERRQ(ierr); 2405 ierr = PetscSectionSetFieldComponents(leafSection, f, numComp);CHKERRQ(ierr); 2406 } 2407 ierr = PetscSectionGetChart(rootSection, &rpStart, &rpEnd);CHKERRQ(ierr); 2408 ierr = PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);CHKERRQ(ierr); 2409 rpEnd = PetscMin(rpEnd,nroots); 2410 rpEnd = PetscMax(rpStart,rpEnd); 2411 /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */ 2412 lsub = (PetscBool)(nroots != rpEnd - rpStart); 2413 ierr = MPIU_Allreduce(&lsub,&sub,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)sf));CHKERRQ(ierr); 2414 if (sub) { 2415 ierr = ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);CHKERRQ(ierr); 2416 ierr = ISGetIndices(selected, &indices);CHKERRQ(ierr); 2417 ierr = PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);CHKERRQ(ierr); 2418 ierr = ISRestoreIndices(selected, &indices);CHKERRQ(ierr); 2419 ierr = ISDestroy(&selected);CHKERRQ(ierr); 2420 } else { 2421 ierr = PetscObjectReference((PetscObject)sf);CHKERRQ(ierr); 2422 embedSF = sf; 2423 } 2424 ierr = PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);CHKERRQ(ierr); 2425 lpEnd++; 2426 2427 ierr = PetscSectionSetChart(leafSection, lpStart, lpEnd);CHKERRQ(ierr); 2428 /* Could fuse these at the cost of a copy and extra allocation */ 2429 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);CHKERRQ(ierr); 2430 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);CHKERRQ(ierr); 2431 for (f = 0; f < numFields; ++f) { 2432 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);CHKERRQ(ierr); 2433 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);CHKERRQ(ierr); 2434 } 2435 if (remoteOffsets) { 2436 ierr = PetscMalloc1(lpEnd - lpStart, remoteOffsets);CHKERRQ(ierr); 2437 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr); 2438 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr); 2439 } 2440 ierr = PetscSFDestroy(&embedSF);CHKERRQ(ierr); 2441 ierr = PetscSectionSetUp(leafSection);CHKERRQ(ierr); 2442 ierr = PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);CHKERRQ(ierr); 2443 PetscFunctionReturn(0); 2444 } 2445 2446 /*@C 2447 PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes 2448 2449 Collective on sf 2450 2451 Input Parameters: 2452 + sf - The SF 2453 . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots) 2454 - leafSection - Data layout of local points for incoming data (this is layout for SF leaves) 2455 2456 Output Parameter: 2457 . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL 2458 2459 Level: developer 2460 2461 .seealso: PetscSFCreate() 2462 @*/ 2463 PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets) 2464 { 2465 PetscSF embedSF; 2466 const PetscInt *indices; 2467 IS selected; 2468 PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0; 2469 PetscErrorCode ierr; 2470 2471 PetscFunctionBegin; 2472 *remoteOffsets = NULL; 2473 ierr = PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);CHKERRQ(ierr); 2474 if (numRoots < 0) PetscFunctionReturn(0); 2475 ierr = PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);CHKERRQ(ierr); 2476 ierr = PetscSectionGetChart(rootSection, &rpStart, &rpEnd);CHKERRQ(ierr); 2477 ierr = PetscSectionGetChart(leafSection, &lpStart, &lpEnd);CHKERRQ(ierr); 2478 ierr = ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);CHKERRQ(ierr); 2479 ierr = ISGetIndices(selected, &indices);CHKERRQ(ierr); 2480 ierr = PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);CHKERRQ(ierr); 2481 ierr = ISRestoreIndices(selected, &indices);CHKERRQ(ierr); 2482 ierr = ISDestroy(&selected);CHKERRQ(ierr); 2483 ierr = PetscCalloc1(lpEnd - lpStart, remoteOffsets);CHKERRQ(ierr); 2484 ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr); 2485 ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr); 2486 ierr = PetscSFDestroy(&embedSF);CHKERRQ(ierr); 2487 ierr = PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);CHKERRQ(ierr); 2488 PetscFunctionReturn(0); 2489 } 2490 2491 /*@C 2492 PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points 2493 2494 Collective on sf 2495 2496 Input Parameters: 2497 + sf - The SF 2498 . rootSection - Data layout of remote points for outgoing data (this is usually the serial section) 2499 . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL 2500 - leafSection - Data layout of local points for incoming data (this is the distributed section) 2501 2502 Output Parameters: 2503 - sectionSF - The new SF 2504 2505 Note: Either rootSection or remoteOffsets can be specified 2506 2507 Level: advanced 2508 2509 .seealso: PetscSFCreate() 2510 @*/ 2511 PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF) 2512 { 2513 MPI_Comm comm; 2514 const PetscInt *localPoints; 2515 const PetscSFNode *remotePoints; 2516 PetscInt lpStart, lpEnd; 2517 PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0; 2518 PetscInt *localIndices; 2519 PetscSFNode *remoteIndices; 2520 PetscInt i, ind; 2521 PetscErrorCode ierr; 2522 2523 PetscFunctionBegin; 2524 PetscValidHeaderSpecific(sf,PETSCSF_CLASSID,1); 2525 PetscValidPointer(rootSection,2); 2526 /* Cannot check PetscValidIntPointer(remoteOffsets,3) because it can be NULL if sf does not reference any points in leafSection */ 2527 PetscValidPointer(leafSection,4); 2528 PetscValidPointer(sectionSF,5); 2529 ierr = PetscObjectGetComm((PetscObject)sf,&comm);CHKERRQ(ierr); 2530 ierr = PetscSFCreate(comm, sectionSF);CHKERRQ(ierr); 2531 ierr = PetscSectionGetChart(leafSection, &lpStart, &lpEnd);CHKERRQ(ierr); 2532 ierr = PetscSectionGetStorageSize(rootSection, &numSectionRoots);CHKERRQ(ierr); 2533 ierr = PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);CHKERRQ(ierr); 2534 if (numRoots < 0) PetscFunctionReturn(0); 2535 ierr = PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);CHKERRQ(ierr); 2536 for (i = 0; i < numPoints; ++i) { 2537 PetscInt localPoint = localPoints ? localPoints[i] : i; 2538 PetscInt dof; 2539 2540 if ((localPoint >= lpStart) && (localPoint < lpEnd)) { 2541 ierr = PetscSectionGetDof(leafSection, localPoint, &dof);CHKERRQ(ierr); 2542 numIndices += dof; 2543 } 2544 } 2545 ierr = PetscMalloc1(numIndices, &localIndices);CHKERRQ(ierr); 2546 ierr = PetscMalloc1(numIndices, &remoteIndices);CHKERRQ(ierr); 2547 /* Create new index graph */ 2548 for (i = 0, ind = 0; i < numPoints; ++i) { 2549 PetscInt localPoint = localPoints ? localPoints[i] : i; 2550 PetscInt rank = remotePoints[i].rank; 2551 2552 if ((localPoint >= lpStart) && (localPoint < lpEnd)) { 2553 PetscInt remoteOffset = remoteOffsets[localPoint-lpStart]; 2554 PetscInt loff, dof, d; 2555 2556 ierr = PetscSectionGetOffset(leafSection, localPoint, &loff);CHKERRQ(ierr); 2557 ierr = PetscSectionGetDof(leafSection, localPoint, &dof);CHKERRQ(ierr); 2558 for (d = 0; d < dof; ++d, ++ind) { 2559 localIndices[ind] = loff+d; 2560 remoteIndices[ind].rank = rank; 2561 remoteIndices[ind].index = remoteOffset+d; 2562 } 2563 } 2564 } 2565 if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices); 2566 ierr = PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);CHKERRQ(ierr); 2567 ierr = PetscSFSetUp(*sectionSF);CHKERRQ(ierr); 2568 ierr = PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);CHKERRQ(ierr); 2569 PetscFunctionReturn(0); 2570 } 2571 2572 /*@ 2573 PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section 2574 2575 Collective on section 2576 2577 Input Parameters: 2578 + section - The PetscSection 2579 . obj - A PetscObject which serves as the key for this index 2580 . clSection - Section giving the size of the closure of each point 2581 - clPoints - IS giving the points in each closure 2582 2583 Note: We compress out closure points with no dofs in this section 2584 2585 Level: advanced 2586 2587 .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex() 2588 @*/ 2589 PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints) 2590 { 2591 PetscErrorCode ierr; 2592 2593 PetscFunctionBegin; 2594 if (section->clObj != obj) {ierr = PetscFree(section->clPerm);CHKERRQ(ierr);ierr = PetscFree(section->clInvPerm);CHKERRQ(ierr);} 2595 section->clObj = obj; 2596 ierr = PetscSectionDestroy(§ion->clSection);CHKERRQ(ierr); 2597 ierr = ISDestroy(§ion->clPoints);CHKERRQ(ierr); 2598 section->clSection = clSection; 2599 section->clPoints = clPoints; 2600 PetscFunctionReturn(0); 2601 } 2602 2603 /*@ 2604 PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section 2605 2606 Collective on section 2607 2608 Input Parameters: 2609 + section - The PetscSection 2610 - obj - A PetscObject which serves as the key for this index 2611 2612 Output Parameters: 2613 + clSection - Section giving the size of the closure of each point 2614 - clPoints - IS giving the points in each closure 2615 2616 Note: We compress out closure points with no dofs in this section 2617 2618 Level: advanced 2619 2620 .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex() 2621 @*/ 2622 PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints) 2623 { 2624 PetscFunctionBegin; 2625 if (section->clObj == obj) { 2626 if (clSection) *clSection = section->clSection; 2627 if (clPoints) *clPoints = section->clPoints; 2628 } else { 2629 if (clSection) *clSection = NULL; 2630 if (clPoints) *clPoints = NULL; 2631 } 2632 PetscFunctionReturn(0); 2633 } 2634 2635 PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm) 2636 { 2637 PetscInt i; 2638 PetscErrorCode ierr; 2639 2640 PetscFunctionBegin; 2641 if (section->clObj != obj) { 2642 ierr = PetscSectionDestroy(§ion->clSection);CHKERRQ(ierr); 2643 ierr = ISDestroy(§ion->clPoints);CHKERRQ(ierr); 2644 } 2645 section->clObj = obj; 2646 ierr = PetscFree(section->clPerm);CHKERRQ(ierr); 2647 ierr = PetscFree(section->clInvPerm);CHKERRQ(ierr); 2648 section->clSize = clSize; 2649 if (mode == PETSC_COPY_VALUES) { 2650 ierr = PetscMalloc1(clSize, §ion->clPerm);CHKERRQ(ierr); 2651 ierr = PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));CHKERRQ(ierr); 2652 ierr = PetscArraycpy(section->clPerm, clPerm, clSize);CHKERRQ(ierr); 2653 } else if (mode == PETSC_OWN_POINTER) { 2654 section->clPerm = clPerm; 2655 } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays"); 2656 ierr = PetscMalloc1(clSize, §ion->clInvPerm);CHKERRQ(ierr); 2657 for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i; 2658 PetscFunctionReturn(0); 2659 } 2660 2661 /*@ 2662 PetscSectionSetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex. 2663 2664 Not Collective 2665 2666 Input Parameters: 2667 + section - The PetscSection 2668 . obj - A PetscObject which serves as the key for this index 2669 - perm - Permutation of the cell dof closure 2670 2671 Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for 2672 other points (like faces). 2673 2674 Level: intermediate 2675 2676 .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode 2677 @*/ 2678 PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm) 2679 { 2680 const PetscInt *clPerm = NULL; 2681 PetscInt clSize = 0; 2682 PetscErrorCode ierr; 2683 2684 PetscFunctionBegin; 2685 if (perm) { 2686 ierr = ISGetLocalSize(perm, &clSize);CHKERRQ(ierr); 2687 ierr = ISGetIndices(perm, &clPerm);CHKERRQ(ierr); 2688 } 2689 ierr = PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);CHKERRQ(ierr); 2690 if (perm) {ierr = ISRestoreIndices(perm, &clPerm);CHKERRQ(ierr);} 2691 PetscFunctionReturn(0); 2692 } 2693 2694 PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[]) 2695 { 2696 PetscFunctionBegin; 2697 if (section->clObj == obj) { 2698 if (size) *size = section->clSize; 2699 if (perm) *perm = section->clPerm; 2700 } else { 2701 if (size) *size = 0; 2702 if (perm) *perm = NULL; 2703 } 2704 PetscFunctionReturn(0); 2705 } 2706 2707 /*@ 2708 PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex. 2709 2710 Not collective 2711 2712 Input Parameters: 2713 + section - The PetscSection 2714 - obj - A PetscObject which serves as the key for this index 2715 2716 Output Parameter: 2717 . perm - The dof closure permutation 2718 2719 Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for 2720 other points (like faces). 2721 2722 The user must destroy the IS that is returned. 2723 2724 Level: intermediate 2725 2726 .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex() 2727 @*/ 2728 PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm) 2729 { 2730 const PetscInt *clPerm; 2731 PetscInt clSize; 2732 PetscErrorCode ierr; 2733 2734 PetscFunctionBegin; 2735 ierr = PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);CHKERRQ(ierr); 2736 ierr = ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);CHKERRQ(ierr); 2737 PetscFunctionReturn(0); 2738 } 2739 2740 PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[]) 2741 { 2742 PetscFunctionBegin; 2743 if (section->clObj == obj) { 2744 if (size) *size = section->clSize; 2745 if (perm) *perm = section->clInvPerm; 2746 } else { 2747 if (size) *size = 0; 2748 if (perm) *perm = NULL; 2749 } 2750 PetscFunctionReturn(0); 2751 } 2752 2753 /*@ 2754 PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex. 2755 2756 Not collective 2757 2758 Input Parameters: 2759 + section - The PetscSection 2760 - obj - A PetscObject which serves as the key for this index 2761 2762 Output Parameters: 2763 + size - The dof closure size 2764 - perm - The dof closure permutation 2765 2766 Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for 2767 other points (like faces). 2768 2769 The user must destroy the IS that is returned. 2770 2771 Level: intermediate 2772 2773 .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex() 2774 @*/ 2775 PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm) 2776 { 2777 const PetscInt *clPerm; 2778 PetscInt clSize; 2779 PetscErrorCode ierr; 2780 2781 PetscFunctionBegin; 2782 ierr = PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);CHKERRQ(ierr); 2783 ierr = ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);CHKERRQ(ierr); 2784 PetscFunctionReturn(0); 2785 } 2786 2787 /*@ 2788 PetscSectionGetField - Get the subsection associated with a single field 2789 2790 Input Parameters: 2791 + s - The PetscSection 2792 - field - The field number 2793 2794 Output Parameter: 2795 . subs - The subsection for the given field 2796 2797 Level: intermediate 2798 2799 .seealso: PetscSectionSetNumFields() 2800 @*/ 2801 PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs) 2802 { 2803 PetscFunctionBegin; 2804 PetscValidHeaderSpecific(s,PETSC_SECTION_CLASSID,1); 2805 PetscValidPointer(subs,3); 2806 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); 2807 *subs = s->field[field]; 2808 PetscFunctionReturn(0); 2809 } 2810 2811 PetscClassId PETSC_SECTION_SYM_CLASSID; 2812 PetscFunctionList PetscSectionSymList = NULL; 2813 2814 /*@ 2815 PetscSectionSymCreate - Creates an empty PetscSectionSym object. 2816 2817 Collective 2818 2819 Input Parameter: 2820 . comm - the MPI communicator 2821 2822 Output Parameter: 2823 . sym - pointer to the new set of symmetries 2824 2825 Level: developer 2826 2827 .seealso: PetscSectionSym, PetscSectionSymDestroy() 2828 @*/ 2829 PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym) 2830 { 2831 PetscErrorCode ierr; 2832 2833 PetscFunctionBegin; 2834 PetscValidPointer(sym,2); 2835 ierr = ISInitializePackage();CHKERRQ(ierr); 2836 ierr = PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);CHKERRQ(ierr); 2837 PetscFunctionReturn(0); 2838 } 2839 2840 /*@C 2841 PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation. 2842 2843 Collective on PetscSectionSym 2844 2845 Input Parameters: 2846 + sym - The section symmetry object 2847 - method - The name of the section symmetry type 2848 2849 Level: developer 2850 2851 .seealso: PetscSectionSymGetType(), PetscSectionSymCreate() 2852 @*/ 2853 PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method) 2854 { 2855 PetscErrorCode (*r)(PetscSectionSym); 2856 PetscBool match; 2857 PetscErrorCode ierr; 2858 2859 PetscFunctionBegin; 2860 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID,1); 2861 ierr = PetscObjectTypeCompare((PetscObject) sym, method, &match);CHKERRQ(ierr); 2862 if (match) PetscFunctionReturn(0); 2863 2864 ierr = PetscFunctionListFind(PetscSectionSymList,method,&r);CHKERRQ(ierr); 2865 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method); 2866 if (sym->ops->destroy) { 2867 ierr = (*sym->ops->destroy)(sym);CHKERRQ(ierr); 2868 sym->ops->destroy = NULL; 2869 } 2870 ierr = (*r)(sym);CHKERRQ(ierr); 2871 ierr = PetscObjectChangeTypeName((PetscObject)sym,method);CHKERRQ(ierr); 2872 PetscFunctionReturn(0); 2873 } 2874 2875 2876 /*@C 2877 PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym. 2878 2879 Not Collective 2880 2881 Input Parameter: 2882 . sym - The section symmetry 2883 2884 Output Parameter: 2885 . type - The index set type name 2886 2887 Level: developer 2888 2889 .seealso: PetscSectionSymSetType(), PetscSectionSymCreate() 2890 @*/ 2891 PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type) 2892 { 2893 PetscFunctionBegin; 2894 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID,1); 2895 PetscValidCharPointer(type,2); 2896 *type = ((PetscObject)sym)->type_name; 2897 PetscFunctionReturn(0); 2898 } 2899 2900 /*@C 2901 PetscSectionSymRegister - Adds a new section symmetry implementation 2902 2903 Not Collective 2904 2905 Input Parameters: 2906 + name - The name of a new user-defined creation routine 2907 - create_func - The creation routine itself 2908 2909 Notes: 2910 PetscSectionSymRegister() may be called multiple times to add several user-defined vectors 2911 2912 Level: developer 2913 2914 .seealso: PetscSectionSymCreate(), PetscSectionSymSetType() 2915 @*/ 2916 PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym)) 2917 { 2918 PetscErrorCode ierr; 2919 2920 PetscFunctionBegin; 2921 ierr = ISInitializePackage();CHKERRQ(ierr); 2922 ierr = PetscFunctionListAdd(&PetscSectionSymList,sname,function);CHKERRQ(ierr); 2923 PetscFunctionReturn(0); 2924 } 2925 2926 /*@ 2927 PetscSectionSymDestroy - Destroys a section symmetry. 2928 2929 Collective on PetscSectionSym 2930 2931 Input Parameters: 2932 . sym - the section symmetry 2933 2934 Level: developer 2935 2936 .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy() 2937 @*/ 2938 PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym) 2939 { 2940 SymWorkLink link,next; 2941 PetscErrorCode ierr; 2942 2943 PetscFunctionBegin; 2944 if (!*sym) PetscFunctionReturn(0); 2945 PetscValidHeaderSpecific((*sym),PETSC_SECTION_SYM_CLASSID,1); 2946 if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; PetscFunctionReturn(0);} 2947 if ((*sym)->ops->destroy) { 2948 ierr = (*(*sym)->ops->destroy)(*sym);CHKERRQ(ierr); 2949 } 2950 if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out"); 2951 for (link=(*sym)->workin; link; link=next) { 2952 next = link->next; 2953 ierr = PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);CHKERRQ(ierr); 2954 ierr = PetscFree(link);CHKERRQ(ierr); 2955 } 2956 (*sym)->workin = NULL; 2957 ierr = PetscHeaderDestroy(sym);CHKERRQ(ierr); 2958 PetscFunctionReturn(0); 2959 } 2960 2961 /*@C 2962 PetscSectionSymView - Displays a section symmetry 2963 2964 Collective on PetscSectionSym 2965 2966 Input Parameters: 2967 + sym - the index set 2968 - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF. 2969 2970 Level: developer 2971 2972 .seealso: PetscViewerASCIIOpen() 2973 @*/ 2974 PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer) 2975 { 2976 PetscErrorCode ierr; 2977 2978 PetscFunctionBegin; 2979 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1); 2980 if (!viewer) { 2981 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);CHKERRQ(ierr); 2982 } 2983 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 2984 PetscCheckSameComm(sym,1,viewer,2); 2985 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);CHKERRQ(ierr); 2986 if (sym->ops->view) { 2987 ierr = (*sym->ops->view)(sym,viewer);CHKERRQ(ierr); 2988 } 2989 PetscFunctionReturn(0); 2990 } 2991 2992 /*@ 2993 PetscSectionSetSym - Set the symmetries for the data referred to by the section 2994 2995 Collective on PetscSection 2996 2997 Input Parameters: 2998 + section - the section describing data layout 2999 - sym - the symmetry describing the affect of orientation on the access of the data 3000 3001 Level: developer 3002 3003 .seealso: PetscSectionGetSym(), PetscSectionSymCreate() 3004 @*/ 3005 PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym) 3006 { 3007 PetscErrorCode ierr; 3008 3009 PetscFunctionBegin; 3010 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3011 ierr = PetscSectionSymDestroy(&(section->sym));CHKERRQ(ierr); 3012 if (sym) { 3013 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,2); 3014 PetscCheckSameComm(section,1,sym,2); 3015 ierr = PetscObjectReference((PetscObject) sym);CHKERRQ(ierr); 3016 } 3017 section->sym = sym; 3018 PetscFunctionReturn(0); 3019 } 3020 3021 /*@ 3022 PetscSectionGetSym - Get the symmetries for the data referred to by the section 3023 3024 Not collective 3025 3026 Input Parameters: 3027 . section - the section describing data layout 3028 3029 Output Parameters: 3030 . sym - the symmetry describing the affect of orientation on the access of the data 3031 3032 Level: developer 3033 3034 .seealso: PetscSectionSetSym(), PetscSectionSymCreate() 3035 @*/ 3036 PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym) 3037 { 3038 PetscFunctionBegin; 3039 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3040 *sym = section->sym; 3041 PetscFunctionReturn(0); 3042 } 3043 3044 /*@ 3045 PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section 3046 3047 Collective on PetscSection 3048 3049 Input Parameters: 3050 + section - the section describing data layout 3051 . field - the field number 3052 - sym - the symmetry describing the affect of orientation on the access of the data 3053 3054 Level: developer 3055 3056 .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate() 3057 @*/ 3058 PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym) 3059 { 3060 PetscErrorCode ierr; 3061 3062 PetscFunctionBegin; 3063 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3064 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); 3065 ierr = PetscSectionSetSym(section->field[field],sym);CHKERRQ(ierr); 3066 PetscFunctionReturn(0); 3067 } 3068 3069 /*@ 3070 PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section 3071 3072 Collective on PetscSection 3073 3074 Input Parameters: 3075 + section - the section describing data layout 3076 - field - the field number 3077 3078 Output Parameters: 3079 . sym - the symmetry describing the affect of orientation on the access of the data 3080 3081 Level: developer 3082 3083 .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate() 3084 @*/ 3085 PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym) 3086 { 3087 PetscFunctionBegin; 3088 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3089 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); 3090 *sym = section->field[field]->sym; 3091 PetscFunctionReturn(0); 3092 } 3093 3094 /*@C 3095 PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations. 3096 3097 Not collective 3098 3099 Input Parameters: 3100 + section - the section 3101 . numPoints - the number of points 3102 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3103 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3104 context, see DMPlexGetConeOrientation()). 3105 3106 Output Parameter: 3107 + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity). 3108 - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all 3109 identity). 3110 3111 Example of usage, gathering dofs into a local array (lArray) from a section array (sArray): 3112 .vb 3113 const PetscInt **perms; 3114 const PetscScalar **rots; 3115 PetscInt lOffset; 3116 3117 PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots); 3118 for (i = 0, lOffset = 0; i < numPoints; i++) { 3119 PetscInt point = points[2*i], dof, sOffset; 3120 const PetscInt *perm = perms ? perms[i] : NULL; 3121 const PetscScalar *rot = rots ? rots[i] : NULL; 3122 3123 PetscSectionGetDof(section,point,&dof); 3124 PetscSectionGetOffset(section,point,&sOffset); 3125 3126 if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}} 3127 else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}} 3128 if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }} 3129 lOffset += dof; 3130 } 3131 PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots); 3132 .ve 3133 3134 Example of usage, adding dofs into a section array (sArray) from a local array (lArray): 3135 .vb 3136 const PetscInt **perms; 3137 const PetscScalar **rots; 3138 PetscInt lOffset; 3139 3140 PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots); 3141 for (i = 0, lOffset = 0; i < numPoints; i++) { 3142 PetscInt point = points[2*i], dof, sOffset; 3143 const PetscInt *perm = perms ? perms[i] : NULL; 3144 const PetscScalar *rot = rots ? rots[i] : NULL; 3145 3146 PetscSectionGetDof(section,point,&dof); 3147 PetscSectionGetOffset(section,point,&sOff); 3148 3149 if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}} 3150 else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}} 3151 offset += dof; 3152 } 3153 PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots); 3154 .ve 3155 3156 Level: developer 3157 3158 .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym() 3159 @*/ 3160 PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3161 { 3162 PetscSectionSym sym; 3163 PetscErrorCode ierr; 3164 3165 PetscFunctionBegin; 3166 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3167 if (numPoints) PetscValidIntPointer(points,3); 3168 if (perms) *perms = NULL; 3169 if (rots) *rots = NULL; 3170 sym = section->sym; 3171 if (sym && (perms || rots)) { 3172 SymWorkLink link; 3173 3174 if (sym->workin) { 3175 link = sym->workin; 3176 sym->workin = sym->workin->next; 3177 } else { 3178 ierr = PetscNewLog(sym,&link);CHKERRQ(ierr); 3179 } 3180 if (numPoints > link->numPoints) { 3181 ierr = PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);CHKERRQ(ierr); 3182 ierr = PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);CHKERRQ(ierr); 3183 link->numPoints = numPoints; 3184 } 3185 link->next = sym->workout; 3186 sym->workout = link; 3187 ierr = PetscArrayzero((PetscInt**)link->perms,numPoints);CHKERRQ(ierr); 3188 ierr = PetscArrayzero((PetscInt**)link->rots,numPoints);CHKERRQ(ierr); 3189 ierr = (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);CHKERRQ(ierr); 3190 if (perms) *perms = link->perms; 3191 if (rots) *rots = link->rots; 3192 } 3193 PetscFunctionReturn(0); 3194 } 3195 3196 /*@C 3197 PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms() 3198 3199 Not collective 3200 3201 Input Parameters: 3202 + section - the section 3203 . numPoints - the number of points 3204 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3205 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3206 context, see DMPlexGetConeOrientation()). 3207 3208 Output Parameter: 3209 + perms - The permutations for the given orientations: set to NULL at conclusion 3210 - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion 3211 3212 Level: developer 3213 3214 .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym() 3215 @*/ 3216 PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3217 { 3218 PetscSectionSym sym; 3219 3220 PetscFunctionBegin; 3221 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3222 sym = section->sym; 3223 if (sym && (perms || rots)) { 3224 SymWorkLink *p,link; 3225 3226 for (p=&sym->workout; (link=*p); p=&link->next) { 3227 if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) { 3228 *p = link->next; 3229 link->next = sym->workin; 3230 sym->workin = link; 3231 if (perms) *perms = NULL; 3232 if (rots) *rots = NULL; 3233 PetscFunctionReturn(0); 3234 } 3235 } 3236 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out"); 3237 } 3238 PetscFunctionReturn(0); 3239 } 3240 3241 /*@C 3242 PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations. 3243 3244 Not collective 3245 3246 Input Parameters: 3247 + section - the section 3248 . field - the field of the section 3249 . numPoints - the number of points 3250 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3251 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3252 context, see DMPlexGetConeOrientation()). 3253 3254 Output Parameter: 3255 + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity). 3256 - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all 3257 identity). 3258 3259 Level: developer 3260 3261 .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms() 3262 @*/ 3263 PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3264 { 3265 PetscErrorCode ierr; 3266 3267 PetscFunctionBegin; 3268 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3269 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); 3270 ierr = PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);CHKERRQ(ierr); 3271 PetscFunctionReturn(0); 3272 } 3273 3274 /*@C 3275 PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms() 3276 3277 Not collective 3278 3279 Input Parameters: 3280 + section - the section 3281 . field - the field number 3282 . numPoints - the number of points 3283 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3284 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3285 context, see DMPlexGetConeOrientation()). 3286 3287 Output Parameter: 3288 + perms - The permutations for the given orientations: set to NULL at conclusion 3289 - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion 3290 3291 Level: developer 3292 3293 .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym() 3294 @*/ 3295 PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3296 { 3297 PetscErrorCode ierr; 3298 3299 PetscFunctionBegin; 3300 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3301 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); 3302 ierr = PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);CHKERRQ(ierr); 3303 PetscFunctionReturn(0); 3304 } 3305 3306 /*@ 3307 PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset 3308 3309 Not collective 3310 3311 Input Parameter: 3312 . s - the global PetscSection 3313 3314 Output Parameters: 3315 . flg - the flag 3316 3317 Level: developer 3318 3319 .seealso: PetscSectionSetChart(), PetscSectionCreate() 3320 @*/ 3321 PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg) 3322 { 3323 PetscFunctionBegin; 3324 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 3325 *flg = s->useFieldOff; 3326 PetscFunctionReturn(0); 3327 } 3328 3329 /*@ 3330 PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset 3331 3332 Not collective 3333 3334 Input Parameters: 3335 + s - the global PetscSection 3336 - flg - the flag 3337 3338 Level: developer 3339 3340 .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate() 3341 @*/ 3342 PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg) 3343 { 3344 PetscFunctionBegin; 3345 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 3346 s->useFieldOff = flg; 3347 PetscFunctionReturn(0); 3348 } 3349 3350 #define PetscSectionExpandPoints_Loop(TYPE) \ 3351 { \ 3352 PetscInt i, n, o0, o1, size; \ 3353 TYPE *a0 = (TYPE*)origArray, *a1; \ 3354 ierr = PetscSectionGetStorageSize(s, &size);CHKERRQ(ierr); \ 3355 ierr = PetscMalloc1(size, &a1);CHKERRQ(ierr); \ 3356 for (i=0; i<npoints; i++) { \ 3357 ierr = PetscSectionGetOffset(origSection, points_[i], &o0);CHKERRQ(ierr); \ 3358 ierr = PetscSectionGetOffset(s, i, &o1);CHKERRQ(ierr); \ 3359 ierr = PetscSectionGetDof(s, i, &n);CHKERRQ(ierr); \ 3360 ierr = PetscMemcpy(&a1[o1], &a0[o0], n*unitsize);CHKERRQ(ierr); \ 3361 } \ 3362 *newArray = (void*)a1; \ 3363 } 3364 3365 /*@ 3366 PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points. 3367 3368 Not collective 3369 3370 Input Parameters: 3371 + origSection - the PetscSection describing the layout of the array 3372 . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL) 3373 . origArray - the array; its size must be equal to the storage size of origSection 3374 - points - IS with points to extract; its indices must lie in the chart of origSection 3375 3376 Output Parameters: 3377 + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs) 3378 - newArray - the array of the extracted DOFs; its size is the storage size of newSection 3379 3380 Level: developer 3381 3382 .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate() 3383 @*/ 3384 PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[]) 3385 { 3386 PetscSection s; 3387 const PetscInt *points_; 3388 PetscInt i, n, npoints, pStart, pEnd; 3389 PetscMPIInt unitsize; 3390 PetscErrorCode ierr; 3391 3392 PetscFunctionBegin; 3393 PetscValidHeaderSpecific(origSection, PETSC_SECTION_CLASSID, 1); 3394 PetscValidPointer(origArray, 3); 3395 PetscValidHeaderSpecific(points, IS_CLASSID, 4); 3396 if (newSection) PetscValidPointer(newSection, 5); 3397 if (newArray) PetscValidPointer(newArray, 6); 3398 ierr = MPI_Type_size(dataType, &unitsize);CHKERRQ(ierr); 3399 ierr = ISGetLocalSize(points, &npoints);CHKERRQ(ierr); 3400 ierr = ISGetIndices(points, &points_);CHKERRQ(ierr); 3401 ierr = PetscSectionGetChart(origSection, &pStart, &pEnd);CHKERRQ(ierr); 3402 ierr = PetscSectionCreate(PETSC_COMM_SELF, &s);CHKERRQ(ierr); 3403 ierr = PetscSectionSetChart(s, 0, npoints);CHKERRQ(ierr); 3404 for (i=0; i<npoints; i++) { 3405 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); 3406 ierr = PetscSectionGetDof(origSection, points_[i], &n);CHKERRQ(ierr); 3407 ierr = PetscSectionSetDof(s, i, n);CHKERRQ(ierr); 3408 } 3409 ierr = PetscSectionSetUp(s);CHKERRQ(ierr); 3410 if (newArray) { 3411 if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);} 3412 else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);} 3413 else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);} 3414 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype"); 3415 } 3416 if (newSection) { 3417 *newSection = s; 3418 } else { 3419 ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); 3420 } 3421 PetscFunctionReturn(0); 3422 } 3423