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