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