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