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