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