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 = NULL; 1694 const char *name; 1695 1696 ierr = PetscSpaceGetHeightSubspace(P, height, &subP);CHKERRQ(ierr); 1697 ierr = PetscDualSpaceGetHeightSubspace(Q, height, &subQ);CHKERRQ(ierr); 1698 if (subQ) { 1699 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);CHKERRQ(ierr); 1700 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 1701 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 1702 ierr = PetscFEGetType(fe, &fetype);CHKERRQ(ierr); 1703 ierr = PetscFESetType(sub, fetype);CHKERRQ(ierr); 1704 ierr = PetscFESetBasisSpace(sub, subP);CHKERRQ(ierr); 1705 ierr = PetscFESetDualSpace(sub, subQ);CHKERRQ(ierr); 1706 ierr = PetscFESetNumComponents(sub, Nc);CHKERRQ(ierr); 1707 ierr = PetscFESetUp(sub);CHKERRQ(ierr); 1708 ierr = PetscFESetQuadrature(sub, subq);CHKERRQ(ierr); 1709 } 1710 fe->subspaces[height-1] = sub; 1711 } 1712 *subfe = fe->subspaces[height-1]; 1713 } else { 1714 *subfe = NULL; 1715 } 1716 PetscFunctionReturn(0); 1717 } 1718 1719 /*@ 1720 PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used 1721 to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more 1722 sparsity). It is also used to create an interpolation between regularly refined meshes. 1723 1724 Collective on fem 1725 1726 Input Parameter: 1727 . fe - The initial PetscFE 1728 1729 Output Parameter: 1730 . feRef - The refined PetscFE 1731 1732 Level: advanced 1733 1734 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 1735 @*/ 1736 PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef) 1737 { 1738 PetscSpace P, Pref; 1739 PetscDualSpace Q, Qref; 1740 DM K, Kref; 1741 PetscQuadrature q, qref; 1742 const PetscReal *v0, *jac; 1743 PetscInt numComp, numSubelements; 1744 PetscInt cStart, cEnd, c; 1745 PetscDualSpace *cellSpaces; 1746 PetscErrorCode ierr; 1747 1748 PetscFunctionBegin; 1749 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1750 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1751 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1752 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr); 1753 /* Create space */ 1754 ierr = PetscObjectReference((PetscObject) P);CHKERRQ(ierr); 1755 Pref = P; 1756 /* Create dual space */ 1757 ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr); 1758 ierr = PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);CHKERRQ(ierr); 1759 ierr = DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);CHKERRQ(ierr); 1760 ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr); 1761 ierr = DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);CHKERRQ(ierr); 1762 ierr = PetscMalloc1(cEnd - cStart, &cellSpaces);CHKERRQ(ierr); 1763 /* TODO: fix for non-uniform refinement */ 1764 for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q; 1765 ierr = PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);CHKERRQ(ierr); 1766 ierr = PetscFree(cellSpaces);CHKERRQ(ierr); 1767 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 1768 ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr); 1769 /* Create element */ 1770 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);CHKERRQ(ierr); 1771 ierr = PetscFESetType(*feRef, PETSCFECOMPOSITE);CHKERRQ(ierr); 1772 ierr = PetscFESetBasisSpace(*feRef, Pref);CHKERRQ(ierr); 1773 ierr = PetscFESetDualSpace(*feRef, Qref);CHKERRQ(ierr); 1774 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 1775 ierr = PetscFESetNumComponents(*feRef, numComp);CHKERRQ(ierr); 1776 ierr = PetscFESetUp(*feRef);CHKERRQ(ierr); 1777 ierr = PetscSpaceDestroy(&Pref);CHKERRQ(ierr); 1778 ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr); 1779 /* Create quadrature */ 1780 ierr = PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);CHKERRQ(ierr); 1781 ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr); 1782 ierr = PetscFESetQuadrature(*feRef, qref);CHKERRQ(ierr); 1783 ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr); 1784 PetscFunctionReturn(0); 1785 } 1786 1787 /*@C 1788 PetscFECreateDefault - Create a PetscFE for basic FEM computation 1789 1790 Collective 1791 1792 Input Parameters: 1793 + comm - The MPI comm 1794 . dim - The spatial dimension 1795 . Nc - The number of components 1796 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1797 . prefix - The options prefix, or NULL 1798 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1799 1800 Output Parameter: 1801 . fem - The PetscFE object 1802 1803 Note: 1804 Each object is SetFromOption() during creation, so that the object may be customized from the command line. 1805 1806 Level: beginner 1807 1808 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1809 @*/ 1810 PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem) 1811 { 1812 PetscQuadrature q, fq; 1813 DM K; 1814 PetscSpace P; 1815 PetscDualSpace Q; 1816 PetscInt order, quadPointsPerEdge; 1817 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1818 PetscErrorCode ierr; 1819 1820 PetscFunctionBegin; 1821 /* Create space */ 1822 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1823 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr); 1824 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1825 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1826 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1827 ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr); 1828 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1829 ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr); 1830 ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr); 1831 /* Create dual space */ 1832 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1833 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1834 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr); 1835 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1836 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1837 ierr = DMDestroy(&K);CHKERRQ(ierr); 1838 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1839 ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr); 1840 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1841 ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr); 1842 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1843 /* Create element */ 1844 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1845 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr); 1846 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1847 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1848 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1849 ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr); 1850 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1851 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1852 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1853 /* Create quadrature (with specified order if given) */ 1854 qorder = qorder >= 0 ? qorder : order; 1855 ierr = PetscObjectOptionsBegin((PetscObject)*fem);CHKERRQ(ierr); 1856 ierr = PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);CHKERRQ(ierr); 1857 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1858 quadPointsPerEdge = PetscMax(qorder + 1,1); 1859 if (isSimplex) { 1860 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1861 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1862 } else { 1863 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1864 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1865 } 1866 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1867 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1868 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1869 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1870 PetscFunctionReturn(0); 1871 } 1872 1873 /*@ 1874 PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k 1875 1876 Collective 1877 1878 Input Parameters: 1879 + comm - The MPI comm 1880 . dim - The spatial dimension 1881 . Nc - The number of components 1882 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1883 . k - The degree k of the space 1884 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1885 1886 Output Parameter: 1887 . fem - The PetscFE object 1888 1889 Level: beginner 1890 1891 Notes: 1892 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. 1893 1894 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1895 @*/ 1896 PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem) 1897 { 1898 PetscQuadrature q, fq; 1899 DM K; 1900 PetscSpace P; 1901 PetscDualSpace Q; 1902 PetscInt quadPointsPerEdge; 1903 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1904 char name[64]; 1905 PetscErrorCode ierr; 1906 1907 PetscFunctionBegin; 1908 /* Create space */ 1909 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1910 ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1911 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1912 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1913 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1914 ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr); 1915 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1916 /* Create dual space */ 1917 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1918 ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1919 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1920 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1921 ierr = DMDestroy(&K);CHKERRQ(ierr); 1922 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1923 ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr); 1924 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1925 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1926 /* Create element */ 1927 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1928 ierr = PetscSNPrintf(name, 64, "P%d", (int) k);CHKERRQ(ierr); 1929 ierr = PetscObjectSetName((PetscObject) *fem, name);CHKERRQ(ierr); 1930 ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr); 1931 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1932 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1933 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1934 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1935 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1936 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1937 /* Create quadrature (with specified order if given) */ 1938 qorder = qorder >= 0 ? qorder : k; 1939 quadPointsPerEdge = PetscMax(qorder + 1,1); 1940 if (isSimplex) { 1941 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1942 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1943 } else { 1944 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1945 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1946 } 1947 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1948 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1949 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1950 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1951 PetscFunctionReturn(0); 1952 } 1953 1954 /*@C 1955 PetscFESetName - Names the FE and its subobjects 1956 1957 Not collective 1958 1959 Input Parameters: 1960 + fe - The PetscFE 1961 - name - The name 1962 1963 Level: intermediate 1964 1965 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1966 @*/ 1967 PetscErrorCode PetscFESetName(PetscFE fe, const char name[]) 1968 { 1969 PetscSpace P; 1970 PetscDualSpace Q; 1971 PetscErrorCode ierr; 1972 1973 PetscFunctionBegin; 1974 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1975 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1976 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 1977 ierr = PetscObjectSetName((PetscObject) P, name);CHKERRQ(ierr); 1978 ierr = PetscObjectSetName((PetscObject) Q, name);CHKERRQ(ierr); 1979 PetscFunctionReturn(0); 1980 } 1981 1982 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[]) 1983 { 1984 PetscInt dOffset = 0, fOffset = 0, f; 1985 PetscErrorCode ierr; 1986 1987 for (f = 0; f < Nf; ++f) { 1988 PetscFE fe; 1989 const PetscInt cdim = T[f]->cdim; 1990 const PetscInt Nq = T[f]->Np; 1991 const PetscInt Nbf = T[f]->Nb; 1992 const PetscInt Ncf = T[f]->Nc; 1993 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 1994 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 1995 PetscInt b, c, d; 1996 1997 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 1998 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 1999 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2000 for (b = 0; b < Nbf; ++b) { 2001 for (c = 0; c < Ncf; ++c) { 2002 const PetscInt cidx = b*Ncf+c; 2003 2004 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2005 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2006 } 2007 } 2008 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2009 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2010 if (u_t) { 2011 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2012 for (b = 0; b < Nbf; ++b) { 2013 for (c = 0; c < Ncf; ++c) { 2014 const PetscInt cidx = b*Ncf+c; 2015 2016 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2017 } 2018 } 2019 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2020 } 2021 fOffset += Ncf; 2022 dOffset += Nbf; 2023 } 2024 return 0; 2025 } 2026 2027 PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_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[]) 2028 { 2029 PetscInt dOffset = 0, fOffset = 0, g; 2030 PetscErrorCode ierr; 2031 2032 for (g = 0; g < 2*Nf-1; ++g) { 2033 if (!T[g/2]) continue; 2034 { 2035 PetscFE fe; 2036 const PetscInt f = g/2; 2037 const PetscInt cdim = T[f]->cdim; 2038 const PetscInt Nq = T[f]->Np; 2039 const PetscInt Nbf = T[f]->Nb; 2040 const PetscInt Ncf = T[f]->Nc; 2041 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2042 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2043 PetscInt b, c, d; 2044 2045 fe = (PetscFE) ds->disc[f]; 2046 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2047 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2048 for (b = 0; b < Nbf; ++b) { 2049 for (c = 0; c < Ncf; ++c) { 2050 const PetscInt cidx = b*Ncf+c; 2051 2052 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2053 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2054 } 2055 } 2056 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2057 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2058 if (u_t) { 2059 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2060 for (b = 0; b < Nbf; ++b) { 2061 for (c = 0; c < Ncf; ++c) { 2062 const PetscInt cidx = b*Ncf+c; 2063 2064 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2065 } 2066 } 2067 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2068 } 2069 fOffset += Ncf; 2070 dOffset += Nbf; 2071 } 2072 } 2073 return 0; 2074 } 2075 2076 PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[]) 2077 { 2078 PetscFE fe; 2079 PetscTabulation Tc; 2080 PetscInt b, c; 2081 PetscErrorCode ierr; 2082 2083 if (!prob) return 0; 2084 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 2085 ierr = PetscFEGetFaceCentroidTabulation(fe, &Tc);CHKERRQ(ierr); 2086 { 2087 const PetscReal *faceBasis = Tc->T[0]; 2088 const PetscInt Nb = Tc->Nb; 2089 const PetscInt Nc = Tc->Nc; 2090 2091 for (c = 0; c < Nc; ++c) {u[c] = 0.0;} 2092 for (b = 0; b < Nb; ++b) { 2093 for (c = 0; c < Nc; ++c) { 2094 const PetscInt cidx = b*Nc+c; 2095 2096 u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx]; 2097 } 2098 } 2099 } 2100 return 0; 2101 } 2102 2103 PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2104 { 2105 const PetscInt dE = T->cdim; /* fegeom->dimEmbed */ 2106 const PetscInt Nq = T->Np; 2107 const PetscInt Nb = T->Nb; 2108 const PetscInt Nc = T->Nc; 2109 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2110 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2111 PetscInt q, b, c, d; 2112 PetscErrorCode ierr; 2113 2114 for (b = 0; b < Nb; ++b) elemVec[b] = 0.0; 2115 for (q = 0; q < Nq; ++q) { 2116 for (b = 0; b < Nb; ++b) { 2117 for (c = 0; c < Nc; ++c) { 2118 const PetscInt bcidx = b*Nc+c; 2119 2120 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2121 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2122 } 2123 } 2124 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2125 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2126 for (b = 0; b < Nb; ++b) { 2127 for (c = 0; c < Nc; ++c) { 2128 const PetscInt bcidx = b*Nc+c; 2129 const PetscInt qcidx = q*Nc+c; 2130 2131 elemVec[b] += tmpBasis[bcidx]*f0[qcidx]; 2132 for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2133 } 2134 } 2135 } 2136 return(0); 2137 } 2138 2139 PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2140 { 2141 const PetscInt dE = T->cdim; 2142 const PetscInt Nq = T->Np; 2143 const PetscInt Nb = T->Nb; 2144 const PetscInt Nc = T->Nc; 2145 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2146 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2147 PetscInt q, b, c, d, s; 2148 PetscErrorCode ierr; 2149 2150 for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0; 2151 for (q = 0; q < Nq; ++q) { 2152 for (b = 0; b < Nb; ++b) { 2153 for (c = 0; c < Nc; ++c) { 2154 const PetscInt bcidx = b*Nc+c; 2155 2156 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2157 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2158 } 2159 } 2160 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2161 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2162 for (s = 0; s < 2; ++s) { 2163 for (b = 0; b < Nb; ++b) { 2164 for (c = 0; c < Nc; ++c) { 2165 const PetscInt bcidx = b*Nc+c; 2166 const PetscInt qcidx = (q*2+s)*Nc+c; 2167 2168 elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx]; 2169 for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2170 } 2171 } 2172 } 2173 } 2174 return(0); 2175 } 2176 2177 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[]) 2178 { 2179 const PetscInt dE = TI->cdim; 2180 const PetscInt NqI = TI->Np; 2181 const PetscInt NbI = TI->Nb; 2182 const PetscInt NcI = TI->Nc; 2183 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2184 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2185 const PetscInt NqJ = TJ->Np; 2186 const PetscInt NbJ = TJ->Nb; 2187 const PetscInt NcJ = TJ->Nc; 2188 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2189 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2190 PetscInt f, fc, g, gc, df, dg; 2191 PetscErrorCode ierr; 2192 2193 for (f = 0; f < NbI; ++f) { 2194 for (fc = 0; fc < NcI; ++fc) { 2195 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2196 2197 tmpBasisI[fidx] = basisI[fidx]; 2198 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2199 } 2200 } 2201 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2202 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2203 for (g = 0; g < NbJ; ++g) { 2204 for (gc = 0; gc < NcJ; ++gc) { 2205 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2206 2207 tmpBasisJ[gidx] = basisJ[gidx]; 2208 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2209 } 2210 } 2211 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2212 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2213 for (f = 0; f < NbI; ++f) { 2214 for (fc = 0; fc < NcI; ++fc) { 2215 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2216 const PetscInt i = offsetI+f; /* Element matrix row */ 2217 for (g = 0; g < NbJ; ++g) { 2218 for (gc = 0; gc < NcJ; ++gc) { 2219 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2220 const PetscInt j = offsetJ+g; /* Element matrix column */ 2221 const PetscInt fOff = eOffset+i*totDim+j; 2222 2223 elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx]; 2224 for (df = 0; df < dE; ++df) { 2225 elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2226 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx]; 2227 for (dg = 0; dg < dE; ++dg) { 2228 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2229 } 2230 } 2231 } 2232 } 2233 } 2234 } 2235 return(0); 2236 } 2237 2238 PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, 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[]) 2239 { 2240 const PetscInt dE = TI->cdim; 2241 const PetscInt NqI = TI->Np; 2242 const PetscInt NbI = TI->Nb; 2243 const PetscInt NcI = TI->Nc; 2244 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2245 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2246 const PetscInt NqJ = TJ->Np; 2247 const PetscInt NbJ = TJ->Nb; 2248 const PetscInt NcJ = TJ->Nc; 2249 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2250 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2251 const PetscInt Ns = isHybridI ? 1 : 2; 2252 const PetscInt Nt = isHybridJ ? 1 : 2; 2253 PetscInt f, fc, g, gc, df, dg, s, t; 2254 PetscErrorCode ierr; 2255 2256 for (f = 0; f < NbI; ++f) { 2257 for (fc = 0; fc < NcI; ++fc) { 2258 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2259 2260 tmpBasisI[fidx] = basisI[fidx]; 2261 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2262 } 2263 } 2264 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2265 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2266 for (g = 0; g < NbJ; ++g) { 2267 for (gc = 0; gc < NcJ; ++gc) { 2268 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2269 2270 tmpBasisJ[gidx] = basisJ[gidx]; 2271 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2272 } 2273 } 2274 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2275 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2276 for (s = 0; s < Ns; ++s) { 2277 for (f = 0; f < NbI; ++f) { 2278 for (fc = 0; fc < NcI; ++fc) { 2279 const PetscInt sc = NcI*s+fc; /* components from each side of the surface */ 2280 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2281 const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */ 2282 for (t = 0; t < Nt; ++t) { 2283 for (g = 0; g < NbJ; ++g) { 2284 for (gc = 0; gc < NcJ; ++gc) { 2285 const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */ 2286 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2287 const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */ 2288 const PetscInt fOff = eOffset+i*totDim+j; 2289 2290 elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx]; 2291 for (df = 0; df < dE; ++df) { 2292 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2293 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx]; 2294 for (dg = 0; dg < dE; ++dg) { 2295 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2296 } 2297 } 2298 } 2299 } 2300 } 2301 } 2302 } 2303 } 2304 return(0); 2305 } 2306 2307 PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom) 2308 { 2309 PetscDualSpace dsp; 2310 DM dm; 2311 PetscQuadrature quadDef; 2312 PetscInt dim, cdim, Nq; 2313 PetscErrorCode ierr; 2314 2315 PetscFunctionBegin; 2316 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); 2317 ierr = PetscDualSpaceGetDM(dsp, &dm);CHKERRQ(ierr); 2318 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2319 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 2320 ierr = PetscFEGetQuadrature(fe, &quadDef);CHKERRQ(ierr); 2321 quad = quad ? quad : quadDef; 2322 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2323 ierr = PetscMalloc1(Nq*cdim, &cgeom->v);CHKERRQ(ierr); 2324 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->J);CHKERRQ(ierr); 2325 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);CHKERRQ(ierr); 2326 ierr = PetscMalloc1(Nq, &cgeom->detJ);CHKERRQ(ierr); 2327 cgeom->dim = dim; 2328 cgeom->dimEmbed = cdim; 2329 cgeom->numCells = 1; 2330 cgeom->numPoints = Nq; 2331 ierr = DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);CHKERRQ(ierr); 2332 PetscFunctionReturn(0); 2333 } 2334 2335 PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom) 2336 { 2337 PetscErrorCode ierr; 2338 2339 PetscFunctionBegin; 2340 ierr = PetscFree(cgeom->v);CHKERRQ(ierr); 2341 ierr = PetscFree(cgeom->J);CHKERRQ(ierr); 2342 ierr = PetscFree(cgeom->invJ);CHKERRQ(ierr); 2343 ierr = PetscFree(cgeom->detJ);CHKERRQ(ierr); 2344 PetscFunctionReturn(0); 2345 } 2346