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