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 /*@ 2516 PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section 2517 2518 Collective on section 2519 2520 Input Parameters: 2521 + section - The PetscSection 2522 . obj - A PetscObject which serves as the key for this index 2523 . clSection - Section giving the size of the closure of each point 2524 - clPoints - IS giving the points in each closure 2525 2526 Note: We compress out closure points with no dofs in this section 2527 2528 Level: advanced 2529 2530 .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex() 2531 @*/ 2532 PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints) 2533 { 2534 PetscErrorCode ierr; 2535 2536 PetscFunctionBegin; 2537 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 2538 PetscValidHeaderSpecific(clSection,PETSC_SECTION_CLASSID,3); 2539 PetscValidHeaderSpecific(clPoints,IS_CLASSID,4); 2540 if (section->clObj != obj) {ierr = PetscSectionResetClosurePermutation(section);CHKERRQ(ierr);} 2541 section->clObj = obj; 2542 ierr = PetscObjectReference((PetscObject)clSection);CHKERRQ(ierr); 2543 ierr = PetscObjectReference((PetscObject)clPoints);CHKERRQ(ierr); 2544 ierr = PetscSectionDestroy(§ion->clSection);CHKERRQ(ierr); 2545 ierr = ISDestroy(§ion->clPoints);CHKERRQ(ierr); 2546 section->clSection = clSection; 2547 section->clPoints = clPoints; 2548 PetscFunctionReturn(0); 2549 } 2550 2551 /*@ 2552 PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section 2553 2554 Collective on section 2555 2556 Input Parameters: 2557 + section - The PetscSection 2558 - obj - A PetscObject which serves as the key for this index 2559 2560 Output Parameters: 2561 + clSection - Section giving the size of the closure of each point 2562 - clPoints - IS giving the points in each closure 2563 2564 Note: We compress out closure points with no dofs in this section 2565 2566 Level: advanced 2567 2568 .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex() 2569 @*/ 2570 PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints) 2571 { 2572 PetscFunctionBegin; 2573 if (section->clObj == obj) { 2574 if (clSection) *clSection = section->clSection; 2575 if (clPoints) *clPoints = section->clPoints; 2576 } else { 2577 if (clSection) *clSection = NULL; 2578 if (clPoints) *clPoints = NULL; 2579 } 2580 PetscFunctionReturn(0); 2581 } 2582 2583 PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm) 2584 { 2585 PetscInt i; 2586 khiter_t iter; 2587 int new_entry; 2588 PetscSectionClosurePermKey key = {depth, clSize}; 2589 PetscSectionClosurePermVal *val; 2590 PetscErrorCode ierr; 2591 2592 PetscFunctionBegin; 2593 if (section->clObj != obj) { 2594 ierr = PetscSectionDestroy(§ion->clSection);CHKERRQ(ierr); 2595 ierr = ISDestroy(§ion->clPoints);CHKERRQ(ierr); 2596 } 2597 section->clObj = obj; 2598 if (!section->clHash) {ierr = PetscClPermCreate(§ion->clHash);CHKERRQ(ierr);} 2599 iter = kh_put(ClPerm, section->clHash, key, &new_entry); 2600 val = &kh_val(section->clHash, iter); 2601 if (!new_entry) { 2602 ierr = PetscFree(val->perm);CHKERRQ(ierr); 2603 ierr = PetscFree(val->invPerm);CHKERRQ(ierr); 2604 } 2605 if (mode == PETSC_COPY_VALUES) { 2606 ierr = PetscMalloc1(clSize, &val->perm);CHKERRQ(ierr); 2607 ierr = PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));CHKERRQ(ierr); 2608 ierr = PetscArraycpy(val->perm, clPerm, clSize);CHKERRQ(ierr); 2609 } else if (mode == PETSC_OWN_POINTER) { 2610 val->perm = clPerm; 2611 } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays"); 2612 ierr = PetscMalloc1(clSize, &val->invPerm);CHKERRQ(ierr); 2613 for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i; 2614 PetscFunctionReturn(0); 2615 } 2616 2617 /*@ 2618 PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex. 2619 2620 Not Collective 2621 2622 Input Parameters: 2623 + section - The PetscSection 2624 . obj - A PetscObject which serves as the key for this index (usually a DM) 2625 . depth - Depth of points on which to apply the given permutation 2626 - perm - Permutation of the cell dof closure 2627 2628 Note: 2629 The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a 2630 mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for 2631 each topology and degree. 2632 2633 This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for 2634 exotic/enriched spaces on mixed topology meshes. 2635 2636 Level: intermediate 2637 2638 .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode 2639 @*/ 2640 PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm) 2641 { 2642 const PetscInt *clPerm = NULL; 2643 PetscInt clSize = 0; 2644 PetscErrorCode ierr; 2645 2646 PetscFunctionBegin; 2647 if (perm) { 2648 ierr = ISGetLocalSize(perm, &clSize);CHKERRQ(ierr); 2649 ierr = ISGetIndices(perm, &clPerm);CHKERRQ(ierr); 2650 } 2651 ierr = PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);CHKERRQ(ierr); 2652 if (perm) {ierr = ISRestoreIndices(perm, &clPerm);CHKERRQ(ierr);} 2653 PetscFunctionReturn(0); 2654 } 2655 2656 PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[]) 2657 { 2658 PetscErrorCode ierr; 2659 2660 PetscFunctionBegin; 2661 if (section->clObj == obj) { 2662 PetscSectionClosurePermKey k = {depth, size}; 2663 PetscSectionClosurePermVal v; 2664 ierr = PetscClPermGet(section->clHash, k, &v);CHKERRQ(ierr); 2665 if (perm) *perm = v.perm; 2666 } else { 2667 if (perm) *perm = NULL; 2668 } 2669 PetscFunctionReturn(0); 2670 } 2671 2672 /*@ 2673 PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex. 2674 2675 Not collective 2676 2677 Input Parameters: 2678 + section - The PetscSection 2679 . obj - A PetscObject which serves as the key for this index (usually a DM) 2680 . depth - Depth stratum on which to obtain closure permutation 2681 - clSize - Closure size to be permuted (e.g., may vary with element topology and degree) 2682 2683 Output Parameter: 2684 . perm - The dof closure permutation 2685 2686 Note: 2687 The user must destroy the IS that is returned. 2688 2689 Level: intermediate 2690 2691 .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex() 2692 @*/ 2693 PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm) 2694 { 2695 const PetscInt *clPerm; 2696 PetscErrorCode ierr; 2697 2698 PetscFunctionBegin; 2699 ierr = PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);CHKERRQ(ierr); 2700 ierr = ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);CHKERRQ(ierr); 2701 PetscFunctionReturn(0); 2702 } 2703 2704 PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[]) 2705 { 2706 PetscErrorCode ierr; 2707 2708 PetscFunctionBegin; 2709 if (section->clObj == obj && section->clHash) { 2710 PetscSectionClosurePermKey k = {depth, size}; 2711 PetscSectionClosurePermVal v; 2712 ierr = PetscClPermGet(section->clHash, k, &v);CHKERRQ(ierr); 2713 if (perm) *perm = v.invPerm; 2714 } else { 2715 if (perm) *perm = NULL; 2716 } 2717 PetscFunctionReturn(0); 2718 } 2719 2720 /*@ 2721 PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex. 2722 2723 Not collective 2724 2725 Input Parameters: 2726 + section - The PetscSection 2727 . obj - A PetscObject which serves as the key for this index (usually a DM) 2728 . depth - Depth stratum on which to obtain closure permutation 2729 - clSize - Closure size to be permuted (e.g., may vary with element topology and degree) 2730 2731 Output Parameters: 2732 . perm - The dof closure permutation 2733 2734 Note: 2735 The user must destroy the IS that is returned. 2736 2737 Level: intermediate 2738 2739 .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex() 2740 @*/ 2741 PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm) 2742 { 2743 const PetscInt *clPerm; 2744 PetscErrorCode ierr; 2745 2746 PetscFunctionBegin; 2747 ierr = PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);CHKERRQ(ierr); 2748 ierr = ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);CHKERRQ(ierr); 2749 PetscFunctionReturn(0); 2750 } 2751 2752 /*@ 2753 PetscSectionGetField - Get the subsection associated with a single field 2754 2755 Input Parameters: 2756 + s - The PetscSection 2757 - field - The field number 2758 2759 Output Parameter: 2760 . subs - The subsection for the given field 2761 2762 Level: intermediate 2763 2764 .seealso: PetscSectionSetNumFields() 2765 @*/ 2766 PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs) 2767 { 2768 PetscFunctionBegin; 2769 PetscValidHeaderSpecific(s,PETSC_SECTION_CLASSID,1); 2770 PetscValidPointer(subs,3); 2771 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); 2772 *subs = s->field[field]; 2773 PetscFunctionReturn(0); 2774 } 2775 2776 PetscClassId PETSC_SECTION_SYM_CLASSID; 2777 PetscFunctionList PetscSectionSymList = NULL; 2778 2779 /*@ 2780 PetscSectionSymCreate - Creates an empty PetscSectionSym object. 2781 2782 Collective 2783 2784 Input Parameter: 2785 . comm - the MPI communicator 2786 2787 Output Parameter: 2788 . sym - pointer to the new set of symmetries 2789 2790 Level: developer 2791 2792 .seealso: PetscSectionSym, PetscSectionSymDestroy() 2793 @*/ 2794 PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym) 2795 { 2796 PetscErrorCode ierr; 2797 2798 PetscFunctionBegin; 2799 PetscValidPointer(sym,2); 2800 ierr = ISInitializePackage();CHKERRQ(ierr); 2801 ierr = PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);CHKERRQ(ierr); 2802 PetscFunctionReturn(0); 2803 } 2804 2805 /*@C 2806 PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation. 2807 2808 Collective on PetscSectionSym 2809 2810 Input Parameters: 2811 + sym - The section symmetry object 2812 - method - The name of the section symmetry type 2813 2814 Level: developer 2815 2816 .seealso: PetscSectionSymGetType(), PetscSectionSymCreate() 2817 @*/ 2818 PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method) 2819 { 2820 PetscErrorCode (*r)(PetscSectionSym); 2821 PetscBool match; 2822 PetscErrorCode ierr; 2823 2824 PetscFunctionBegin; 2825 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID,1); 2826 ierr = PetscObjectTypeCompare((PetscObject) sym, method, &match);CHKERRQ(ierr); 2827 if (match) PetscFunctionReturn(0); 2828 2829 ierr = PetscFunctionListFind(PetscSectionSymList,method,&r);CHKERRQ(ierr); 2830 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method); 2831 if (sym->ops->destroy) { 2832 ierr = (*sym->ops->destroy)(sym);CHKERRQ(ierr); 2833 sym->ops->destroy = NULL; 2834 } 2835 ierr = (*r)(sym);CHKERRQ(ierr); 2836 ierr = PetscObjectChangeTypeName((PetscObject)sym,method);CHKERRQ(ierr); 2837 PetscFunctionReturn(0); 2838 } 2839 2840 2841 /*@C 2842 PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym. 2843 2844 Not Collective 2845 2846 Input Parameter: 2847 . sym - The section symmetry 2848 2849 Output Parameter: 2850 . type - The index set type name 2851 2852 Level: developer 2853 2854 .seealso: PetscSectionSymSetType(), PetscSectionSymCreate() 2855 @*/ 2856 PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type) 2857 { 2858 PetscFunctionBegin; 2859 PetscValidHeaderSpecific(sym, PETSC_SECTION_SYM_CLASSID,1); 2860 PetscValidCharPointer(type,2); 2861 *type = ((PetscObject)sym)->type_name; 2862 PetscFunctionReturn(0); 2863 } 2864 2865 /*@C 2866 PetscSectionSymRegister - Adds a new section symmetry implementation 2867 2868 Not Collective 2869 2870 Input Parameters: 2871 + name - The name of a new user-defined creation routine 2872 - create_func - The creation routine itself 2873 2874 Notes: 2875 PetscSectionSymRegister() may be called multiple times to add several user-defined vectors 2876 2877 Level: developer 2878 2879 .seealso: PetscSectionSymCreate(), PetscSectionSymSetType() 2880 @*/ 2881 PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym)) 2882 { 2883 PetscErrorCode ierr; 2884 2885 PetscFunctionBegin; 2886 ierr = ISInitializePackage();CHKERRQ(ierr); 2887 ierr = PetscFunctionListAdd(&PetscSectionSymList,sname,function);CHKERRQ(ierr); 2888 PetscFunctionReturn(0); 2889 } 2890 2891 /*@ 2892 PetscSectionSymDestroy - Destroys a section symmetry. 2893 2894 Collective on PetscSectionSym 2895 2896 Input Parameters: 2897 . sym - the section symmetry 2898 2899 Level: developer 2900 2901 .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy() 2902 @*/ 2903 PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym) 2904 { 2905 SymWorkLink link,next; 2906 PetscErrorCode ierr; 2907 2908 PetscFunctionBegin; 2909 if (!*sym) PetscFunctionReturn(0); 2910 PetscValidHeaderSpecific((*sym),PETSC_SECTION_SYM_CLASSID,1); 2911 if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; PetscFunctionReturn(0);} 2912 if ((*sym)->ops->destroy) { 2913 ierr = (*(*sym)->ops->destroy)(*sym);CHKERRQ(ierr); 2914 } 2915 if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out"); 2916 for (link=(*sym)->workin; link; link=next) { 2917 next = link->next; 2918 ierr = PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);CHKERRQ(ierr); 2919 ierr = PetscFree(link);CHKERRQ(ierr); 2920 } 2921 (*sym)->workin = NULL; 2922 ierr = PetscHeaderDestroy(sym);CHKERRQ(ierr); 2923 PetscFunctionReturn(0); 2924 } 2925 2926 /*@C 2927 PetscSectionSymView - Displays a section symmetry 2928 2929 Collective on PetscSectionSym 2930 2931 Input Parameters: 2932 + sym - the index set 2933 - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF. 2934 2935 Level: developer 2936 2937 .seealso: PetscViewerASCIIOpen() 2938 @*/ 2939 PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer) 2940 { 2941 PetscErrorCode ierr; 2942 2943 PetscFunctionBegin; 2944 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,1); 2945 if (!viewer) { 2946 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);CHKERRQ(ierr); 2947 } 2948 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 2949 PetscCheckSameComm(sym,1,viewer,2); 2950 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);CHKERRQ(ierr); 2951 if (sym->ops->view) { 2952 ierr = (*sym->ops->view)(sym,viewer);CHKERRQ(ierr); 2953 } 2954 PetscFunctionReturn(0); 2955 } 2956 2957 /*@ 2958 PetscSectionSetSym - Set the symmetries for the data referred to by the section 2959 2960 Collective on PetscSection 2961 2962 Input Parameters: 2963 + section - the section describing data layout 2964 - sym - the symmetry describing the affect of orientation on the access of the data 2965 2966 Level: developer 2967 2968 .seealso: PetscSectionGetSym(), PetscSectionSymCreate() 2969 @*/ 2970 PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym) 2971 { 2972 PetscErrorCode ierr; 2973 2974 PetscFunctionBegin; 2975 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 2976 ierr = PetscSectionSymDestroy(&(section->sym));CHKERRQ(ierr); 2977 if (sym) { 2978 PetscValidHeaderSpecific(sym,PETSC_SECTION_SYM_CLASSID,2); 2979 PetscCheckSameComm(section,1,sym,2); 2980 ierr = PetscObjectReference((PetscObject) sym);CHKERRQ(ierr); 2981 } 2982 section->sym = sym; 2983 PetscFunctionReturn(0); 2984 } 2985 2986 /*@ 2987 PetscSectionGetSym - Get the symmetries for the data referred to by the section 2988 2989 Not collective 2990 2991 Input Parameters: 2992 . section - the section describing data layout 2993 2994 Output Parameters: 2995 . sym - the symmetry describing the affect of orientation on the access of the data 2996 2997 Level: developer 2998 2999 .seealso: PetscSectionSetSym(), PetscSectionSymCreate() 3000 @*/ 3001 PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym) 3002 { 3003 PetscFunctionBegin; 3004 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3005 *sym = section->sym; 3006 PetscFunctionReturn(0); 3007 } 3008 3009 /*@ 3010 PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section 3011 3012 Collective on PetscSection 3013 3014 Input Parameters: 3015 + section - the section describing data layout 3016 . field - the field number 3017 - sym - the symmetry describing the affect of orientation on the access of the data 3018 3019 Level: developer 3020 3021 .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate() 3022 @*/ 3023 PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym) 3024 { 3025 PetscErrorCode ierr; 3026 3027 PetscFunctionBegin; 3028 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3029 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); 3030 ierr = PetscSectionSetSym(section->field[field],sym);CHKERRQ(ierr); 3031 PetscFunctionReturn(0); 3032 } 3033 3034 /*@ 3035 PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section 3036 3037 Collective on PetscSection 3038 3039 Input Parameters: 3040 + section - the section describing data layout 3041 - field - the field number 3042 3043 Output Parameters: 3044 . sym - the symmetry describing the affect of orientation on the access of the data 3045 3046 Level: developer 3047 3048 .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate() 3049 @*/ 3050 PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym) 3051 { 3052 PetscFunctionBegin; 3053 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3054 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); 3055 *sym = section->field[field]->sym; 3056 PetscFunctionReturn(0); 3057 } 3058 3059 /*@C 3060 PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations. 3061 3062 Not collective 3063 3064 Input Parameters: 3065 + section - the section 3066 . numPoints - the number of points 3067 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3068 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3069 context, see DMPlexGetConeOrientation()). 3070 3071 Output Parameter: 3072 + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity). 3073 - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all 3074 identity). 3075 3076 Example of usage, gathering dofs into a local array (lArray) from a section array (sArray): 3077 .vb 3078 const PetscInt **perms; 3079 const PetscScalar **rots; 3080 PetscInt lOffset; 3081 3082 PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots); 3083 for (i = 0, lOffset = 0; i < numPoints; i++) { 3084 PetscInt point = points[2*i], dof, sOffset; 3085 const PetscInt *perm = perms ? perms[i] : NULL; 3086 const PetscScalar *rot = rots ? rots[i] : NULL; 3087 3088 PetscSectionGetDof(section,point,&dof); 3089 PetscSectionGetOffset(section,point,&sOffset); 3090 3091 if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}} 3092 else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}} 3093 if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }} 3094 lOffset += dof; 3095 } 3096 PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots); 3097 .ve 3098 3099 Example of usage, adding dofs into a section array (sArray) from a local array (lArray): 3100 .vb 3101 const PetscInt **perms; 3102 const PetscScalar **rots; 3103 PetscInt lOffset; 3104 3105 PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots); 3106 for (i = 0, lOffset = 0; i < numPoints; i++) { 3107 PetscInt point = points[2*i], dof, sOffset; 3108 const PetscInt *perm = perms ? perms[i] : NULL; 3109 const PetscScalar *rot = rots ? rots[i] : NULL; 3110 3111 PetscSectionGetDof(section,point,&dof); 3112 PetscSectionGetOffset(section,point,&sOff); 3113 3114 if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}} 3115 else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}} 3116 offset += dof; 3117 } 3118 PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots); 3119 .ve 3120 3121 Level: developer 3122 3123 .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym() 3124 @*/ 3125 PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3126 { 3127 PetscSectionSym sym; 3128 PetscErrorCode ierr; 3129 3130 PetscFunctionBegin; 3131 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3132 if (numPoints) PetscValidIntPointer(points,3); 3133 if (perms) *perms = NULL; 3134 if (rots) *rots = NULL; 3135 sym = section->sym; 3136 if (sym && (perms || rots)) { 3137 SymWorkLink link; 3138 3139 if (sym->workin) { 3140 link = sym->workin; 3141 sym->workin = sym->workin->next; 3142 } else { 3143 ierr = PetscNewLog(sym,&link);CHKERRQ(ierr); 3144 } 3145 if (numPoints > link->numPoints) { 3146 ierr = PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);CHKERRQ(ierr); 3147 ierr = PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);CHKERRQ(ierr); 3148 link->numPoints = numPoints; 3149 } 3150 link->next = sym->workout; 3151 sym->workout = link; 3152 ierr = PetscArrayzero((PetscInt**)link->perms,numPoints);CHKERRQ(ierr); 3153 ierr = PetscArrayzero((PetscInt**)link->rots,numPoints);CHKERRQ(ierr); 3154 ierr = (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);CHKERRQ(ierr); 3155 if (perms) *perms = link->perms; 3156 if (rots) *rots = link->rots; 3157 } 3158 PetscFunctionReturn(0); 3159 } 3160 3161 /*@C 3162 PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms() 3163 3164 Not collective 3165 3166 Input Parameters: 3167 + section - the section 3168 . numPoints - the number of points 3169 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3170 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3171 context, see DMPlexGetConeOrientation()). 3172 3173 Output Parameter: 3174 + perms - The permutations for the given orientations: set to NULL at conclusion 3175 - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion 3176 3177 Level: developer 3178 3179 .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym() 3180 @*/ 3181 PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3182 { 3183 PetscSectionSym sym; 3184 3185 PetscFunctionBegin; 3186 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3187 sym = section->sym; 3188 if (sym && (perms || rots)) { 3189 SymWorkLink *p,link; 3190 3191 for (p=&sym->workout; (link=*p); p=&link->next) { 3192 if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) { 3193 *p = link->next; 3194 link->next = sym->workin; 3195 sym->workin = link; 3196 if (perms) *perms = NULL; 3197 if (rots) *rots = NULL; 3198 PetscFunctionReturn(0); 3199 } 3200 } 3201 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out"); 3202 } 3203 PetscFunctionReturn(0); 3204 } 3205 3206 /*@C 3207 PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations. 3208 3209 Not collective 3210 3211 Input Parameters: 3212 + section - the section 3213 . field - the field of the section 3214 . numPoints - the number of points 3215 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3216 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3217 context, see DMPlexGetConeOrientation()). 3218 3219 Output Parameter: 3220 + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity). 3221 - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all 3222 identity). 3223 3224 Level: developer 3225 3226 .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms() 3227 @*/ 3228 PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3229 { 3230 PetscErrorCode ierr; 3231 3232 PetscFunctionBegin; 3233 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3234 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); 3235 ierr = PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);CHKERRQ(ierr); 3236 PetscFunctionReturn(0); 3237 } 3238 3239 /*@C 3240 PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms() 3241 3242 Not collective 3243 3244 Input Parameters: 3245 + section - the section 3246 . field - the field number 3247 . numPoints - the number of points 3248 - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an 3249 arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that 3250 context, see DMPlexGetConeOrientation()). 3251 3252 Output Parameter: 3253 + perms - The permutations for the given orientations: set to NULL at conclusion 3254 - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion 3255 3256 Level: developer 3257 3258 .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym() 3259 @*/ 3260 PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots) 3261 { 3262 PetscErrorCode ierr; 3263 3264 PetscFunctionBegin; 3265 PetscValidHeaderSpecific(section,PETSC_SECTION_CLASSID,1); 3266 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); 3267 ierr = PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);CHKERRQ(ierr); 3268 PetscFunctionReturn(0); 3269 } 3270 3271 /*@ 3272 PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset 3273 3274 Not collective 3275 3276 Input Parameter: 3277 . s - the global PetscSection 3278 3279 Output Parameters: 3280 . flg - the flag 3281 3282 Level: developer 3283 3284 .seealso: PetscSectionSetChart(), PetscSectionCreate() 3285 @*/ 3286 PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg) 3287 { 3288 PetscFunctionBegin; 3289 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 3290 *flg = s->useFieldOff; 3291 PetscFunctionReturn(0); 3292 } 3293 3294 /*@ 3295 PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset 3296 3297 Not collective 3298 3299 Input Parameters: 3300 + s - the global PetscSection 3301 - flg - the flag 3302 3303 Level: developer 3304 3305 .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate() 3306 @*/ 3307 PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg) 3308 { 3309 PetscFunctionBegin; 3310 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1); 3311 s->useFieldOff = flg; 3312 PetscFunctionReturn(0); 3313 } 3314 3315 #define PetscSectionExpandPoints_Loop(TYPE) \ 3316 { \ 3317 PetscInt i, n, o0, o1, size; \ 3318 TYPE *a0 = (TYPE*)origArray, *a1; \ 3319 ierr = PetscSectionGetStorageSize(s, &size);CHKERRQ(ierr); \ 3320 ierr = PetscMalloc1(size, &a1);CHKERRQ(ierr); \ 3321 for (i=0; i<npoints; i++) { \ 3322 ierr = PetscSectionGetOffset(origSection, points_[i], &o0);CHKERRQ(ierr); \ 3323 ierr = PetscSectionGetOffset(s, i, &o1);CHKERRQ(ierr); \ 3324 ierr = PetscSectionGetDof(s, i, &n);CHKERRQ(ierr); \ 3325 ierr = PetscMemcpy(&a1[o1], &a0[o0], n*unitsize);CHKERRQ(ierr); \ 3326 } \ 3327 *newArray = (void*)a1; \ 3328 } 3329 3330 /*@ 3331 PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points. 3332 3333 Not collective 3334 3335 Input Parameters: 3336 + origSection - the PetscSection describing the layout of the array 3337 . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL) 3338 . origArray - the array; its size must be equal to the storage size of origSection 3339 - points - IS with points to extract; its indices must lie in the chart of origSection 3340 3341 Output Parameters: 3342 + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs) 3343 - newArray - the array of the extracted DOFs; its size is the storage size of newSection 3344 3345 Level: developer 3346 3347 .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate() 3348 @*/ 3349 PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[]) 3350 { 3351 PetscSection s; 3352 const PetscInt *points_; 3353 PetscInt i, n, npoints, pStart, pEnd; 3354 PetscMPIInt unitsize; 3355 PetscErrorCode ierr; 3356 3357 PetscFunctionBegin; 3358 PetscValidHeaderSpecific(origSection, PETSC_SECTION_CLASSID, 1); 3359 PetscValidPointer(origArray, 3); 3360 PetscValidHeaderSpecific(points, IS_CLASSID, 4); 3361 if (newSection) PetscValidPointer(newSection, 5); 3362 if (newArray) PetscValidPointer(newArray, 6); 3363 ierr = MPI_Type_size(dataType, &unitsize);CHKERRQ(ierr); 3364 ierr = ISGetLocalSize(points, &npoints);CHKERRQ(ierr); 3365 ierr = ISGetIndices(points, &points_);CHKERRQ(ierr); 3366 ierr = PetscSectionGetChart(origSection, &pStart, &pEnd);CHKERRQ(ierr); 3367 ierr = PetscSectionCreate(PETSC_COMM_SELF, &s);CHKERRQ(ierr); 3368 ierr = PetscSectionSetChart(s, 0, npoints);CHKERRQ(ierr); 3369 for (i=0; i<npoints; i++) { 3370 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); 3371 ierr = PetscSectionGetDof(origSection, points_[i], &n);CHKERRQ(ierr); 3372 ierr = PetscSectionSetDof(s, i, n);CHKERRQ(ierr); 3373 } 3374 ierr = PetscSectionSetUp(s);CHKERRQ(ierr); 3375 if (newArray) { 3376 if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);} 3377 else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);} 3378 else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);} 3379 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype"); 3380 } 3381 if (newSection) { 3382 *newSection = s; 3383 } else { 3384 ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); 3385 } 3386 PetscFunctionReturn(0); 3387 } 3388