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