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 PetscFEGetDimension - Get the dimension of the finite element space on a cell 1161 1162 Not collective 1163 1164 Input Parameter: 1165 . fe - The PetscFE 1166 1167 Output Parameter: 1168 . dim - The dimension 1169 1170 Level: intermediate 1171 1172 .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension() 1173 @*/ 1174 PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) 1175 { 1176 PetscErrorCode ierr; 1177 1178 PetscFunctionBegin; 1179 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1180 PetscValidPointer(dim, 2); 1181 if (fem->ops->getdimension) {ierr = (*fem->ops->getdimension)(fem, dim);CHKERRQ(ierr);} 1182 PetscFunctionReturn(0); 1183 } 1184 1185 /*@C 1186 PetscFEPushforward - Map the reference element function to real space 1187 1188 Input Parameters: 1189 + fe - The PetscFE 1190 . fegeom - The cell geometry 1191 . Nv - The number of function values 1192 - vals - The function values 1193 1194 Output Parameter: 1195 . vals - The transformed function values 1196 1197 Level: advanced 1198 1199 Note: This just forwards the call onto PetscDualSpacePushforward(). 1200 1201 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1202 1203 .seealso: PetscDualSpacePushforward() 1204 @*/ 1205 PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1206 { 1207 PetscErrorCode ierr; 1208 1209 PetscFunctionBeginHot; 1210 ierr = PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1211 PetscFunctionReturn(0); 1212 } 1213 1214 /*@C 1215 PetscFEPushforwardGradient - Map the reference element function gradient to real space 1216 1217 Input Parameters: 1218 + fe - The PetscFE 1219 . fegeom - The cell geometry 1220 . Nv - The number of function gradient values 1221 - vals - The function gradient values 1222 1223 Output Parameter: 1224 . vals - The transformed function gradient values 1225 1226 Level: advanced 1227 1228 Note: This just forwards the call onto PetscDualSpacePushforwardGradient(). 1229 1230 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1231 1232 .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward() 1233 @*/ 1234 PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1235 { 1236 PetscErrorCode ierr; 1237 1238 PetscFunctionBeginHot; 1239 ierr = PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1240 PetscFunctionReturn(0); 1241 } 1242 1243 /*@C 1244 PetscFEPushforwardHessian - Map the reference element function Hessian to real space 1245 1246 Input Parameters: 1247 + fe - The PetscFE 1248 . fegeom - The cell geometry 1249 . Nv - The number of function Hessian values 1250 - vals - The function Hessian values 1251 1252 Output Parameter: 1253 . vals - The transformed function Hessian values 1254 1255 Level: advanced 1256 1257 Note: This just forwards the call onto PetscDualSpacePushforwardHessian(). 1258 1259 Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. 1260 1261 .seealso: PetscFEPushforward(), PetscDualSpacePushforwardHessian(), PetscDualSpacePushforward() 1262 @*/ 1263 PetscErrorCode PetscFEPushforwardHessian(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) 1264 { 1265 PetscErrorCode ierr; 1266 1267 PetscFunctionBeginHot; 1268 ierr = PetscDualSpacePushforwardHessian(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); 1269 PetscFunctionReturn(0); 1270 } 1271 1272 /* 1273 Purpose: Compute element vector for chunk of elements 1274 1275 Input: 1276 Sizes: 1277 Ne: number of elements 1278 Nf: number of fields 1279 PetscFE 1280 dim: spatial dimension 1281 Nb: number of basis functions 1282 Nc: number of field components 1283 PetscQuadrature 1284 Nq: number of quadrature points 1285 1286 Geometry: 1287 PetscFEGeom[Ne] possibly *Nq 1288 PetscReal v0s[dim] 1289 PetscReal n[dim] 1290 PetscReal jacobians[dim*dim] 1291 PetscReal jacobianInverses[dim*dim] 1292 PetscReal jacobianDeterminants 1293 FEM: 1294 PetscFE 1295 PetscQuadrature 1296 PetscReal quadPoints[Nq*dim] 1297 PetscReal quadWeights[Nq] 1298 PetscReal basis[Nq*Nb*Nc] 1299 PetscReal basisDer[Nq*Nb*Nc*dim] 1300 PetscScalar coefficients[Ne*Nb*Nc] 1301 PetscScalar elemVec[Ne*Nb*Nc] 1302 1303 Problem: 1304 PetscInt f: the active field 1305 f0, f1 1306 1307 Work Space: 1308 PetscFE 1309 PetscScalar f0[Nq*dim]; 1310 PetscScalar f1[Nq*dim*dim]; 1311 PetscScalar u[Nc]; 1312 PetscScalar gradU[Nc*dim]; 1313 PetscReal x[dim]; 1314 PetscScalar realSpaceDer[dim]; 1315 1316 Purpose: Compute element vector for N_cb batches of elements 1317 1318 Input: 1319 Sizes: 1320 N_cb: Number of serial cell batches 1321 1322 Geometry: 1323 PetscReal v0s[Ne*dim] 1324 PetscReal jacobians[Ne*dim*dim] possibly *Nq 1325 PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq 1326 PetscReal jacobianDeterminants[Ne] possibly *Nq 1327 FEM: 1328 static PetscReal quadPoints[Nq*dim] 1329 static PetscReal quadWeights[Nq] 1330 static PetscReal basis[Nq*Nb*Nc] 1331 static PetscReal basisDer[Nq*Nb*Nc*dim] 1332 PetscScalar coefficients[Ne*Nb*Nc] 1333 PetscScalar elemVec[Ne*Nb*Nc] 1334 1335 ex62.c: 1336 PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], 1337 const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], 1338 void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), 1339 void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) 1340 1341 ex52.c: 1342 PetscErrorCode IntegrateLaplacianBatchCPU(PetscInt Ne, PetscInt Nb, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user) 1343 PetscErrorCode IntegrateElasticityBatchCPU(PetscInt Ne, PetscInt Nb, PetscInt Ncomp, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user) 1344 1345 ex52_integrateElement.cu 1346 __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) 1347 1348 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], 1349 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1350 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1351 1352 ex52_integrateElementOpenCL.c: 1353 PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], 1354 const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], 1355 PetscLogEvent event, PetscInt debug, PetscInt pde_op) 1356 1357 __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) 1358 */ 1359 1360 /*@C 1361 PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration 1362 1363 Not collective 1364 1365 Input Parameters: 1366 + prob - The PetscDS specifying the discretizations and continuum functions 1367 . field - The field being integrated 1368 . Ne - The number of elements in the chunk 1369 . cgeom - The cell geometry for each cell in the chunk 1370 . coefficients - The array of FEM basis coefficients for the elements 1371 . probAux - The PetscDS specifying the auxiliary discretizations 1372 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1373 1374 Output Parameter: 1375 . integral - the integral for this field 1376 1377 Level: intermediate 1378 1379 .seealso: PetscFEIntegrateResidual() 1380 @*/ 1381 PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 1382 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1383 { 1384 PetscFE fe; 1385 PetscErrorCode ierr; 1386 1387 PetscFunctionBegin; 1388 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1389 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1390 if (fe->ops->integrate) {ierr = (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} 1391 PetscFunctionReturn(0); 1392 } 1393 1394 /*@C 1395 PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration 1396 1397 Not collective 1398 1399 Input Parameters: 1400 + prob - The PetscDS specifying the discretizations and continuum functions 1401 . field - The field being integrated 1402 . obj_func - The function to be integrated 1403 . Ne - The number of elements in the chunk 1404 . fgeom - The face geometry for each face in the chunk 1405 . coefficients - The array of FEM basis coefficients for the elements 1406 . probAux - The PetscDS specifying the auxiliary discretizations 1407 - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1408 1409 Output Parameter: 1410 . integral - the integral for this field 1411 1412 Level: intermediate 1413 1414 .seealso: PetscFEIntegrateResidual() 1415 @*/ 1416 PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field, 1417 void (*obj_func)(PetscInt, PetscInt, PetscInt, 1418 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1419 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1420 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 1421 PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 1422 { 1423 PetscFE fe; 1424 PetscErrorCode ierr; 1425 1426 PetscFunctionBegin; 1427 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1428 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1429 if (fe->ops->integratebd) {ierr = (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} 1430 PetscFunctionReturn(0); 1431 } 1432 1433 /*@C 1434 PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration 1435 1436 Not collective 1437 1438 Input Parameters: 1439 + ds - The PetscDS specifying the discretizations and continuum functions 1440 . key - The (label+value, field) being integrated 1441 . Ne - The number of elements in the chunk 1442 . cgeom - The cell geometry for each cell in the chunk 1443 . coefficients - The array of FEM basis coefficients for the elements 1444 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1445 . probAux - The PetscDS specifying the auxiliary discretizations 1446 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1447 - t - The time 1448 1449 Output Parameter: 1450 . elemVec - the element residual vectors from each element 1451 1452 Note: 1453 $ Loop over batch of elements (e): 1454 $ Loop over quadrature points (q): 1455 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q 1456 $ Call f_0 and f_1 1457 $ Loop over element vector entries (f,fc --> i): 1458 $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) 1459 1460 Level: intermediate 1461 1462 .seealso: PetscFEIntegrateResidual() 1463 @*/ 1464 PetscErrorCode PetscFEIntegrateResidual(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, 1465 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1466 { 1467 PetscFE fe; 1468 PetscErrorCode ierr; 1469 1470 PetscFunctionBeginHot; 1471 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1472 ierr = PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);CHKERRQ(ierr); 1473 if (fe->ops->integrateresidual) {ierr = (*fe->ops->integrateresidual)(ds, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1474 PetscFunctionReturn(0); 1475 } 1476 1477 /*@C 1478 PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary 1479 1480 Not collective 1481 1482 Input Parameters: 1483 + ds - The PetscDS specifying the discretizations and continuum functions 1484 . wf - The PetscWeakForm object holding the pointwise functions 1485 . key - The (label+value, field) being integrated 1486 . Ne - The number of elements in the chunk 1487 . fgeom - The face geometry for each cell in the chunk 1488 . coefficients - The array of FEM basis coefficients for the elements 1489 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1490 . probAux - The PetscDS specifying the auxiliary discretizations 1491 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1492 - t - The time 1493 1494 Output Parameter: 1495 . elemVec - the element residual vectors from each element 1496 1497 Level: intermediate 1498 1499 .seealso: PetscFEIntegrateResidual() 1500 @*/ 1501 PetscErrorCode PetscFEIntegrateBdResidual(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1502 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1503 { 1504 PetscFE fe; 1505 PetscErrorCode ierr; 1506 1507 PetscFunctionBegin; 1508 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1509 ierr = PetscDSGetDiscretization(ds, key.field, (PetscObject *) &fe);CHKERRQ(ierr); 1510 if (fe->ops->integratebdresidual) {ierr = (*fe->ops->integratebdresidual)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1511 PetscFunctionReturn(0); 1512 } 1513 1514 /*@C 1515 PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration 1516 1517 Not collective 1518 1519 Input Parameters: 1520 + prob - The PetscDS specifying the discretizations and continuum functions 1521 . key - The (label+value, field) being integrated 1522 . Ne - The number of elements in the chunk 1523 . fgeom - The face geometry for each cell in the chunk 1524 . coefficients - The array of FEM basis coefficients for the elements 1525 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1526 . probAux - The PetscDS specifying the auxiliary discretizations 1527 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1528 - t - The time 1529 1530 Output Parameter 1531 . elemVec - the element residual vectors from each element 1532 1533 Level: developer 1534 1535 .seealso: PetscFEIntegrateResidual() 1536 @*/ 1537 PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1538 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 1539 { 1540 PetscFE fe; 1541 PetscErrorCode ierr; 1542 1543 PetscFunctionBegin; 1544 PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); 1545 ierr = PetscDSGetDiscretization(prob, key.field, (PetscObject *) &fe);CHKERRQ(ierr); 1546 if (fe->ops->integratehybridresidual) {ierr = (*fe->ops->integratehybridresidual)(prob, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} 1547 PetscFunctionReturn(0); 1548 } 1549 1550 /*@C 1551 PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration 1552 1553 Not collective 1554 1555 Input Parameters: 1556 + ds - The PetscDS specifying the discretizations and continuum functions 1557 . jtype - The type of matrix pointwise functions that should be used 1558 . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1559 . Ne - The number of elements in the chunk 1560 . cgeom - The cell geometry for each cell in the chunk 1561 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1562 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1563 . probAux - The PetscDS specifying the auxiliary discretizations 1564 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1565 . t - The time 1566 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1567 1568 Output Parameter: 1569 . elemMat - the element matrices for the Jacobian from each element 1570 1571 Note: 1572 $ Loop over batch of elements (e): 1573 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1574 $ Loop over quadrature points (q): 1575 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1576 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1577 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1578 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1579 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1580 Level: intermediate 1581 1582 .seealso: PetscFEIntegrateResidual() 1583 @*/ 1584 PetscErrorCode PetscFEIntegrateJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, 1585 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1586 { 1587 PetscFE fe; 1588 PetscInt Nf; 1589 PetscErrorCode ierr; 1590 1591 PetscFunctionBegin; 1592 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1593 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1594 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1595 if (fe->ops->integratejacobian) {ierr = (*fe->ops->integratejacobian)(ds, jtype, key, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1596 PetscFunctionReturn(0); 1597 } 1598 1599 /*@C 1600 PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration 1601 1602 Not collective 1603 1604 Input Parameters: 1605 + ds - The PetscDS specifying the discretizations and continuum functions 1606 . wf - The PetscWeakForm holding the pointwise functions 1607 . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1608 . Ne - The number of elements in the chunk 1609 . fgeom - The face geometry for each cell in the chunk 1610 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1611 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1612 . probAux - The PetscDS specifying the auxiliary discretizations 1613 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1614 . t - The time 1615 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1616 1617 Output Parameter: 1618 . elemMat - the element matrices for the Jacobian from each element 1619 1620 Note: 1621 $ Loop over batch of elements (e): 1622 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1623 $ Loop over quadrature points (q): 1624 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1625 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1626 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1627 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1628 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1629 Level: intermediate 1630 1631 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1632 @*/ 1633 PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1634 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1635 { 1636 PetscFE fe; 1637 PetscInt Nf; 1638 PetscErrorCode ierr; 1639 1640 PetscFunctionBegin; 1641 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1642 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1643 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1644 if (fe->ops->integratebdjacobian) {ierr = (*fe->ops->integratebdjacobian)(ds, wf, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1645 PetscFunctionReturn(0); 1646 } 1647 1648 /*@C 1649 PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration 1650 1651 Not collective 1652 1653 Input Parameters: 1654 + ds - The PetscDS specifying the discretizations and continuum functions 1655 . jtype - The type of matrix pointwise functions that should be used 1656 . key - The (label+value, fieldI*Nf + fieldJ) being integrated 1657 . Ne - The number of elements in the chunk 1658 . fgeom - The face geometry for each cell in the chunk 1659 . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point 1660 . coefficients_t - The array of FEM basis time derivative coefficients for the elements 1661 . probAux - The PetscDS specifying the auxiliary discretizations 1662 . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements 1663 . t - The time 1664 - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) 1665 1666 Output Parameter 1667 . elemMat - the element matrices for the Jacobian from each element 1668 1669 Note: 1670 $ Loop over batch of elements (e): 1671 $ Loop over element matrix entries (f,fc,g,gc --> i,j): 1672 $ Loop over quadrature points (q): 1673 $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) 1674 $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) 1675 $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1676 $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) 1677 $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) 1678 Level: developer 1679 1680 .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() 1681 @*/ 1682 PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, 1683 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 1684 { 1685 PetscFE fe; 1686 PetscInt Nf; 1687 PetscErrorCode ierr; 1688 1689 PetscFunctionBegin; 1690 PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1); 1691 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 1692 ierr = PetscDSGetDiscretization(ds, key.field / Nf, (PetscObject *) &fe);CHKERRQ(ierr); 1693 if (fe->ops->integratehybridjacobian) {ierr = (*fe->ops->integratehybridjacobian)(ds, jtype, key, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} 1694 PetscFunctionReturn(0); 1695 } 1696 1697 /*@ 1698 PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height 1699 1700 Input Parameters: 1701 + fe - The finite element space 1702 - height - The height of the Plex point 1703 1704 Output Parameter: 1705 . subfe - The subspace of this FE space 1706 1707 Note: For example, if we want the subspace of this space for a face, we would choose height = 1. 1708 1709 Level: advanced 1710 1711 .seealso: PetscFECreateDefault() 1712 @*/ 1713 PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe) 1714 { 1715 PetscSpace P, subP; 1716 PetscDualSpace Q, subQ; 1717 PetscQuadrature subq; 1718 PetscFEType fetype; 1719 PetscInt dim, Nc; 1720 PetscErrorCode ierr; 1721 1722 PetscFunctionBegin; 1723 PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1); 1724 PetscValidPointer(subfe, 3); 1725 if (height == 0) { 1726 *subfe = fe; 1727 PetscFunctionReturn(0); 1728 } 1729 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1730 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1731 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1732 ierr = PetscFEGetFaceQuadrature(fe, &subq);CHKERRQ(ierr); 1733 ierr = PetscDualSpaceGetDimension(Q, &dim);CHKERRQ(ierr); 1734 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); 1735 if (!fe->subspaces) {ierr = PetscCalloc1(dim, &fe->subspaces);CHKERRQ(ierr);} 1736 if (height <= dim) { 1737 if (!fe->subspaces[height-1]) { 1738 PetscFE sub = NULL; 1739 const char *name; 1740 1741 ierr = PetscSpaceGetHeightSubspace(P, height, &subP);CHKERRQ(ierr); 1742 ierr = PetscDualSpaceGetHeightSubspace(Q, height, &subQ);CHKERRQ(ierr); 1743 if (subQ) { 1744 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);CHKERRQ(ierr); 1745 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 1746 ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); 1747 ierr = PetscFEGetType(fe, &fetype);CHKERRQ(ierr); 1748 ierr = PetscFESetType(sub, fetype);CHKERRQ(ierr); 1749 ierr = PetscFESetBasisSpace(sub, subP);CHKERRQ(ierr); 1750 ierr = PetscFESetDualSpace(sub, subQ);CHKERRQ(ierr); 1751 ierr = PetscFESetNumComponents(sub, Nc);CHKERRQ(ierr); 1752 ierr = PetscFESetUp(sub);CHKERRQ(ierr); 1753 ierr = PetscFESetQuadrature(sub, subq);CHKERRQ(ierr); 1754 } 1755 fe->subspaces[height-1] = sub; 1756 } 1757 *subfe = fe->subspaces[height-1]; 1758 } else { 1759 *subfe = NULL; 1760 } 1761 PetscFunctionReturn(0); 1762 } 1763 1764 /*@ 1765 PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used 1766 to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more 1767 sparsity). It is also used to create an interpolation between regularly refined meshes. 1768 1769 Collective on fem 1770 1771 Input Parameter: 1772 . fe - The initial PetscFE 1773 1774 Output Parameter: 1775 . feRef - The refined PetscFE 1776 1777 Level: advanced 1778 1779 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 1780 @*/ 1781 PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef) 1782 { 1783 PetscSpace P, Pref; 1784 PetscDualSpace Q, Qref; 1785 DM K, Kref; 1786 PetscQuadrature q, qref; 1787 const PetscReal *v0, *jac; 1788 PetscInt numComp, numSubelements; 1789 PetscInt cStart, cEnd, c; 1790 PetscDualSpace *cellSpaces; 1791 PetscErrorCode ierr; 1792 1793 PetscFunctionBegin; 1794 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 1795 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1796 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1797 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr); 1798 /* Create space */ 1799 ierr = PetscObjectReference((PetscObject) P);CHKERRQ(ierr); 1800 Pref = P; 1801 /* Create dual space */ 1802 ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr); 1803 ierr = PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);CHKERRQ(ierr); 1804 ierr = DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);CHKERRQ(ierr); 1805 ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr); 1806 ierr = DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);CHKERRQ(ierr); 1807 ierr = PetscMalloc1(cEnd - cStart, &cellSpaces);CHKERRQ(ierr); 1808 /* TODO: fix for non-uniform refinement */ 1809 for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q; 1810 ierr = PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);CHKERRQ(ierr); 1811 ierr = PetscFree(cellSpaces);CHKERRQ(ierr); 1812 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 1813 ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr); 1814 /* Create element */ 1815 ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);CHKERRQ(ierr); 1816 ierr = PetscFESetType(*feRef, PETSCFECOMPOSITE);CHKERRQ(ierr); 1817 ierr = PetscFESetBasisSpace(*feRef, Pref);CHKERRQ(ierr); 1818 ierr = PetscFESetDualSpace(*feRef, Qref);CHKERRQ(ierr); 1819 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 1820 ierr = PetscFESetNumComponents(*feRef, numComp);CHKERRQ(ierr); 1821 ierr = PetscFESetUp(*feRef);CHKERRQ(ierr); 1822 ierr = PetscSpaceDestroy(&Pref);CHKERRQ(ierr); 1823 ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr); 1824 /* Create quadrature */ 1825 ierr = PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);CHKERRQ(ierr); 1826 ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr); 1827 ierr = PetscFESetQuadrature(*feRef, qref);CHKERRQ(ierr); 1828 ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr); 1829 PetscFunctionReturn(0); 1830 } 1831 1832 /*@C 1833 PetscFECreateDefault - Create a PetscFE for basic FEM computation 1834 1835 Collective 1836 1837 Input Parameters: 1838 + comm - The MPI comm 1839 . dim - The spatial dimension 1840 . Nc - The number of components 1841 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1842 . prefix - The options prefix, or NULL 1843 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1844 1845 Output Parameter: 1846 . fem - The PetscFE object 1847 1848 Note: 1849 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. 1850 1851 Level: beginner 1852 1853 .seealso: PetscSpaceSetFromOptions(), PetscDualSpaceSetFromOptions(), PetscFESetFromOptions(), PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1854 @*/ 1855 PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem) 1856 { 1857 PetscQuadrature q, fq; 1858 DM K; 1859 PetscSpace P; 1860 PetscDualSpace Q; 1861 PetscInt order, quadPointsPerEdge; 1862 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1863 PetscErrorCode ierr; 1864 1865 PetscFunctionBegin; 1866 /* Create space */ 1867 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1868 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr); 1869 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1870 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1871 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1872 ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr); 1873 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1874 ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr); 1875 ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr); 1876 /* Create dual space */ 1877 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1878 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1879 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr); 1880 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1881 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1882 ierr = DMDestroy(&K);CHKERRQ(ierr); 1883 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1884 ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr); 1885 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1886 ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr); 1887 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1888 /* Create element */ 1889 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1890 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr); 1891 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1892 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1893 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1894 ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr); 1895 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1896 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1897 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1898 /* Create quadrature (with specified order if given) */ 1899 qorder = qorder >= 0 ? qorder : order; 1900 ierr = PetscObjectOptionsBegin((PetscObject)*fem);CHKERRQ(ierr); 1901 ierr = PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);CHKERRQ(ierr); 1902 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1903 quadPointsPerEdge = PetscMax(qorder + 1,1); 1904 if (isSimplex) { 1905 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1906 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1907 } else { 1908 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1909 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1910 } 1911 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1912 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1913 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1914 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1915 PetscFunctionReturn(0); 1916 } 1917 1918 /*@ 1919 PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k 1920 1921 Collective 1922 1923 Input Parameters: 1924 + comm - The MPI comm 1925 . dim - The spatial dimension 1926 . Nc - The number of components 1927 . isSimplex - Flag for simplex reference cell, otherwise its a tensor product 1928 . k - The degree k of the space 1929 - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree 1930 1931 Output Parameter: 1932 . fem - The PetscFE object 1933 1934 Level: beginner 1935 1936 Notes: 1937 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. 1938 1939 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 1940 @*/ 1941 PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem) 1942 { 1943 PetscQuadrature q, fq; 1944 DM K; 1945 PetscSpace P; 1946 PetscDualSpace Q; 1947 PetscInt quadPointsPerEdge; 1948 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1949 char name[64]; 1950 PetscErrorCode ierr; 1951 1952 PetscFunctionBegin; 1953 /* Create space */ 1954 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1955 ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1956 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 1957 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1958 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1959 ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr); 1960 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1961 /* Create dual space */ 1962 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1963 ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1964 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1965 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1966 ierr = DMDestroy(&K);CHKERRQ(ierr); 1967 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1968 ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr); 1969 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 1970 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1971 /* Create finite element */ 1972 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1973 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1974 ierr = PetscObjectSetName((PetscObject) *fem, name);CHKERRQ(ierr); 1975 ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr); 1976 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1977 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1978 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1979 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1980 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1981 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1982 /* Create quadrature (with specified order if given) */ 1983 qorder = qorder >= 0 ? qorder : k; 1984 quadPointsPerEdge = PetscMax(qorder + 1,1); 1985 if (isSimplex) { 1986 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1987 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1988 } else { 1989 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 1990 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 1991 } 1992 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1993 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1994 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1995 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1996 /* Set finite element name */ 1997 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1998 ierr = PetscFESetName(*fem, name);CHKERRQ(ierr); 1999 PetscFunctionReturn(0); 2000 } 2001 2002 /*@C 2003 PetscFESetName - Names the FE and its subobjects 2004 2005 Not collective 2006 2007 Input Parameters: 2008 + fe - The PetscFE 2009 - name - The name 2010 2011 Level: intermediate 2012 2013 .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() 2014 @*/ 2015 PetscErrorCode PetscFESetName(PetscFE fe, const char name[]) 2016 { 2017 PetscSpace P; 2018 PetscDualSpace Q; 2019 PetscErrorCode ierr; 2020 2021 PetscFunctionBegin; 2022 ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); 2023 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2024 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 2025 ierr = PetscObjectSetName((PetscObject) P, name);CHKERRQ(ierr); 2026 ierr = PetscObjectSetName((PetscObject) Q, name);CHKERRQ(ierr); 2027 PetscFunctionReturn(0); 2028 } 2029 2030 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[]) 2031 { 2032 PetscInt dOffset = 0, fOffset = 0, f, g; 2033 PetscErrorCode ierr; 2034 2035 for (f = 0; f < Nf; ++f) { 2036 PetscFE fe; 2037 const PetscInt k = ds->jetDegree[f]; 2038 const PetscInt cdim = T[f]->cdim; 2039 const PetscInt Nq = T[f]->Np; 2040 const PetscInt Nbf = T[f]->Nb; 2041 const PetscInt Ncf = T[f]->Nc; 2042 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2043 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2044 const PetscReal *Hq = k > 1 ? &T[f]->T[2][(r*Nq+q)*Nbf*Ncf*cdim*cdim] : NULL; 2045 PetscInt hOffset = 0, b, c, d; 2046 2047 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 2048 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2049 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2050 for (b = 0; b < Nbf; ++b) { 2051 for (c = 0; c < Ncf; ++c) { 2052 const PetscInt cidx = b*Ncf+c; 2053 2054 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2055 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2056 } 2057 } 2058 if (k > 1) { 2059 for (g = 0; g < Nf; ++g) hOffset += T[g]->Nc*cdim; 2060 for (d = 0; d < cdim*cdim*Ncf; ++d) u_x[hOffset+fOffset*cdim*cdim+d] = 0.0; 2061 for (b = 0; b < Nbf; ++b) { 2062 for (c = 0; c < Ncf; ++c) { 2063 const PetscInt cidx = b*Ncf+c; 2064 2065 for (d = 0; d < cdim*cdim; ++d) u_x[hOffset+(fOffset+c)*cdim*cdim+d] += Hq[cidx*cdim*cdim+d]*coefficients[dOffset+b]; 2066 } 2067 } 2068 ierr = PetscFEPushforwardHessian(fe, fegeom, 1, &u_x[hOffset+fOffset*cdim*cdim]);CHKERRQ(ierr); 2069 } 2070 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2071 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2072 if (u_t) { 2073 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2074 for (b = 0; b < Nbf; ++b) { 2075 for (c = 0; c < Ncf; ++c) { 2076 const PetscInt cidx = b*Ncf+c; 2077 2078 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2079 } 2080 } 2081 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2082 } 2083 fOffset += Ncf; 2084 dOffset += Nbf; 2085 } 2086 return 0; 2087 } 2088 2089 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[]) 2090 { 2091 PetscInt dOffset = 0, fOffset = 0, g; 2092 PetscErrorCode ierr; 2093 2094 for (g = 0; g < 2*Nf-1; ++g) { 2095 if (!T[g/2]) continue; 2096 { 2097 PetscFE fe; 2098 const PetscInt f = g/2; 2099 const PetscInt cdim = T[f]->cdim; 2100 const PetscInt Nq = T[f]->Np; 2101 const PetscInt Nbf = T[f]->Nb; 2102 const PetscInt Ncf = T[f]->Nc; 2103 const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; 2104 const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; 2105 PetscInt b, c, d; 2106 2107 fe = (PetscFE) ds->disc[f]; 2108 for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; 2109 for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; 2110 for (b = 0; b < Nbf; ++b) { 2111 for (c = 0; c < Ncf; ++c) { 2112 const PetscInt cidx = b*Ncf+c; 2113 2114 u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; 2115 for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; 2116 } 2117 } 2118 ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); 2119 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); 2120 if (u_t) { 2121 for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; 2122 for (b = 0; b < Nbf; ++b) { 2123 for (c = 0; c < Ncf; ++c) { 2124 const PetscInt cidx = b*Ncf+c; 2125 2126 u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; 2127 } 2128 } 2129 ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); 2130 } 2131 fOffset += Ncf; 2132 dOffset += Nbf; 2133 } 2134 } 2135 return 0; 2136 } 2137 2138 PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[]) 2139 { 2140 PetscFE fe; 2141 PetscTabulation Tc; 2142 PetscInt b, c; 2143 PetscErrorCode ierr; 2144 2145 if (!prob) return 0; 2146 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 2147 ierr = PetscFEGetFaceCentroidTabulation(fe, &Tc);CHKERRQ(ierr); 2148 { 2149 const PetscReal *faceBasis = Tc->T[0]; 2150 const PetscInt Nb = Tc->Nb; 2151 const PetscInt Nc = Tc->Nc; 2152 2153 for (c = 0; c < Nc; ++c) {u[c] = 0.0;} 2154 for (b = 0; b < Nb; ++b) { 2155 for (c = 0; c < Nc; ++c) { 2156 u[c] += coefficients[b] * faceBasis[(faceLoc*Nb + b)*Nc + c]; 2157 } 2158 } 2159 } 2160 return 0; 2161 } 2162 2163 PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscInt e, PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2164 { 2165 PetscFEGeom pgeom; 2166 const PetscInt dEt = T->cdim; 2167 const PetscInt dE = fegeom->dimEmbed; 2168 const PetscInt Nq = T->Np; 2169 const PetscInt Nb = T->Nb; 2170 const PetscInt Nc = T->Nc; 2171 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2172 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dEt]; 2173 PetscInt q, b, c, d; 2174 PetscErrorCode ierr; 2175 2176 for (b = 0; b < Nb; ++b) elemVec[b] = 0.0; 2177 for (q = 0; q < Nq; ++q) { 2178 for (b = 0; b < Nb; ++b) { 2179 for (c = 0; c < Nc; ++c) { 2180 const PetscInt bcidx = b*Nc+c; 2181 2182 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2183 for (d = 0; d < dEt; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dEt+bcidx*dEt+d]; 2184 } 2185 } 2186 ierr = PetscFEGeomGetCellPoint(fegeom, e, q, &pgeom);CHKERRQ(ierr); 2187 ierr = PetscFEPushforward(fe, &pgeom, Nb, tmpBasis);CHKERRQ(ierr); 2188 ierr = PetscFEPushforwardGradient(fe, &pgeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2189 for (b = 0; b < Nb; ++b) { 2190 for (c = 0; c < Nc; ++c) { 2191 const PetscInt bcidx = b*Nc+c; 2192 const PetscInt qcidx = q*Nc+c; 2193 2194 elemVec[b] += tmpBasis[bcidx]*f0[qcidx]; 2195 for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2196 } 2197 } 2198 } 2199 return(0); 2200 } 2201 2202 PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) 2203 { 2204 const PetscInt dE = T->cdim; 2205 const PetscInt Nq = T->Np; 2206 const PetscInt Nb = T->Nb; 2207 const PetscInt Nc = T->Nc; 2208 const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; 2209 const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; 2210 PetscInt q, b, c, d, s; 2211 PetscErrorCode ierr; 2212 2213 for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0; 2214 for (q = 0; q < Nq; ++q) { 2215 for (b = 0; b < Nb; ++b) { 2216 for (c = 0; c < Nc; ++c) { 2217 const PetscInt bcidx = b*Nc+c; 2218 2219 tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; 2220 for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; 2221 } 2222 } 2223 ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); 2224 ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); 2225 for (s = 0; s < 2; ++s) { 2226 for (b = 0; b < Nb; ++b) { 2227 for (c = 0; c < Nc; ++c) { 2228 const PetscInt bcidx = b*Nc+c; 2229 const PetscInt qcidx = (q*2+s)*Nc+c; 2230 2231 elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx]; 2232 for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; 2233 } 2234 } 2235 } 2236 } 2237 return(0); 2238 } 2239 2240 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[]) 2241 { 2242 const PetscInt dE = TI->cdim; 2243 const PetscInt NqI = TI->Np; 2244 const PetscInt NbI = TI->Nb; 2245 const PetscInt NcI = TI->Nc; 2246 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2247 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2248 const PetscInt NqJ = TJ->Np; 2249 const PetscInt NbJ = TJ->Nb; 2250 const PetscInt NcJ = TJ->Nc; 2251 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2252 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2253 PetscInt f, fc, g, gc, df, dg; 2254 PetscErrorCode ierr; 2255 2256 for (f = 0; f < NbI; ++f) { 2257 for (fc = 0; fc < NcI; ++fc) { 2258 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2259 2260 tmpBasisI[fidx] = basisI[fidx]; 2261 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2262 } 2263 } 2264 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2265 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2266 for (g = 0; g < NbJ; ++g) { 2267 for (gc = 0; gc < NcJ; ++gc) { 2268 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2269 2270 tmpBasisJ[gidx] = basisJ[gidx]; 2271 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2272 } 2273 } 2274 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2275 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2276 for (f = 0; f < NbI; ++f) { 2277 for (fc = 0; fc < NcI; ++fc) { 2278 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2279 const PetscInt i = offsetI+f; /* Element matrix row */ 2280 for (g = 0; g < NbJ; ++g) { 2281 for (gc = 0; gc < NcJ; ++gc) { 2282 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2283 const PetscInt j = offsetJ+g; /* Element matrix column */ 2284 const PetscInt fOff = eOffset+i*totDim+j; 2285 2286 elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx]; 2287 for (df = 0; df < dE; ++df) { 2288 elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2289 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx]; 2290 for (dg = 0; dg < dE; ++dg) { 2291 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2292 } 2293 } 2294 } 2295 } 2296 } 2297 } 2298 return(0); 2299 } 2300 2301 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[]) 2302 { 2303 const PetscInt dE = TI->cdim; 2304 const PetscInt NqI = TI->Np; 2305 const PetscInt NbI = TI->Nb; 2306 const PetscInt NcI = TI->Nc; 2307 const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; 2308 const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; 2309 const PetscInt NqJ = TJ->Np; 2310 const PetscInt NbJ = TJ->Nb; 2311 const PetscInt NcJ = TJ->Nc; 2312 const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; 2313 const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; 2314 const PetscInt Ns = isHybridI ? 1 : 2; 2315 const PetscInt Nt = isHybridJ ? 1 : 2; 2316 PetscInt f, fc, g, gc, df, dg, s, t; 2317 PetscErrorCode ierr; 2318 2319 for (f = 0; f < NbI; ++f) { 2320 for (fc = 0; fc < NcI; ++fc) { 2321 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2322 2323 tmpBasisI[fidx] = basisI[fidx]; 2324 for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; 2325 } 2326 } 2327 ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); 2328 ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); 2329 for (g = 0; g < NbJ; ++g) { 2330 for (gc = 0; gc < NcJ; ++gc) { 2331 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2332 2333 tmpBasisJ[gidx] = basisJ[gidx]; 2334 for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; 2335 } 2336 } 2337 ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); 2338 ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); 2339 for (s = 0; s < Ns; ++s) { 2340 for (f = 0; f < NbI; ++f) { 2341 for (fc = 0; fc < NcI; ++fc) { 2342 const PetscInt sc = NcI*s+fc; /* components from each side of the surface */ 2343 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 2344 const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */ 2345 for (t = 0; t < Nt; ++t) { 2346 for (g = 0; g < NbJ; ++g) { 2347 for (gc = 0; gc < NcJ; ++gc) { 2348 const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */ 2349 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 2350 const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */ 2351 const PetscInt fOff = eOffset+i*totDim+j; 2352 2353 elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx]; 2354 for (df = 0; df < dE; ++df) { 2355 elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; 2356 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx]; 2357 for (dg = 0; dg < dE; ++dg) { 2358 elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; 2359 } 2360 } 2361 } 2362 } 2363 } 2364 } 2365 } 2366 } 2367 return(0); 2368 } 2369 2370 PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom) 2371 { 2372 PetscDualSpace dsp; 2373 DM dm; 2374 PetscQuadrature quadDef; 2375 PetscInt dim, cdim, Nq; 2376 PetscErrorCode ierr; 2377 2378 PetscFunctionBegin; 2379 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); 2380 ierr = PetscDualSpaceGetDM(dsp, &dm);CHKERRQ(ierr); 2381 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2382 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 2383 ierr = PetscFEGetQuadrature(fe, &quadDef);CHKERRQ(ierr); 2384 quad = quad ? quad : quadDef; 2385 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2386 ierr = PetscMalloc1(Nq*cdim, &cgeom->v);CHKERRQ(ierr); 2387 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->J);CHKERRQ(ierr); 2388 ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);CHKERRQ(ierr); 2389 ierr = PetscMalloc1(Nq, &cgeom->detJ);CHKERRQ(ierr); 2390 cgeom->dim = dim; 2391 cgeom->dimEmbed = cdim; 2392 cgeom->numCells = 1; 2393 cgeom->numPoints = Nq; 2394 ierr = DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);CHKERRQ(ierr); 2395 PetscFunctionReturn(0); 2396 } 2397 2398 PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom) 2399 { 2400 PetscErrorCode ierr; 2401 2402 PetscFunctionBegin; 2403 ierr = PetscFree(cgeom->v);CHKERRQ(ierr); 2404 ierr = PetscFree(cgeom->J);CHKERRQ(ierr); 2405 ierr = PetscFree(cgeom->invJ);CHKERRQ(ierr); 2406 ierr = PetscFree(cgeom->detJ);CHKERRQ(ierr); 2407 PetscFunctionReturn(0); 2408 } 2409