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