1 /* Basis Jet Tabulation 2 3 We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We 4 follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can 5 be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis 6 as a prime basis. 7 8 \psi_i = \sum_k \alpha_{ki} \phi_k 9 10 Our nodal basis is defined in terms of the dual basis $n_j$ 11 12 n_j \cdot \psi_i = \delta_{ji} 13 14 and we may act on the first equation to obtain 15 16 n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k 17 \delta_{ji} = \sum_k \alpha_{ki} V_{jk} 18 I = V \alpha 19 20 so the coefficients of the nodal basis in the prime basis are 21 22 \alpha = V^{-1} 23 24 We will define the dual basis vectors $n_j$ using a quadrature rule. 25 26 Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials 27 (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can 28 be implemented exactly as in FIAT using functionals $L_j$. 29 30 I will have to count the degrees correctly for the Legendre product when we are on simplices. 31 32 We will have three objects: 33 - Space, P: this just need point evaluation I think 34 - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q 35 - FEM: This keeps {P, P', Q} 36 */ 37 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 38 #include <petscdmplex.h> 39 40 PetscBool FEcite = PETSC_FALSE; 41 const char FECitation[] = "@article{kirby2004,\n" 42 " title = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n" 43 " journal = {ACM Transactions on Mathematical Software},\n" 44 " author = {Robert C. Kirby},\n" 45 " volume = {30},\n" 46 " number = {4},\n" 47 " pages = {502--516},\n" 48 " doi = {10.1145/1039813.1039820},\n" 49 " year = {2004}\n}\n"; 50 51 PetscClassId PETSCFE_CLASSID = 0; 52 53 PetscLogEvent PETSCFE_SetUp; 54 55 PetscFunctionList PetscFEList = NULL; 56 PetscBool PetscFERegisterAllCalled = PETSC_FALSE; 57 58 /*@C 59 PetscFERegister - Adds a new PetscFE implementation 60 61 Not Collective 62 63 Input Parameters: 64 + name - The name of a new user-defined creation routine 65 - create_func - The creation routine itself 66 67 Notes: 68 PetscFERegister() may be called multiple times to add several user-defined PetscFEs 69 70 Sample usage: 71 .vb 72 PetscFERegister("my_fe", MyPetscFECreate); 73 .ve 74 75 Then, your PetscFE type can be chosen with the procedural interface via 76 .vb 77 PetscFECreate(MPI_Comm, PetscFE *); 78 PetscFESetType(PetscFE, "my_fe"); 79 .ve 80 or at runtime via the option 81 .vb 82 -petscfe_type my_fe 83 .ve 84 85 Level: advanced 86 87 .seealso: PetscFERegisterAll(), PetscFERegisterDestroy() 88 89 @*/ 90 PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE)) 91 { 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 ierr = PetscFunctionListAdd(&PetscFEList, sname, function);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 99 /*@C 100 PetscFESetType - Builds a particular PetscFE 101 102 Collective on fem 103 104 Input Parameters: 105 + fem - The PetscFE object 106 - name - The kind of FEM space 107 108 Options Database Key: 109 . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types 110 111 Level: intermediate 112 113 .seealso: PetscFEGetType(), PetscFECreate() 114 @*/ 115 PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name) 116 { 117 PetscErrorCode (*r)(PetscFE); 118 PetscBool match; 119 PetscErrorCode ierr; 120 121 PetscFunctionBegin; 122 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 123 ierr = PetscObjectTypeCompare((PetscObject) fem, name, &match);CHKERRQ(ierr); 124 if (match) PetscFunctionReturn(0); 125 126 if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);} 127 ierr = PetscFunctionListFind(PetscFEList, name, &r);CHKERRQ(ierr); 128 if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name); 129 130 if (fem->ops->destroy) { 131 ierr = (*fem->ops->destroy)(fem);CHKERRQ(ierr); 132 fem->ops->destroy = NULL; 133 } 134 ierr = (*r)(fem);CHKERRQ(ierr); 135 ierr = PetscObjectChangeTypeName((PetscObject) fem, name);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 } 138 139 /*@C 140 PetscFEGetType - Gets the PetscFE type name (as a string) from the object. 141 142 Not Collective 143 144 Input Parameter: 145 . fem - The PetscFE 146 147 Output Parameter: 148 . name - The PetscFE type name 149 150 Level: intermediate 151 152 .seealso: PetscFESetType(), PetscFECreate() 153 @*/ 154 PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name) 155 { 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 160 PetscValidPointer(name, 2); 161 if (!PetscFERegisterAllCalled) { 162 ierr = PetscFERegisterAll();CHKERRQ(ierr); 163 } 164 *name = ((PetscObject) fem)->type_name; 165 PetscFunctionReturn(0); 166 } 167 168 /*@C 169 PetscFEViewFromOptions - View from Options 170 171 Collective on PetscFE 172 173 Input Parameters: 174 + A - the PetscFE object 175 . obj - Optional object 176 - name - command line option 177 178 Level: intermediate 179 .seealso: PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate() 180 @*/ 181 PetscErrorCode PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[]) 182 { 183 PetscErrorCode ierr; 184 185 PetscFunctionBegin; 186 PetscValidHeaderSpecific(A,PETSCFE_CLASSID,1); 187 ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 /*@C 192 PetscFEView - Views a PetscFE 193 194 Collective on fem 195 196 Input Parameter: 197 + fem - the PetscFE object to view 198 - viewer - the viewer 199 200 Level: beginner 201 202 .seealso PetscFEDestroy() 203 @*/ 204 PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer) 205 { 206 PetscBool iascii; 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 211 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 212 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);CHKERRQ(ierr);} 213 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);CHKERRQ(ierr); 214 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 215 if (fem->ops->view) {ierr = (*fem->ops->view)(fem, viewer);CHKERRQ(ierr);} 216 PetscFunctionReturn(0); 217 } 218 219 /*@ 220 PetscFESetFromOptions - sets parameters in a PetscFE from the options database 221 222 Collective on fem 223 224 Input Parameter: 225 . fem - the PetscFE object to set options for 226 227 Options Database: 228 + -petscfe_num_blocks - the number of cell blocks to integrate concurrently 229 - -petscfe_num_batches - the number of cell batches to integrate serially 230 231 Level: intermediate 232 233 .seealso PetscFEView() 234 @*/ 235 PetscErrorCode PetscFESetFromOptions(PetscFE fem) 236 { 237 const char *defaultType; 238 char name[256]; 239 PetscBool flg; 240 PetscErrorCode ierr; 241 242 PetscFunctionBegin; 243 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 244 if (!((PetscObject) fem)->type_name) { 245 defaultType = PETSCFEBASIC; 246 } else { 247 defaultType = ((PetscObject) fem)->type_name; 248 } 249 if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);} 250 251 ierr = PetscObjectOptionsBegin((PetscObject) fem);CHKERRQ(ierr); 252 ierr = PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);CHKERRQ(ierr); 253 if (flg) { 254 ierr = PetscFESetType(fem, name);CHKERRQ(ierr); 255 } else if (!((PetscObject) fem)->type_name) { 256 ierr = PetscFESetType(fem, defaultType);CHKERRQ(ierr); 257 } 258 ierr = PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);CHKERRQ(ierr); 259 ierr = PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);CHKERRQ(ierr); 260 if (fem->ops->setfromoptions) { 261 ierr = (*fem->ops->setfromoptions)(PetscOptionsObject,fem);CHKERRQ(ierr); 262 } 263 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 264 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);CHKERRQ(ierr); 265 ierr = PetscOptionsEnd();CHKERRQ(ierr); 266 ierr = PetscFEViewFromOptions(fem, NULL, "-petscfe_view");CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 /*@C 271 PetscFESetUp - Construct data structures for the PetscFE 272 273 Collective on fem 274 275 Input Parameter: 276 . fem - the PetscFE object to setup 277 278 Level: intermediate 279 280 .seealso PetscFEView(), PetscFEDestroy() 281 @*/ 282 PetscErrorCode PetscFESetUp(PetscFE fem) 283 { 284 PetscErrorCode ierr; 285 286 PetscFunctionBegin; 287 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 288 if (fem->setupcalled) PetscFunctionReturn(0); 289 ierr = PetscLogEventBegin(PETSCFE_SetUp, fem, 0, 0, 0);CHKERRQ(ierr); 290 fem->setupcalled = PETSC_TRUE; 291 if (fem->ops->setup) {ierr = (*fem->ops->setup)(fem);CHKERRQ(ierr);} 292 ierr = PetscLogEventEnd(PETSCFE_SetUp, fem, 0, 0, 0);CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 /*@ 297 PetscFEDestroy - Destroys a PetscFE object 298 299 Collective on fem 300 301 Input Parameter: 302 . fem - the PetscFE object to destroy 303 304 Level: beginner 305 306 .seealso PetscFEView() 307 @*/ 308 PetscErrorCode PetscFEDestroy(PetscFE *fem) 309 { 310 PetscErrorCode ierr; 311 312 PetscFunctionBegin; 313 if (!*fem) PetscFunctionReturn(0); 314 PetscValidHeaderSpecific((*fem), PETSCFE_CLASSID, 1); 315 316 if (--((PetscObject)(*fem))->refct > 0) {*fem = NULL; PetscFunctionReturn(0);} 317 ((PetscObject) (*fem))->refct = 0; 318 319 if ((*fem)->subspaces) { 320 PetscInt dim, d; 321 322 ierr = PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);CHKERRQ(ierr); 323 for (d = 0; d < dim; ++d) {ierr = PetscFEDestroy(&(*fem)->subspaces[d]);CHKERRQ(ierr);} 324 } 325 ierr = PetscFree((*fem)->subspaces);CHKERRQ(ierr); 326 ierr = PetscFree((*fem)->invV);CHKERRQ(ierr); 327 ierr = PetscTabulationDestroy(&(*fem)->T);CHKERRQ(ierr); 328 ierr = PetscTabulationDestroy(&(*fem)->Tf);CHKERRQ(ierr); 329 ierr = PetscTabulationDestroy(&(*fem)->Tc);CHKERRQ(ierr); 330 ierr = PetscSpaceDestroy(&(*fem)->basisSpace);CHKERRQ(ierr); 331 ierr = PetscDualSpaceDestroy(&(*fem)->dualSpace);CHKERRQ(ierr); 332 ierr = PetscQuadratureDestroy(&(*fem)->quadrature);CHKERRQ(ierr); 333 ierr = PetscQuadratureDestroy(&(*fem)->faceQuadrature);CHKERRQ(ierr); 334 #ifdef PETSC_HAVE_LIBCEED 335 ierr = CeedBasisDestroy(&(*fem)->ceedBasis);CHKERRQ(ierr); 336 ierr = CeedDestroy(&(*fem)->ceed);CHKERRQ(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 if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc); 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 if ((qNc != 1) && (Nc != qNc)) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_SIZ, "FE components %D != Quadrature components %D and non-scalar quadrature", Nc, qNc); 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 Parameter: 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 if (fem->T && k > fem->T->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->T->K); 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 Parameter: 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 derviatives 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 if (fem->Tf && k > fem->Tf->K) SETERRQ2(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_OUTOFRANGE, "Requested %D derivatives, but only tabulated %D", k, fem->Tf->K); 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 if (T->K != (!cdim ? 0 : K)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation K %D must match requested K %D", T->K, !cdim ? 0 : K); 1040 if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb); 1041 if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc); 1042 if (T->cdim != cdim) SETERRQ2(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 /*@ 1161 PetscFEGetDimension - Get the dimension of the finite element space on a cell 1162 1163 Not collective 1164 1165 Input Parameter: 1166 . fe - The PetscFE 1167 1168 Output Parameter: 1169 . dim - The dimension 1170 1171 Level: intermediate 1172 1173 .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension() 1174 @*/ 1175 PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) 1176 { 1177 PetscErrorCode ierr; 1178 1179 PetscFunctionBegin; 1180 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1181 PetscValidPointer(dim, 2); 1182 if (fem->ops->getdimension) {ierr = (*fem->ops->getdimension)(fem, dim);CHKERRQ(ierr);} 1183 PetscFunctionReturn(0); 1184 } 1185 1186 /*@C 1187 PetscFEPushforward - Map the reference element function to real space 1188 1189 Input Parameters: 1190 + fe - The PetscFE 1191 . fegeom - The cell geometry 1192 . Nv - The number of function values 1193 - vals - The function values 1194 1195 Output Parameter: 1196 . vals - The transformed function values 1197 1198 Level: advanced 1199 1200 Note: This just forwards the call onto PetscDualSpacePushforward(). 1201 1202 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1203 1204 .seealso: PetscDualSpacePushforward() 1205 @*/ 1206 PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1207 { 1208 PetscErrorCode ierr; 1209 1210 PetscFunctionBeginHot; 1211 ierr = PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1212 PetscFunctionReturn(0); 1213 } 1214 1215 /*@C 1216 PetscFEPushforwardGradient - Map the reference element function gradient to real space 1217 1218 Input Parameters: 1219 + fe - The PetscFE 1220 . fegeom - The cell geometry 1221 . Nv - The number of function gradient values 1222 - vals - The function gradient values 1223 1224 Output Parameter: 1225 . vals - The transformed function gradient values 1226 1227 Level: advanced 1228 1229 Note: This just forwards the call onto PetscDualSpacePushforwardGradient(). 1230 1231 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1232 1233 .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward() 1234 @*/ 1235 PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1236 { 1237 PetscErrorCode ierr; 1238 1239 PetscFunctionBeginHot; 1240 ierr = PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1241 PetscFunctionReturn(0); 1242 } 1243 1244 /*@C 1245 PetscFEPushforwardHessian - Map the reference element function Hessian to real space 1246 1247 Input Parameters: 1248 + fe - The PetscFE 1249 . fegeom - The cell geometry 1250 . Nv - The number of function Hessian values 1251 - vals - The function Hessian values 1252 1253 Output Parameter: 1254 . vals - The transformed function Hessian values 1255 1256 Level: advanced 1257 1258 Note: This just forwards the call onto PetscDualSpacePushforwardHessian(). 1259 1260 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1261 1262 .seealso: PetscFEPushforward(), PetscDualSpacePushforwardHessian(), PetscDualSpacePushforward() 1263 @*/ 1264 PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1265 { 1266 PetscErrorCode ierr; 1267 1268 PetscFunctionBeginHot; 1269 ierr = PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1270 PetscFunctionReturn(0); 1271 } 1272 1273 /* 1274 Purpose: Compute element vector for chunk of elements 1275 1276 Input: 1277 Sizes: 1278 Ne: number of elements 1279 Nf: number of fields 1280 PetscFE 1281 dim: spatial dimension 1282 Nb: number of basis functions 1283 Nc: number of field components 1284 PetscQuadrature 1285 Nq: number of quadrature points 1286 1287 Geometry: 1288 PetscFEGeom[Ne] possibly *Nq 1289 PetscReal v0s[dim] 1290 PetscReal n[dim] 1291 PetscReal jacobians[dim*dim] 1292 PetscReal jacobianInverses[dim*dim] 1293 PetscReal jacobianDeterminants 1294 FEM: 1295 PetscFE 1296 PetscQuadrature 1297 PetscReal quadPoints[Nq*dim] 1298 PetscReal quadWeights[Nq] 1299 PetscReal basis[Nq*Nb*Nc] 1300 PetscReal basisDer[Nq*Nb*Nc*dim] 1301 PetscScalar coefficients[Ne*Nb*Nc] 1302 PetscScalar elemVec[Ne*Nb*Nc] 1303 1304 Problem: 1305 PetscInt f: the active field 1306 f0, f1 1307 1308 Work Space: 1309 PetscFE 1310 PetscScalar f0[Nq*dim]; 1311 PetscScalar f1[Nq*dim*dim]; 1312 PetscScalar u[Nc]; 1313 PetscScalar gradU[Nc*dim]; 1314 PetscReal x[dim]; 1315 PetscScalar realSpaceDer[dim]; 1316 1317 Purpose: Compute element vector for N_cb batches of elements 1318 1319 Input: 1320 Sizes: 1321 N_cb: Number of serial cell batches 1322 1323 Geometry: 1324 PetscReal v0s[Ne*dim] 1325 PetscReal jacobians[Ne*dim*dim] possibly *Nq 1326 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 1327 PetscReal jacobianDeterminants[Ne] possibly *Nq 1328 FEM: 1329 static PetscReal quadPoints[Nq*dim] 1330 static PetscReal quadWeights[Nq] 1331 static PetscReal basis[Nq*Nb*Nc] 1332 static PetscReal basisDer[Nq*Nb*Nc*dim] 1333 PetscScalar coefficients[Ne*Nb*Nc] 1334 PetscScalar elemVec[Ne*Nb*Nc] 1335 1336 ex62.c: 1337 PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], 1338 const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], 1339 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), 1340 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) 1341 1342 ex52.c: 1343 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) 1344 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) 1345 1346 ex52_integrateElement.cu 1347 __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) 1348 1349 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], 1350 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1351 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1352 1353 ex52_integrateElementOpenCL.c: 1354 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], 1355 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1356 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1357 1358 __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) 1359 */ 1360 1361 /*@C 1362 PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration 1363 1364 Not collective 1365 1366 Input Parameters: 1367 + prob - The PetscDS specifying the discretizations and continuum functions 1368 . field - The field being integrated 1369 . Ne - The number of elements in the chunk 1370 . cgeom - The cell geometry for each cell in the chunk 1371 . coefficients - The array of FEM basis coefficients for the elements 1372 . probAux - The PetscDS specifying the auxiliary discretizations 1373 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1374 1375 Output Parameter: 1376 . integral - the integral for this field 1377 1378 Level: intermediate 1379 1380 .seealso: PetscFEIntegrateResidual() 1381 @*/ 1382 PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 1383 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1384 { 1385 PetscFE fe; 1386 PetscErrorCode ierr; 1387 1388 PetscFunctionBegin; 1389 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1390 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1391 if (fe->ops->integrate) {ierr = (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} 1392 PetscFunctionReturn(0); 1393 } 1394 1395 /*@C 1396 PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration 1397 1398 Not collective 1399 1400 Input Parameters: 1401 + prob - The PetscDS specifying the discretizations and continuum functions 1402 . field - The field being integrated 1403 . obj_func - The function to be integrated 1404 . Ne - The number of elements in the chunk 1405 . fgeom - The face geometry for each face in the chunk 1406 . coefficients - The array of FEM basis coefficients for the elements 1407 . probAux - The PetscDS specifying the auxiliary discretizations 1408 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1409 1410 Output Parameter: 1411 . integral - the integral for this field 1412 1413 Level: intermediate 1414 1415 .seealso: PetscFEIntegrateResidual() 1416 @*/ 1417 PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field, 1418 void (*obj_func)(PetscInt, PetscInt, PetscInt, 1419 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1420 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1421 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 1422 PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1423 { 1424 PetscFE fe; 1425 PetscErrorCode ierr; 1426 1427 PetscFunctionBegin; 1428 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1429 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1430 if (fe->ops->integratebd) {ierr = (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} 1431 PetscFunctionReturn(0); 1432 } 1433 1434 /*@C 1435 PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration 1436 1437 Not collective 1438 1439 Input Parameters: 1440 + ds - The PetscDS specifying the discretizations and continuum functions 1441 . key - The (label+value, field) being integrated 1442 . Ne - The number of elements in the chunk 1443 . cgeom - The cell geometry for each cell in the chunk 1444 . coefficients - The array of FEM basis coefficients for the elements 1445 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1446 . probAux - The PetscDS specifying the auxiliary discretizations 1447 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1448 - t - The time 1449 1450 Output Parameter: 1451 . elemVec - the element residual vectors from each element 1452 1453 Note: 1454 $ Loop over batch of elements (e): 1455 $ Loop over quadrature points (q): 1456 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q 1457 $ Call f_0 and f_1 1458 $ Loop over element vector entries (f,fc --> i): 1459 $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) 1460 1461 Level: intermediate 1462 1463 .seealso: PetscFEIntegrateResidual() 1464 @*/ 1465 PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *cgeom, 1466 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1467 { 1468 PetscFE fe; 1469 PetscErrorCode ierr; 1470 1471 PetscFunctionBeginHot; 1472 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1473 ierr = PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);CHKERRQ(ierr); 1474 if (fe->ops->integrateresidual) {ierr = (*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1475 PetscFunctionReturn(0); 1476 } 1477 1478 /*@C 1479 PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary 1480 1481 Not collective 1482 1483 Input Parameters: 1484 + ds - The PetscDS specifying the discretizations and continuum functions 1485 . wf - The PetscWeakForm object holding the pointwise functions 1486 . key - The (label+value, field) being integrated 1487 . Ne - The number of elements in the chunk 1488 . fgeom - The face geometry for each cell in the chunk 1489 . coefficients - The array of FEM basis coefficients for the elements 1490 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1491 . probAux - The PetscDS specifying the auxiliary discretizations 1492 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1493 - t - The time 1494 1495 Output Parameter: 1496 . elemVec - the element residual vectors from each element 1497 1498 Level: intermediate 1499 1500 .seealso: PetscFEIntegrateResidual() 1501 @*/ 1502 PetscErrorCode PetscFEIntegrateBdResidual(PetscDS ds, PetscWeakForm wf, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1503 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1504 { 1505 PetscFE fe; 1506 PetscErrorCode ierr; 1507 1508 PetscFunctionBegin; 1509 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1510 ierr = PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);CHKERRQ(ierr); 1511 if (fe->ops->integratebdresidual) {ierr = (*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1512 PetscFunctionReturn(0); 1513 } 1514 1515 /*@C 1516 PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration 1517 1518 Not collective 1519 1520 Input Parameters: 1521 + prob - The PetscDS specifying the discretizations and continuum functions 1522 . key - The (label+value, field) being integrated 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, PetscHashFormKey key, 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, 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 . Ne - The number of elements in the chunk 1561 . cgeom - The cell geometry for each cell in the chunk 1562 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1563 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1564 . probAux - The PetscDS specifying the auxiliary discretizations 1565 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1566 . t - The time 1567 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1568 1569 Output Parameter: 1570 . elemMat - the element matrices for the Jacobian from each element 1571 1572 Note: 1573 $ Loop over batch of elements (e): 1574 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1575 $ Loop over quadrature points (q): 1576 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1577 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1578 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1579 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1580 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1581 Level: intermediate 1582 1583 .seealso: PetscFEIntegrateResidual() 1584 @*/ 1585 PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *cgeom, 1586 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1587 { 1588 PetscFE fe; 1589 PetscInt Nf; 1590 PetscErrorCode ierr; 1591 1592 PetscFunctionBegin; 1593 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1594 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1595 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1596 if (fe->ops->integratejacobian) {ierr = (*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1597 PetscFunctionReturn(0); 1598 } 1599 1600 /*@C 1601 PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration 1602 1603 Not collective 1604 1605 Input Parameters: 1606 + ds - The PetscDS specifying the discretizations and continuum functions 1607 . wf - The PetscWeakForm holding the pointwise functions 1608 . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1609 . Ne - The number of elements in the chunk 1610 . fgeom - The face geometry for each cell in the chunk 1611 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1612 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1613 . probAux - The PetscDS specifying the auxiliary discretizations 1614 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1615 . t - The time 1616 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1617 1618 Output Parameter: 1619 . elemMat - the element matrices for the Jacobian from each element 1620 1621 Note: 1622 $ Loop over batch of elements (e): 1623 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1624 $ Loop over quadrature points (q): 1625 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1626 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1627 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1628 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1629 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1630 Level: intermediate 1631 1632 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1633 @*/ 1634 PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1635 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1636 { 1637 PetscFE fe; 1638 PetscInt Nf; 1639 PetscErrorCode ierr; 1640 1641 PetscFunctionBegin; 1642 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1643 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1644 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1645 if (fe->ops->integratebdjacobian) {ierr = (*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1646 PetscFunctionReturn(0); 1647 } 1648 1649 /*@C 1650 PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration 1651 1652 Not collective 1653 1654 Input Parameters: 1655 + ds - The PetscDS specifying the discretizations and continuum functions 1656 . jtype - The type of matrix pointwise functions that should be used 1657 . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1658 . Ne - The number of elements in the chunk 1659 . fgeom - The face geometry for each cell in the chunk 1660 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1661 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1662 . probAux - The PetscDS specifying the auxiliary discretizations 1663 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1664 . t - The time 1665 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1666 1667 Output Parameter 1668 . elemMat - the element matrices for the Jacobian from each element 1669 1670 Note: 1671 $ Loop over batch of elements (e): 1672 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1673 $ Loop over quadrature points (q): 1674 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1675 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1676 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1677 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1678 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1679 Level: developer 1680 1681 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1682 @*/ 1683 PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1684 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1685 { 1686 PetscFE fe; 1687 PetscInt Nf; 1688 PetscErrorCode ierr; 1689 1690 PetscFunctionBegin; 1691 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1692 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1693 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1694 if (fe->ops->integratehybridjacobian) {ierr = (*fe->ops->integratehybridjacobian)(ds, jtype, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1695 PetscFunctionReturn(0); 1696 } 1697 1698 /*@ 1699 PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height 1700 1701 Input Parameters: 1702 + fe - The finite element space 1703 - height - The height of the Plex point 1704 1705 Output Parameter: 1706 . subfe - The subspace of this FE space 1707 1708 Note: For example, if we want the subspace of this space for a face, we would choose height = 1. 1709 1710 Level: advanced 1711 1712 .seealso: PetscFECreateDefault() 1713 @*/ 1714 PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe) 1715 { 1716 PetscSpace P, subP; 1717 PetscDualSpace Q, subQ; 1718 PetscQuadrature subq; 1719 PetscFEType fetype; 1720 PetscInt dim, Nc; 1721 PetscErrorCode ierr; 1722 1723 PetscFunctionBegin; 1724 PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1); 1725 PetscValidPointer(subfe, 3); 1726 if (height == 0) { 1727 *subfe = fe; 1728 PetscFunctionReturn(0); 1729 } 1730 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1731 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1732 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1733 ierr = PetscFEGetFaceQuadrature(fe, &subq);CHKERRQ(ierr); 1734 ierr = PetscDualSpaceGetDimension(Q, &dim);CHKERRQ(ierr); 1735 if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim); 1736 if (!fe->subspaces) {ierr = PetscCalloc1(dim, &fe->subspaces);CHKERRQ(ierr);} 1737 if (height <= dim) { 1738 if (!fe->subspaces[height-1]) { 1739 PetscFE sub = NULL; 1740 const char *name; 1741 1742 ierr = PetscSpaceGetHeightSubspace(P, height, &subP);CHKERRQ(ierr); 1743 ierr = PetscDualSpaceGetHeightSubspace(Q, height, &subQ);CHKERRQ(ierr); 1744 if (subQ) { 1745 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);CHKERRQ(ierr); 1746 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 1747 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 1748 ierr = PetscFEGetType(fe, &fetype);CHKERRQ(ierr); 1749 ierr = PetscFESetType(sub, fetype);CHKERRQ(ierr); 1750 ierr = PetscFESetBasisSpace(sub, subP);CHKERRQ(ierr); 1751 ierr = PetscFESetDualSpace(sub, subQ);CHKERRQ(ierr); 1752 ierr = PetscFESetNumComponents(sub, Nc);CHKERRQ(ierr); 1753 ierr = PetscFESetUp(sub);CHKERRQ(ierr); 1754 ierr = PetscFESetQuadrature(sub, subq);CHKERRQ(ierr); 1755 } 1756 fe->subspaces[height-1] = sub; 1757 } 1758 *subfe = fe->subspaces[height-1]; 1759 } else { 1760 *subfe = NULL; 1761 } 1762 PetscFunctionReturn(0); 1763 } 1764 1765 /*@ 1766 PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used 1767 to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more 1768 sparsity). It is also used to create an interpolation between regularly refined meshes. 1769 1770 Collective on fem 1771 1772 Input Parameter: 1773 . fe - The initial PetscFE 1774 1775 Output Parameter: 1776 . feRef - The refined PetscFE 1777 1778 Level: advanced 1779 1780 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 1781 @*/ 1782 PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef) 1783 { 1784 PetscSpace P, Pref; 1785 PetscDualSpace Q, Qref; 1786 DM K, Kref; 1787 PetscQuadrature q, qref; 1788 const PetscReal *v0, *jac; 1789 PetscInt numComp, numSubelements; 1790 PetscInt cStart, cEnd, c; 1791 PetscDualSpace *cellSpaces; 1792 PetscErrorCode ierr; 1793 1794 PetscFunctionBegin; 1795 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1796 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1797 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1798 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr); 1799 /* Create space */ 1800 ierr = PetscObjectReference((PetscObject) P);CHKERRQ(ierr); 1801 Pref = P; 1802 /* Create dual space */ 1803 ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr); 1804 ierr = PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);CHKERRQ(ierr); 1805 ierr = DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);CHKERRQ(ierr); 1806 ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr); 1807 ierr = DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);CHKERRQ(ierr); 1808 ierr = PetscMalloc1(cEnd - cStart, &cellSpaces);CHKERRQ(ierr); 1809 /* TODO: fix for non-uniform refinement */ 1810 for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q; 1811 ierr = PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);CHKERRQ(ierr); 1812 ierr = PetscFree(cellSpaces);CHKERRQ(ierr); 1813 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 1814 ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr); 1815 /* Create element */ 1816 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);CHKERRQ(ierr); 1817 ierr = PetscFESetType(*feRef, PETSCFECOMPOSITE);CHKERRQ(ierr); 1818 ierr = PetscFESetBasisSpace(*feRef, Pref);CHKERRQ(ierr); 1819 ierr = PetscFESetDualSpace(*feRef, Qref);CHKERRQ(ierr); 1820 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 1821 ierr = PetscFESetNumComponents(*feRef, numComp);CHKERRQ(ierr); 1822 ierr = PetscFESetUp(*feRef);CHKERRQ(ierr); 1823 ierr = PetscSpaceDestroy(&Pref);CHKERRQ(ierr); 1824 ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr); 1825 /* Create quadrature */ 1826 ierr = PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);CHKERRQ(ierr); 1827 ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr); 1828 ierr = PetscFESetQuadrature(*feRef, qref);CHKERRQ(ierr); 1829 ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr); 1830 PetscFunctionReturn(0); 1831 } 1832 1833 /*@C 1834 PetscFECreateDefault - Create a PetscFE for basic FEM computation 1835 1836 Collective 1837 1838 Input Parameters: 1839 + comm - The MPI comm 1840 . dim - The spatial dimension 1841 . Nc - The number of components 1842 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1843 . prefix - The options prefix, or NULL 1844 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1845 1846 Output Parameter: 1847 . fem - The PetscFE object 1848 1849 Note: 1850 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. 1851 1852 Level: beginner 1853 1854 .seealso: PetscSpaceSetFromOptions(), PetscDualSpaceSetFromOptions(), PetscFESetFromOptions(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1855 @*/ 1856 PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem) 1857 { 1858 PetscQuadrature q, fq; 1859 DM K; 1860 PetscSpace P; 1861 PetscDualSpace Q; 1862 PetscInt order, quadPointsPerEdge; 1863 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1864 PetscErrorCode ierr; 1865 1866 PetscFunctionBegin; 1867 /* Create space */ 1868 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1869 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr); 1870 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1871 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1872 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1873 ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr); 1874 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1875 ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr); 1876 ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr); 1877 /* Create dual space */ 1878 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1879 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1880 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr); 1881 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1882 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1883 ierr = DMDestroy(&K);CHKERRQ(ierr); 1884 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1885 ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr); 1886 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1887 ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr); 1888 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1889 /* Create element */ 1890 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1891 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr); 1892 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1893 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1894 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1895 ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr); 1896 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1897 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1898 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1899 /* Create quadrature (with specified order if given) */ 1900 qorder = qorder >= 0 ? qorder : order; 1901 ierr = PetscObjectOptionsBegin((PetscObject)*fem);CHKERRQ(ierr); 1902 ierr = PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);CHKERRQ(ierr); 1903 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1904 quadPointsPerEdge = PetscMax(qorder + 1,1); 1905 if (isSimplex) { 1906 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1907 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1908 } else { 1909 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1910 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1911 } 1912 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1913 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1914 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1915 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1916 PetscFunctionReturn(0); 1917 } 1918 1919 /*@ 1920 PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k 1921 1922 Collective 1923 1924 Input Parameters: 1925 + comm - The MPI comm 1926 . dim - The spatial dimension 1927 . Nc - The number of components 1928 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1929 . k - The degree k of the space 1930 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1931 1932 Output Parameter: 1933 . fem - The PetscFE object 1934 1935 Level: beginner 1936 1937 Notes: 1938 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. 1939 1940 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1941 @*/ 1942 PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem) 1943 { 1944 PetscQuadrature q, fq; 1945 DM K; 1946 PetscSpace P; 1947 PetscDualSpace Q; 1948 PetscInt quadPointsPerEdge; 1949 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1950 char name[64]; 1951 PetscErrorCode ierr; 1952 1953 PetscFunctionBegin; 1954 /* Create space */ 1955 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1956 ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1957 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1958 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1959 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1960 ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr); 1961 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1962 /* Create dual space */ 1963 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1964 ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1965 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1966 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1967 ierr = DMDestroy(&K);CHKERRQ(ierr); 1968 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1969 ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr); 1970 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1971 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1972 /* Create finite element */ 1973 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1974 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1975 ierr = PetscObjectSetName((PetscObject) *fem, name);CHKERRQ(ierr); 1976 ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr); 1977 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1978 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1979 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1980 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1981 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1982 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1983 /* Create quadrature (with specified order if given) */ 1984 qorder = qorder >= 0 ? qorder : k; 1985 quadPointsPerEdge = PetscMax(qorder + 1,1); 1986 if (isSimplex) { 1987 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1988 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1989 } else { 1990 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1991 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1992 } 1993 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1994 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1995 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1996 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1997 /* Set finite element name */ 1998 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1999 ierr = PetscFESetName(*fem, name);CHKERRQ(ierr); 2000 PetscFunctionReturn(0); 2001 } 2002 2003 /*@C 2004 PetscFESetName - Names the FE and its subobjects 2005 2006 Not collective 2007 2008 Input Parameters: 2009 + fe - The PetscFE 2010 - name - The name 2011 2012 Level: intermediate 2013 2014 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 2015 @*/ 2016 PetscErrorCode PetscFESetName(PetscFE fe, const char name[]) 2017 { 2018 PetscSpace P; 2019 PetscDualSpace Q; 2020 PetscErrorCode ierr; 2021 2022 PetscFunctionBegin; 2023 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 2024 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2025 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 2026 ierr = PetscObjectSetName((PetscObject) P, name);CHKERRQ(ierr); 2027 ierr = PetscObjectSetName((PetscObject) Q, name);CHKERRQ(ierr); 2028 PetscFunctionReturn(0); 2029 } 2030 2031 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[]) 2032 { 2033 PetscInt dOffset = 0, fOffset = 0, f, g; 2034 PetscErrorCode ierr; 2035 2036 for (f = 0; f < Nf; ++f) { 2037 PetscFE fe; 2038 const PetscInt k = ds->jetDegree[f]; 2039 const PetscInt cdim = T[f]->cdim; 2040 const PetscInt Nq = T[f]->Np; 2041 const PetscInt Nbf = T[f]->Nb; 2042 const PetscInt Ncf = T[f]->Nc; 2043 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2044 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2045 const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r*Nq+q)*Nbf*Ncf*cdim*cdim] : NULL; 2046 PetscInt hOffset = 0, b, c, d; 2047 2048 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 2049 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2050 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2051 for (b = 0; b < Nbf; ++b) { 2052 for (c = 0; c < Ncf; ++c) { 2053 const PetscInt cidx = b*Ncf+c; 2054 2055 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2056 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2057 } 2058 } 2059 if (k > 1) { 2060 for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc*cdim; 2061 for (d = 0; d < cdim*cdim*Ncf; ++d) u_x[hOffset+fOffset*cdim*cdim+d] = 0.0; 2062 for (b = 0; b < Nbf; ++b) { 2063 for (c = 0; c < Ncf; ++c) { 2064 const PetscInt cidx = b*Ncf+c; 2065 2066 for (d = 0; d < cdim*cdim; ++d) u_x[hOffset+(fOffset+c)*cdim*cdim+d] += Hq[cidx*cdim*cdim+d]*coefficients[dOffset+b]; 2067 } 2068 } 2069 ierr = PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset+fOffset*cdim*cdim]);CHKERRQ(ierr); 2070 } 2071 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2072 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2073 if (u_t) { 2074 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2075 for (b = 0; b < Nbf; ++b) { 2076 for (c = 0; c < Ncf; ++c) { 2077 const PetscInt cidx = b*Ncf+c; 2078 2079 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2080 } 2081 } 2082 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2083 } 2084 fOffset += Ncf; 2085 dOffset += Nbf; 2086 } 2087 return 0; 2088 } 2089 2090 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[]) 2091 { 2092 PetscInt dOffset = 0, fOffset = 0, g; 2093 PetscErrorCode ierr; 2094 2095 for (g = 0; g < 2*Nf-1; ++g) { 2096 if (!T[g/2]) continue; 2097 { 2098 PetscFE fe; 2099 const PetscInt f = g/2; 2100 const PetscInt cdim = T[f]->cdim; 2101 const PetscInt Nq = T[f]->Np; 2102 const PetscInt Nbf = T[f]->Nb; 2103 const PetscInt Ncf = T[f]->Nc; 2104 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2105 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2106 PetscInt b, c, d; 2107 2108 fe = (PetscFE) ds->disc[f]; 2109 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2110 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2111 for (b = 0; b < Nbf; ++b) { 2112 for (c = 0; c < Ncf; ++c) { 2113 const PetscInt cidx = b*Ncf+c; 2114 2115 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2116 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2117 } 2118 } 2119 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2120 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2121 if (u_t) { 2122 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2123 for (b = 0; b < Nbf; ++b) { 2124 for (c = 0; c < Ncf; ++c) { 2125 const PetscInt cidx = b*Ncf+c; 2126 2127 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2128 } 2129 } 2130 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2131 } 2132 fOffset += Ncf; 2133 dOffset += Nbf; 2134 } 2135 } 2136 return 0; 2137 } 2138 2139 PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[]) 2140 { 2141 PetscFE fe; 2142 PetscTabulation Tc; 2143 PetscInt b, c; 2144 PetscErrorCode ierr; 2145 2146 if (!prob) return 0; 2147 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 2148 ierr = PetscFEGetFaceCentroidTabulation(fe, &Tc);CHKERRQ(ierr); 2149 { 2150 const PetscReal *faceBasis = Tc->T[0]; 2151 const PetscInt Nb = Tc->Nb; 2152 const PetscInt Nc = Tc->Nc; 2153 2154 for (c = 0; c < Nc; ++c) {u[c] = 0.0;} 2155 for (b = 0; b < Nb; ++b) { 2156 for (c = 0; c < Nc; ++c) { 2157 u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c]; 2158 } 2159 } 2160 } 2161 return 0; 2162 } 2163 2164 PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2165 { 2166 const PetscInt dE = T->cdim; /* fegeom->dimEmbed */ 2167 const PetscInt Nq = T->Np; 2168 const PetscInt Nb = T->Nb; 2169 const PetscInt Nc = T->Nc; 2170 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2171 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2172 PetscInt q, b, c, d; 2173 PetscErrorCode ierr; 2174 2175 for (b = 0; b < Nb; ++b) elemVec[b] = 0.0; 2176 for (q = 0; q < Nq; ++q) { 2177 for (b = 0; b < Nb; ++b) { 2178 for (c = 0; c < Nc; ++c) { 2179 const PetscInt bcidx = b*Nc+c; 2180 2181 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2182 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2183 } 2184 } 2185 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2186 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2187 for (b = 0; b < Nb; ++b) { 2188 for (c = 0; c < Nc; ++c) { 2189 const PetscInt bcidx = b*Nc+c; 2190 const PetscInt qcidx = q*Nc+c; 2191 2192 elemVec[b] += tmpBasis[bcidx]*f0[qcidx]; 2193 for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2194 } 2195 } 2196 } 2197 return(0); 2198 } 2199 2200 PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2201 { 2202 const PetscInt dE = T->cdim; 2203 const PetscInt Nq = T->Np; 2204 const PetscInt Nb = T->Nb; 2205 const PetscInt Nc = T->Nc; 2206 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2207 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2208 PetscInt q, b, c, d, s; 2209 PetscErrorCode ierr; 2210 2211 for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0; 2212 for (q = 0; q < Nq; ++q) { 2213 for (b = 0; b < Nb; ++b) { 2214 for (c = 0; c < Nc; ++c) { 2215 const PetscInt bcidx = b*Nc+c; 2216 2217 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2218 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2219 } 2220 } 2221 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2222 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2223 for (s = 0; s < 2; ++s) { 2224 for (b = 0; b < Nb; ++b) { 2225 for (c = 0; c < Nc; ++c) { 2226 const PetscInt bcidx = b*Nc+c; 2227 const PetscInt qcidx = (q*2+s)*Nc+c; 2228 2229 elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx]; 2230 for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2231 } 2232 } 2233 } 2234 } 2235 return(0); 2236 } 2237 2238 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[]) 2239 { 2240 const PetscInt dE = TI->cdim; 2241 const PetscInt NqI = TI->Np; 2242 const PetscInt NbI = TI->Nb; 2243 const PetscInt NcI = TI->Nc; 2244 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2245 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2246 const PetscInt NqJ = TJ->Np; 2247 const PetscInt NbJ = TJ->Nb; 2248 const PetscInt NcJ = TJ->Nc; 2249 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2250 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2251 PetscInt f, fc, g, gc, df, dg; 2252 PetscErrorCode ierr; 2253 2254 for (f = 0; f < NbI; ++f) { 2255 for (fc = 0; fc < NcI; ++fc) { 2256 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2257 2258 tmpBasisI[fidx] = basisI[fidx]; 2259 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2260 } 2261 } 2262 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2263 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2264 for (g = 0; g < NbJ; ++g) { 2265 for (gc = 0; gc < NcJ; ++gc) { 2266 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2267 2268 tmpBasisJ[gidx] = basisJ[gidx]; 2269 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2270 } 2271 } 2272 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2273 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2274 for (f = 0; f < NbI; ++f) { 2275 for (fc = 0; fc < NcI; ++fc) { 2276 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2277 const PetscInt i = offsetI+f; /* Element matrix row */ 2278 for (g = 0; g < NbJ; ++g) { 2279 for (gc = 0; gc < NcJ; ++gc) { 2280 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2281 const PetscInt j = offsetJ+g; /* Element matrix column */ 2282 const PetscInt fOff = eOffset+i*totDim+j; 2283 2284 elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx]; 2285 for (df = 0; df < dE; ++df) { 2286 elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2287 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx]; 2288 for (dg = 0; dg < dE; ++dg) { 2289 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2290 } 2291 } 2292 } 2293 } 2294 } 2295 } 2296 return(0); 2297 } 2298 2299 PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE feI, PetscBool isHybridI, PetscFE feJ, PetscBool isHybridJ, PetscInt r, PetscInt q, PetscTabulation TI, PetscScalar tmpBasisI[], PetscScalar tmpBasisDerI[], PetscTabulation TJ, PetscScalar tmpBasisJ[], PetscScalar tmpBasisDerJ[], PetscFEGeom *fegeom, const PetscScalar g0[], const PetscScalar g1[], const PetscScalar g2[], const PetscScalar g3[], PetscInt eOffset, PetscInt totDim, PetscInt offsetI, PetscInt offsetJ, PetscScalar elemMat[]) 2300 { 2301 const PetscInt dE = TI->cdim; 2302 const PetscInt NqI = TI->Np; 2303 const PetscInt NbI = TI->Nb; 2304 const PetscInt NcI = TI->Nc; 2305 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2306 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2307 const PetscInt NqJ = TJ->Np; 2308 const PetscInt NbJ = TJ->Nb; 2309 const PetscInt NcJ = TJ->Nc; 2310 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2311 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2312 const PetscInt Ns = isHybridI ? 1 : 2; 2313 const PetscInt Nt = isHybridJ ? 1 : 2; 2314 PetscInt f, fc, g, gc, df, dg, s, t; 2315 PetscErrorCode ierr; 2316 2317 for (f = 0; f < NbI; ++f) { 2318 for (fc = 0; fc < NcI; ++fc) { 2319 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2320 2321 tmpBasisI[fidx] = basisI[fidx]; 2322 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2323 } 2324 } 2325 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2326 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2327 for (g = 0; g < NbJ; ++g) { 2328 for (gc = 0; gc < NcJ; ++gc) { 2329 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2330 2331 tmpBasisJ[gidx] = basisJ[gidx]; 2332 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2333 } 2334 } 2335 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2336 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2337 for (s = 0; s < Ns; ++s) { 2338 for (f = 0; f < NbI; ++f) { 2339 for (fc = 0; fc < NcI; ++fc) { 2340 const PetscInt sc = NcI*s+fc; /* components from each side of the surface */ 2341 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2342 const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */ 2343 for (t = 0; t < Nt; ++t) { 2344 for (g = 0; g < NbJ; ++g) { 2345 for (gc = 0; gc < NcJ; ++gc) { 2346 const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */ 2347 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2348 const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */ 2349 const PetscInt fOff = eOffset+i*totDim+j; 2350 2351 elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx]; 2352 for (df = 0; df < dE; ++df) { 2353 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2354 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx]; 2355 for (dg = 0; dg < dE; ++dg) { 2356 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2357 } 2358 } 2359 } 2360 } 2361 } 2362 } 2363 } 2364 } 2365 return(0); 2366 } 2367 2368 PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom) 2369 { 2370 PetscDualSpace dsp; 2371 DM dm; 2372 PetscQuadrature quadDef; 2373 PetscInt dim, cdim, Nq; 2374 PetscErrorCode ierr; 2375 2376 PetscFunctionBegin; 2377 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); 2378 ierr = PetscDualSpaceGetDM(dsp, &dm);CHKERRQ(ierr); 2379 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2380 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 2381 ierr = PetscFEGetQuadrature(fe, &quadDef);CHKERRQ(ierr); 2382 quad = quad ? quad : quadDef; 2383 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2384 ierr = PetscMalloc1(Nq*cdim, &cgeom->v);CHKERRQ(ierr); 2385 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->J);CHKERRQ(ierr); 2386 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);CHKERRQ(ierr); 2387 ierr = PetscMalloc1(Nq, &cgeom->detJ);CHKERRQ(ierr); 2388 cgeom->dim = dim; 2389 cgeom->dimEmbed = cdim; 2390 cgeom->numCells = 1; 2391 cgeom->numPoints = Nq; 2392 ierr = DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);CHKERRQ(ierr); 2393 PetscFunctionReturn(0); 2394 } 2395 2396 PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom) 2397 { 2398 PetscErrorCode ierr; 2399 2400 PetscFunctionBegin; 2401 ierr = PetscFree(cgeom->v);CHKERRQ(ierr); 2402 ierr = PetscFree(cgeom->J);CHKERRQ(ierr); 2403 ierr = PetscFree(cgeom->invJ);CHKERRQ(ierr); 2404 ierr = PetscFree(cgeom->detJ);CHKERRQ(ierr); 2405 PetscFunctionReturn(0); 2406 } 2407