1 /* Basis Jet Tabulation 2 3 We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We 4 follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can 5 be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis 6 as a prime basis. 7 8 \psi_i = \sum_k \alpha_{ki} \phi_k 9 10 Our nodal basis is defined in terms of the dual basis $n_j$ 11 12 n_j \cdot \psi_i = \delta_{ji} 13 14 and we may act on the first equation to obtain 15 16 n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k 17 \delta_{ji} = \sum_k \alpha_{ki} V_{jk} 18 I = V \alpha 19 20 so the coefficients of the nodal basis in the prime basis are 21 22 \alpha = V^{-1} 23 24 We will define the dual basis vectors $n_j$ using a quadrature rule. 25 26 Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials 27 (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can 28 be implemented exactly as in FIAT using functionals $L_j$. 29 30 I will have to count the degrees correctly for the Legendre product when we are on simplices. 31 32 We will have three objects: 33 - Space, P: this just need point evaluation I think 34 - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q 35 - FEM: This keeps {P, P', Q} 36 */ 37 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 38 #include <petscdmplex.h> 39 40 PetscBool FEcite = PETSC_FALSE; 41 const char FECitation[] = "@article{kirby2004,\n" 42 " title = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n" 43 " journal = {ACM Transactions on Mathematical Software},\n" 44 " author = {Robert C. Kirby},\n" 45 " volume = {30},\n" 46 " number = {4},\n" 47 " pages = {502--516},\n" 48 " doi = {10.1145/1039813.1039820},\n" 49 " year = {2004}\n}\n"; 50 51 PetscClassId PETSCFE_CLASSID = 0; 52 53 PetscFunctionList PetscFEList = NULL; 54 PetscBool PetscFERegisterAllCalled = PETSC_FALSE; 55 56 /*@C 57 PetscFERegister - Adds a new PetscFE implementation 58 59 Not Collective 60 61 Input Parameters: 62 + name - The name of a new user-defined creation routine 63 - create_func - The creation routine itself 64 65 Notes: 66 PetscFERegister() may be called multiple times to add several user-defined PetscFEs 67 68 Sample usage: 69 .vb 70 PetscFERegister("my_fe", MyPetscFECreate); 71 .ve 72 73 Then, your PetscFE type can be chosen with the procedural interface via 74 .vb 75 PetscFECreate(MPI_Comm, PetscFE *); 76 PetscFESetType(PetscFE, "my_fe"); 77 .ve 78 or at runtime via the option 79 .vb 80 -petscfe_type my_fe 81 .ve 82 83 Level: advanced 84 85 .seealso: PetscFERegisterAll(), PetscFERegisterDestroy() 86 87 @*/ 88 PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE)) 89 { 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 ierr = PetscFunctionListAdd(&PetscFEList, sname, function);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 97 /*@C 98 PetscFESetType - Builds a particular PetscFE 99 100 Collective on fem 101 102 Input Parameters: 103 + fem - The PetscFE object 104 - name - The kind of FEM space 105 106 Options Database Key: 107 . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types 108 109 Level: intermediate 110 111 .seealso: PetscFEGetType(), PetscFECreate() 112 @*/ 113 PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name) 114 { 115 PetscErrorCode (*r)(PetscFE); 116 PetscBool match; 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 121 ierr = PetscObjectTypeCompare((PetscObject) fem, name, &match);CHKERRQ(ierr); 122 if (match) PetscFunctionReturn(0); 123 124 if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);} 125 ierr = PetscFunctionListFind(PetscFEList, name, &r);CHKERRQ(ierr); 126 if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name); 127 128 if (fem->ops->destroy) { 129 ierr = (*fem->ops->destroy)(fem);CHKERRQ(ierr); 130 fem->ops->destroy = NULL; 131 } 132 ierr = (*r)(fem);CHKERRQ(ierr); 133 ierr = PetscObjectChangeTypeName((PetscObject) fem, name);CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 /*@C 138 PetscFEGetType - Gets the PetscFE type name (as a string) from the object. 139 140 Not Collective 141 142 Input Parameter: 143 . fem - The PetscFE 144 145 Output Parameter: 146 . name - The PetscFE type name 147 148 Level: intermediate 149 150 .seealso: PetscFESetType(), PetscFECreate() 151 @*/ 152 PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name) 153 { 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 158 PetscValidPointer(name, 2); 159 if (!PetscFERegisterAllCalled) { 160 ierr = PetscFERegisterAll();CHKERRQ(ierr); 161 } 162 *name = ((PetscObject) fem)->type_name; 163 PetscFunctionReturn(0); 164 } 165 166 /*@C 167 PetscFEViewFromOptions - View from Options 168 169 Collective on PetscFE 170 171 Input Parameters: 172 + A - the PetscFE object 173 . obj - Optional object 174 - name - command line option 175 176 Level: intermediate 177 .seealso: PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate() 178 @*/ 179 PetscErrorCode PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[]) 180 { 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(A,PETSCFE_CLASSID,1); 185 ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 /*@C 190 PetscFEView - Views a PetscFE 191 192 Collective on fem 193 194 Input Parameter: 195 + fem - the PetscFE object to view 196 - viewer - the viewer 197 198 Level: beginner 199 200 .seealso PetscFEDestroy() 201 @*/ 202 PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer) 203 { 204 PetscBool iascii; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 209 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 210 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);CHKERRQ(ierr);} 211 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);CHKERRQ(ierr); 212 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 213 if (fem->ops->view) {ierr = (*fem->ops->view)(fem, viewer);CHKERRQ(ierr);} 214 PetscFunctionReturn(0); 215 } 216 217 /*@ 218 PetscFESetFromOptions - sets parameters in a PetscFE from the options database 219 220 Collective on fem 221 222 Input Parameter: 223 . fem - the PetscFE object to set options for 224 225 Options Database: 226 + -petscfe_num_blocks - the number of cell blocks to integrate concurrently 227 - -petscfe_num_batches - the number of cell batches to integrate serially 228 229 Level: intermediate 230 231 .seealso PetscFEView() 232 @*/ 233 PetscErrorCode PetscFESetFromOptions(PetscFE fem) 234 { 235 const char *defaultType; 236 char name[256]; 237 PetscBool flg; 238 PetscErrorCode ierr; 239 240 PetscFunctionBegin; 241 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 242 if (!((PetscObject) fem)->type_name) { 243 defaultType = PETSCFEBASIC; 244 } else { 245 defaultType = ((PetscObject) fem)->type_name; 246 } 247 if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);} 248 249 ierr = PetscObjectOptionsBegin((PetscObject) fem);CHKERRQ(ierr); 250 ierr = PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);CHKERRQ(ierr); 251 if (flg) { 252 ierr = PetscFESetType(fem, name);CHKERRQ(ierr); 253 } else if (!((PetscObject) fem)->type_name) { 254 ierr = PetscFESetType(fem, defaultType);CHKERRQ(ierr); 255 } 256 ierr = PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);CHKERRQ(ierr); 257 ierr = PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);CHKERRQ(ierr); 258 if (fem->ops->setfromoptions) { 259 ierr = (*fem->ops->setfromoptions)(PetscOptionsObject,fem);CHKERRQ(ierr); 260 } 261 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 262 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);CHKERRQ(ierr); 263 ierr = PetscOptionsEnd();CHKERRQ(ierr); 264 ierr = PetscFEViewFromOptions(fem, NULL, "-petscfe_view");CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 268 /*@C 269 PetscFESetUp - Construct data structures for the PetscFE 270 271 Collective on fem 272 273 Input Parameter: 274 . fem - the PetscFE object to setup 275 276 Level: intermediate 277 278 .seealso PetscFEView(), PetscFEDestroy() 279 @*/ 280 PetscErrorCode PetscFESetUp(PetscFE fem) 281 { 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 286 if (fem->setupcalled) PetscFunctionReturn(0); 287 fem->setupcalled = PETSC_TRUE; 288 if (fem->ops->setup) {ierr = (*fem->ops->setup)(fem);CHKERRQ(ierr);} 289 PetscFunctionReturn(0); 290 } 291 292 /*@ 293 PetscFEDestroy - Destroys a PetscFE object 294 295 Collective on fem 296 297 Input Parameter: 298 . fem - the PetscFE object to destroy 299 300 Level: beginner 301 302 .seealso PetscFEView() 303 @*/ 304 PetscErrorCode PetscFEDestroy(PetscFE *fem) 305 { 306 PetscErrorCode ierr; 307 308 PetscFunctionBegin; 309 if (!*fem) PetscFunctionReturn(0); 310 PetscValidHeaderSpecific((*fem), PETSCFE_CLASSID, 1); 311 312 if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; PetscFunctionReturn(0);} 313 ((PetscObject) (*fem))->refct = 0; 314 315 if ((*fem)->subspaces) { 316 PetscInt dim, d; 317 318 ierr = PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);CHKERRQ(ierr); 319 for (d = 0; d < dim; ++d) {ierr = PetscFEDestroy(&(*fem)->subspaces[d]);CHKERRQ(ierr);} 320 } 321 ierr = PetscFree((*fem)->subspaces);CHKERRQ(ierr); 322 ierr = PetscFree((*fem)->invV);CHKERRQ(ierr); 323 ierr = PetscTabulationDestroy(&(*fem)->T);CHKERRQ(ierr); 324 ierr = PetscTabulationDestroy(&(*fem)->Tf);CHKERRQ(ierr); 325 ierr = PetscTabulationDestroy(&(*fem)->Tc);CHKERRQ(ierr); 326 ierr = PetscSpaceDestroy(&(*fem)->basisSpace);CHKERRQ(ierr); 327 ierr = PetscDualSpaceDestroy(&(*fem)->dualSpace);CHKERRQ(ierr); 328 ierr = PetscQuadratureDestroy(&(*fem)->quadrature);CHKERRQ(ierr); 329 ierr = PetscQuadratureDestroy(&(*fem)->faceQuadrature);CHKERRQ(ierr); 330 331 if ((*fem)->ops->destroy) {ierr = (*(*fem)->ops->destroy)(*fem);CHKERRQ(ierr);} 332 ierr = PetscHeaderDestroy(fem);CHKERRQ(ierr); 333 PetscFunctionReturn(0); 334 } 335 336 /*@ 337 PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType(). 338 339 Collective 340 341 Input Parameter: 342 . comm - The communicator for the PetscFE object 343 344 Output Parameter: 345 . fem - The PetscFE object 346 347 Level: beginner 348 349 .seealso: PetscFESetType(), PETSCFEGALERKIN 350 @*/ 351 PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem) 352 { 353 PetscFE f; 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 PetscValidPointer(fem, 2); 358 ierr = PetscCitationsRegister(FECitation,&FEcite);CHKERRQ(ierr); 359 *fem = NULL; 360 ierr = PetscFEInitializePackage();CHKERRQ(ierr); 361 362 ierr = PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);CHKERRQ(ierr); 363 364 f->basisSpace = NULL; 365 f->dualSpace = NULL; 366 f->numComponents = 1; 367 f->subspaces = NULL; 368 f->invV = NULL; 369 f->T = NULL; 370 f->Tf = NULL; 371 f->Tc = NULL; 372 ierr = PetscArrayzero(&f->quadrature, 1);CHKERRQ(ierr); 373 ierr = PetscArrayzero(&f->faceQuadrature, 1);CHKERRQ(ierr); 374 f->blockSize = 0; 375 f->numBlocks = 1; 376 f->batchSize = 0; 377 f->numBatches = 1; 378 379 *fem = f; 380 PetscFunctionReturn(0); 381 } 382 383 /*@ 384 PetscFEGetSpatialDimension - Returns the spatial dimension of the element 385 386 Not collective 387 388 Input Parameter: 389 . fem - The PetscFE object 390 391 Output Parameter: 392 . dim - The spatial dimension 393 394 Level: intermediate 395 396 .seealso: PetscFECreate() 397 @*/ 398 PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim) 399 { 400 DM dm; 401 PetscErrorCode ierr; 402 403 PetscFunctionBegin; 404 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 405 PetscValidPointer(dim, 2); 406 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 407 ierr = DMGetDimension(dm, dim);CHKERRQ(ierr); 408 PetscFunctionReturn(0); 409 } 410 411 /*@ 412 PetscFESetNumComponents - Sets the number of components in the element 413 414 Not collective 415 416 Input Parameters: 417 + fem - The PetscFE object 418 - comp - The number of field components 419 420 Level: intermediate 421 422 .seealso: PetscFECreate() 423 @*/ 424 PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp) 425 { 426 PetscFunctionBegin; 427 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 428 fem->numComponents = comp; 429 PetscFunctionReturn(0); 430 } 431 432 /*@ 433 PetscFEGetNumComponents - Returns the number of components in the element 434 435 Not collective 436 437 Input Parameter: 438 . fem - The PetscFE object 439 440 Output Parameter: 441 . comp - The number of field components 442 443 Level: intermediate 444 445 .seealso: PetscFECreate() 446 @*/ 447 PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp) 448 { 449 PetscFunctionBegin; 450 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 451 PetscValidPointer(comp, 2); 452 *comp = fem->numComponents; 453 PetscFunctionReturn(0); 454 } 455 456 /*@ 457 PetscFESetTileSizes - Sets the tile sizes for evaluation 458 459 Not collective 460 461 Input Parameters: 462 + fem - The PetscFE object 463 . blockSize - The number of elements in a block 464 . numBlocks - The number of blocks in a batch 465 . batchSize - The number of elements in a batch 466 - numBatches - The number of batches in a chunk 467 468 Level: intermediate 469 470 .seealso: PetscFECreate() 471 @*/ 472 PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches) 473 { 474 PetscFunctionBegin; 475 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 476 fem->blockSize = blockSize; 477 fem->numBlocks = numBlocks; 478 fem->batchSize = batchSize; 479 fem->numBatches = numBatches; 480 PetscFunctionReturn(0); 481 } 482 483 /*@ 484 PetscFEGetTileSizes - Returns the tile sizes for evaluation 485 486 Not collective 487 488 Input Parameter: 489 . fem - The PetscFE object 490 491 Output Parameters: 492 + blockSize - The number of elements in a block 493 . numBlocks - The number of blocks in a batch 494 . batchSize - The number of elements in a batch 495 - numBatches - The number of batches in a chunk 496 497 Level: intermediate 498 499 .seealso: PetscFECreate() 500 @*/ 501 PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches) 502 { 503 PetscFunctionBegin; 504 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 505 if (blockSize) PetscValidPointer(blockSize, 2); 506 if (numBlocks) PetscValidPointer(numBlocks, 3); 507 if (batchSize) PetscValidPointer(batchSize, 4); 508 if (numBatches) PetscValidPointer(numBatches, 5); 509 if (blockSize) *blockSize = fem->blockSize; 510 if (numBlocks) *numBlocks = fem->numBlocks; 511 if (batchSize) *batchSize = fem->batchSize; 512 if (numBatches) *numBatches = fem->numBatches; 513 PetscFunctionReturn(0); 514 } 515 516 /*@ 517 PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution 518 519 Not collective 520 521 Input Parameter: 522 . fem - The PetscFE object 523 524 Output Parameter: 525 . sp - The PetscSpace object 526 527 Level: intermediate 528 529 .seealso: PetscFECreate() 530 @*/ 531 PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp) 532 { 533 PetscFunctionBegin; 534 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 535 PetscValidPointer(sp, 2); 536 *sp = fem->basisSpace; 537 PetscFunctionReturn(0); 538 } 539 540 /*@ 541 PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution 542 543 Not collective 544 545 Input Parameters: 546 + fem - The PetscFE object 547 - sp - The PetscSpace object 548 549 Level: intermediate 550 551 .seealso: PetscFECreate() 552 @*/ 553 PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp) 554 { 555 PetscErrorCode ierr; 556 557 PetscFunctionBegin; 558 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 559 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 2); 560 ierr = PetscSpaceDestroy(&fem->basisSpace);CHKERRQ(ierr); 561 fem->basisSpace = sp; 562 ierr = PetscObjectReference((PetscObject) fem->basisSpace);CHKERRQ(ierr); 563 PetscFunctionReturn(0); 564 } 565 566 /*@ 567 PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product 568 569 Not collective 570 571 Input Parameter: 572 . fem - The PetscFE object 573 574 Output Parameter: 575 . sp - The PetscDualSpace object 576 577 Level: intermediate 578 579 .seealso: PetscFECreate() 580 @*/ 581 PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp) 582 { 583 PetscFunctionBegin; 584 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 585 PetscValidPointer(sp, 2); 586 *sp = fem->dualSpace; 587 PetscFunctionReturn(0); 588 } 589 590 /*@ 591 PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product 592 593 Not collective 594 595 Input Parameters: 596 + fem - The PetscFE object 597 - sp - The PetscDualSpace object 598 599 Level: intermediate 600 601 .seealso: PetscFECreate() 602 @*/ 603 PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp) 604 { 605 PetscErrorCode ierr; 606 607 PetscFunctionBegin; 608 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 609 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 2); 610 ierr = PetscDualSpaceDestroy(&fem->dualSpace);CHKERRQ(ierr); 611 fem->dualSpace = sp; 612 ierr = PetscObjectReference((PetscObject) fem->dualSpace);CHKERRQ(ierr); 613 PetscFunctionReturn(0); 614 } 615 616 /*@ 617 PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products 618 619 Not collective 620 621 Input Parameter: 622 . fem - The PetscFE object 623 624 Output Parameter: 625 . q - The PetscQuadrature object 626 627 Level: intermediate 628 629 .seealso: PetscFECreate() 630 @*/ 631 PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q) 632 { 633 PetscFunctionBegin; 634 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 635 PetscValidPointer(q, 2); 636 *q = fem->quadrature; 637 PetscFunctionReturn(0); 638 } 639 640 /*@ 641 PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products 642 643 Not collective 644 645 Input Parameters: 646 + fem - The PetscFE object 647 - q - The PetscQuadrature object 648 649 Level: intermediate 650 651 .seealso: PetscFECreate() 652 @*/ 653 PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q) 654 { 655 PetscInt Nc, qNc; 656 PetscErrorCode ierr; 657 658 PetscFunctionBegin; 659 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 660 if (q == fem->quadrature) PetscFunctionReturn(0); 661 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 662 ierr = PetscQuadratureGetNumComponents(q, &qNc);CHKERRQ(ierr); 663 if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc); 664 ierr = PetscTabulationDestroy(&fem->T);CHKERRQ(ierr); 665 ierr = PetscTabulationDestroy(&fem->Tc);CHKERRQ(ierr); 666 ierr = PetscObjectReference((PetscObject) q);CHKERRQ(ierr); 667 ierr = PetscQuadratureDestroy(&fem->quadrature);CHKERRQ(ierr); 668 fem->quadrature = q; 669 PetscFunctionReturn(0); 670 } 671 672 /*@ 673 PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces 674 675 Not collective 676 677 Input Parameter: 678 . fem - The PetscFE object 679 680 Output Parameter: 681 . q - The PetscQuadrature object 682 683 Level: intermediate 684 685 .seealso: PetscFECreate() 686 @*/ 687 PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q) 688 { 689 PetscFunctionBegin; 690 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 691 PetscValidPointer(q, 2); 692 *q = fem->faceQuadrature; 693 PetscFunctionReturn(0); 694 } 695 696 /*@ 697 PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces 698 699 Not collective 700 701 Input Parameters: 702 + fem - The PetscFE object 703 - q - The PetscQuadrature object 704 705 Level: intermediate 706 707 .seealso: PetscFECreate() 708 @*/ 709 PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q) 710 { 711 PetscInt Nc, qNc; 712 PetscErrorCode ierr; 713 714 PetscFunctionBegin; 715 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 716 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 717 ierr = PetscQuadratureGetNumComponents(q, &qNc);CHKERRQ(ierr); 718 if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc); 719 ierr = PetscTabulationDestroy(&fem->Tf);CHKERRQ(ierr); 720 ierr = PetscQuadratureDestroy(&fem->faceQuadrature);CHKERRQ(ierr); 721 fem->faceQuadrature = q; 722 ierr = PetscObjectReference((PetscObject) q);CHKERRQ(ierr); 723 PetscFunctionReturn(0); 724 } 725 726 /*@ 727 PetscFECopyQuadrature - Copy both volumetric and surface quadrature 728 729 Not collective 730 731 Input Parameters: 732 + sfe - The PetscFE source for the quadratures 733 - tfe - The PetscFE target for the quadratures 734 735 Level: intermediate 736 737 .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature() 738 @*/ 739 PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe) 740 { 741 PetscQuadrature q; 742 PetscErrorCode ierr; 743 744 PetscFunctionBegin; 745 PetscValidHeaderSpecific(sfe, PETSCFE_CLASSID, 1); 746 PetscValidHeaderSpecific(tfe, PETSCFE_CLASSID, 2); 747 ierr = PetscFEGetQuadrature(sfe, &q);CHKERRQ(ierr); 748 ierr = PetscFESetQuadrature(tfe, q);CHKERRQ(ierr); 749 ierr = PetscFEGetFaceQuadrature(sfe, &q);CHKERRQ(ierr); 750 ierr = PetscFESetFaceQuadrature(tfe, q);CHKERRQ(ierr); 751 PetscFunctionReturn(0); 752 } 753 754 /*@C 755 PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension 756 757 Not collective 758 759 Input Parameter: 760 . fem - The PetscFE object 761 762 Output Parameter: 763 . numDof - Array with the number of dofs per dimension 764 765 Level: intermediate 766 767 .seealso: PetscFECreate() 768 @*/ 769 PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof) 770 { 771 PetscErrorCode ierr; 772 773 PetscFunctionBegin; 774 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 775 PetscValidPointer(numDof, 2); 776 ierr = PetscDualSpaceGetNumDof(fem->dualSpace, numDof);CHKERRQ(ierr); 777 PetscFunctionReturn(0); 778 } 779 780 /*@C 781 PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell 782 783 Not collective 784 785 Input Parameter: 786 . fem - The PetscFE object 787 788 Output Parameter: 789 . T - The basis function values and derivatives at quadrature points 790 791 Note: 792 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c 793 $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d 794 $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e 795 796 Level: intermediate 797 798 .seealso: PetscFECreateTabulation(), PetscTabulationDestroy() 799 @*/ 800 PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscTabulation *T) 801 { 802 PetscInt npoints; 803 const PetscReal *points; 804 PetscErrorCode ierr; 805 806 PetscFunctionBegin; 807 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 808 PetscValidPointer(T, 2); 809 ierr = PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr); 810 if (!fem->T) {ierr = PetscFECreateTabulation(fem, 1, npoints, points, 1, &fem->T);CHKERRQ(ierr);} 811 *T = fem->T; 812 PetscFunctionReturn(0); 813 } 814 815 /*@C 816 PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell 817 818 Not collective 819 820 Input Parameter: 821 . fem - The PetscFE object 822 823 Output Parameters: 824 . Tf - The basis function values and derviatives at face quadrature points 825 826 Note: 827 $ T->T[0] = Bf[((f*Nq + q)*pdim + i)*Nc + c] is the value at point f,q for basis function i and component c 828 $ T->T[1] = Df[(((f*Nq + q)*pdim + i)*Nc + c)*dim + d] is the derivative value at point f,q for basis function i, component c, in direction d 829 $ T->T[2] = Hf[((((f*Nq + q)*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point f,q for basis function i, component c, in directions d and e 830 831 Level: intermediate 832 833 .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy() 834 @*/ 835 PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscTabulation *Tf) 836 { 837 PetscErrorCode ierr; 838 839 PetscFunctionBegin; 840 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 841 PetscValidPointer(Tf, 2); 842 if (!fem->Tf) { 843 const PetscReal xi0[3] = {-1., -1., -1.}; 844 PetscReal v0[3], J[9], detJ; 845 PetscQuadrature fq; 846 PetscDualSpace sp; 847 DM dm; 848 const PetscInt *faces; 849 PetscInt dim, numFaces, f, npoints, q; 850 const PetscReal *points; 851 PetscReal *facePoints; 852 853 ierr = PetscFEGetDualSpace(fem, &sp);CHKERRQ(ierr); 854 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 855 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 856 ierr = DMPlexGetConeSize(dm, 0, &numFaces);CHKERRQ(ierr); 857 ierr = DMPlexGetCone(dm, 0, &faces);CHKERRQ(ierr); 858 ierr = PetscFEGetFaceQuadrature(fem, &fq);CHKERRQ(ierr); 859 if (fq) { 860 ierr = PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr); 861 ierr = PetscMalloc1(numFaces*npoints*dim, &facePoints);CHKERRQ(ierr); 862 for (f = 0; f < numFaces; ++f) { 863 ierr = DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 864 for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]); 865 } 866 ierr = PetscFECreateTabulation(fem, numFaces, npoints, facePoints, 1, &fem->Tf);CHKERRQ(ierr); 867 ierr = PetscFree(facePoints);CHKERRQ(ierr); 868 } 869 } 870 *Tf = fem->Tf; 871 PetscFunctionReturn(0); 872 } 873 874 /*@C 875 PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points 876 877 Not collective 878 879 Input Parameter: 880 . fem - The PetscFE object 881 882 Output Parameters: 883 . Tc - The basis function values at face centroid points 884 885 Note: 886 $ T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c 887 888 Level: intermediate 889 890 .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy() 891 @*/ 892 PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc) 893 { 894 PetscErrorCode ierr; 895 896 PetscFunctionBegin; 897 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 898 PetscValidPointer(Tc, 2); 899 if (!fem->Tc) { 900 PetscDualSpace sp; 901 DM dm; 902 const PetscInt *cone; 903 PetscReal *centroids; 904 PetscInt dim, numFaces, f; 905 906 ierr = PetscFEGetDualSpace(fem, &sp);CHKERRQ(ierr); 907 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 908 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 909 ierr = DMPlexGetConeSize(dm, 0, &numFaces);CHKERRQ(ierr); 910 ierr = DMPlexGetCone(dm, 0, &cone);CHKERRQ(ierr); 911 ierr = PetscMalloc1(numFaces*dim, ¢roids);CHKERRQ(ierr); 912 for (f = 0; f < numFaces; ++f) {ierr = DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);CHKERRQ(ierr);} 913 ierr = PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);CHKERRQ(ierr); 914 ierr = PetscFree(centroids);CHKERRQ(ierr); 915 } 916 *Tc = fem->Tc; 917 PetscFunctionReturn(0); 918 } 919 920 /*@C 921 PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided. 922 923 Not collective 924 925 Input Parameters: 926 + fem - The PetscFE object 927 . nrepl - The number of replicas 928 . npoints - The number of tabulation points in a replica 929 . points - The tabulation point coordinates 930 - K - The number of derivatives calculated 931 932 Output Parameter: 933 . T - The basis function values and derivatives at tabulation points 934 935 Note: 936 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c 937 $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d 938 $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e 939 940 Level: intermediate 941 942 .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy() 943 @*/ 944 PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T) 945 { 946 DM dm; 947 PetscDualSpace Q; 948 PetscInt Nb; /* Dimension of FE space P */ 949 PetscInt Nc; /* Field components */ 950 PetscInt cdim; /* Reference coordinate dimension */ 951 PetscInt k; 952 PetscErrorCode ierr; 953 954 PetscFunctionBegin; 955 if (!npoints || !fem->dualSpace || K < 0) { 956 *T = NULL; 957 PetscFunctionReturn(0); 958 } 959 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 960 PetscValidPointer(points, 4); 961 PetscValidPointer(T, 6); 962 ierr = PetscFEGetDualSpace(fem, &Q);CHKERRQ(ierr); 963 ierr = PetscDualSpaceGetDM(Q, &dm);CHKERRQ(ierr); 964 ierr = DMGetDimension(dm, &cdim);CHKERRQ(ierr); 965 ierr = PetscDualSpaceGetDimension(Q, &Nb);CHKERRQ(ierr); 966 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 967 ierr = PetscMalloc1(1, T);CHKERRQ(ierr); 968 (*T)->K = !cdim ? 0 : K; 969 (*T)->Nr = nrepl; 970 (*T)->Np = npoints; 971 (*T)->Nb = Nb; 972 (*T)->Nc = Nc; 973 (*T)->cdim = cdim; 974 ierr = PetscMalloc1((*T)->K+1, &(*T)->T);CHKERRQ(ierr); 975 for (k = 0; k <= (*T)->K; ++k) { 976 ierr = PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);CHKERRQ(ierr); 977 } 978 ierr = (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);CHKERRQ(ierr); 979 PetscFunctionReturn(0); 980 } 981 982 /*@C 983 PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided. 984 985 Not collective 986 987 Input Parameters: 988 + fem - The PetscFE object 989 . npoints - The number of tabulation points 990 . points - The tabulation point coordinates 991 . K - The number of derivatives calculated 992 - T - An existing tabulation object with enough allocated space 993 994 Output Parameter: 995 . T - The basis function values and derivatives at tabulation points 996 997 Note: 998 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c 999 $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d 1000 $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e 1001 1002 Level: intermediate 1003 1004 .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy() 1005 @*/ 1006 PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T) 1007 { 1008 PetscErrorCode ierr; 1009 1010 PetscFunctionBeginHot; 1011 if (!npoints || !fem->dualSpace || K < 0) PetscFunctionReturn(0); 1012 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1013 PetscValidPointer(points, 3); 1014 PetscValidPointer(T, 5); 1015 if (PetscDefined(USE_DEBUG)) { 1016 DM dm; 1017 PetscDualSpace Q; 1018 PetscInt Nb; /* Dimension of FE space P */ 1019 PetscInt Nc; /* Field components */ 1020 PetscInt cdim; /* Reference coordinate dimension */ 1021 1022 ierr = PetscFEGetDualSpace(fem, &Q);CHKERRQ(ierr); 1023 ierr = PetscDualSpaceGetDM(Q, &dm);CHKERRQ(ierr); 1024 ierr = DMGetDimension(dm, &cdim);CHKERRQ(ierr); 1025 ierr = PetscDualSpaceGetDimension(Q, &Nb);CHKERRQ(ierr); 1026 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 1027 if (T->K != (!cdim ? 0 : K)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %D must match requested K %D", T->K, !cdim ? 0 : K); 1028 if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb); 1029 if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc); 1030 if (T->cdim != cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim); 1031 } 1032 T->Nr = 1; 1033 T->Np = npoints; 1034 ierr = (*fem->ops->createtabulation)(fem, npoints, points, K, T);CHKERRQ(ierr); 1035 PetscFunctionReturn(0); 1036 } 1037 1038 /*@C 1039 PetscTabulationDestroy - Frees memory from the associated tabulation. 1040 1041 Not collective 1042 1043 Input Parameter: 1044 . T - The tabulation 1045 1046 Level: intermediate 1047 1048 .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation() 1049 @*/ 1050 PetscErrorCode PetscTabulationDestroy(PetscTabulation *T) 1051 { 1052 PetscInt k; 1053 PetscErrorCode ierr; 1054 1055 PetscFunctionBegin; 1056 PetscValidPointer(T, 1); 1057 if (!T || !(*T)) PetscFunctionReturn(0); 1058 for (k = 0; k <= (*T)->K; ++k) {ierr = PetscFree((*T)->T[k]);CHKERRQ(ierr);} 1059 ierr = PetscFree((*T)->T);CHKERRQ(ierr); 1060 ierr = PetscFree(*T);CHKERRQ(ierr); 1061 *T = NULL; 1062 PetscFunctionReturn(0); 1063 } 1064 1065 PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE) 1066 { 1067 PetscSpace bsp, bsubsp; 1068 PetscDualSpace dsp, dsubsp; 1069 PetscInt dim, depth, numComp, i, j, coneSize, order; 1070 PetscFEType type; 1071 DM dm; 1072 DMLabel label; 1073 PetscReal *xi, *v, *J, detJ; 1074 const char *name; 1075 PetscQuadrature origin, fullQuad, subQuad; 1076 PetscErrorCode ierr; 1077 1078 PetscFunctionBegin; 1079 PetscValidHeaderSpecific(fe,PETSCFE_CLASSID,1); 1080 PetscValidPointer(trFE,3); 1081 ierr = PetscFEGetBasisSpace(fe,&bsp);CHKERRQ(ierr); 1082 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 1083 ierr = PetscDualSpaceGetDM(dsp,&dm);CHKERRQ(ierr); 1084 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 1085 ierr = DMPlexGetDepthLabel(dm,&label);CHKERRQ(ierr); 1086 ierr = DMLabelGetValue(label,refPoint,&depth);CHKERRQ(ierr); 1087 ierr = PetscCalloc1(depth,&xi);CHKERRQ(ierr); 1088 ierr = PetscMalloc1(dim,&v);CHKERRQ(ierr); 1089 ierr = PetscMalloc1(dim*dim,&J);CHKERRQ(ierr); 1090 for (i = 0; i < depth; i++) xi[i] = 0.; 1091 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,&origin);CHKERRQ(ierr); 1092 ierr = PetscQuadratureSetData(origin,depth,0,1,xi,NULL);CHKERRQ(ierr); 1093 ierr = DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);CHKERRQ(ierr); 1094 /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */ 1095 for (i = 1; i < dim; i++) { 1096 for (j = 0; j < depth; j++) { 1097 J[i * depth + j] = J[i * dim + j]; 1098 } 1099 } 1100 ierr = PetscQuadratureDestroy(&origin);CHKERRQ(ierr); 1101 ierr = PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);CHKERRQ(ierr); 1102 ierr = PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);CHKERRQ(ierr); 1103 ierr = PetscSpaceSetUp(bsubsp);CHKERRQ(ierr); 1104 ierr = PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);CHKERRQ(ierr); 1105 ierr = PetscFEGetType(fe,&type);CHKERRQ(ierr); 1106 ierr = PetscFESetType(*trFE,type);CHKERRQ(ierr); 1107 ierr = PetscFEGetNumComponents(fe,&numComp);CHKERRQ(ierr); 1108 ierr = PetscFESetNumComponents(*trFE,numComp);CHKERRQ(ierr); 1109 ierr = PetscFESetBasisSpace(*trFE,bsubsp);CHKERRQ(ierr); 1110 ierr = PetscFESetDualSpace(*trFE,dsubsp);CHKERRQ(ierr); 1111 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 1112 if (name) {ierr = PetscFESetName(*trFE, name);CHKERRQ(ierr);} 1113 ierr = PetscFEGetQuadrature(fe,&fullQuad);CHKERRQ(ierr); 1114 ierr = PetscQuadratureGetOrder(fullQuad,&order);CHKERRQ(ierr); 1115 ierr = DMPlexGetConeSize(dm,refPoint,&coneSize);CHKERRQ(ierr); 1116 if (coneSize == 2 * depth) { 1117 ierr = PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);CHKERRQ(ierr); 1118 } else { 1119 ierr = PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);CHKERRQ(ierr); 1120 } 1121 ierr = PetscFESetQuadrature(*trFE,subQuad);CHKERRQ(ierr); 1122 ierr = PetscFESetUp(*trFE);CHKERRQ(ierr); 1123 ierr = PetscQuadratureDestroy(&subQuad);CHKERRQ(ierr); 1124 ierr = PetscSpaceDestroy(&bsubsp);CHKERRQ(ierr); 1125 PetscFunctionReturn(0); 1126 } 1127 1128 PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE) 1129 { 1130 PetscInt hStart, hEnd; 1131 PetscDualSpace dsp; 1132 DM dm; 1133 PetscErrorCode ierr; 1134 1135 PetscFunctionBegin; 1136 PetscValidHeaderSpecific(fe,PETSCFE_CLASSID,1); 1137 PetscValidPointer(trFE,3); 1138 *trFE = NULL; 1139 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 1140 ierr = PetscDualSpaceGetDM(dsp,&dm);CHKERRQ(ierr); 1141 ierr = DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);CHKERRQ(ierr); 1142 if (hEnd <= hStart) PetscFunctionReturn(0); 1143 ierr = PetscFECreatePointTrace(fe,hStart,trFE);CHKERRQ(ierr); 1144 PetscFunctionReturn(0); 1145 } 1146 1147 1148 /*@ 1149 PetscFEGetDimension - Get the dimension of the finite element space on a cell 1150 1151 Not collective 1152 1153 Input Parameter: 1154 . fe - The PetscFE 1155 1156 Output Parameter: 1157 . dim - The dimension 1158 1159 Level: intermediate 1160 1161 .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension() 1162 @*/ 1163 PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) 1164 { 1165 PetscErrorCode ierr; 1166 1167 PetscFunctionBegin; 1168 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1169 PetscValidPointer(dim, 2); 1170 if (fem->ops->getdimension) {ierr = (*fem->ops->getdimension)(fem, dim);CHKERRQ(ierr);} 1171 PetscFunctionReturn(0); 1172 } 1173 1174 /*@C 1175 PetscFEPushforward - Map the reference element function to real space 1176 1177 Input Parameters: 1178 + fe - The PetscFE 1179 . fegeom - The cell geometry 1180 . Nv - The number of function values 1181 - vals - The function values 1182 1183 Output Parameter: 1184 . vals - The transformed function values 1185 1186 Level: advanced 1187 1188 Note: This just forwards the call onto PetscDualSpacePushforward(). 1189 1190 Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1191 1192 .seealso: PetscDualSpacePushforward() 1193 @*/ 1194 PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1195 { 1196 PetscErrorCode ierr; 1197 1198 PetscFunctionBeginHot; 1199 ierr = PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1200 PetscFunctionReturn(0); 1201 } 1202 1203 /*@C 1204 PetscFEPushforwardGradient - Map the reference element function gradient to real space 1205 1206 Input Parameters: 1207 + fe - The PetscFE 1208 . fegeom - The cell geometry 1209 . Nv - The number of function gradient values 1210 - vals - The function gradient values 1211 1212 Output Parameter: 1213 . vals - The transformed function gradient values 1214 1215 Level: advanced 1216 1217 Note: This just forwards the call onto PetscDualSpacePushforwardGradient(). 1218 1219 Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1220 1221 .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward() 1222 @*/ 1223 PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1224 { 1225 PetscErrorCode ierr; 1226 1227 PetscFunctionBeginHot; 1228 ierr = PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1229 PetscFunctionReturn(0); 1230 } 1231 1232 /* 1233 Purpose: Compute element vector for chunk of elements 1234 1235 Input: 1236 Sizes: 1237 Ne: number of elements 1238 Nf: number of fields 1239 PetscFE 1240 dim: spatial dimension 1241 Nb: number of basis functions 1242 Nc: number of field components 1243 PetscQuadrature 1244 Nq: number of quadrature points 1245 1246 Geometry: 1247 PetscFEGeom[Ne] possibly *Nq 1248 PetscReal v0s[dim] 1249 PetscReal n[dim] 1250 PetscReal jacobians[dim*dim] 1251 PetscReal jacobianInverses[dim*dim] 1252 PetscReal jacobianDeterminants 1253 FEM: 1254 PetscFE 1255 PetscQuadrature 1256 PetscReal quadPoints[Nq*dim] 1257 PetscReal quadWeights[Nq] 1258 PetscReal basis[Nq*Nb*Nc] 1259 PetscReal basisDer[Nq*Nb*Nc*dim] 1260 PetscScalar coefficients[Ne*Nb*Nc] 1261 PetscScalar elemVec[Ne*Nb*Nc] 1262 1263 Problem: 1264 PetscInt f: the active field 1265 f0, f1 1266 1267 Work Space: 1268 PetscFE 1269 PetscScalar f0[Nq*dim]; 1270 PetscScalar f1[Nq*dim*dim]; 1271 PetscScalar u[Nc]; 1272 PetscScalar gradU[Nc*dim]; 1273 PetscReal x[dim]; 1274 PetscScalar realSpaceDer[dim]; 1275 1276 Purpose: Compute element vector for N_cb batches of elements 1277 1278 Input: 1279 Sizes: 1280 N_cb: Number of serial cell batches 1281 1282 Geometry: 1283 PetscReal v0s[Ne*dim] 1284 PetscReal jacobians[Ne*dim*dim] possibly *Nq 1285 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 1286 PetscReal jacobianDeterminants[Ne] possibly *Nq 1287 FEM: 1288 static PetscReal quadPoints[Nq*dim] 1289 static PetscReal quadWeights[Nq] 1290 static PetscReal basis[Nq*Nb*Nc] 1291 static PetscReal basisDer[Nq*Nb*Nc*dim] 1292 PetscScalar coefficients[Ne*Nb*Nc] 1293 PetscScalar elemVec[Ne*Nb*Nc] 1294 1295 ex62.c: 1296 PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], 1297 const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], 1298 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), 1299 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) 1300 1301 ex52.c: 1302 PetscErrorCode IntegrateLaplacianBatchCPU(PetscInt Ne, PetscInt Nb, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user) 1303 PetscErrorCode IntegrateElasticityBatchCPU(PetscInt Ne, PetscInt Nb, PetscInt Ncomp, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user) 1304 1305 ex52_integrateElement.cu 1306 __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) 1307 1308 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], 1309 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1310 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1311 1312 ex52_integrateElementOpenCL.c: 1313 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], 1314 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1315 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1316 1317 __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) 1318 */ 1319 1320 /*@C 1321 PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration 1322 1323 Not collective 1324 1325 Input Parameters: 1326 + prob - The PetscDS specifying the discretizations and continuum functions 1327 . field - The field being integrated 1328 . Ne - The number of elements in the chunk 1329 . cgeom - The cell geometry for each cell in the chunk 1330 . coefficients - The array of FEM basis coefficients for the elements 1331 . probAux - The PetscDS specifying the auxiliary discretizations 1332 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1333 1334 Output Parameter: 1335 . integral - the integral for this field 1336 1337 Level: intermediate 1338 1339 .seealso: PetscFEIntegrateResidual() 1340 @*/ 1341 PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 1342 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1343 { 1344 PetscFE fe; 1345 PetscErrorCode ierr; 1346 1347 PetscFunctionBegin; 1348 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1349 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1350 if (fe->ops->integrate) {ierr = (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} 1351 PetscFunctionReturn(0); 1352 } 1353 1354 /*@C 1355 PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration 1356 1357 Not collective 1358 1359 Input Parameters: 1360 + prob - The PetscDS specifying the discretizations and continuum functions 1361 . field - The field being integrated 1362 . obj_func - The function to be integrated 1363 . Ne - The number of elements in the chunk 1364 . fgeom - The face geometry for each face in the chunk 1365 . coefficients - The array of FEM basis coefficients for the elements 1366 . probAux - The PetscDS specifying the auxiliary discretizations 1367 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1368 1369 Output Parameter: 1370 . integral - the integral for this field 1371 1372 Level: intermediate 1373 1374 .seealso: PetscFEIntegrateResidual() 1375 @*/ 1376 PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field, 1377 void (*obj_func)(PetscInt, PetscInt, PetscInt, 1378 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1379 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1380 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 1381 PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1382 { 1383 PetscFE fe; 1384 PetscErrorCode ierr; 1385 1386 PetscFunctionBegin; 1387 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1388 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1389 if (fe->ops->integratebd) {ierr = (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} 1390 PetscFunctionReturn(0); 1391 } 1392 1393 /*@C 1394 PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration 1395 1396 Not collective 1397 1398 Input Parameters: 1399 + prob - The PetscDS specifying the discretizations and continuum functions 1400 . field - The field being integrated 1401 . Ne - The number of elements in the chunk 1402 . cgeom - The cell geometry for each cell in the chunk 1403 . coefficients - The array of FEM basis coefficients for the elements 1404 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1405 . probAux - The PetscDS specifying the auxiliary discretizations 1406 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1407 - t - The time 1408 1409 Output Parameter: 1410 . elemVec - the element residual vectors from each element 1411 1412 Note: 1413 $ Loop over batch of elements (e): 1414 $ Loop over quadrature points (q): 1415 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q 1416 $ Call f_0 and f_1 1417 $ Loop over element vector entries (f,fc --> i): 1418 $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) 1419 1420 Level: intermediate 1421 1422 .seealso: PetscFEIntegrateResidual() 1423 @*/ 1424 PetscErrorCode PetscFEIntegrateResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 1425 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1426 { 1427 PetscFE fe; 1428 PetscErrorCode ierr; 1429 1430 PetscFunctionBegin; 1431 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1432 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1433 if (fe->ops->integrateresidual) {ierr = (*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1434 PetscFunctionReturn(0); 1435 } 1436 1437 /*@C 1438 PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary 1439 1440 Not collective 1441 1442 Input Parameters: 1443 + prob - The PetscDS specifying the discretizations and continuum functions 1444 . field - The field being integrated 1445 . Ne - The number of elements in the chunk 1446 . fgeom - The face geometry for each cell in the chunk 1447 . coefficients - The array of FEM basis coefficients for the elements 1448 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1449 . probAux - The PetscDS specifying the auxiliary discretizations 1450 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1451 - t - The time 1452 1453 Output Parameter: 1454 . elemVec - the element residual vectors from each element 1455 1456 Level: intermediate 1457 1458 .seealso: PetscFEIntegrateResidual() 1459 @*/ 1460 PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, 1461 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1462 { 1463 PetscFE fe; 1464 PetscErrorCode ierr; 1465 1466 PetscFunctionBegin; 1467 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1468 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1469 if (fe->ops->integratebdresidual) {ierr = (*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1470 PetscFunctionReturn(0); 1471 } 1472 1473 /*@C 1474 PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration 1475 1476 Not collective 1477 1478 Input Parameters: 1479 + prob - The PetscDS specifying the discretizations and continuum functions 1480 . field - The field being integrated 1481 . Ne - The number of elements in the chunk 1482 . fgeom - The face geometry for each cell in the chunk 1483 . coefficients - The array of FEM basis coefficients for the elements 1484 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1485 . probAux - The PetscDS specifying the auxiliary discretizations 1486 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1487 - t - The time 1488 1489 Output Parameter 1490 . elemVec - the element residual vectors from each element 1491 1492 Level: developer 1493 1494 .seealso: PetscFEIntegrateResidual() 1495 @*/ 1496 PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, 1497 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1498 { 1499 PetscFE fe; 1500 PetscErrorCode ierr; 1501 1502 PetscFunctionBegin; 1503 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1504 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1505 if (fe->ops->integratehybridresidual) {ierr = (*fe->ops->integratehybridresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1506 PetscFunctionReturn(0); 1507 } 1508 1509 /*@C 1510 PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration 1511 1512 Not collective 1513 1514 Input Parameters: 1515 + prob - The PetscDS specifying the discretizations and continuum functions 1516 . jtype - The type of matrix pointwise functions that should be used 1517 . fieldI - The test field being integrated 1518 . fieldJ - The basis field being integrated 1519 . Ne - The number of elements in the chunk 1520 . cgeom - The cell geometry for each cell in the chunk 1521 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1522 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1523 . probAux - The PetscDS specifying the auxiliary discretizations 1524 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1525 . t - The time 1526 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1527 1528 Output Parameter: 1529 . elemMat - the element matrices for the Jacobian from each element 1530 1531 Note: 1532 $ Loop over batch of elements (e): 1533 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1534 $ Loop over quadrature points (q): 1535 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1536 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1537 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1538 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1539 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1540 Level: intermediate 1541 1542 .seealso: PetscFEIntegrateResidual() 1543 @*/ 1544 PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom, 1545 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1546 { 1547 PetscFE fe; 1548 PetscErrorCode ierr; 1549 1550 PetscFunctionBegin; 1551 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1552 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1553 if (fe->ops->integratejacobian) {ierr = (*fe->ops->integratejacobian)(prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1554 PetscFunctionReturn(0); 1555 } 1556 1557 /*@C 1558 PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration 1559 1560 Not collective 1561 1562 Input Parameters: 1563 + prob - The PetscDS specifying the discretizations and continuum functions 1564 . fieldI - The test field being integrated 1565 . fieldJ - The basis field being integrated 1566 . Ne - The number of elements in the chunk 1567 . fgeom - The face geometry for each cell in the chunk 1568 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1569 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1570 . probAux - The PetscDS specifying the auxiliary discretizations 1571 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1572 . t - The time 1573 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1574 1575 Output Parameter: 1576 . elemMat - the element matrices for the Jacobian from each element 1577 1578 Note: 1579 $ Loop over batch of elements (e): 1580 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1581 $ Loop over quadrature points (q): 1582 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1583 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1584 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1585 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1586 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1587 Level: intermediate 1588 1589 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1590 @*/ 1591 PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, 1592 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1593 { 1594 PetscFE fe; 1595 PetscErrorCode ierr; 1596 1597 PetscFunctionBegin; 1598 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1599 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1600 if (fe->ops->integratebdjacobian) {ierr = (*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1601 PetscFunctionReturn(0); 1602 } 1603 1604 /*@C 1605 PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration 1606 1607 Not collective 1608 1609 Input Parameters: 1610 . prob - The PetscDS specifying the discretizations and continuum functions 1611 . jtype - The type of matrix pointwise functions that should be used 1612 . fieldI - The test field being integrated 1613 . fieldJ - The basis field being integrated 1614 . Ne - The number of elements in the chunk 1615 . fgeom - The face geometry for each cell in the chunk 1616 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1617 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1618 . probAux - The PetscDS specifying the auxiliary discretizations 1619 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1620 . t - The time 1621 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1622 1623 Output Parameter 1624 . elemMat - the element matrices for the Jacobian from each element 1625 1626 Note: 1627 $ Loop over batch of elements (e): 1628 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1629 $ Loop over quadrature points (q): 1630 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1631 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1632 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1633 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1634 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1635 Level: developer 1636 1637 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1638 @*/ 1639 PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, 1640 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1641 { 1642 PetscFE fe; 1643 PetscErrorCode ierr; 1644 1645 PetscFunctionBegin; 1646 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1647 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1648 if (fe->ops->integratehybridjacobian) {ierr = (*fe->ops->integratehybridjacobian)(prob, jtype, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1649 PetscFunctionReturn(0); 1650 } 1651 1652 /*@ 1653 PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height 1654 1655 Input Parameters: 1656 + fe - The finite element space 1657 - height - The height of the Plex point 1658 1659 Output Parameter: 1660 . subfe - The subspace of this FE space 1661 1662 Note: For example, if we want the subspace of this space for a face, we would choose height = 1. 1663 1664 Level: advanced 1665 1666 .seealso: PetscFECreateDefault() 1667 @*/ 1668 PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe) 1669 { 1670 PetscSpace P, subP; 1671 PetscDualSpace Q, subQ; 1672 PetscQuadrature subq; 1673 PetscFEType fetype; 1674 PetscInt dim, Nc; 1675 PetscErrorCode ierr; 1676 1677 PetscFunctionBegin; 1678 PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1); 1679 PetscValidPointer(subfe, 3); 1680 if (height == 0) { 1681 *subfe = fe; 1682 PetscFunctionReturn(0); 1683 } 1684 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1685 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1686 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1687 ierr = PetscFEGetFaceQuadrature(fe, &subq);CHKERRQ(ierr); 1688 ierr = PetscDualSpaceGetDimension(Q, &dim);CHKERRQ(ierr); 1689 if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);} 1690 if (!fe->subspaces) {ierr = PetscCalloc1(dim, &fe->subspaces);CHKERRQ(ierr);} 1691 if (height <= dim) { 1692 if (!fe->subspaces[height-1]) { 1693 PetscFE sub; 1694 const char *name; 1695 1696 ierr = PetscSpaceGetHeightSubspace(P, height, &subP);CHKERRQ(ierr); 1697 ierr = PetscDualSpaceGetHeightSubspace(Q, height, &subQ);CHKERRQ(ierr); 1698 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);CHKERRQ(ierr); 1699 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 1700 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 1701 ierr = PetscFEGetType(fe, &fetype);CHKERRQ(ierr); 1702 ierr = PetscFESetType(sub, fetype);CHKERRQ(ierr); 1703 ierr = PetscFESetBasisSpace(sub, subP);CHKERRQ(ierr); 1704 ierr = PetscFESetDualSpace(sub, subQ);CHKERRQ(ierr); 1705 ierr = PetscFESetNumComponents(sub, Nc);CHKERRQ(ierr); 1706 ierr = PetscFESetUp(sub);CHKERRQ(ierr); 1707 ierr = PetscFESetQuadrature(sub, subq);CHKERRQ(ierr); 1708 fe->subspaces[height-1] = sub; 1709 } 1710 *subfe = fe->subspaces[height-1]; 1711 } else { 1712 *subfe = NULL; 1713 } 1714 PetscFunctionReturn(0); 1715 } 1716 1717 /*@ 1718 PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used 1719 to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more 1720 sparsity). It is also used to create an interpolation between regularly refined meshes. 1721 1722 Collective on fem 1723 1724 Input Parameter: 1725 . fe - The initial PetscFE 1726 1727 Output Parameter: 1728 . feRef - The refined PetscFE 1729 1730 Level: advanced 1731 1732 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 1733 @*/ 1734 PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef) 1735 { 1736 PetscSpace P, Pref; 1737 PetscDualSpace Q, Qref; 1738 DM K, Kref; 1739 PetscQuadrature q, qref; 1740 const PetscReal *v0, *jac; 1741 PetscInt numComp, numSubelements; 1742 PetscInt cStart, cEnd, c; 1743 PetscDualSpace *cellSpaces; 1744 PetscErrorCode ierr; 1745 1746 PetscFunctionBegin; 1747 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1748 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1749 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1750 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr); 1751 /* Create space */ 1752 ierr = PetscObjectReference((PetscObject) P);CHKERRQ(ierr); 1753 Pref = P; 1754 /* Create dual space */ 1755 ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr); 1756 ierr = PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);CHKERRQ(ierr); 1757 ierr = DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);CHKERRQ(ierr); 1758 ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr); 1759 ierr = DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);CHKERRQ(ierr); 1760 ierr = PetscMalloc1(cEnd - cStart, &cellSpaces);CHKERRQ(ierr); 1761 /* TODO: fix for non-uniform refinement */ 1762 for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q; 1763 ierr = PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);CHKERRQ(ierr); 1764 ierr = PetscFree(cellSpaces);CHKERRQ(ierr); 1765 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 1766 ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr); 1767 /* Create element */ 1768 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);CHKERRQ(ierr); 1769 ierr = PetscFESetType(*feRef, PETSCFECOMPOSITE);CHKERRQ(ierr); 1770 ierr = PetscFESetBasisSpace(*feRef, Pref);CHKERRQ(ierr); 1771 ierr = PetscFESetDualSpace(*feRef, Qref);CHKERRQ(ierr); 1772 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 1773 ierr = PetscFESetNumComponents(*feRef, numComp);CHKERRQ(ierr); 1774 ierr = PetscFESetUp(*feRef);CHKERRQ(ierr); 1775 ierr = PetscSpaceDestroy(&Pref);CHKERRQ(ierr); 1776 ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr); 1777 /* Create quadrature */ 1778 ierr = PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);CHKERRQ(ierr); 1779 ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr); 1780 ierr = PetscFESetQuadrature(*feRef, qref);CHKERRQ(ierr); 1781 ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr); 1782 PetscFunctionReturn(0); 1783 } 1784 1785 /*@C 1786 PetscFECreateDefault - Create a PetscFE for basic FEM computation 1787 1788 Collective 1789 1790 Input Parameters: 1791 + comm - The MPI comm 1792 . dim - The spatial dimension 1793 . Nc - The number of components 1794 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1795 . prefix - The options prefix, or NULL 1796 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1797 1798 Output Parameter: 1799 . fem - The PetscFE object 1800 1801 Note: 1802 Each object is SetFromOption() during creation, so that the object may be customized from the command line. 1803 1804 Level: beginner 1805 1806 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1807 @*/ 1808 PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem) 1809 { 1810 PetscQuadrature q, fq; 1811 DM K; 1812 PetscSpace P; 1813 PetscDualSpace Q; 1814 PetscInt order, quadPointsPerEdge; 1815 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1816 PetscErrorCode ierr; 1817 1818 PetscFunctionBegin; 1819 /* Create space */ 1820 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1821 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr); 1822 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1823 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1824 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1825 ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr); 1826 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1827 ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr); 1828 ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr); 1829 /* Create dual space */ 1830 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1831 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1832 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr); 1833 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1834 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1835 ierr = DMDestroy(&K);CHKERRQ(ierr); 1836 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1837 ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr); 1838 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1839 ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr); 1840 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1841 /* Create element */ 1842 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1843 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr); 1844 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1845 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1846 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1847 ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr); 1848 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1849 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1850 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1851 /* Create quadrature (with specified order if given) */ 1852 qorder = qorder >= 0 ? qorder : order; 1853 ierr = PetscObjectOptionsBegin((PetscObject)*fem);CHKERRQ(ierr); 1854 ierr = PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);CHKERRQ(ierr); 1855 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1856 quadPointsPerEdge = PetscMax(qorder + 1,1); 1857 if (isSimplex) { 1858 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1859 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1860 } else { 1861 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1862 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1863 } 1864 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1865 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1866 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1867 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1868 PetscFunctionReturn(0); 1869 } 1870 1871 /*@ 1872 PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k 1873 1874 Collective 1875 1876 Input Parameters: 1877 + comm - The MPI comm 1878 . dim - The spatial dimension 1879 . Nc - The number of components 1880 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1881 . k - The degree k of the space 1882 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1883 1884 Output Parameter: 1885 . fem - The PetscFE object 1886 1887 Level: beginner 1888 1889 Notes: 1890 For simplices, this element is the space of maximum polynomial degree k, otherwise it is a tensor product of 1D polynomials, each with maximal degree k. 1891 1892 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1893 @*/ 1894 PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem) 1895 { 1896 PetscQuadrature q, fq; 1897 DM K; 1898 PetscSpace P; 1899 PetscDualSpace Q; 1900 PetscInt quadPointsPerEdge; 1901 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1902 char name[64]; 1903 PetscErrorCode ierr; 1904 1905 PetscFunctionBegin; 1906 /* Create space */ 1907 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1908 ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1909 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1910 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1911 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1912 ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr); 1913 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1914 /* Create dual space */ 1915 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1916 ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1917 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1918 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1919 ierr = DMDestroy(&K);CHKERRQ(ierr); 1920 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1921 ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr); 1922 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1923 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1924 /* Create element */ 1925 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1926 ierr = PetscSNPrintf(name, 64, "P%d", (int) k);CHKERRQ(ierr); 1927 ierr = PetscObjectSetName((PetscObject) *fem, name);CHKERRQ(ierr); 1928 ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr); 1929 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1930 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1931 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1932 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1933 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1934 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1935 /* Create quadrature (with specified order if given) */ 1936 qorder = qorder >= 0 ? qorder : k; 1937 quadPointsPerEdge = PetscMax(qorder + 1,1); 1938 if (isSimplex) { 1939 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1940 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1941 } else { 1942 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1943 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1944 } 1945 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1946 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1947 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1948 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1949 PetscFunctionReturn(0); 1950 } 1951 1952 /*@C 1953 PetscFESetName - Names the FE and its subobjects 1954 1955 Not collective 1956 1957 Input Parameters: 1958 + fe - The PetscFE 1959 - name - The name 1960 1961 Level: intermediate 1962 1963 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1964 @*/ 1965 PetscErrorCode PetscFESetName(PetscFE fe, const char name[]) 1966 { 1967 PetscSpace P; 1968 PetscDualSpace Q; 1969 PetscErrorCode ierr; 1970 1971 PetscFunctionBegin; 1972 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1973 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1974 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 1975 ierr = PetscObjectSetName((PetscObject) P, name);CHKERRQ(ierr); 1976 ierr = PetscObjectSetName((PetscObject) Q, name);CHKERRQ(ierr); 1977 PetscFunctionReturn(0); 1978 } 1979 1980 PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS ds, PetscInt Nf, PetscInt r, PetscInt q, PetscTabulation T[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[]) 1981 { 1982 PetscInt dOffset = 0, fOffset = 0, f; 1983 PetscErrorCode ierr; 1984 1985 for (f = 0; f < Nf; ++f) { 1986 PetscFE fe; 1987 const PetscInt cdim = T[f]->cdim; 1988 const PetscInt Nq = T[f]->Np; 1989 const PetscInt Nbf = T[f]->Nb; 1990 const PetscInt Ncf = T[f]->Nc; 1991 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 1992 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 1993 PetscInt b, c, d; 1994 1995 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 1996 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 1997 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 1998 for (b = 0; b < Nbf; ++b) { 1999 for (c = 0; c < Ncf; ++c) { 2000 const PetscInt cidx = b*Ncf+c; 2001 2002 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2003 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2004 } 2005 } 2006 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2007 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2008 if (u_t) { 2009 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2010 for (b = 0; b < Nbf; ++b) { 2011 for (c = 0; c < Ncf; ++c) { 2012 const PetscInt cidx = b*Ncf+c; 2013 2014 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2015 } 2016 } 2017 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2018 } 2019 fOffset += Ncf; 2020 dOffset += Nbf; 2021 } 2022 return 0; 2023 } 2024 2025 PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS ds, PetscInt dim, PetscInt Nf, const PetscInt Nb[], const PetscInt Nc[], PetscInt q, PetscReal *basisField[], PetscReal *basisFieldDer[], PetscFEGeom *fegeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[]) 2026 { 2027 const PetscInt dE = fegeom->dimEmbed; 2028 PetscInt dOffset = 0, fOffset = 0, g; 2029 PetscErrorCode ierr; 2030 2031 for (g = 0; g < 2*Nf-1; ++g) { 2032 PetscFE fe; 2033 const PetscInt f = g/2; 2034 const PetscInt Nbf = Nb[f], Ncf = Nc[f]; 2035 const PetscReal *Bq = &basisField[f][q*Nbf*Ncf]; 2036 const PetscReal *Dq = &basisFieldDer[f][q*Nbf*Ncf*dim]; 2037 PetscInt b, c, d; 2038 2039 fe = (PetscFE) ds->disc[f]; 2040 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2041 for (d = 0; d < dim*Ncf; ++d) u_x[fOffset*dim+d] = 0.0; 2042 for (b = 0; b < Nbf; ++b) { 2043 for (c = 0; c < Ncf; ++c) { 2044 const PetscInt cidx = b*Ncf+c; 2045 2046 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2047 for (d = 0; d < dim; ++d) u_x[(fOffset+c)*dE+d] += Dq[cidx*dim+d]*coefficients[dOffset+b]; 2048 } 2049 } 2050 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2051 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*dim]);CHKERRQ(ierr); 2052 if (u_t) { 2053 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2054 for (b = 0; b < Nbf; ++b) { 2055 for (c = 0; c < Ncf; ++c) { 2056 const PetscInt cidx = b*Ncf+c; 2057 2058 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2059 } 2060 } 2061 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2062 } 2063 fOffset += Ncf; 2064 dOffset += Nbf; 2065 } 2066 return 0; 2067 } 2068 2069 PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[]) 2070 { 2071 PetscFE fe; 2072 PetscTabulation Tc; 2073 PetscInt b, c; 2074 PetscErrorCode ierr; 2075 2076 if (!prob) return 0; 2077 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 2078 ierr = PetscFEGetFaceCentroidTabulation(fe, &Tc);CHKERRQ(ierr); 2079 { 2080 const PetscReal *faceBasis = Tc->T[0]; 2081 const PetscInt Nb = Tc->Nb; 2082 const PetscInt Nc = Tc->Nc; 2083 2084 for (c = 0; c < Nc; ++c) {u[c] = 0.0;} 2085 for (b = 0; b < Nb; ++b) { 2086 for (c = 0; c < Nc; ++c) { 2087 const PetscInt cidx = b*Nc+c; 2088 2089 u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx]; 2090 } 2091 } 2092 } 2093 return 0; 2094 } 2095 2096 PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2097 { 2098 const PetscInt dE = T->cdim; /* fegeom->dimEmbed */ 2099 const PetscInt Nq = T->Np; 2100 const PetscInt Nb = T->Nb; 2101 const PetscInt Nc = T->Nc; 2102 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2103 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dim]; 2104 PetscInt q, b, c, d; 2105 PetscErrorCode ierr; 2106 2107 for (b = 0; b < Nb; ++b) elemVec[b] = 0.0; 2108 for (q = 0; q < Nq; ++q) { 2109 for (b = 0; b < Nb; ++b) { 2110 for (c = 0; c < Nc; ++c) { 2111 const PetscInt bcidx = b*Nc+c; 2112 2113 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2114 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2115 } 2116 } 2117 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2118 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2119 for (b = 0; b < Nb; ++b) { 2120 for (c = 0; c < Nc; ++c) { 2121 const PetscInt bcidx = b*Nc+c; 2122 const PetscInt qcidx = q*Nc+c; 2123 2124 elemVec[b] += tmpBasis[bcidx]*f0[qcidx]; 2125 for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2126 } 2127 } 2128 } 2129 return(0); 2130 } 2131 2132 PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2133 { 2134 const PetscInt dE = T->cdim; 2135 const PetscInt Nq = T->Np; 2136 const PetscInt Nb = T->Nb; 2137 const PetscInt Nc = T->Nc; 2138 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2139 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2140 PetscInt q, b, c, d, s; 2141 PetscErrorCode ierr; 2142 2143 for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0; 2144 for (q = 0; q < Nq; ++q) { 2145 for (b = 0; b < Nb; ++b) { 2146 for (c = 0; c < Nc; ++c) { 2147 const PetscInt bcidx = b*Nc+c; 2148 2149 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2150 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2151 } 2152 } 2153 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2154 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2155 for (s = 0; s < 2; ++s) { 2156 for (b = 0; b < Nb; ++b) { 2157 for (c = 0; c < Nc; ++c) { 2158 const PetscInt bcidx = b*Nc+c; 2159 const PetscInt qcidx = (q*2+s)*Nc+c; 2160 2161 elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx]; 2162 for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2163 } 2164 } 2165 } 2166 } 2167 return(0); 2168 } 2169 2170 PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE feI, PetscFE feJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[]) 2171 { 2172 const PetscInt dE = TI->cdim; 2173 const PetscInt NqI = TI->Np; 2174 const PetscInt NbI = TI->Nb; 2175 const PetscInt NcI = TI->Nc; 2176 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2177 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dim]; 2178 const PetscInt NqJ = TJ->Np; 2179 const PetscInt NbJ = TJ->Nb; 2180 const PetscInt NcJ = TJ->Nc; 2181 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2182 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dim]; 2183 PetscInt f, fc, g, gc, df, dg; 2184 PetscErrorCode ierr; 2185 2186 for (f = 0; f < NbI; ++f) { 2187 for (fc = 0; fc < NcI; ++fc) { 2188 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2189 2190 tmpBasisI[fidx] = basisI[fidx]; 2191 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2192 } 2193 } 2194 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2195 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2196 for (g = 0; g < NbJ; ++g) { 2197 for (gc = 0; gc < NcJ; ++gc) { 2198 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2199 2200 tmpBasisJ[gidx] = basisJ[gidx]; 2201 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2202 } 2203 } 2204 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2205 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2206 for (f = 0; f < NbI; ++f) { 2207 for (fc = 0; fc < NcI; ++fc) { 2208 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2209 const PetscInt i = offsetI+f; /* Element matrix row */ 2210 for (g = 0; g < NbJ; ++g) { 2211 for (gc = 0; gc < NcJ; ++gc) { 2212 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2213 const PetscInt j = offsetJ+g; /* Element matrix column */ 2214 const PetscInt fOff = eOffset+i*totDim+j; 2215 2216 elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx]; 2217 for (df = 0; df < dE; ++df) { 2218 elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2219 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx]; 2220 for (dg = 0; dg < dE; ++dg) { 2221 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2222 } 2223 } 2224 } 2225 } 2226 } 2227 } 2228 return(0); 2229 } 2230 2231 PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt dim, PetscInt NbI, PetscInt NcI, const PetscReal basisI[], const PetscReal basisDerI[], PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscInt NbJ, PetscInt NcJ, const PetscReal basisJ[], const PetscReal basisDerJ[], PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[]) 2232 { 2233 const PetscInt dE = fegeom->dimEmbed; 2234 const PetscInt Ns = isHybridI ? 1 : 2; 2235 const PetscInt Nt = isHybridJ ? 1 : 2; 2236 PetscInt f, fc, g, gc, df, dg, s, t; 2237 PetscErrorCode ierr; 2238 2239 for (f = 0; f < NbI; ++f) { 2240 for (fc = 0; fc < NcI; ++fc) { 2241 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2242 2243 tmpBasisI[fidx] = basisI[fidx]; 2244 for (df = 0; df < dim; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dim+df]; 2245 } 2246 } 2247 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2248 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2249 for (g = 0; g < NbJ; ++g) { 2250 for (gc = 0; gc < NcJ; ++gc) { 2251 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2252 2253 tmpBasisJ[gidx] = basisJ[gidx]; 2254 for (dg = 0; dg < dim; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dim+dg]; 2255 } 2256 } 2257 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2258 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2259 for (s = 0; s < Ns; ++s) { 2260 for (f = 0; f < NbI; ++f) { 2261 for (fc = 0; fc < NcI; ++fc) { 2262 const PetscInt sc = NcI*s+fc; /* components from each side of the surface */ 2263 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2264 const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */ 2265 for (t = 0; t < Nt; ++t) { 2266 for (g = 0; g < NbJ; ++g) { 2267 for (gc = 0; gc < NcJ; ++gc) { 2268 const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */ 2269 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2270 const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */ 2271 const PetscInt fOff = eOffset+i*totDim+j; 2272 2273 elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx]; 2274 for (df = 0; df < dE; ++df) { 2275 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2276 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx]; 2277 for (dg = 0; dg < dE; ++dg) { 2278 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2279 } 2280 } 2281 } 2282 } 2283 } 2284 } 2285 } 2286 } 2287 return(0); 2288 } 2289 2290 PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom) 2291 { 2292 PetscDualSpace dsp; 2293 DM dm; 2294 PetscQuadrature quadDef; 2295 PetscInt dim, cdim, Nq; 2296 PetscErrorCode ierr; 2297 2298 PetscFunctionBegin; 2299 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); 2300 ierr = PetscDualSpaceGetDM(dsp, &dm);CHKERRQ(ierr); 2301 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2302 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 2303 ierr = PetscFEGetQuadrature(fe, &quadDef);CHKERRQ(ierr); 2304 quad = quad ? quad : quadDef; 2305 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2306 ierr = PetscMalloc1(Nq*cdim, &cgeom->v);CHKERRQ(ierr); 2307 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->J);CHKERRQ(ierr); 2308 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);CHKERRQ(ierr); 2309 ierr = PetscMalloc1(Nq, &cgeom->detJ);CHKERRQ(ierr); 2310 cgeom->dim = dim; 2311 cgeom->dimEmbed = cdim; 2312 cgeom->numCells = 1; 2313 cgeom->numPoints = Nq; 2314 ierr = DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);CHKERRQ(ierr); 2315 PetscFunctionReturn(0); 2316 } 2317 2318 PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom) 2319 { 2320 PetscErrorCode ierr; 2321 2322 PetscFunctionBegin; 2323 ierr = PetscFree(cgeom->v);CHKERRQ(ierr); 2324 ierr = PetscFree(cgeom->J);CHKERRQ(ierr); 2325 ierr = PetscFree(cgeom->invJ);CHKERRQ(ierr); 2326 ierr = PetscFree(cgeom->detJ);CHKERRQ(ierr); 2327 PetscFunctionReturn(0); 2328 } 2329