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