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