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