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 PetscFELimitDegree - Copy a `PetscFE` but limit the degree to be in the given range 2177 2178 Collective 2179 2180 Input Parameters: 2181 + fe - The `PetscFE` 2182 . minDegree - The minimum degree, or `PETSC_DETERMINE` for no limit 2183 - maxDegree - The maximum degree, or `PETSC_DETERMINE` for no limit 2184 2185 Output Parameter: 2186 . newfe - The `PetscFE` object 2187 2188 Level: advanced 2189 2190 Note: 2191 This currently only works for Lagrange elements. 2192 2193 .seealso: `PetscFECreateLagrange()`, `PetscFECreateDefault()`, `PetscFECreateByCell()`, `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()` 2194 @*/ 2195 PetscErrorCode PetscFELimitDegree(PetscFE fe, PetscInt minDegree, PetscInt maxDegree, PetscFE *newfe) 2196 { 2197 PetscDualSpace Q; 2198 PetscBool islag, issum; 2199 PetscInt oldk = 0, k; 2200 2201 PetscFunctionBegin; 2202 PetscCall(PetscFEGetDualSpace(fe, &Q)); 2203 PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACELAGRANGE, &islag)); 2204 PetscCall(PetscObjectTypeCompare((PetscObject)Q, PETSCDUALSPACESUM, &issum)); 2205 if (islag) { 2206 PetscCall(PetscDualSpaceGetOrder(Q, &oldk)); 2207 } else if (issum) { 2208 PetscDualSpace subQ; 2209 2210 PetscCall(PetscDualSpaceSumGetSubspace(Q, 0, &subQ)); 2211 PetscCall(PetscDualSpaceGetOrder(subQ, &oldk)); 2212 } else { 2213 PetscCall(PetscObjectReference((PetscObject)fe)); 2214 *newfe = fe; 2215 PetscFunctionReturn(PETSC_SUCCESS); 2216 } 2217 k = oldk; 2218 if (minDegree >= 0) k = PetscMax(k, minDegree); 2219 if (maxDegree >= 0) k = PetscMin(k, maxDegree); 2220 if (k != oldk) { 2221 DM K; 2222 PetscSpace P; 2223 PetscQuadrature q; 2224 DMPolytopeType ct; 2225 PetscInt dim, Nc; 2226 2227 PetscCall(PetscFEGetBasisSpace(fe, &P)); 2228 PetscCall(PetscSpaceGetNumVariables(P, &dim)); 2229 PetscCall(PetscSpaceGetNumComponents(P, &Nc)); 2230 PetscCall(PetscDualSpaceGetDM(Q, &K)); 2231 PetscCall(DMPlexGetCellType(K, 0, &ct)); 2232 PetscCall(PetscFECreateLagrangeByCell(PetscObjectComm((PetscObject)fe), dim, Nc, ct, k, PETSC_DETERMINE, newfe)); 2233 PetscCall(PetscFEGetQuadrature(fe, &q)); 2234 PetscCall(PetscFESetQuadrature(*newfe, q)); 2235 } else { 2236 PetscCall(PetscObjectReference((PetscObject)fe)); 2237 *newfe = fe; 2238 } 2239 PetscFunctionReturn(PETSC_SUCCESS); 2240 } 2241 2242 /*@ 2243 PetscFESetName - Names the `PetscFE` and its subobjects 2244 2245 Not Collective 2246 2247 Input Parameters: 2248 + fe - The `PetscFE` 2249 - name - The name 2250 2251 Level: intermediate 2252 2253 .seealso: `PetscFECreate()`, `PetscSpaceCreate()`, `PetscDualSpaceCreate()` 2254 @*/ 2255 PetscErrorCode PetscFESetName(PetscFE fe, const char name[]) 2256 { 2257 PetscSpace P; 2258 PetscDualSpace Q; 2259 2260 PetscFunctionBegin; 2261 PetscCall(PetscFEGetBasisSpace(fe, &P)); 2262 PetscCall(PetscFEGetDualSpace(fe, &Q)); 2263 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 2264 PetscCall(PetscObjectSetName((PetscObject)P, name)); 2265 PetscCall(PetscObjectSetName((PetscObject)Q, name)); 2266 PetscFunctionReturn(PETSC_SUCCESS); 2267 } 2268 2269 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[]) 2270 { 2271 PetscInt dOffset = 0, fOffset = 0, f, g; 2272 2273 for (f = 0; f < Nf; ++f) { 2274 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); 2275 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); 2276 PetscFE fe; 2277 const PetscInt k = ds->jetDegree[f]; 2278 const PetscInt cdim = T[f]->cdim; 2279 const PetscInt dE = fegeom->dimEmbed; 2280 const PetscInt Nq = T[f]->Np; 2281 const PetscInt Nbf = T[f]->Nb; 2282 const PetscInt Ncf = T[f]->Nc; 2283 const PetscReal *Bq = &T[f]->T[0][(r * Nq + q) * Nbf * Ncf]; 2284 const PetscReal *Dq = &T[f]->T[1][(r * Nq + q) * Nbf * Ncf * cdim]; 2285 const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r * Nq + q) * Nbf * Ncf * cdim * cdim] : NULL; 2286 PetscInt hOffset = 0, b, c, d; 2287 2288 PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe)); 2289 for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0; 2290 for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0; 2291 for (b = 0; b < Nbf; ++b) { 2292 for (c = 0; c < Ncf; ++c) { 2293 const PetscInt cidx = b * Ncf + c; 2294 2295 u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b]; 2296 for (d = 0; d < cdim; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * cdim + d] * coefficients[dOffset + b]; 2297 } 2298 } 2299 if (k > 1) { 2300 for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc * dE; 2301 for (d = 0; d < dE * dE * Ncf; ++d) u_x[hOffset + fOffset * dE * dE + d] = 0.0; 2302 for (b = 0; b < Nbf; ++b) { 2303 for (c = 0; c < Ncf; ++c) { 2304 const PetscInt cidx = b * Ncf + c; 2305 2306 for (d = 0; d < cdim * cdim; ++d) u_x[hOffset + (fOffset + c) * dE * dE + d] += Hq[cidx * cdim * cdim + d] * coefficients[dOffset + b]; 2307 } 2308 } 2309 PetscCall(PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset + fOffset * dE * dE])); 2310 } 2311 PetscCall(PetscFEPushforward(fe, fegeom, 1, &u[fOffset])); 2312 PetscCall(PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE])); 2313 if (u_t) { 2314 for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0; 2315 for (b = 0; b < Nbf; ++b) { 2316 for (c = 0; c < Ncf; ++c) { 2317 const PetscInt cidx = b * Ncf + c; 2318 2319 u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b]; 2320 } 2321 } 2322 PetscCall(PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset])); 2323 } 2324 fOffset += Ncf; 2325 dOffset += Nbf; 2326 } 2327 return PETSC_SUCCESS; 2328 } 2329 2330 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[]) 2331 { 2332 PetscInt dOffset = 0, fOffset = 0, f, g; 2333 2334 /* f is the field number in the DS, g is the field number in u[] */ 2335 for (f = 0, g = 0; f < Nf; ++f) { 2336 PetscBool isCohesive; 2337 PetscInt Ns, s; 2338 2339 if (!Tab[f]) continue; 2340 PetscCall(PetscDSGetCohesive(ds, f, &isCohesive)); 2341 Ns = isCohesive ? 1 : 2; 2342 { 2343 PetscTabulation T = isCohesive ? Tab[f] : Tabf[f]; 2344 PetscFE fe = (PetscFE)ds->disc[f]; 2345 const PetscInt dEt = T->cdim; 2346 const PetscInt dE = fegeom->dimEmbed; 2347 const PetscInt Nq = T->Np; 2348 const PetscInt Nbf = T->Nb; 2349 const PetscInt Ncf = T->Nc; 2350 2351 for (s = 0; s < Ns; ++s, ++g) { 2352 const PetscInt r = isCohesive ? rc : rf[s]; 2353 const PetscInt q = isCohesive ? qc : qf[s]; 2354 const PetscReal *Bq = &T->T[0][(r * Nq + q) * Nbf * Ncf]; 2355 const PetscReal *Dq = &T->T[1][(r * Nq + q) * Nbf * Ncf * dEt]; 2356 PetscInt b, c, d; 2357 2358 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); 2359 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); 2360 for (c = 0; c < Ncf; ++c) u[fOffset + c] = 0.0; 2361 for (d = 0; d < dE * Ncf; ++d) u_x[fOffset * dE + d] = 0.0; 2362 for (b = 0; b < Nbf; ++b) { 2363 for (c = 0; c < Ncf; ++c) { 2364 const PetscInt cidx = b * Ncf + c; 2365 2366 u[fOffset + c] += Bq[cidx] * coefficients[dOffset + b]; 2367 for (d = 0; d < dEt; ++d) u_x[(fOffset + c) * dE + d] += Dq[cidx * dEt + d] * coefficients[dOffset + b]; 2368 } 2369 } 2370 PetscCall(PetscFEPushforward(fe, fegeom, 1, &u[fOffset])); 2371 PetscCall(PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset * dE])); 2372 if (u_t) { 2373 for (c = 0; c < Ncf; ++c) u_t[fOffset + c] = 0.0; 2374 for (b = 0; b < Nbf; ++b) { 2375 for (c = 0; c < Ncf; ++c) { 2376 const PetscInt cidx = b * Ncf + c; 2377 2378 u_t[fOffset + c] += Bq[cidx] * coefficients_t[dOffset + b]; 2379 } 2380 } 2381 PetscCall(PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset])); 2382 } 2383 fOffset += Ncf; 2384 dOffset += Nbf; 2385 } 2386 } 2387 } 2388 return PETSC_SUCCESS; 2389 } 2390 2391 PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[]) 2392 { 2393 PetscFE fe; 2394 PetscTabulation Tc; 2395 PetscInt b, c; 2396 2397 if (!prob) return PETSC_SUCCESS; 2398 PetscCall(PetscDSGetDiscretization(prob, field, (PetscObject *)&fe)); 2399 PetscCall(PetscFEGetFaceCentroidTabulation(fe, &Tc)); 2400 { 2401 const PetscReal *faceBasis = Tc->T[0]; 2402 const PetscInt Nb = Tc->Nb; 2403 const PetscInt Nc = Tc->Nc; 2404 2405 for (c = 0; c < Nc; ++c) u[c] = 0.0; 2406 for (b = 0; b < Nb; ++b) { 2407 for (c = 0; c < Nc; ++c) u[c] += coefficients[b] * faceBasis[(faceLoc * Nb + b) * Nc + c]; 2408 } 2409 } 2410 return PETSC_SUCCESS; 2411 } 2412 2413 PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2414 { 2415 PetscFEGeom pgeom; 2416 const PetscInt dEt = T->cdim; 2417 const PetscInt dE = fegeom->dimEmbed; 2418 const PetscInt Nq = T->Np; 2419 const PetscInt Nb = T->Nb; 2420 const PetscInt Nc = T->Nc; 2421 const PetscReal *basis = &T->T[0][r * Nq * Nb * Nc]; 2422 const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dEt]; 2423 PetscInt q, b, c, d; 2424 2425 for (q = 0; q < Nq; ++q) { 2426 for (b = 0; b < Nb; ++b) { 2427 for (c = 0; c < Nc; ++c) { 2428 const PetscInt bcidx = b * Nc + c; 2429 2430 tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx]; 2431 for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dEt + bcidx * dEt + d]; 2432 for (d = dEt; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = 0.0; 2433 } 2434 } 2435 PetscCall(PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom)); 2436 PetscCall(PetscFEPushforward(fe, &pgeom, Nb, tmpBasis)); 2437 PetscCall(PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer)); 2438 for (b = 0; b < Nb; ++b) { 2439 for (c = 0; c < Nc; ++c) { 2440 const PetscInt bcidx = b * Nc + c; 2441 const PetscInt qcidx = q * Nc + c; 2442 2443 elemVec[b] += tmpBasis[bcidx] * f0[qcidx]; 2444 for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d]; 2445 } 2446 } 2447 } 2448 return PETSC_SUCCESS; 2449 } 2450 2451 PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscInt side, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2452 { 2453 const PetscInt dE = T->cdim; 2454 const PetscInt Nq = T->Np; 2455 const PetscInt Nb = T->Nb; 2456 const PetscInt Nc = T->Nc; 2457 const PetscReal *basis = &T->T[0][r * Nq * Nb * Nc]; 2458 const PetscReal *basisDer = &T->T[1][r * Nq * Nb * Nc * dE]; 2459 2460 for (PetscInt q = 0; q < Nq; ++q) { 2461 for (PetscInt b = 0; b < Nb; ++b) { 2462 for (PetscInt c = 0; c < Nc; ++c) { 2463 const PetscInt bcidx = b * Nc + c; 2464 2465 tmpBasis[bcidx] = basis[q * Nb * Nc + bcidx]; 2466 for (PetscInt d = 0; d < dE; ++d) tmpBasisDer[bcidx * dE + d] = basisDer[q * Nb * Nc * dE + bcidx * dE + d]; 2467 } 2468 } 2469 PetscCall(PetscFEPushforward(fe, fegeom, Nb, tmpBasis)); 2470 // TODO This is currently broken since we do not pull the geometry down to the lower dimension 2471 // PetscCall(PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer)); 2472 if (side == 2) { 2473 // Integrating over whole cohesive cell, so insert for both sides 2474 for (PetscInt s = 0; s < 2; ++s) { 2475 for (PetscInt b = 0; b < Nb; ++b) { 2476 for (PetscInt c = 0; c < Nc; ++c) { 2477 const PetscInt bcidx = b * Nc + c; 2478 const PetscInt qcidx = (q * 2 + s) * Nc + c; 2479 2480 elemVec[Nb * s + b] += tmpBasis[bcidx] * f0[qcidx]; 2481 for (PetscInt d = 0; d < dE; ++d) elemVec[Nb * s + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d]; 2482 } 2483 } 2484 } 2485 } else { 2486 // Integrating over endcaps of cohesive cell, so insert for correct side 2487 for (PetscInt b = 0; b < Nb; ++b) { 2488 for (PetscInt c = 0; c < Nc; ++c) { 2489 const PetscInt bcidx = b * Nc + c; 2490 const PetscInt qcidx = q * Nc + c; 2491 2492 elemVec[Nb * side + b] += tmpBasis[bcidx] * f0[qcidx]; 2493 for (PetscInt d = 0; d < dE; ++d) elemVec[Nb * side + b] += tmpBasisDer[bcidx * dE + d] * f1[qcidx * dE + d]; 2494 } 2495 } 2496 } 2497 } 2498 return PETSC_SUCCESS; 2499 } 2500 2501 #define petsc_elemmat_kernel_g1(_NbI, _NcI, _NbJ, _NcJ, _dE) \ 2502 do { \ 2503 for (PetscInt fc = 0; fc < (_NcI); ++fc) { \ 2504 for (PetscInt gc = 0; gc < (_NcJ); ++gc) { \ 2505 const PetscScalar *G = g1 + (fc * (_NcJ) + gc) * _dE; \ 2506 for (PetscInt f = 0; f < (_NbI); ++f) { \ 2507 const PetscScalar tBIv = tmpBasisI[f * (_NcI) + fc]; \ 2508 for (PetscInt g = 0; g < (_NbJ); ++g) { \ 2509 const PetscScalar *tBDJ = tmpBasisDerJ + (g * (_NcJ) + gc) * (_dE); \ 2510 PetscScalar s = 0.0; \ 2511 for (PetscInt df = 0; df < _dE; ++df) { s += G[df] * tBDJ[df]; } \ 2512 elemMat[(offsetI + f) * totDim + (offsetJ + g)] += s * tBIv; \ 2513 } \ 2514 } \ 2515 } \ 2516 } \ 2517 } while (0) 2518 2519 #define petsc_elemmat_kernel_g2(_NbI, _NcI, _NbJ, _NcJ, _dE) \ 2520 do { \ 2521 for (PetscInt gc = 0; gc < (_NcJ); ++gc) { \ 2522 for (PetscInt fc = 0; fc < (_NcI); ++fc) { \ 2523 const PetscScalar *G = g2 + (fc * (_NcJ) + gc) * _dE; \ 2524 for (PetscInt g = 0; g < (_NbJ); ++g) { \ 2525 const PetscScalar tBJv = tmpBasisJ[g * (_NcJ) + gc]; \ 2526 for (PetscInt f = 0; f < (_NbI); ++f) { \ 2527 const PetscScalar *tBDI = tmpBasisDerI + (f * (_NcI) + fc) * (_dE); \ 2528 PetscScalar s = 0.0; \ 2529 for (PetscInt df = 0; df < _dE; ++df) { s += tBDI[df] * G[df]; } \ 2530 elemMat[(offsetI + f) * totDim + (offsetJ + g)] += s * tBJv; \ 2531 } \ 2532 } \ 2533 } \ 2534 } \ 2535 } while (0) 2536 2537 #define petsc_elemmat_kernel_g3(_NbI, _NcI, _NbJ, _NcJ, _dE) \ 2538 do { \ 2539 for (PetscInt fc = 0; fc < (_NcI); ++fc) { \ 2540 for (PetscInt gc = 0; gc < (_NcJ); ++gc) { \ 2541 const PetscScalar *G = g3 + (fc * (_NcJ) + gc) * (_dE) * (_dE); \ 2542 for (PetscInt f = 0; f < (_NbI); ++f) { \ 2543 const PetscScalar *tBDI = tmpBasisDerI + (f * (_NcI) + fc) * (_dE); \ 2544 for (PetscInt g = 0; g < (_NbJ); ++g) { \ 2545 PetscScalar s = 0.0; \ 2546 const PetscScalar *tBDJ = tmpBasisDerJ + (g * (_NcJ) + gc) * (_dE); \ 2547 for (PetscInt df = 0; df < (_dE); ++df) { \ 2548 for (PetscInt dg = 0; dg < (_dE); ++dg) { s += tBDI[df] * G[df * (_dE) + dg] * tBDJ[dg]; } \ 2549 } \ 2550 elemMat[(offsetI + f) * totDim + (offsetJ + g)] += s; \ 2551 } \ 2552 } \ 2553 } \ 2554 } \ 2555 } while (0) 2556 2557 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 totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[]) 2558 { 2559 const PetscInt cdim = TI->cdim; 2560 const PetscInt dE = fegeom->dimEmbed; 2561 const PetscInt NqI = TI->Np; 2562 const PetscInt NbI = TI->Nb; 2563 const PetscInt NcI = TI->Nc; 2564 const PetscReal *basisI = &TI->T[0][(r * NqI + q) * NbI * NcI]; 2565 const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * cdim]; 2566 const PetscInt NqJ = TJ->Np; 2567 const PetscInt NbJ = TJ->Nb; 2568 const PetscInt NcJ = TJ->Nc; 2569 const PetscReal *basisJ = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ]; 2570 const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * cdim]; 2571 2572 for (PetscInt f = 0; f < NbI; ++f) { 2573 for (PetscInt fc = 0; fc < NcI; ++fc) { 2574 const PetscInt fidx = f * NcI + fc; /* Test function basis index */ 2575 2576 tmpBasisI[fidx] = basisI[fidx]; 2577 for (PetscInt df = 0; df < cdim; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * cdim + df]; 2578 } 2579 } 2580 PetscCall(PetscFEPushforward(feI, fegeom, NbI, tmpBasisI)); 2581 PetscCall(PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI)); 2582 if (feI != feJ) { 2583 for (PetscInt g = 0; g < NbJ; ++g) { 2584 for (PetscInt gc = 0; gc < NcJ; ++gc) { 2585 const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */ 2586 2587 tmpBasisJ[gidx] = basisJ[gidx]; 2588 for (PetscInt dg = 0; dg < cdim; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * cdim + dg]; 2589 } 2590 } 2591 PetscCall(PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ)); 2592 PetscCall(PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ)); 2593 } else { 2594 tmpBasisJ = tmpBasisI; 2595 tmpBasisDerJ = tmpBasisDerI; 2596 } 2597 if (PetscUnlikely(g0)) { 2598 for (PetscInt f = 0; f < NbI; ++f) { 2599 const PetscInt i = offsetI + f; /* Element matrix row */ 2600 2601 for (PetscInt fc = 0; fc < NcI; ++fc) { 2602 const PetscScalar bI = tmpBasisI[f * NcI + fc]; /* Test function basis value */ 2603 2604 for (PetscInt g = 0; g < NbJ; ++g) { 2605 const PetscInt j = offsetJ + g; /* Element matrix column */ 2606 const PetscInt fOff = i * totDim + j; 2607 2608 for (PetscInt gc = 0; gc < NcJ; ++gc) { elemMat[fOff] += bI * g0[fc * NcJ + gc] * tmpBasisJ[g * NcJ + gc]; } 2609 } 2610 } 2611 } 2612 } 2613 if (PetscUnlikely(g1)) { 2614 #if 1 2615 if (dE == 2) { 2616 petsc_elemmat_kernel_g1(NbI, NcI, NbJ, NcJ, 2); 2617 } else if (dE == 3) { 2618 petsc_elemmat_kernel_g1(NbI, NcI, NbJ, NcJ, 3); 2619 } else { 2620 petsc_elemmat_kernel_g1(NbI, NcI, NbJ, NcJ, dE); 2621 } 2622 #else 2623 for (PetscInt f = 0; f < NbI; ++f) { 2624 const PetscInt i = offsetI + f; /* Element matrix row */ 2625 2626 for (PetscInt fc = 0; fc < NcI; ++fc) { 2627 const PetscScalar bI = tmpBasisI[f * NcI + fc]; /* Test function basis value */ 2628 2629 for (PetscInt g = 0; g < NbJ; ++g) { 2630 const PetscInt j = offsetJ + g; /* Element matrix column */ 2631 const PetscInt fOff = i * totDim + j; 2632 2633 for (PetscInt gc = 0; gc < NcJ; ++gc) { 2634 const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */ 2635 2636 for (PetscInt df = 0; df < dE; ++df) { elemMat[fOff] += bI * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df]; } 2637 } 2638 } 2639 } 2640 } 2641 #endif 2642 } 2643 if (PetscUnlikely(g2)) { 2644 #if 1 2645 if (dE == 2) { 2646 petsc_elemmat_kernel_g2(NbI, NcI, NbJ, NcJ, 2); 2647 } else if (dE == 3) { 2648 petsc_elemmat_kernel_g2(NbI, NcI, NbJ, NcJ, 3); 2649 } else { 2650 petsc_elemmat_kernel_g2(NbI, NcI, NbJ, NcJ, dE); 2651 } 2652 #else 2653 for (PetscInt g = 0; g < NbJ; ++g) { 2654 const PetscInt j = offsetJ + g; /* Element matrix column */ 2655 2656 for (PetscInt gc = 0; gc < NcJ; ++gc) { 2657 const PetscScalar bJ = tmpBasisJ[g * NcJ + gc]; /* Trial function basis value */ 2658 2659 for (PetscInt f = 0; f < NbI; ++f) { 2660 const PetscInt i = offsetI + f; /* Element matrix row */ 2661 const PetscInt fOff = i * totDim + j; 2662 2663 for (PetscInt fc = 0; fc < NcI; ++fc) { 2664 const PetscInt fidx = f * NcI + fc; /* Test function basis index */ 2665 2666 for (PetscInt df = 0; df < dE; ++df) { elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * bJ; } 2667 } 2668 } 2669 } 2670 } 2671 #endif 2672 } 2673 if (PetscUnlikely(g3)) { 2674 #if 1 2675 if (dE == 2) { 2676 petsc_elemmat_kernel_g3(NbI, NcI, NbJ, NcJ, 2); 2677 } else if (dE == 3) { 2678 petsc_elemmat_kernel_g3(NbI, NcI, NbJ, NcJ, 3); 2679 } else { 2680 petsc_elemmat_kernel_g3(NbI, NcI, NbJ, NcJ, dE); 2681 } 2682 #else 2683 for (PetscInt f = 0; f < NbI; ++f) { 2684 const PetscInt i = offsetI + f; /* Element matrix row */ 2685 2686 for (PetscInt fc = 0; fc < NcI; ++fc) { 2687 const PetscInt fidx = f * NcI + fc; /* Test function basis index */ 2688 2689 for (PetscInt g = 0; g < NbJ; ++g) { 2690 const PetscInt j = offsetJ + g; /* Element matrix column */ 2691 const PetscInt fOff = i * totDim + j; 2692 2693 for (PetscInt gc = 0; gc < NcJ; ++gc) { 2694 const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */ 2695 2696 for (PetscInt df = 0; df < dE; ++df) { 2697 for (PetscInt dg = 0; dg < dE; ++dg) { elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg]; } 2698 } 2699 } 2700 } 2701 } 2702 } 2703 #endif 2704 } 2705 return PETSC_SUCCESS; 2706 } 2707 2708 #undef petsc_elemmat_kernel_g1 2709 #undef petsc_elemmat_kernel_g2 2710 #undef petsc_elemmat_kernel_g3 2711 2712 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[]) 2713 { 2714 const PetscInt dE = TI->cdim; 2715 const PetscInt NqI = TI->Np; 2716 const PetscInt NbI = TI->Nb; 2717 const PetscInt NcI = TI->Nc; 2718 const PetscReal *basisI = &TI->T[0][(r * NqI + q) * NbI * NcI]; 2719 const PetscReal *basisDerI = &TI->T[1][(r * NqI + q) * NbI * NcI * dE]; 2720 const PetscInt NqJ = TJ->Np; 2721 const PetscInt NbJ = TJ->Nb; 2722 const PetscInt NcJ = TJ->Nc; 2723 const PetscReal *basisJ = &TJ->T[0][(r * NqJ + q) * NbJ * NcJ]; 2724 const PetscReal *basisDerJ = &TJ->T[1][(r * NqJ + q) * NbJ * NcJ * dE]; 2725 const PetscInt so = isHybridI ? 0 : s; 2726 const PetscInt to = isHybridJ ? 0 : t; 2727 PetscInt f, fc, g, gc, df, dg; 2728 2729 for (f = 0; f < NbI; ++f) { 2730 for (fc = 0; fc < NcI; ++fc) { 2731 const PetscInt fidx = f * NcI + fc; /* Test function basis index */ 2732 2733 tmpBasisI[fidx] = basisI[fidx]; 2734 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx * dE + df] = basisDerI[fidx * dE + df]; 2735 } 2736 } 2737 PetscCall(PetscFEPushforward(feI, fegeom, NbI, tmpBasisI)); 2738 PetscCall(PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI)); 2739 for (g = 0; g < NbJ; ++g) { 2740 for (gc = 0; gc < NcJ; ++gc) { 2741 const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */ 2742 2743 tmpBasisJ[gidx] = basisJ[gidx]; 2744 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx * dE + dg] = basisDerJ[gidx * dE + dg]; 2745 } 2746 } 2747 PetscCall(PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ)); 2748 // TODO This is currently broken since we do not pull the geometry down to the lower dimension 2749 // PetscCall(PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ)); 2750 for (f = 0; f < NbI; ++f) { 2751 for (fc = 0; fc < NcI; ++fc) { 2752 const PetscInt fidx = f * NcI + fc; /* Test function basis index */ 2753 const PetscInt i = offsetI + NbI * so + f; /* Element matrix row */ 2754 for (g = 0; g < NbJ; ++g) { 2755 for (gc = 0; gc < NcJ; ++gc) { 2756 const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */ 2757 const PetscInt j = offsetJ + NbJ * to + g; /* Element matrix column */ 2758 const PetscInt fOff = eOffset + i * totDim + j; 2759 2760 elemMat[fOff] += tmpBasisI[fidx] * g0[fc * NcJ + gc] * tmpBasisJ[gidx]; 2761 for (df = 0; df < dE; ++df) { 2762 elemMat[fOff] += tmpBasisI[fidx] * g1[(fc * NcJ + gc) * dE + df] * tmpBasisDerJ[gidx * dE + df]; 2763 elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g2[(fc * NcJ + gc) * dE + df] * tmpBasisJ[gidx]; 2764 for (dg = 0; dg < dE; ++dg) elemMat[fOff] += tmpBasisDerI[fidx * dE + df] * g3[((fc * NcJ + gc) * dE + df) * dE + dg] * tmpBasisDerJ[gidx * dE + dg]; 2765 } 2766 } 2767 } 2768 } 2769 } 2770 return PETSC_SUCCESS; 2771 } 2772 2773 PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom) 2774 { 2775 PetscDualSpace dsp; 2776 DM dm; 2777 PetscQuadrature quadDef; 2778 PetscInt dim, cdim, Nq; 2779 2780 PetscFunctionBegin; 2781 PetscCall(PetscFEGetDualSpace(fe, &dsp)); 2782 PetscCall(PetscDualSpaceGetDM(dsp, &dm)); 2783 PetscCall(DMGetDimension(dm, &dim)); 2784 PetscCall(DMGetCoordinateDim(dm, &cdim)); 2785 PetscCall(PetscFEGetQuadrature(fe, &quadDef)); 2786 quad = quad ? quad : quadDef; 2787 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL)); 2788 PetscCall(PetscMalloc1(Nq * cdim, &cgeom->v)); 2789 PetscCall(PetscMalloc1(Nq * cdim * cdim, &cgeom->J)); 2790 PetscCall(PetscMalloc1(Nq * cdim * cdim, &cgeom->invJ)); 2791 PetscCall(PetscMalloc1(Nq, &cgeom->detJ)); 2792 cgeom->dim = dim; 2793 cgeom->dimEmbed = cdim; 2794 cgeom->numCells = 1; 2795 cgeom->numPoints = Nq; 2796 PetscCall(DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ)); 2797 PetscFunctionReturn(PETSC_SUCCESS); 2798 } 2799 2800 PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom) 2801 { 2802 PetscFunctionBegin; 2803 PetscCall(PetscFree(cgeom->v)); 2804 PetscCall(PetscFree(cgeom->J)); 2805 PetscCall(PetscFree(cgeom->invJ)); 2806 PetscCall(PetscFree(cgeom->detJ)); 2807 PetscFunctionReturn(PETSC_SUCCESS); 2808 } 2809 2810 #if 0 2811 PetscErrorCode PetscFEUpdateElementMat_Internal_SparseIndices(PetscTabulation TI, PetscTabulation TJ, PetscInt dimEmbed, const PetscInt g0[], const PetscInt g1[], const PetscInt g2[], const PetscInt g3[], PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscInt *n_g0, PetscInt **g0_idxs_out, PetscInt *n_g1, PetscInt **g1_idxs_out, PetscInt *n_g2, PetscInt **g2_idxs_out, PetscInt *n_g3, PetscInt **g3_idxs_out) 2812 { 2813 const PetscInt dE = dimEmbed; 2814 const PetscInt NbI = TI->Nb; 2815 const PetscInt NcI = TI->Nc; 2816 const PetscInt NbJ = TJ->Nb; 2817 const PetscInt NcJ = TJ->Nc; 2818 PetscBool has_g0 = g0 ? PETSC_TRUE : PETSC_FALSE; 2819 PetscBool has_g1 = g1 ? PETSC_TRUE : PETSC_FALSE; 2820 PetscBool has_g2 = g2 ? PETSC_TRUE : PETSC_FALSE; 2821 PetscBool has_g3 = g3 ? PETSC_TRUE : PETSC_FALSE; 2822 PetscInt *g0_idxs = NULL, *g1_idxs = NULL, *g2_idxs = NULL, *g3_idxs = NULL; 2823 PetscInt g0_i, g1_i, g2_i, g3_i; 2824 2825 PetscFunctionBegin; 2826 g0_i = g1_i = g2_i = g3_i = 0; 2827 if (has_g0) 2828 for (PetscInt i = 0; i < NcI * NcJ; i++) 2829 if (g0[i]) g0_i += NbI * NbJ; 2830 if (has_g1) 2831 for (PetscInt i = 0; i < NcI * NcJ * dE; i++) 2832 if (g1[i]) g1_i += NbI * NbJ; 2833 if (has_g2) 2834 for (PetscInt i = 0; i < NcI * NcJ * dE; i++) 2835 if (g2[i]) g2_i += NbI * NbJ; 2836 if (has_g3) 2837 for (PetscInt i = 0; i < NcI * NcJ * dE * dE; i++) 2838 if (g3[i]) g3_i += NbI * NbJ; 2839 if (g0_i == NbI * NbJ * NcI * NcJ) g0_i = 0; 2840 if (g1_i == NbI * NbJ * NcI * NcJ * dE) g1_i = 0; 2841 if (g2_i == NbI * NbJ * NcI * NcJ * dE) g2_i = 0; 2842 if (g3_i == NbI * NbJ * NcI * NcJ * dE * dE) g3_i = 0; 2843 has_g0 = g0_i ? PETSC_TRUE : PETSC_FALSE; 2844 has_g1 = g1_i ? PETSC_TRUE : PETSC_FALSE; 2845 has_g2 = g2_i ? PETSC_TRUE : PETSC_FALSE; 2846 has_g3 = g3_i ? PETSC_TRUE : PETSC_FALSE; 2847 if (has_g0) PetscCall(PetscMalloc1(4 * g0_i, &g0_idxs)); 2848 if (has_g1) PetscCall(PetscMalloc1(4 * g1_i, &g1_idxs)); 2849 if (has_g2) PetscCall(PetscMalloc1(4 * g2_i, &g2_idxs)); 2850 if (has_g3) PetscCall(PetscMalloc1(4 * g3_i, &g3_idxs)); 2851 g0_i = g1_i = g2_i = g3_i = 0; 2852 2853 for (PetscInt f = 0; f < NbI; ++f) { 2854 const PetscInt i = offsetI + f; /* Element matrix row */ 2855 for (PetscInt fc = 0; fc < NcI; ++fc) { 2856 const PetscInt fidx = f * NcI + fc; /* Test function basis index */ 2857 2858 for (PetscInt g = 0; g < NbJ; ++g) { 2859 const PetscInt j = offsetJ + g; /* Element matrix column */ 2860 const PetscInt fOff = i * totDim + j; 2861 for (PetscInt gc = 0; gc < NcJ; ++gc) { 2862 const PetscInt gidx = g * NcJ + gc; /* Trial function basis index */ 2863 2864 if (has_g0) { 2865 if (g0[fc * NcJ + gc]) { 2866 g0_idxs[4 * g0_i + 0] = fidx; 2867 g0_idxs[4 * g0_i + 1] = fc * NcJ + gc; 2868 g0_idxs[4 * g0_i + 2] = gidx; 2869 g0_idxs[4 * g0_i + 3] = fOff; 2870 g0_i++; 2871 } 2872 } 2873 2874 for (PetscInt df = 0; df < dE; ++df) { 2875 if (has_g1) { 2876 if (g1[(fc * NcJ + gc) * dE + df]) { 2877 g1_idxs[4 * g1_i + 0] = fidx; 2878 g1_idxs[4 * g1_i + 1] = (fc * NcJ + gc) * dE + df; 2879 g1_idxs[4 * g1_i + 2] = gidx * dE + df; 2880 g1_idxs[4 * g1_i + 3] = fOff; 2881 g1_i++; 2882 } 2883 } 2884 if (has_g2) { 2885 if (g2[(fc * NcJ + gc) * dE + df]) { 2886 g2_idxs[4 * g2_i + 0] = fidx * dE + df; 2887 g2_idxs[4 * g2_i + 1] = (fc * NcJ + gc) * dE + df; 2888 g2_idxs[4 * g2_i + 2] = gidx; 2889 g2_idxs[4 * g2_i + 3] = fOff; 2890 g2_i++; 2891 } 2892 } 2893 if (has_g3) { 2894 for (PetscInt dg = 0; dg < dE; ++dg) { 2895 if (g3[((fc * NcJ + gc) * dE + df) * dE + dg]) { 2896 g3_idxs[4 * g3_i + 0] = fidx * dE + df; 2897 g3_idxs[4 * g3_i + 1] = ((fc * NcJ + gc) * dE + df) * dE + dg; 2898 g3_idxs[4 * g3_i + 2] = gidx * dE + dg; 2899 g3_idxs[4 * g3_i + 3] = fOff; 2900 g3_i++; 2901 } 2902 } 2903 } 2904 } 2905 } 2906 } 2907 } 2908 } 2909 *n_g0 = g0_i; 2910 *n_g1 = g1_i; 2911 *n_g2 = g2_i; 2912 *n_g3 = g3_i; 2913 2914 *g0_idxs_out = g0_idxs; 2915 *g1_idxs_out = g1_idxs; 2916 *g2_idxs_out = g2_idxs; 2917 *g3_idxs_out = g3_idxs; 2918 PetscFunctionReturn(PETSC_SUCCESS); 2919 } 2920 2921 //example HOW TO USE 2922 for (PetscInt i = 0; i < g0_sparse_n; i++) { 2923 PetscInt bM = g0_sparse_idxs[4 * i + 0]; 2924 PetscInt bN = g0_sparse_idxs[4 * i + 1]; 2925 PetscInt bK = g0_sparse_idxs[4 * i + 2]; 2926 PetscInt bO = g0_sparse_idxs[4 * i + 3]; 2927 elemMat[bO] += tmpBasisI[bM] * g0[bN] * tmpBasisJ[bK]; 2928 } 2929 #endif 2930