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