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 + prob - The PetscDS specifying the discretizations and continuum functions 1437 . field - The 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 prob, PetscInt field, 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 PetscFunctionBegin; 1468 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1469 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1470 if (fe->ops->integrateresidual) {ierr = (*fe->ops->integrateresidual)(prob, field, 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 + prob - The PetscDS specifying the discretizations and continuum functions 1481 . field - The field being integrated 1482 . Ne - The number of elements in the chunk 1483 . fgeom - The face geometry for each cell in the chunk 1484 . coefficients - The array of FEM basis coefficients for the elements 1485 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1486 . probAux - The PetscDS specifying the auxiliary discretizations 1487 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1488 - t - The time 1489 1490 Output Parameter: 1491 . elemVec - the element residual vectors from each element 1492 1493 Level: intermediate 1494 1495 .seealso: PetscFEIntegrateResidual() 1496 @*/ 1497 PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, 1498 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1499 { 1500 PetscFE fe; 1501 PetscErrorCode ierr; 1502 1503 PetscFunctionBegin; 1504 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1505 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1506 if (fe->ops->integratebdresidual) {ierr = (*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1507 PetscFunctionReturn(0); 1508 } 1509 1510 /*@C 1511 PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration 1512 1513 Not collective 1514 1515 Input Parameters: 1516 + prob - The PetscDS specifying the discretizations and continuum functions 1517 . field - The field being integrated 1518 . Ne - The number of elements in the chunk 1519 . fgeom - The face geometry for each cell in the chunk 1520 . coefficients - The array of FEM basis coefficients for the elements 1521 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1522 . probAux - The PetscDS specifying the auxiliary discretizations 1523 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1524 - t - The time 1525 1526 Output Parameter 1527 . elemVec - the element residual vectors from each element 1528 1529 Level: developer 1530 1531 .seealso: PetscFEIntegrateResidual() 1532 @*/ 1533 PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, 1534 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1535 { 1536 PetscFE fe; 1537 PetscErrorCode ierr; 1538 1539 PetscFunctionBegin; 1540 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1541 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1542 if (fe->ops->integratehybridresidual) {ierr = (*fe->ops->integratehybridresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1543 PetscFunctionReturn(0); 1544 } 1545 1546 /*@C 1547 PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration 1548 1549 Not collective 1550 1551 Input Parameters: 1552 + prob - The PetscDS specifying the discretizations and continuum functions 1553 . jtype - The type of matrix pointwise functions that should be used 1554 . fieldI - The test field being integrated 1555 . fieldJ - The basis field 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 prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, 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 PetscErrorCode ierr; 1586 1587 PetscFunctionBegin; 1588 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1589 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1590 if (fe->ops->integratejacobian) {ierr = (*fe->ops->integratejacobian)(prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1591 PetscFunctionReturn(0); 1592 } 1593 1594 /*@C 1595 PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration 1596 1597 Not collective 1598 1599 Input Parameters: 1600 + prob - The PetscDS specifying the discretizations and continuum functions 1601 . fieldI - The test field being integrated 1602 . fieldJ - The basis field being integrated 1603 . Ne - The number of elements in the chunk 1604 . fgeom - The face geometry for each cell in the chunk 1605 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1606 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1607 . probAux - The PetscDS specifying the auxiliary discretizations 1608 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1609 . t - The time 1610 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1611 1612 Output Parameter: 1613 . elemMat - the element matrices for the Jacobian from each element 1614 1615 Note: 1616 $ Loop over batch of elements (e): 1617 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1618 $ Loop over quadrature points (q): 1619 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1620 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1621 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1622 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1623 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1624 Level: intermediate 1625 1626 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1627 @*/ 1628 PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, 1629 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1630 { 1631 PetscFE fe; 1632 PetscErrorCode ierr; 1633 1634 PetscFunctionBegin; 1635 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1636 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1637 if (fe->ops->integratebdjacobian) {ierr = (*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1638 PetscFunctionReturn(0); 1639 } 1640 1641 /*@C 1642 PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration 1643 1644 Not collective 1645 1646 Input Parameters: 1647 + prob - The PetscDS specifying the discretizations and continuum functions 1648 . jtype - The type of matrix pointwise functions that should be used 1649 . fieldI - The test field being integrated 1650 . fieldJ - The basis field being integrated 1651 . Ne - The number of elements in the chunk 1652 . fgeom - The face geometry for each cell in the chunk 1653 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1654 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1655 . probAux - The PetscDS specifying the auxiliary discretizations 1656 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1657 . t - The time 1658 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1659 1660 Output Parameter 1661 . elemMat - the element matrices for the Jacobian from each element 1662 1663 Note: 1664 $ Loop over batch of elements (e): 1665 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1666 $ Loop over quadrature points (q): 1667 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1668 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1669 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1670 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1671 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1672 Level: developer 1673 1674 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1675 @*/ 1676 PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, 1677 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1678 { 1679 PetscFE fe; 1680 PetscErrorCode ierr; 1681 1682 PetscFunctionBegin; 1683 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1684 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1685 if (fe->ops->integratehybridjacobian) {ierr = (*fe->ops->integratehybridjacobian)(prob, jtype, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1686 PetscFunctionReturn(0); 1687 } 1688 1689 /*@ 1690 PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height 1691 1692 Input Parameters: 1693 + fe - The finite element space 1694 - height - The height of the Plex point 1695 1696 Output Parameter: 1697 . subfe - The subspace of this FE space 1698 1699 Note: For example, if we want the subspace of this space for a face, we would choose height = 1. 1700 1701 Level: advanced 1702 1703 .seealso: PetscFECreateDefault() 1704 @*/ 1705 PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe) 1706 { 1707 PetscSpace P, subP; 1708 PetscDualSpace Q, subQ; 1709 PetscQuadrature subq; 1710 PetscFEType fetype; 1711 PetscInt dim, Nc; 1712 PetscErrorCode ierr; 1713 1714 PetscFunctionBegin; 1715 PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1); 1716 PetscValidPointer(subfe, 3); 1717 if (height == 0) { 1718 *subfe = fe; 1719 PetscFunctionReturn(0); 1720 } 1721 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1722 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1723 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1724 ierr = PetscFEGetFaceQuadrature(fe, &subq);CHKERRQ(ierr); 1725 ierr = PetscDualSpaceGetDimension(Q, &dim);CHKERRQ(ierr); 1726 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); 1727 if (!fe->subspaces) {ierr = PetscCalloc1(dim, &fe->subspaces);CHKERRQ(ierr);} 1728 if (height <= dim) { 1729 if (!fe->subspaces[height-1]) { 1730 PetscFE sub = NULL; 1731 const char *name; 1732 1733 ierr = PetscSpaceGetHeightSubspace(P, height, &subP);CHKERRQ(ierr); 1734 ierr = PetscDualSpaceGetHeightSubspace(Q, height, &subQ);CHKERRQ(ierr); 1735 if (subQ) { 1736 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);CHKERRQ(ierr); 1737 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 1738 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 1739 ierr = PetscFEGetType(fe, &fetype);CHKERRQ(ierr); 1740 ierr = PetscFESetType(sub, fetype);CHKERRQ(ierr); 1741 ierr = PetscFESetBasisSpace(sub, subP);CHKERRQ(ierr); 1742 ierr = PetscFESetDualSpace(sub, subQ);CHKERRQ(ierr); 1743 ierr = PetscFESetNumComponents(sub, Nc);CHKERRQ(ierr); 1744 ierr = PetscFESetUp(sub);CHKERRQ(ierr); 1745 ierr = PetscFESetQuadrature(sub, subq);CHKERRQ(ierr); 1746 } 1747 fe->subspaces[height-1] = sub; 1748 } 1749 *subfe = fe->subspaces[height-1]; 1750 } else { 1751 *subfe = NULL; 1752 } 1753 PetscFunctionReturn(0); 1754 } 1755 1756 /*@ 1757 PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used 1758 to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more 1759 sparsity). It is also used to create an interpolation between regularly refined meshes. 1760 1761 Collective on fem 1762 1763 Input Parameter: 1764 . fe - The initial PetscFE 1765 1766 Output Parameter: 1767 . feRef - The refined PetscFE 1768 1769 Level: advanced 1770 1771 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 1772 @*/ 1773 PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef) 1774 { 1775 PetscSpace P, Pref; 1776 PetscDualSpace Q, Qref; 1777 DM K, Kref; 1778 PetscQuadrature q, qref; 1779 const PetscReal *v0, *jac; 1780 PetscInt numComp, numSubelements; 1781 PetscInt cStart, cEnd, c; 1782 PetscDualSpace *cellSpaces; 1783 PetscErrorCode ierr; 1784 1785 PetscFunctionBegin; 1786 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1787 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1788 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1789 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr); 1790 /* Create space */ 1791 ierr = PetscObjectReference((PetscObject) P);CHKERRQ(ierr); 1792 Pref = P; 1793 /* Create dual space */ 1794 ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr); 1795 ierr = PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);CHKERRQ(ierr); 1796 ierr = DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);CHKERRQ(ierr); 1797 ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr); 1798 ierr = DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);CHKERRQ(ierr); 1799 ierr = PetscMalloc1(cEnd - cStart, &cellSpaces);CHKERRQ(ierr); 1800 /* TODO: fix for non-uniform refinement */ 1801 for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q; 1802 ierr = PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);CHKERRQ(ierr); 1803 ierr = PetscFree(cellSpaces);CHKERRQ(ierr); 1804 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 1805 ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr); 1806 /* Create element */ 1807 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);CHKERRQ(ierr); 1808 ierr = PetscFESetType(*feRef, PETSCFECOMPOSITE);CHKERRQ(ierr); 1809 ierr = PetscFESetBasisSpace(*feRef, Pref);CHKERRQ(ierr); 1810 ierr = PetscFESetDualSpace(*feRef, Qref);CHKERRQ(ierr); 1811 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 1812 ierr = PetscFESetNumComponents(*feRef, numComp);CHKERRQ(ierr); 1813 ierr = PetscFESetUp(*feRef);CHKERRQ(ierr); 1814 ierr = PetscSpaceDestroy(&Pref);CHKERRQ(ierr); 1815 ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr); 1816 /* Create quadrature */ 1817 ierr = PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);CHKERRQ(ierr); 1818 ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr); 1819 ierr = PetscFESetQuadrature(*feRef, qref);CHKERRQ(ierr); 1820 ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr); 1821 PetscFunctionReturn(0); 1822 } 1823 1824 /*@C 1825 PetscFECreateDefault - Create a PetscFE for basic FEM computation 1826 1827 Collective 1828 1829 Input Parameters: 1830 + comm - The MPI comm 1831 . dim - The spatial dimension 1832 . Nc - The number of components 1833 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1834 . prefix - The options prefix, or NULL 1835 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1836 1837 Output Parameter: 1838 . fem - The PetscFE object 1839 1840 Note: 1841 Each object is SetFromOption() during creation, so that the object may be customized from the command line. 1842 1843 Level: beginner 1844 1845 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1846 @*/ 1847 PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem) 1848 { 1849 PetscQuadrature q, fq; 1850 DM K; 1851 PetscSpace P; 1852 PetscDualSpace Q; 1853 PetscInt order, quadPointsPerEdge; 1854 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1855 PetscErrorCode ierr; 1856 1857 PetscFunctionBegin; 1858 /* Create space */ 1859 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1860 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr); 1861 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1862 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1863 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1864 ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr); 1865 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1866 ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr); 1867 ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr); 1868 /* Create dual space */ 1869 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1870 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1871 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr); 1872 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1873 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1874 ierr = DMDestroy(&K);CHKERRQ(ierr); 1875 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1876 ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr); 1877 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1878 ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr); 1879 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1880 /* Create element */ 1881 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1882 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr); 1883 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1884 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1885 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1886 ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr); 1887 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1888 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1889 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1890 /* Create quadrature (with specified order if given) */ 1891 qorder = qorder >= 0 ? qorder : order; 1892 ierr = PetscObjectOptionsBegin((PetscObject)*fem);CHKERRQ(ierr); 1893 ierr = PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);CHKERRQ(ierr); 1894 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1895 quadPointsPerEdge = PetscMax(qorder + 1,1); 1896 if (isSimplex) { 1897 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1898 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1899 } else { 1900 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1901 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1902 } 1903 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1904 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1905 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1906 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909 1910 /*@ 1911 PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k 1912 1913 Collective 1914 1915 Input Parameters: 1916 + comm - The MPI comm 1917 . dim - The spatial dimension 1918 . Nc - The number of components 1919 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1920 . k - The degree k of the space 1921 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1922 1923 Output Parameter: 1924 . fem - The PetscFE object 1925 1926 Level: beginner 1927 1928 Notes: 1929 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. 1930 1931 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1932 @*/ 1933 PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem) 1934 { 1935 PetscQuadrature q, fq; 1936 DM K; 1937 PetscSpace P; 1938 PetscDualSpace Q; 1939 PetscInt quadPointsPerEdge; 1940 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1941 char name[64]; 1942 PetscErrorCode ierr; 1943 1944 PetscFunctionBegin; 1945 /* Create space */ 1946 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1947 ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1948 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1949 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1950 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1951 ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr); 1952 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1953 /* Create dual space */ 1954 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1955 ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1956 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1957 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1958 ierr = DMDestroy(&K);CHKERRQ(ierr); 1959 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1960 ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr); 1961 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1962 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1963 /* Create finite element */ 1964 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1965 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1966 ierr = PetscObjectSetName((PetscObject) *fem, name);CHKERRQ(ierr); 1967 ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr); 1968 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1969 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1970 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1971 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1972 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1973 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1974 /* Create quadrature (with specified order if given) */ 1975 qorder = qorder >= 0 ? qorder : k; 1976 quadPointsPerEdge = PetscMax(qorder + 1,1); 1977 if (isSimplex) { 1978 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1979 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1980 } else { 1981 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1982 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1983 } 1984 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1985 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1986 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1987 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1988 /* Set finite element name */ 1989 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1990 ierr = PetscFESetName(*fem, name);CHKERRQ(ierr); 1991 PetscFunctionReturn(0); 1992 } 1993 1994 /*@C 1995 PetscFESetName - Names the FE and its subobjects 1996 1997 Not collective 1998 1999 Input Parameters: 2000 + fe - The PetscFE 2001 - name - The name 2002 2003 Level: intermediate 2004 2005 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 2006 @*/ 2007 PetscErrorCode PetscFESetName(PetscFE fe, const char name[]) 2008 { 2009 PetscSpace P; 2010 PetscDualSpace Q; 2011 PetscErrorCode ierr; 2012 2013 PetscFunctionBegin; 2014 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 2015 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2016 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 2017 ierr = PetscObjectSetName((PetscObject) P, name);CHKERRQ(ierr); 2018 ierr = PetscObjectSetName((PetscObject) Q, name);CHKERRQ(ierr); 2019 PetscFunctionReturn(0); 2020 } 2021 2022 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[]) 2023 { 2024 PetscInt dOffset = 0, fOffset = 0, f, g; 2025 PetscErrorCode ierr; 2026 2027 for (f = 0; f < Nf; ++f) { 2028 PetscFE fe; 2029 const PetscInt k = ds->jetDegree[f]; 2030 const PetscInt cdim = T[f]->cdim; 2031 const PetscInt Nq = T[f]->Np; 2032 const PetscInt Nbf = T[f]->Nb; 2033 const PetscInt Ncf = T[f]->Nc; 2034 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2035 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2036 const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r*Nq+q)*Nbf*Ncf*cdim*cdim] : NULL; 2037 PetscInt hOffset = 0, b, c, d; 2038 2039 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 2040 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2041 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2042 for (b = 0; b < Nbf; ++b) { 2043 for (c = 0; c < Ncf; ++c) { 2044 const PetscInt cidx = b*Ncf+c; 2045 2046 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2047 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2048 } 2049 } 2050 if (k > 1) { 2051 for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc*cdim; 2052 for (d = 0; d < cdim*cdim*Ncf; ++d) u_x[hOffset+fOffset*cdim*cdim+d] = 0.0; 2053 for (b = 0; b < Nbf; ++b) { 2054 for (c = 0; c < Ncf; ++c) { 2055 const PetscInt cidx = b*Ncf+c; 2056 2057 for (d = 0; d < cdim*cdim; ++d) u_x[hOffset+(fOffset+c)*cdim*cdim+d] += Hq[cidx*cdim*cdim+d]*coefficients[dOffset+b]; 2058 } 2059 } 2060 ierr = PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset+fOffset*cdim*cdim]);CHKERRQ(ierr); 2061 } 2062 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2063 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2064 if (u_t) { 2065 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2066 for (b = 0; b < Nbf; ++b) { 2067 for (c = 0; c < Ncf; ++c) { 2068 const PetscInt cidx = b*Ncf+c; 2069 2070 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2071 } 2072 } 2073 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2074 } 2075 fOffset += Ncf; 2076 dOffset += Nbf; 2077 } 2078 return 0; 2079 } 2080 2081 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[]) 2082 { 2083 PetscInt dOffset = 0, fOffset = 0, g; 2084 PetscErrorCode ierr; 2085 2086 for (g = 0; g < 2*Nf-1; ++g) { 2087 if (!T[g/2]) continue; 2088 { 2089 PetscFE fe; 2090 const PetscInt f = g/2; 2091 const PetscInt cdim = T[f]->cdim; 2092 const PetscInt Nq = T[f]->Np; 2093 const PetscInt Nbf = T[f]->Nb; 2094 const PetscInt Ncf = T[f]->Nc; 2095 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2096 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2097 PetscInt b, c, d; 2098 2099 fe = (PetscFE) ds->disc[f]; 2100 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2101 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2102 for (b = 0; b < Nbf; ++b) { 2103 for (c = 0; c < Ncf; ++c) { 2104 const PetscInt cidx = b*Ncf+c; 2105 2106 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2107 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2108 } 2109 } 2110 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2111 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2112 if (u_t) { 2113 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2114 for (b = 0; b < Nbf; ++b) { 2115 for (c = 0; c < Ncf; ++c) { 2116 const PetscInt cidx = b*Ncf+c; 2117 2118 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2119 } 2120 } 2121 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2122 } 2123 fOffset += Ncf; 2124 dOffset += Nbf; 2125 } 2126 } 2127 return 0; 2128 } 2129 2130 PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[]) 2131 { 2132 PetscFE fe; 2133 PetscTabulation Tc; 2134 PetscInt b, c; 2135 PetscErrorCode ierr; 2136 2137 if (!prob) return 0; 2138 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 2139 ierr = PetscFEGetFaceCentroidTabulation(fe, &Tc);CHKERRQ(ierr); 2140 { 2141 const PetscReal *faceBasis = Tc->T[0]; 2142 const PetscInt Nb = Tc->Nb; 2143 const PetscInt Nc = Tc->Nc; 2144 2145 for (c = 0; c < Nc; ++c) {u[c] = 0.0;} 2146 for (b = 0; b < Nb; ++b) { 2147 for (c = 0; c < Nc; ++c) { 2148 u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c]; 2149 } 2150 } 2151 } 2152 return 0; 2153 } 2154 2155 PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2156 { 2157 const PetscInt dE = T->cdim; /* fegeom->dimEmbed */ 2158 const PetscInt Nq = T->Np; 2159 const PetscInt Nb = T->Nb; 2160 const PetscInt Nc = T->Nc; 2161 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2162 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2163 PetscInt q, b, c, d; 2164 PetscErrorCode ierr; 2165 2166 for (b = 0; b < Nb; ++b) elemVec[b] = 0.0; 2167 for (q = 0; q < Nq; ++q) { 2168 for (b = 0; b < Nb; ++b) { 2169 for (c = 0; c < Nc; ++c) { 2170 const PetscInt bcidx = b*Nc+c; 2171 2172 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2173 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2174 } 2175 } 2176 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2177 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2178 for (b = 0; b < Nb; ++b) { 2179 for (c = 0; c < Nc; ++c) { 2180 const PetscInt bcidx = b*Nc+c; 2181 const PetscInt qcidx = q*Nc+c; 2182 2183 elemVec[b] += tmpBasis[bcidx]*f0[qcidx]; 2184 for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2185 } 2186 } 2187 } 2188 return(0); 2189 } 2190 2191 PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2192 { 2193 const PetscInt dE = T->cdim; 2194 const PetscInt Nq = T->Np; 2195 const PetscInt Nb = T->Nb; 2196 const PetscInt Nc = T->Nc; 2197 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2198 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2199 PetscInt q, b, c, d, s; 2200 PetscErrorCode ierr; 2201 2202 for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0; 2203 for (q = 0; q < Nq; ++q) { 2204 for (b = 0; b < Nb; ++b) { 2205 for (c = 0; c < Nc; ++c) { 2206 const PetscInt bcidx = b*Nc+c; 2207 2208 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2209 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2210 } 2211 } 2212 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2213 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2214 for (s = 0; s < 2; ++s) { 2215 for (b = 0; b < Nb; ++b) { 2216 for (c = 0; c < Nc; ++c) { 2217 const PetscInt bcidx = b*Nc+c; 2218 const PetscInt qcidx = (q*2+s)*Nc+c; 2219 2220 elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx]; 2221 for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2222 } 2223 } 2224 } 2225 } 2226 return(0); 2227 } 2228 2229 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[]) 2230 { 2231 const PetscInt dE = TI->cdim; 2232 const PetscInt NqI = TI->Np; 2233 const PetscInt NbI = TI->Nb; 2234 const PetscInt NcI = TI->Nc; 2235 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2236 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2237 const PetscInt NqJ = TJ->Np; 2238 const PetscInt NbJ = TJ->Nb; 2239 const PetscInt NcJ = TJ->Nc; 2240 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2241 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2242 PetscInt f, fc, g, gc, df, dg; 2243 PetscErrorCode ierr; 2244 2245 for (f = 0; f < NbI; ++f) { 2246 for (fc = 0; fc < NcI; ++fc) { 2247 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2248 2249 tmpBasisI[fidx] = basisI[fidx]; 2250 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2251 } 2252 } 2253 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2254 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2255 for (g = 0; g < NbJ; ++g) { 2256 for (gc = 0; gc < NcJ; ++gc) { 2257 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2258 2259 tmpBasisJ[gidx] = basisJ[gidx]; 2260 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2261 } 2262 } 2263 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2264 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2265 for (f = 0; f < NbI; ++f) { 2266 for (fc = 0; fc < NcI; ++fc) { 2267 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2268 const PetscInt i = offsetI+f; /* Element matrix row */ 2269 for (g = 0; g < NbJ; ++g) { 2270 for (gc = 0; gc < NcJ; ++gc) { 2271 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2272 const PetscInt j = offsetJ+g; /* Element matrix column */ 2273 const PetscInt fOff = eOffset+i*totDim+j; 2274 2275 elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx]; 2276 for (df = 0; df < dE; ++df) { 2277 elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2278 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx]; 2279 for (dg = 0; dg < dE; ++dg) { 2280 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2281 } 2282 } 2283 } 2284 } 2285 } 2286 } 2287 return(0); 2288 } 2289 2290 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[]) 2291 { 2292 const PetscInt dE = TI->cdim; 2293 const PetscInt NqI = TI->Np; 2294 const PetscInt NbI = TI->Nb; 2295 const PetscInt NcI = TI->Nc; 2296 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2297 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2298 const PetscInt NqJ = TJ->Np; 2299 const PetscInt NbJ = TJ->Nb; 2300 const PetscInt NcJ = TJ->Nc; 2301 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2302 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2303 const PetscInt Ns = isHybridI ? 1 : 2; 2304 const PetscInt Nt = isHybridJ ? 1 : 2; 2305 PetscInt f, fc, g, gc, df, dg, s, t; 2306 PetscErrorCode ierr; 2307 2308 for (f = 0; f < NbI; ++f) { 2309 for (fc = 0; fc < NcI; ++fc) { 2310 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2311 2312 tmpBasisI[fidx] = basisI[fidx]; 2313 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2314 } 2315 } 2316 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2317 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2318 for (g = 0; g < NbJ; ++g) { 2319 for (gc = 0; gc < NcJ; ++gc) { 2320 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2321 2322 tmpBasisJ[gidx] = basisJ[gidx]; 2323 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2324 } 2325 } 2326 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2327 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2328 for (s = 0; s < Ns; ++s) { 2329 for (f = 0; f < NbI; ++f) { 2330 for (fc = 0; fc < NcI; ++fc) { 2331 const PetscInt sc = NcI*s+fc; /* components from each side of the surface */ 2332 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2333 const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */ 2334 for (t = 0; t < Nt; ++t) { 2335 for (g = 0; g < NbJ; ++g) { 2336 for (gc = 0; gc < NcJ; ++gc) { 2337 const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */ 2338 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2339 const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */ 2340 const PetscInt fOff = eOffset+i*totDim+j; 2341 2342 elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx]; 2343 for (df = 0; df < dE; ++df) { 2344 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2345 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx]; 2346 for (dg = 0; dg < dE; ++dg) { 2347 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2348 } 2349 } 2350 } 2351 } 2352 } 2353 } 2354 } 2355 } 2356 return(0); 2357 } 2358 2359 PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom) 2360 { 2361 PetscDualSpace dsp; 2362 DM dm; 2363 PetscQuadrature quadDef; 2364 PetscInt dim, cdim, Nq; 2365 PetscErrorCode ierr; 2366 2367 PetscFunctionBegin; 2368 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); 2369 ierr = PetscDualSpaceGetDM(dsp, &dm);CHKERRQ(ierr); 2370 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2371 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 2372 ierr = PetscFEGetQuadrature(fe, &quadDef);CHKERRQ(ierr); 2373 quad = quad ? quad : quadDef; 2374 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2375 ierr = PetscMalloc1(Nq*cdim, &cgeom->v);CHKERRQ(ierr); 2376 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->J);CHKERRQ(ierr); 2377 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);CHKERRQ(ierr); 2378 ierr = PetscMalloc1(Nq, &cgeom->detJ);CHKERRQ(ierr); 2379 cgeom->dim = dim; 2380 cgeom->dimEmbed = cdim; 2381 cgeom->numCells = 1; 2382 cgeom->numPoints = Nq; 2383 ierr = DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);CHKERRQ(ierr); 2384 PetscFunctionReturn(0); 2385 } 2386 2387 PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom) 2388 { 2389 PetscErrorCode ierr; 2390 2391 PetscFunctionBegin; 2392 ierr = PetscFree(cgeom->v);CHKERRQ(ierr); 2393 ierr = PetscFree(cgeom->J);CHKERRQ(ierr); 2394 ierr = PetscFree(cgeom->invJ);CHKERRQ(ierr); 2395 ierr = PetscFree(cgeom->detJ);CHKERRQ(ierr); 2396 PetscFunctionReturn(0); 2397 } 2398