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