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 - k - The highest derivative we need to tabulate, very often 1 792 793 Output Parameter: 794 . T - The basis function values and derivatives at quadrature points 795 796 Note: 797 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c 798 $ 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 799 $ 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 800 801 Level: intermediate 802 803 .seealso: PetscFECreateTabulation(), PetscTabulationDestroy() 804 @*/ 805 PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscInt k, PetscTabulation *T) 806 { 807 PetscInt npoints; 808 const PetscReal *points; 809 PetscErrorCode ierr; 810 811 PetscFunctionBegin; 812 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 813 PetscValidPointer(T, 2); 814 ierr = PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr); 815 if (!fem->T) {ierr = PetscFECreateTabulation(fem, 1, npoints, points, k, &fem->T);CHKERRQ(ierr);} 816 if (fem->T && k > fem->T->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->T->K); 817 *T = fem->T; 818 PetscFunctionReturn(0); 819 } 820 821 /*@C 822 PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell 823 824 Not collective 825 826 Input Parameter: 827 + fem - The PetscFE object 828 - k - The highest derivative we need to tabulate, very often 1 829 830 Output Parameters: 831 . Tf - The basis function values and derviatives at face quadrature points 832 833 Note: 834 $ 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 835 $ 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 836 $ 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 837 838 Level: intermediate 839 840 .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy() 841 @*/ 842 PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscInt k, PetscTabulation *Tf) 843 { 844 PetscErrorCode ierr; 845 846 PetscFunctionBegin; 847 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 848 PetscValidPointer(Tf, 2); 849 if (!fem->Tf) { 850 const PetscReal xi0[3] = {-1., -1., -1.}; 851 PetscReal v0[3], J[9], detJ; 852 PetscQuadrature fq; 853 PetscDualSpace sp; 854 DM dm; 855 const PetscInt *faces; 856 PetscInt dim, numFaces, f, npoints, q; 857 const PetscReal *points; 858 PetscReal *facePoints; 859 860 ierr = PetscFEGetDualSpace(fem, &sp);CHKERRQ(ierr); 861 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 862 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 863 ierr = DMPlexGetConeSize(dm, 0, &numFaces);CHKERRQ(ierr); 864 ierr = DMPlexGetCone(dm, 0, &faces);CHKERRQ(ierr); 865 ierr = PetscFEGetFaceQuadrature(fem, &fq);CHKERRQ(ierr); 866 if (fq) { 867 ierr = PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr); 868 ierr = PetscMalloc1(numFaces*npoints*dim, &facePoints);CHKERRQ(ierr); 869 for (f = 0; f < numFaces; ++f) { 870 ierr = DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 871 for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]); 872 } 873 ierr = PetscFECreateTabulation(fem, numFaces, npoints, facePoints, k, &fem->Tf);CHKERRQ(ierr); 874 ierr = PetscFree(facePoints);CHKERRQ(ierr); 875 } 876 } 877 if (fem->Tf && k > fem->Tf->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->Tf->K); 878 *Tf = fem->Tf; 879 PetscFunctionReturn(0); 880 } 881 882 /*@C 883 PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points 884 885 Not collective 886 887 Input Parameter: 888 . fem - The PetscFE object 889 890 Output Parameters: 891 . Tc - The basis function values at face centroid points 892 893 Note: 894 $ T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c 895 896 Level: intermediate 897 898 .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy() 899 @*/ 900 PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc) 901 { 902 PetscErrorCode ierr; 903 904 PetscFunctionBegin; 905 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 906 PetscValidPointer(Tc, 2); 907 if (!fem->Tc) { 908 PetscDualSpace sp; 909 DM dm; 910 const PetscInt *cone; 911 PetscReal *centroids; 912 PetscInt dim, numFaces, f; 913 914 ierr = PetscFEGetDualSpace(fem, &sp);CHKERRQ(ierr); 915 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 916 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 917 ierr = DMPlexGetConeSize(dm, 0, &numFaces);CHKERRQ(ierr); 918 ierr = DMPlexGetCone(dm, 0, &cone);CHKERRQ(ierr); 919 ierr = PetscMalloc1(numFaces*dim, ¢roids);CHKERRQ(ierr); 920 for (f = 0; f < numFaces; ++f) {ierr = DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);CHKERRQ(ierr);} 921 ierr = PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);CHKERRQ(ierr); 922 ierr = PetscFree(centroids);CHKERRQ(ierr); 923 } 924 *Tc = fem->Tc; 925 PetscFunctionReturn(0); 926 } 927 928 /*@C 929 PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided. 930 931 Not collective 932 933 Input Parameters: 934 + fem - The PetscFE object 935 . nrepl - The number of replicas 936 . npoints - The number of tabulation points in a replica 937 . points - The tabulation point coordinates 938 - K - The number of derivatives calculated 939 940 Output Parameter: 941 . T - The basis function values and derivatives at tabulation points 942 943 Note: 944 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c 945 $ 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 946 $ 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 947 948 Level: intermediate 949 950 .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy() 951 @*/ 952 PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T) 953 { 954 DM dm; 955 PetscDualSpace Q; 956 PetscInt Nb; /* Dimension of FE space P */ 957 PetscInt Nc; /* Field components */ 958 PetscInt cdim; /* Reference coordinate dimension */ 959 PetscInt k; 960 PetscErrorCode ierr; 961 962 PetscFunctionBegin; 963 if (!npoints || !fem->dualSpace || K < 0) { 964 *T = NULL; 965 PetscFunctionReturn(0); 966 } 967 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 968 PetscValidPointer(points, 4); 969 PetscValidPointer(T, 6); 970 ierr = PetscFEGetDualSpace(fem, &Q);CHKERRQ(ierr); 971 ierr = PetscDualSpaceGetDM(Q, &dm);CHKERRQ(ierr); 972 ierr = DMGetDimension(dm, &cdim);CHKERRQ(ierr); 973 ierr = PetscDualSpaceGetDimension(Q, &Nb);CHKERRQ(ierr); 974 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 975 ierr = PetscMalloc1(1, T);CHKERRQ(ierr); 976 (*T)->K = !cdim ? 0 : K; 977 (*T)->Nr = nrepl; 978 (*T)->Np = npoints; 979 (*T)->Nb = Nb; 980 (*T)->Nc = Nc; 981 (*T)->cdim = cdim; 982 ierr = PetscMalloc1((*T)->K+1, &(*T)->T);CHKERRQ(ierr); 983 for (k = 0; k <= (*T)->K; ++k) { 984 ierr = PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);CHKERRQ(ierr); 985 } 986 ierr = (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);CHKERRQ(ierr); 987 PetscFunctionReturn(0); 988 } 989 990 /*@C 991 PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided. 992 993 Not collective 994 995 Input Parameters: 996 + fem - The PetscFE object 997 . npoints - The number of tabulation points 998 . points - The tabulation point coordinates 999 . K - The number of derivatives calculated 1000 - T - An existing tabulation object with enough allocated space 1001 1002 Output Parameter: 1003 . T - The basis function values and derivatives at tabulation points 1004 1005 Note: 1006 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c 1007 $ 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 1008 $ 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 1009 1010 Level: intermediate 1011 1012 .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy() 1013 @*/ 1014 PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T) 1015 { 1016 PetscErrorCode ierr; 1017 1018 PetscFunctionBeginHot; 1019 if (!npoints || !fem->dualSpace || K < 0) PetscFunctionReturn(0); 1020 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1021 PetscValidPointer(points, 3); 1022 PetscValidPointer(T, 5); 1023 if (PetscDefined(USE_DEBUG)) { 1024 DM dm; 1025 PetscDualSpace Q; 1026 PetscInt Nb; /* Dimension of FE space P */ 1027 PetscInt Nc; /* Field components */ 1028 PetscInt cdim; /* Reference coordinate dimension */ 1029 1030 ierr = PetscFEGetDualSpace(fem, &Q);CHKERRQ(ierr); 1031 ierr = PetscDualSpaceGetDM(Q, &dm);CHKERRQ(ierr); 1032 ierr = DMGetDimension(dm, &cdim);CHKERRQ(ierr); 1033 ierr = PetscDualSpaceGetDimension(Q, &Nb);CHKERRQ(ierr); 1034 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 1035 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); 1036 if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb); 1037 if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc); 1038 if (T->cdim != cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim); 1039 } 1040 T->Nr = 1; 1041 T->Np = npoints; 1042 ierr = (*fem->ops->createtabulation)(fem, npoints, points, K, T);CHKERRQ(ierr); 1043 PetscFunctionReturn(0); 1044 } 1045 1046 /*@C 1047 PetscTabulationDestroy - Frees memory from the associated tabulation. 1048 1049 Not collective 1050 1051 Input Parameter: 1052 . T - The tabulation 1053 1054 Level: intermediate 1055 1056 .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation() 1057 @*/ 1058 PetscErrorCode PetscTabulationDestroy(PetscTabulation *T) 1059 { 1060 PetscInt k; 1061 PetscErrorCode ierr; 1062 1063 PetscFunctionBegin; 1064 PetscValidPointer(T, 1); 1065 if (!T || !(*T)) PetscFunctionReturn(0); 1066 for (k = 0; k <= (*T)->K; ++k) {ierr = PetscFree((*T)->T[k]);CHKERRQ(ierr);} 1067 ierr = PetscFree((*T)->T);CHKERRQ(ierr); 1068 ierr = PetscFree(*T);CHKERRQ(ierr); 1069 *T = NULL; 1070 PetscFunctionReturn(0); 1071 } 1072 1073 PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE) 1074 { 1075 PetscSpace bsp, bsubsp; 1076 PetscDualSpace dsp, dsubsp; 1077 PetscInt dim, depth, numComp, i, j, coneSize, order; 1078 PetscFEType type; 1079 DM dm; 1080 DMLabel label; 1081 PetscReal *xi, *v, *J, detJ; 1082 const char *name; 1083 PetscQuadrature origin, fullQuad, subQuad; 1084 PetscErrorCode ierr; 1085 1086 PetscFunctionBegin; 1087 PetscValidHeaderSpecific(fe,PETSCFE_CLASSID,1); 1088 PetscValidPointer(trFE,3); 1089 ierr = PetscFEGetBasisSpace(fe,&bsp);CHKERRQ(ierr); 1090 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 1091 ierr = PetscDualSpaceGetDM(dsp,&dm);CHKERRQ(ierr); 1092 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 1093 ierr = DMPlexGetDepthLabel(dm,&label);CHKERRQ(ierr); 1094 ierr = DMLabelGetValue(label,refPoint,&depth);CHKERRQ(ierr); 1095 ierr = PetscCalloc1(depth,&xi);CHKERRQ(ierr); 1096 ierr = PetscMalloc1(dim,&v);CHKERRQ(ierr); 1097 ierr = PetscMalloc1(dim*dim,&J);CHKERRQ(ierr); 1098 for (i = 0; i < depth; i++) xi[i] = 0.; 1099 ierr = PetscQuadratureCreate(PETSC_COMM_SELF,&origin);CHKERRQ(ierr); 1100 ierr = PetscQuadratureSetData(origin,depth,0,1,xi,NULL);CHKERRQ(ierr); 1101 ierr = DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);CHKERRQ(ierr); 1102 /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */ 1103 for (i = 1; i < dim; i++) { 1104 for (j = 0; j < depth; j++) { 1105 J[i * depth + j] = J[i * dim + j]; 1106 } 1107 } 1108 ierr = PetscQuadratureDestroy(&origin);CHKERRQ(ierr); 1109 ierr = PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);CHKERRQ(ierr); 1110 ierr = PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);CHKERRQ(ierr); 1111 ierr = PetscSpaceSetUp(bsubsp);CHKERRQ(ierr); 1112 ierr = PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);CHKERRQ(ierr); 1113 ierr = PetscFEGetType(fe,&type);CHKERRQ(ierr); 1114 ierr = PetscFESetType(*trFE,type);CHKERRQ(ierr); 1115 ierr = PetscFEGetNumComponents(fe,&numComp);CHKERRQ(ierr); 1116 ierr = PetscFESetNumComponents(*trFE,numComp);CHKERRQ(ierr); 1117 ierr = PetscFESetBasisSpace(*trFE,bsubsp);CHKERRQ(ierr); 1118 ierr = PetscFESetDualSpace(*trFE,dsubsp);CHKERRQ(ierr); 1119 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 1120 if (name) {ierr = PetscFESetName(*trFE, name);CHKERRQ(ierr);} 1121 ierr = PetscFEGetQuadrature(fe,&fullQuad);CHKERRQ(ierr); 1122 ierr = PetscQuadratureGetOrder(fullQuad,&order);CHKERRQ(ierr); 1123 ierr = DMPlexGetConeSize(dm,refPoint,&coneSize);CHKERRQ(ierr); 1124 if (coneSize == 2 * depth) { 1125 ierr = PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);CHKERRQ(ierr); 1126 } else { 1127 ierr = PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);CHKERRQ(ierr); 1128 } 1129 ierr = PetscFESetQuadrature(*trFE,subQuad);CHKERRQ(ierr); 1130 ierr = PetscFESetUp(*trFE);CHKERRQ(ierr); 1131 ierr = PetscQuadratureDestroy(&subQuad);CHKERRQ(ierr); 1132 ierr = PetscSpaceDestroy(&bsubsp);CHKERRQ(ierr); 1133 PetscFunctionReturn(0); 1134 } 1135 1136 PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE) 1137 { 1138 PetscInt hStart, hEnd; 1139 PetscDualSpace dsp; 1140 DM dm; 1141 PetscErrorCode ierr; 1142 1143 PetscFunctionBegin; 1144 PetscValidHeaderSpecific(fe,PETSCFE_CLASSID,1); 1145 PetscValidPointer(trFE,3); 1146 *trFE = NULL; 1147 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 1148 ierr = PetscDualSpaceGetDM(dsp,&dm);CHKERRQ(ierr); 1149 ierr = DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);CHKERRQ(ierr); 1150 if (hEnd <= hStart) PetscFunctionReturn(0); 1151 ierr = PetscFECreatePointTrace(fe,hStart,trFE);CHKERRQ(ierr); 1152 PetscFunctionReturn(0); 1153 } 1154 1155 1156 /*@ 1157 PetscFEGetDimension - Get the dimension of the finite element space on a cell 1158 1159 Not collective 1160 1161 Input Parameter: 1162 . fe - The PetscFE 1163 1164 Output Parameter: 1165 . dim - The dimension 1166 1167 Level: intermediate 1168 1169 .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension() 1170 @*/ 1171 PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) 1172 { 1173 PetscErrorCode ierr; 1174 1175 PetscFunctionBegin; 1176 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1177 PetscValidPointer(dim, 2); 1178 if (fem->ops->getdimension) {ierr = (*fem->ops->getdimension)(fem, dim);CHKERRQ(ierr);} 1179 PetscFunctionReturn(0); 1180 } 1181 1182 /*@C 1183 PetscFEPushforward - Map the reference element function to real space 1184 1185 Input Parameters: 1186 + fe - The PetscFE 1187 . fegeom - The cell geometry 1188 . Nv - The number of function values 1189 - vals - The function values 1190 1191 Output Parameter: 1192 . vals - The transformed function values 1193 1194 Level: advanced 1195 1196 Note: This just forwards the call onto PetscDualSpacePushforward(). 1197 1198 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1199 1200 .seealso: PetscDualSpacePushforward() 1201 @*/ 1202 PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1203 { 1204 PetscErrorCode ierr; 1205 1206 PetscFunctionBeginHot; 1207 ierr = PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1208 PetscFunctionReturn(0); 1209 } 1210 1211 /*@C 1212 PetscFEPushforwardGradient - Map the reference element function gradient to real space 1213 1214 Input Parameters: 1215 + fe - The PetscFE 1216 . fegeom - The cell geometry 1217 . Nv - The number of function gradient values 1218 - vals - The function gradient values 1219 1220 Output Parameter: 1221 . vals - The transformed function gradient values 1222 1223 Level: advanced 1224 1225 Note: This just forwards the call onto PetscDualSpacePushforwardGradient(). 1226 1227 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1228 1229 .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward() 1230 @*/ 1231 PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1232 { 1233 PetscErrorCode ierr; 1234 1235 PetscFunctionBeginHot; 1236 ierr = PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1237 PetscFunctionReturn(0); 1238 } 1239 1240 /*@C 1241 PetscFEPushforwardHessian - Map the reference element function Hessian to real space 1242 1243 Input Parameters: 1244 + fe - The PetscFE 1245 . fegeom - The cell geometry 1246 . Nv - The number of function Hessian values 1247 - vals - The function Hessian values 1248 1249 Output Parameter: 1250 . vals - The transformed function Hessian values 1251 1252 Level: advanced 1253 1254 Note: This just forwards the call onto PetscDualSpacePushforwardHessian(). 1255 1256 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1257 1258 .seealso: PetscFEPushforward(), PetscDualSpacePushforwardHessian(), PetscDualSpacePushforward() 1259 @*/ 1260 PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1261 { 1262 PetscErrorCode ierr; 1263 1264 PetscFunctionBeginHot; 1265 ierr = PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1266 PetscFunctionReturn(0); 1267 } 1268 1269 /* 1270 Purpose: Compute element vector for chunk of elements 1271 1272 Input: 1273 Sizes: 1274 Ne: number of elements 1275 Nf: number of fields 1276 PetscFE 1277 dim: spatial dimension 1278 Nb: number of basis functions 1279 Nc: number of field components 1280 PetscQuadrature 1281 Nq: number of quadrature points 1282 1283 Geometry: 1284 PetscFEGeom[Ne] possibly *Nq 1285 PetscReal v0s[dim] 1286 PetscReal n[dim] 1287 PetscReal jacobians[dim*dim] 1288 PetscReal jacobianInverses[dim*dim] 1289 PetscReal jacobianDeterminants 1290 FEM: 1291 PetscFE 1292 PetscQuadrature 1293 PetscReal quadPoints[Nq*dim] 1294 PetscReal quadWeights[Nq] 1295 PetscReal basis[Nq*Nb*Nc] 1296 PetscReal basisDer[Nq*Nb*Nc*dim] 1297 PetscScalar coefficients[Ne*Nb*Nc] 1298 PetscScalar elemVec[Ne*Nb*Nc] 1299 1300 Problem: 1301 PetscInt f: the active field 1302 f0, f1 1303 1304 Work Space: 1305 PetscFE 1306 PetscScalar f0[Nq*dim]; 1307 PetscScalar f1[Nq*dim*dim]; 1308 PetscScalar u[Nc]; 1309 PetscScalar gradU[Nc*dim]; 1310 PetscReal x[dim]; 1311 PetscScalar realSpaceDer[dim]; 1312 1313 Purpose: Compute element vector for N_cb batches of elements 1314 1315 Input: 1316 Sizes: 1317 N_cb: Number of serial cell batches 1318 1319 Geometry: 1320 PetscReal v0s[Ne*dim] 1321 PetscReal jacobians[Ne*dim*dim] possibly *Nq 1322 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 1323 PetscReal jacobianDeterminants[Ne] possibly *Nq 1324 FEM: 1325 static PetscReal quadPoints[Nq*dim] 1326 static PetscReal quadWeights[Nq] 1327 static PetscReal basis[Nq*Nb*Nc] 1328 static PetscReal basisDer[Nq*Nb*Nc*dim] 1329 PetscScalar coefficients[Ne*Nb*Nc] 1330 PetscScalar elemVec[Ne*Nb*Nc] 1331 1332 ex62.c: 1333 PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], 1334 const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], 1335 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), 1336 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) 1337 1338 ex52.c: 1339 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) 1340 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) 1341 1342 ex52_integrateElement.cu 1343 __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) 1344 1345 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], 1346 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1347 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1348 1349 ex52_integrateElementOpenCL.c: 1350 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], 1351 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1352 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1353 1354 __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) 1355 */ 1356 1357 /*@C 1358 PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration 1359 1360 Not collective 1361 1362 Input Parameters: 1363 + prob - The PetscDS specifying the discretizations and continuum functions 1364 . field - The field being integrated 1365 . Ne - The number of elements in the chunk 1366 . cgeom - The cell geometry for each cell in the chunk 1367 . coefficients - The array of FEM basis coefficients for the elements 1368 . probAux - The PetscDS specifying the auxiliary discretizations 1369 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1370 1371 Output Parameter: 1372 . integral - the integral for this field 1373 1374 Level: intermediate 1375 1376 .seealso: PetscFEIntegrateResidual() 1377 @*/ 1378 PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 1379 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1380 { 1381 PetscFE fe; 1382 PetscErrorCode ierr; 1383 1384 PetscFunctionBegin; 1385 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1386 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1387 if (fe->ops->integrate) {ierr = (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} 1388 PetscFunctionReturn(0); 1389 } 1390 1391 /*@C 1392 PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration 1393 1394 Not collective 1395 1396 Input Parameters: 1397 + prob - The PetscDS specifying the discretizations and continuum functions 1398 . field - The field being integrated 1399 . obj_func - The function to be integrated 1400 . Ne - The number of elements in the chunk 1401 . fgeom - The face geometry for each face in the chunk 1402 . coefficients - The array of FEM basis coefficients for the elements 1403 . probAux - The PetscDS specifying the auxiliary discretizations 1404 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1405 1406 Output Parameter: 1407 . integral - the integral for this field 1408 1409 Level: intermediate 1410 1411 .seealso: PetscFEIntegrateResidual() 1412 @*/ 1413 PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field, 1414 void (*obj_func)(PetscInt, PetscInt, PetscInt, 1415 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1416 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1417 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 1418 PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1419 { 1420 PetscFE fe; 1421 PetscErrorCode ierr; 1422 1423 PetscFunctionBegin; 1424 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1425 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1426 if (fe->ops->integratebd) {ierr = (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} 1427 PetscFunctionReturn(0); 1428 } 1429 1430 /*@C 1431 PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration 1432 1433 Not collective 1434 1435 Input Parameters: 1436 + ds - The PetscDS specifying the discretizations and continuum functions 1437 . key - The (label+value, field) being integrated 1438 . Ne - The number of elements in the chunk 1439 . cgeom - The cell geometry for each cell in the chunk 1440 . coefficients - The array of FEM basis coefficients for the elements 1441 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1442 . probAux - The PetscDS specifying the auxiliary discretizations 1443 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1444 - t - The time 1445 1446 Output Parameter: 1447 . elemVec - the element residual vectors from each element 1448 1449 Note: 1450 $ Loop over batch of elements (e): 1451 $ Loop over quadrature points (q): 1452 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q 1453 $ Call f_0 and f_1 1454 $ Loop over element vector entries (f,fc --> i): 1455 $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) 1456 1457 Level: intermediate 1458 1459 .seealso: PetscFEIntegrateResidual() 1460 @*/ 1461 PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *cgeom, 1462 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1463 { 1464 PetscFE fe; 1465 PetscErrorCode ierr; 1466 1467 PetscFunctionBeginHot; 1468 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1469 ierr = PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);CHKERRQ(ierr); 1470 if (fe->ops->integrateresidual) {ierr = (*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1471 PetscFunctionReturn(0); 1472 } 1473 1474 /*@C 1475 PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary 1476 1477 Not collective 1478 1479 Input Parameters: 1480 + ds - The PetscDS specifying the discretizations and continuum functions 1481 . wf - The PetscWeakForm object holding the pointwise functions 1482 . key - The (label+value, field) being integrated 1483 . Ne - The number of elements in the chunk 1484 . fgeom - The face geometry for each cell in the chunk 1485 . coefficients - The array of FEM basis coefficients for the elements 1486 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1487 . probAux - The PetscDS specifying the auxiliary discretizations 1488 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1489 - t - The time 1490 1491 Output Parameter: 1492 . elemVec - the element residual vectors from each element 1493 1494 Level: intermediate 1495 1496 .seealso: PetscFEIntegrateResidual() 1497 @*/ 1498 PetscErrorCode PetscFEIntegrateBdResidual(PetscDS ds, PetscWeakForm wf, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1499 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1500 { 1501 PetscFE fe; 1502 PetscErrorCode ierr; 1503 1504 PetscFunctionBegin; 1505 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1506 ierr = PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);CHKERRQ(ierr); 1507 if (fe->ops->integratebdresidual) {ierr = (*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1508 PetscFunctionReturn(0); 1509 } 1510 1511 /*@C 1512 PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration 1513 1514 Not collective 1515 1516 Input Parameters: 1517 + prob - The PetscDS specifying the discretizations and continuum functions 1518 . key - The (label+value, field) being integrated 1519 . Ne - The number of elements in the chunk 1520 . fgeom - The face geometry for each cell in the chunk 1521 . coefficients - The array of FEM basis coefficients for the elements 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 1527 Output Parameter 1528 . elemVec - the element residual vectors from each element 1529 1530 Level: developer 1531 1532 .seealso: PetscFEIntegrateResidual() 1533 @*/ 1534 PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1535 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1536 { 1537 PetscFE fe; 1538 PetscErrorCode ierr; 1539 1540 PetscFunctionBegin; 1541 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1542 ierr = PetscDSGetDiscretization(prob, key.field, (PetscObject *) &fe);CHKERRQ(ierr); 1543 if (fe->ops->integratehybridresidual) {ierr = (*fe->ops->integratehybridresidual)(prob, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1544 PetscFunctionReturn(0); 1545 } 1546 1547 /*@C 1548 PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration 1549 1550 Not collective 1551 1552 Input Parameters: 1553 + ds - The PetscDS specifying the discretizations and continuum functions 1554 . jtype - The type of matrix pointwise functions that should be used 1555 . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1556 . Ne - The number of elements in the chunk 1557 . cgeom - The cell geometry for each cell in the chunk 1558 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1559 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1560 . probAux - The PetscDS specifying the auxiliary discretizations 1561 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1562 . t - The time 1563 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1564 1565 Output Parameter: 1566 . elemMat - the element matrices for the Jacobian from each element 1567 1568 Note: 1569 $ Loop over batch of elements (e): 1570 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1571 $ Loop over quadrature points (q): 1572 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1573 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1574 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1575 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1576 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1577 Level: intermediate 1578 1579 .seealso: PetscFEIntegrateResidual() 1580 @*/ 1581 PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *cgeom, 1582 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1583 { 1584 PetscFE fe; 1585 PetscInt Nf; 1586 PetscErrorCode ierr; 1587 1588 PetscFunctionBegin; 1589 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1590 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1591 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1592 if (fe->ops->integratejacobian) {ierr = (*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1593 PetscFunctionReturn(0); 1594 } 1595 1596 /*@C 1597 PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration 1598 1599 Not collective 1600 1601 Input Parameters: 1602 + ds - The PetscDS specifying the discretizations and continuum functions 1603 . wf - The PetscWeakForm holding the pointwise functions 1604 . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1605 . Ne - The number of elements in the chunk 1606 . fgeom - The face geometry for each cell in the chunk 1607 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1608 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1609 . probAux - The PetscDS specifying the auxiliary discretizations 1610 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1611 . t - The time 1612 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1613 1614 Output Parameter: 1615 . elemMat - the element matrices for the Jacobian from each element 1616 1617 Note: 1618 $ Loop over batch of elements (e): 1619 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1620 $ Loop over quadrature points (q): 1621 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1622 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1623 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1624 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1625 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1626 Level: intermediate 1627 1628 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1629 @*/ 1630 PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1631 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1632 { 1633 PetscFE fe; 1634 PetscInt Nf; 1635 PetscErrorCode ierr; 1636 1637 PetscFunctionBegin; 1638 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1639 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1640 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1641 if (fe->ops->integratebdjacobian) {ierr = (*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1642 PetscFunctionReturn(0); 1643 } 1644 1645 /*@C 1646 PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration 1647 1648 Not collective 1649 1650 Input Parameters: 1651 + ds - The PetscDS specifying the discretizations and continuum functions 1652 . jtype - The type of matrix pointwise functions that should be used 1653 . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1654 . Ne - The number of elements in the chunk 1655 . fgeom - The face geometry for each cell in the chunk 1656 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1657 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1658 . probAux - The PetscDS specifying the auxiliary discretizations 1659 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1660 . t - The time 1661 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1662 1663 Output Parameter 1664 . elemMat - the element matrices for the Jacobian from each element 1665 1666 Note: 1667 $ Loop over batch of elements (e): 1668 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1669 $ Loop over quadrature points (q): 1670 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1671 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1672 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1673 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1674 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1675 Level: developer 1676 1677 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1678 @*/ 1679 PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1680 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1681 { 1682 PetscFE fe; 1683 PetscInt Nf; 1684 PetscErrorCode ierr; 1685 1686 PetscFunctionBegin; 1687 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1688 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1689 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1690 if (fe->ops->integratehybridjacobian) {ierr = (*fe->ops->integratehybridjacobian)(ds, jtype, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1691 PetscFunctionReturn(0); 1692 } 1693 1694 /*@ 1695 PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height 1696 1697 Input Parameters: 1698 + fe - The finite element space 1699 - height - The height of the Plex point 1700 1701 Output Parameter: 1702 . subfe - The subspace of this FE space 1703 1704 Note: For example, if we want the subspace of this space for a face, we would choose height = 1. 1705 1706 Level: advanced 1707 1708 .seealso: PetscFECreateDefault() 1709 @*/ 1710 PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe) 1711 { 1712 PetscSpace P, subP; 1713 PetscDualSpace Q, subQ; 1714 PetscQuadrature subq; 1715 PetscFEType fetype; 1716 PetscInt dim, Nc; 1717 PetscErrorCode ierr; 1718 1719 PetscFunctionBegin; 1720 PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1); 1721 PetscValidPointer(subfe, 3); 1722 if (height == 0) { 1723 *subfe = fe; 1724 PetscFunctionReturn(0); 1725 } 1726 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1727 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1728 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1729 ierr = PetscFEGetFaceQuadrature(fe, &subq);CHKERRQ(ierr); 1730 ierr = PetscDualSpaceGetDimension(Q, &dim);CHKERRQ(ierr); 1731 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); 1732 if (!fe->subspaces) {ierr = PetscCalloc1(dim, &fe->subspaces);CHKERRQ(ierr);} 1733 if (height <= dim) { 1734 if (!fe->subspaces[height-1]) { 1735 PetscFE sub = NULL; 1736 const char *name; 1737 1738 ierr = PetscSpaceGetHeightSubspace(P, height, &subP);CHKERRQ(ierr); 1739 ierr = PetscDualSpaceGetHeightSubspace(Q, height, &subQ);CHKERRQ(ierr); 1740 if (subQ) { 1741 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);CHKERRQ(ierr); 1742 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 1743 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 1744 ierr = PetscFEGetType(fe, &fetype);CHKERRQ(ierr); 1745 ierr = PetscFESetType(sub, fetype);CHKERRQ(ierr); 1746 ierr = PetscFESetBasisSpace(sub, subP);CHKERRQ(ierr); 1747 ierr = PetscFESetDualSpace(sub, subQ);CHKERRQ(ierr); 1748 ierr = PetscFESetNumComponents(sub, Nc);CHKERRQ(ierr); 1749 ierr = PetscFESetUp(sub);CHKERRQ(ierr); 1750 ierr = PetscFESetQuadrature(sub, subq);CHKERRQ(ierr); 1751 } 1752 fe->subspaces[height-1] = sub; 1753 } 1754 *subfe = fe->subspaces[height-1]; 1755 } else { 1756 *subfe = NULL; 1757 } 1758 PetscFunctionReturn(0); 1759 } 1760 1761 /*@ 1762 PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used 1763 to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more 1764 sparsity). It is also used to create an interpolation between regularly refined meshes. 1765 1766 Collective on fem 1767 1768 Input Parameter: 1769 . fe - The initial PetscFE 1770 1771 Output Parameter: 1772 . feRef - The refined PetscFE 1773 1774 Level: advanced 1775 1776 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 1777 @*/ 1778 PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef) 1779 { 1780 PetscSpace P, Pref; 1781 PetscDualSpace Q, Qref; 1782 DM K, Kref; 1783 PetscQuadrature q, qref; 1784 const PetscReal *v0, *jac; 1785 PetscInt numComp, numSubelements; 1786 PetscInt cStart, cEnd, c; 1787 PetscDualSpace *cellSpaces; 1788 PetscErrorCode ierr; 1789 1790 PetscFunctionBegin; 1791 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1792 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1793 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1794 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr); 1795 /* Create space */ 1796 ierr = PetscObjectReference((PetscObject) P);CHKERRQ(ierr); 1797 Pref = P; 1798 /* Create dual space */ 1799 ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr); 1800 ierr = PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);CHKERRQ(ierr); 1801 ierr = DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);CHKERRQ(ierr); 1802 ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr); 1803 ierr = DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);CHKERRQ(ierr); 1804 ierr = PetscMalloc1(cEnd - cStart, &cellSpaces);CHKERRQ(ierr); 1805 /* TODO: fix for non-uniform refinement */ 1806 for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q; 1807 ierr = PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);CHKERRQ(ierr); 1808 ierr = PetscFree(cellSpaces);CHKERRQ(ierr); 1809 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 1810 ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr); 1811 /* Create element */ 1812 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);CHKERRQ(ierr); 1813 ierr = PetscFESetType(*feRef, PETSCFECOMPOSITE);CHKERRQ(ierr); 1814 ierr = PetscFESetBasisSpace(*feRef, Pref);CHKERRQ(ierr); 1815 ierr = PetscFESetDualSpace(*feRef, Qref);CHKERRQ(ierr); 1816 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 1817 ierr = PetscFESetNumComponents(*feRef, numComp);CHKERRQ(ierr); 1818 ierr = PetscFESetUp(*feRef);CHKERRQ(ierr); 1819 ierr = PetscSpaceDestroy(&Pref);CHKERRQ(ierr); 1820 ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr); 1821 /* Create quadrature */ 1822 ierr = PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);CHKERRQ(ierr); 1823 ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr); 1824 ierr = PetscFESetQuadrature(*feRef, qref);CHKERRQ(ierr); 1825 ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr); 1826 PetscFunctionReturn(0); 1827 } 1828 1829 /*@C 1830 PetscFECreateDefault - Create a PetscFE for basic FEM computation 1831 1832 Collective 1833 1834 Input Parameters: 1835 + comm - The MPI comm 1836 . dim - The spatial dimension 1837 . Nc - The number of components 1838 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1839 . prefix - The options prefix, or NULL 1840 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1841 1842 Output Parameter: 1843 . fem - The PetscFE object 1844 1845 Note: 1846 Each subobject is SetFromOption() during creation, so that the object may be customized from the command line, using the prefix specified above. See the links below for the particular options available. 1847 1848 Level: beginner 1849 1850 .seealso: PetscSpaceSetFromOptions(), PetscDualSpaceSetFromOptions(), PetscFESetFromOptions(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1851 @*/ 1852 PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem) 1853 { 1854 PetscQuadrature q, fq; 1855 DM K; 1856 PetscSpace P; 1857 PetscDualSpace Q; 1858 PetscInt order, quadPointsPerEdge; 1859 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1860 PetscErrorCode ierr; 1861 1862 PetscFunctionBegin; 1863 /* Create space */ 1864 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1865 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr); 1866 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1867 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1868 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1869 ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr); 1870 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1871 ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr); 1872 ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr); 1873 /* Create dual space */ 1874 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1875 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1876 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr); 1877 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1878 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1879 ierr = DMDestroy(&K);CHKERRQ(ierr); 1880 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1881 ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr); 1882 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1883 ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr); 1884 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1885 /* Create element */ 1886 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1887 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr); 1888 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1889 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1890 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1891 ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr); 1892 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1893 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1894 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1895 /* Create quadrature (with specified order if given) */ 1896 qorder = qorder >= 0 ? qorder : order; 1897 ierr = PetscObjectOptionsBegin((PetscObject)*fem);CHKERRQ(ierr); 1898 ierr = PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);CHKERRQ(ierr); 1899 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1900 quadPointsPerEdge = PetscMax(qorder + 1,1); 1901 if (isSimplex) { 1902 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1903 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1904 } else { 1905 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1906 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1907 } 1908 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1909 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1910 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1911 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1912 PetscFunctionReturn(0); 1913 } 1914 1915 /*@ 1916 PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k 1917 1918 Collective 1919 1920 Input Parameters: 1921 + comm - The MPI comm 1922 . dim - The spatial dimension 1923 . Nc - The number of components 1924 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1925 . k - The degree k of the space 1926 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1927 1928 Output Parameter: 1929 . fem - The PetscFE object 1930 1931 Level: beginner 1932 1933 Notes: 1934 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. 1935 1936 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1937 @*/ 1938 PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem) 1939 { 1940 PetscQuadrature q, fq; 1941 DM K; 1942 PetscSpace P; 1943 PetscDualSpace Q; 1944 PetscInt quadPointsPerEdge; 1945 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1946 char name[64]; 1947 PetscErrorCode ierr; 1948 1949 PetscFunctionBegin; 1950 /* Create space */ 1951 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1952 ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1953 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1954 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1955 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1956 ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr); 1957 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1958 /* Create dual space */ 1959 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1960 ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1961 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1962 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1963 ierr = DMDestroy(&K);CHKERRQ(ierr); 1964 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1965 ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr); 1966 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1967 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1968 /* Create finite element */ 1969 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1970 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1971 ierr = PetscObjectSetName((PetscObject) *fem, name);CHKERRQ(ierr); 1972 ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr); 1973 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1974 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1975 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1976 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1977 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1978 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1979 /* Create quadrature (with specified order if given) */ 1980 qorder = qorder >= 0 ? qorder : k; 1981 quadPointsPerEdge = PetscMax(qorder + 1,1); 1982 if (isSimplex) { 1983 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1984 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1985 } else { 1986 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1987 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1988 } 1989 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1990 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1991 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1992 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1993 /* Set finite element name */ 1994 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1995 ierr = PetscFESetName(*fem, name);CHKERRQ(ierr); 1996 PetscFunctionReturn(0); 1997 } 1998 1999 /*@C 2000 PetscFESetName - Names the FE and its subobjects 2001 2002 Not collective 2003 2004 Input Parameters: 2005 + fe - The PetscFE 2006 - name - The name 2007 2008 Level: intermediate 2009 2010 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 2011 @*/ 2012 PetscErrorCode PetscFESetName(PetscFE fe, const char name[]) 2013 { 2014 PetscSpace P; 2015 PetscDualSpace Q; 2016 PetscErrorCode ierr; 2017 2018 PetscFunctionBegin; 2019 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 2020 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2021 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 2022 ierr = PetscObjectSetName((PetscObject) P, name);CHKERRQ(ierr); 2023 ierr = PetscObjectSetName((PetscObject) Q, name);CHKERRQ(ierr); 2024 PetscFunctionReturn(0); 2025 } 2026 2027 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[]) 2028 { 2029 PetscInt dOffset = 0, fOffset = 0, f, g; 2030 PetscErrorCode ierr; 2031 2032 for (f = 0; f < Nf; ++f) { 2033 PetscFE fe; 2034 const PetscInt k = ds->jetDegree[f]; 2035 const PetscInt cdim = T[f]->cdim; 2036 const PetscInt Nq = T[f]->Np; 2037 const PetscInt Nbf = T[f]->Nb; 2038 const PetscInt Ncf = T[f]->Nc; 2039 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2040 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2041 const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r*Nq+q)*Nbf*Ncf*cdim*cdim] : NULL; 2042 PetscInt hOffset = 0, b, c, d; 2043 2044 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 2045 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2046 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2047 for (b = 0; b < Nbf; ++b) { 2048 for (c = 0; c < Ncf; ++c) { 2049 const PetscInt cidx = b*Ncf+c; 2050 2051 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2052 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2053 } 2054 } 2055 if (k > 1) { 2056 for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc*cdim; 2057 for (d = 0; d < cdim*cdim*Ncf; ++d) u_x[hOffset+fOffset*cdim*cdim+d] = 0.0; 2058 for (b = 0; b < Nbf; ++b) { 2059 for (c = 0; c < Ncf; ++c) { 2060 const PetscInt cidx = b*Ncf+c; 2061 2062 for (d = 0; d < cdim*cdim; ++d) u_x[hOffset+(fOffset+c)*cdim*cdim+d] += Hq[cidx*cdim*cdim+d]*coefficients[dOffset+b]; 2063 } 2064 } 2065 ierr = PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset+fOffset*cdim*cdim]);CHKERRQ(ierr); 2066 } 2067 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2068 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2069 if (u_t) { 2070 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2071 for (b = 0; b < Nbf; ++b) { 2072 for (c = 0; c < Ncf; ++c) { 2073 const PetscInt cidx = b*Ncf+c; 2074 2075 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2076 } 2077 } 2078 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2079 } 2080 fOffset += Ncf; 2081 dOffset += Nbf; 2082 } 2083 return 0; 2084 } 2085 2086 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[]) 2087 { 2088 PetscInt dOffset = 0, fOffset = 0, g; 2089 PetscErrorCode ierr; 2090 2091 for (g = 0; g < 2*Nf-1; ++g) { 2092 if (!T[g/2]) continue; 2093 { 2094 PetscFE fe; 2095 const PetscInt f = g/2; 2096 const PetscInt cdim = T[f]->cdim; 2097 const PetscInt Nq = T[f]->Np; 2098 const PetscInt Nbf = T[f]->Nb; 2099 const PetscInt Ncf = T[f]->Nc; 2100 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2101 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2102 PetscInt b, c, d; 2103 2104 fe = (PetscFE) ds->disc[f]; 2105 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2106 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2107 for (b = 0; b < Nbf; ++b) { 2108 for (c = 0; c < Ncf; ++c) { 2109 const PetscInt cidx = b*Ncf+c; 2110 2111 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2112 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2113 } 2114 } 2115 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2116 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2117 if (u_t) { 2118 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2119 for (b = 0; b < Nbf; ++b) { 2120 for (c = 0; c < Ncf; ++c) { 2121 const PetscInt cidx = b*Ncf+c; 2122 2123 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2124 } 2125 } 2126 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2127 } 2128 fOffset += Ncf; 2129 dOffset += Nbf; 2130 } 2131 } 2132 return 0; 2133 } 2134 2135 PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[]) 2136 { 2137 PetscFE fe; 2138 PetscTabulation Tc; 2139 PetscInt b, c; 2140 PetscErrorCode ierr; 2141 2142 if (!prob) return 0; 2143 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 2144 ierr = PetscFEGetFaceCentroidTabulation(fe, &Tc);CHKERRQ(ierr); 2145 { 2146 const PetscReal *faceBasis = Tc->T[0]; 2147 const PetscInt Nb = Tc->Nb; 2148 const PetscInt Nc = Tc->Nc; 2149 2150 for (c = 0; c < Nc; ++c) {u[c] = 0.0;} 2151 for (b = 0; b < Nb; ++b) { 2152 for (c = 0; c < Nc; ++c) { 2153 u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c]; 2154 } 2155 } 2156 } 2157 return 0; 2158 } 2159 2160 PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2161 { 2162 const PetscInt dE = T->cdim; /* fegeom->dimEmbed */ 2163 const PetscInt Nq = T->Np; 2164 const PetscInt Nb = T->Nb; 2165 const PetscInt Nc = T->Nc; 2166 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2167 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2168 PetscInt q, b, c, d; 2169 PetscErrorCode ierr; 2170 2171 for (b = 0; b < Nb; ++b) elemVec[b] = 0.0; 2172 for (q = 0; q < Nq; ++q) { 2173 for (b = 0; b < Nb; ++b) { 2174 for (c = 0; c < Nc; ++c) { 2175 const PetscInt bcidx = b*Nc+c; 2176 2177 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2178 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2179 } 2180 } 2181 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2182 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2183 for (b = 0; b < Nb; ++b) { 2184 for (c = 0; c < Nc; ++c) { 2185 const PetscInt bcidx = b*Nc+c; 2186 const PetscInt qcidx = q*Nc+c; 2187 2188 elemVec[b] += tmpBasis[bcidx]*f0[qcidx]; 2189 for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2190 } 2191 } 2192 } 2193 return(0); 2194 } 2195 2196 PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2197 { 2198 const PetscInt dE = T->cdim; 2199 const PetscInt Nq = T->Np; 2200 const PetscInt Nb = T->Nb; 2201 const PetscInt Nc = T->Nc; 2202 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2203 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2204 PetscInt q, b, c, d, s; 2205 PetscErrorCode ierr; 2206 2207 for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0; 2208 for (q = 0; q < Nq; ++q) { 2209 for (b = 0; b < Nb; ++b) { 2210 for (c = 0; c < Nc; ++c) { 2211 const PetscInt bcidx = b*Nc+c; 2212 2213 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2214 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2215 } 2216 } 2217 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2218 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2219 for (s = 0; s < 2; ++s) { 2220 for (b = 0; b < Nb; ++b) { 2221 for (c = 0; c < Nc; ++c) { 2222 const PetscInt bcidx = b*Nc+c; 2223 const PetscInt qcidx = (q*2+s)*Nc+c; 2224 2225 elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx]; 2226 for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2227 } 2228 } 2229 } 2230 } 2231 return(0); 2232 } 2233 2234 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[]) 2235 { 2236 const PetscInt dE = TI->cdim; 2237 const PetscInt NqI = TI->Np; 2238 const PetscInt NbI = TI->Nb; 2239 const PetscInt NcI = TI->Nc; 2240 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2241 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2242 const PetscInt NqJ = TJ->Np; 2243 const PetscInt NbJ = TJ->Nb; 2244 const PetscInt NcJ = TJ->Nc; 2245 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2246 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2247 PetscInt f, fc, g, gc, df, dg; 2248 PetscErrorCode ierr; 2249 2250 for (f = 0; f < NbI; ++f) { 2251 for (fc = 0; fc < NcI; ++fc) { 2252 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2253 2254 tmpBasisI[fidx] = basisI[fidx]; 2255 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2256 } 2257 } 2258 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2259 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2260 for (g = 0; g < NbJ; ++g) { 2261 for (gc = 0; gc < NcJ; ++gc) { 2262 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2263 2264 tmpBasisJ[gidx] = basisJ[gidx]; 2265 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2266 } 2267 } 2268 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2269 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2270 for (f = 0; f < NbI; ++f) { 2271 for (fc = 0; fc < NcI; ++fc) { 2272 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2273 const PetscInt i = offsetI+f; /* Element matrix row */ 2274 for (g = 0; g < NbJ; ++g) { 2275 for (gc = 0; gc < NcJ; ++gc) { 2276 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2277 const PetscInt j = offsetJ+g; /* Element matrix column */ 2278 const PetscInt fOff = eOffset+i*totDim+j; 2279 2280 elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx]; 2281 for (df = 0; df < dE; ++df) { 2282 elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2283 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx]; 2284 for (dg = 0; dg < dE; ++dg) { 2285 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2286 } 2287 } 2288 } 2289 } 2290 } 2291 } 2292 return(0); 2293 } 2294 2295 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[]) 2296 { 2297 const PetscInt dE = TI->cdim; 2298 const PetscInt NqI = TI->Np; 2299 const PetscInt NbI = TI->Nb; 2300 const PetscInt NcI = TI->Nc; 2301 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2302 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2303 const PetscInt NqJ = TJ->Np; 2304 const PetscInt NbJ = TJ->Nb; 2305 const PetscInt NcJ = TJ->Nc; 2306 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2307 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2308 const PetscInt Ns = isHybridI ? 1 : 2; 2309 const PetscInt Nt = isHybridJ ? 1 : 2; 2310 PetscInt f, fc, g, gc, df, dg, s, t; 2311 PetscErrorCode ierr; 2312 2313 for (f = 0; f < NbI; ++f) { 2314 for (fc = 0; fc < NcI; ++fc) { 2315 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2316 2317 tmpBasisI[fidx] = basisI[fidx]; 2318 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2319 } 2320 } 2321 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2322 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2323 for (g = 0; g < NbJ; ++g) { 2324 for (gc = 0; gc < NcJ; ++gc) { 2325 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2326 2327 tmpBasisJ[gidx] = basisJ[gidx]; 2328 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2329 } 2330 } 2331 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2332 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2333 for (s = 0; s < Ns; ++s) { 2334 for (f = 0; f < NbI; ++f) { 2335 for (fc = 0; fc < NcI; ++fc) { 2336 const PetscInt sc = NcI*s+fc; /* components from each side of the surface */ 2337 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2338 const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */ 2339 for (t = 0; t < Nt; ++t) { 2340 for (g = 0; g < NbJ; ++g) { 2341 for (gc = 0; gc < NcJ; ++gc) { 2342 const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */ 2343 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2344 const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */ 2345 const PetscInt fOff = eOffset+i*totDim+j; 2346 2347 elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx]; 2348 for (df = 0; df < dE; ++df) { 2349 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2350 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx]; 2351 for (dg = 0; dg < dE; ++dg) { 2352 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2353 } 2354 } 2355 } 2356 } 2357 } 2358 } 2359 } 2360 } 2361 return(0); 2362 } 2363 2364 PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom) 2365 { 2366 PetscDualSpace dsp; 2367 DM dm; 2368 PetscQuadrature quadDef; 2369 PetscInt dim, cdim, Nq; 2370 PetscErrorCode ierr; 2371 2372 PetscFunctionBegin; 2373 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); 2374 ierr = PetscDualSpaceGetDM(dsp, &dm);CHKERRQ(ierr); 2375 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2376 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 2377 ierr = PetscFEGetQuadrature(fe, &quadDef);CHKERRQ(ierr); 2378 quad = quad ? quad : quadDef; 2379 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2380 ierr = PetscMalloc1(Nq*cdim, &cgeom->v);CHKERRQ(ierr); 2381 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->J);CHKERRQ(ierr); 2382 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);CHKERRQ(ierr); 2383 ierr = PetscMalloc1(Nq, &cgeom->detJ);CHKERRQ(ierr); 2384 cgeom->dim = dim; 2385 cgeom->dimEmbed = cdim; 2386 cgeom->numCells = 1; 2387 cgeom->numPoints = Nq; 2388 ierr = DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);CHKERRQ(ierr); 2389 PetscFunctionReturn(0); 2390 } 2391 2392 PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom) 2393 { 2394 PetscErrorCode ierr; 2395 2396 PetscFunctionBegin; 2397 ierr = PetscFree(cgeom->v);CHKERRQ(ierr); 2398 ierr = PetscFree(cgeom->J);CHKERRQ(ierr); 2399 ierr = PetscFree(cgeom->invJ);CHKERRQ(ierr); 2400 ierr = PetscFree(cgeom->detJ);CHKERRQ(ierr); 2401 PetscFunctionReturn(0); 2402 } 2403