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