/* Basis Jet Tabulation We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis as a prime basis. \psi_i = \sum_k \alpha_{ki} \phi_k Our nodal basis is defined in terms of the dual basis $n_j$ n_j \cdot \psi_i = \delta_{ji} and we may act on the first equation to obtain n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k \delta_{ji} = \sum_k \alpha_{ki} V_{jk} I = V \alpha so the coefficients of the nodal basis in the prime basis are \alpha = V^{-1} We will define the dual basis vectors $n_j$ using a quadrature rule. Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can be implemented exactly as in FIAT using functionals $L_j$. I will have to count the degrees correctly for the Legendre product when we are on simplices. We will have three objects: - Space, P: this just need point evaluation I think - 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 - FEM: This keeps {P, P', Q} */ #include /*I "petscfe.h" I*/ #include PetscBool FEcite = PETSC_FALSE; const char FECitation[] = "@article{kirby2004,\n" " title = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n" " journal = {ACM Transactions on Mathematical Software},\n" " author = {Robert C. Kirby},\n" " volume = {30},\n" " number = {4},\n" " pages = {502--516},\n" " doi = {10.1145/1039813.1039820},\n" " year = {2004}\n}\n"; PetscClassId PETSCFE_CLASSID = 0; PetscFunctionList PetscFEList = NULL; PetscBool PetscFERegisterAllCalled = PETSC_FALSE; /*@C PetscFERegister - Adds a new PetscFE implementation Not Collective Input Parameters: + name - The name of a new user-defined creation routine - create_func - The creation routine itself Notes: PetscFERegister() may be called multiple times to add several user-defined PetscFEs Sample usage: .vb PetscFERegister("my_fe", MyPetscFECreate); .ve Then, your PetscFE type can be chosen with the procedural interface via .vb PetscFECreate(MPI_Comm, PetscFE *); PetscFESetType(PetscFE, "my_fe"); .ve or at runtime via the option .vb -petscfe_type my_fe .ve Level: advanced .seealso: PetscFERegisterAll(), PetscFERegisterDestroy() @*/ PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE)) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFunctionListAdd(&PetscFEList, sname, function);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFESetType - Builds a particular PetscFE Collective on fem Input Parameters: + fem - The PetscFE object - name - The kind of FEM space Options Database Key: . -petscfe_type - Sets the PetscFE type; use -help for a list of available types Level: intermediate .seealso: PetscFEGetType(), PetscFECreate() @*/ PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name) { PetscErrorCode (*r)(PetscFE); PetscBool match; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); ierr = PetscObjectTypeCompare((PetscObject) fem, name, &match);CHKERRQ(ierr); if (match) PetscFunctionReturn(0); if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);} ierr = PetscFunctionListFind(PetscFEList, name, &r);CHKERRQ(ierr); if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name); if (fem->ops->destroy) { ierr = (*fem->ops->destroy)(fem);CHKERRQ(ierr); fem->ops->destroy = NULL; } ierr = (*r)(fem);CHKERRQ(ierr); ierr = PetscObjectChangeTypeName((PetscObject) fem, name);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFEGetType - Gets the PetscFE type name (as a string) from the object. Not Collective Input Parameter: . fem - The PetscFE Output Parameter: . name - The PetscFE type name Level: intermediate .seealso: PetscFESetType(), PetscFECreate() @*/ PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(name, 2); if (!PetscFERegisterAllCalled) { ierr = PetscFERegisterAll();CHKERRQ(ierr); } *name = ((PetscObject) fem)->type_name; PetscFunctionReturn(0); } /*@C PetscFEViewFromOptions - View from Options Collective on PetscFE Input Parameters: + A - the PetscFE object . obj - Optional object - name - command line option Level: intermediate .seealso: PetscFE(), PetscFEView(), PetscObjectViewFromOptions(), PetscFECreate() @*/ PetscErrorCode PetscFEViewFromOptions(PetscFE A,PetscObject obj,const char name[]) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(A,PETSCFE_CLASSID,1); ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFEView - Views a PetscFE Collective on fem Input Parameter: + fem - the PetscFE object to view - viewer - the viewer Level: beginner .seealso PetscFEDestroy() @*/ PetscErrorCode PetscFEView(PetscFE fem, PetscViewer viewer) { PetscBool iascii; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &viewer);CHKERRQ(ierr);} ierr = PetscObjectPrintClassNamePrefixType((PetscObject)fem, viewer);CHKERRQ(ierr); ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); if (fem->ops->view) {ierr = (*fem->ops->view)(fem, viewer);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@ PetscFESetFromOptions - sets parameters in a PetscFE from the options database Collective on fem Input Parameter: . fem - the PetscFE object to set options for Options Database: + -petscfe_num_blocks - the number of cell blocks to integrate concurrently - -petscfe_num_batches - the number of cell batches to integrate serially Level: intermediate .seealso PetscFEView() @*/ PetscErrorCode PetscFESetFromOptions(PetscFE fem) { const char *defaultType; char name[256]; PetscBool flg; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); if (!((PetscObject) fem)->type_name) { defaultType = PETSCFEBASIC; } else { defaultType = ((PetscObject) fem)->type_name; } if (!PetscFERegisterAllCalled) {ierr = PetscFERegisterAll();CHKERRQ(ierr);} ierr = PetscObjectOptionsBegin((PetscObject) fem);CHKERRQ(ierr); ierr = PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);CHKERRQ(ierr); if (flg) { ierr = PetscFESetType(fem, name);CHKERRQ(ierr); } else if (!((PetscObject) fem)->type_name) { ierr = PetscFESetType(fem, defaultType);CHKERRQ(ierr); } ierr = PetscOptionsBoundedInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL,1);CHKERRQ(ierr); ierr = PetscOptionsBoundedInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL,1);CHKERRQ(ierr); if (fem->ops->setfromoptions) { ierr = (*fem->ops->setfromoptions)(PetscOptionsObject,fem);CHKERRQ(ierr); } /* process any options handlers added with PetscObjectAddOptionsHandler() */ ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);CHKERRQ(ierr); ierr = PetscOptionsEnd();CHKERRQ(ierr); ierr = PetscFEViewFromOptions(fem, NULL, "-petscfe_view");CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFESetUp - Construct data structures for the PetscFE Collective on fem Input Parameter: . fem - the PetscFE object to setup Level: intermediate .seealso PetscFEView(), PetscFEDestroy() @*/ PetscErrorCode PetscFESetUp(PetscFE fem) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); if (fem->setupcalled) PetscFunctionReturn(0); fem->setupcalled = PETSC_TRUE; if (fem->ops->setup) {ierr = (*fem->ops->setup)(fem);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@ PetscFEDestroy - Destroys a PetscFE object Collective on fem Input Parameter: . fem - the PetscFE object to destroy Level: beginner .seealso PetscFEView() @*/ PetscErrorCode PetscFEDestroy(PetscFE *fem) { PetscErrorCode ierr; PetscFunctionBegin; if (!*fem) PetscFunctionReturn(0); PetscValidHeaderSpecific((*fem), PETSCFE_CLASSID, 1); if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; PetscFunctionReturn(0);} ((PetscObject) (*fem))->refct = 0; if ((*fem)->subspaces) { PetscInt dim, d; ierr = PetscDualSpaceGetDimension((*fem)->dualSpace, &dim);CHKERRQ(ierr); for (d = 0; d < dim; ++d) {ierr = PetscFEDestroy(&(*fem)->subspaces[d]);CHKERRQ(ierr);} } ierr = PetscFree((*fem)->subspaces);CHKERRQ(ierr); ierr = PetscFree((*fem)->invV);CHKERRQ(ierr); ierr = PetscTabulationDestroy(&(*fem)->T);CHKERRQ(ierr); ierr = PetscTabulationDestroy(&(*fem)->Tf);CHKERRQ(ierr); ierr = PetscTabulationDestroy(&(*fem)->Tc);CHKERRQ(ierr); ierr = PetscSpaceDestroy(&(*fem)->basisSpace);CHKERRQ(ierr); ierr = PetscDualSpaceDestroy(&(*fem)->dualSpace);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&(*fem)->quadrature);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&(*fem)->faceQuadrature);CHKERRQ(ierr); if ((*fem)->ops->destroy) {ierr = (*(*fem)->ops->destroy)(*fem);CHKERRQ(ierr);} ierr = PetscHeaderDestroy(fem);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType(). Collective Input Parameter: . comm - The communicator for the PetscFE object Output Parameter: . fem - The PetscFE object Level: beginner .seealso: PetscFESetType(), PETSCFEGALERKIN @*/ PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem) { PetscFE f; PetscErrorCode ierr; PetscFunctionBegin; PetscValidPointer(fem, 2); ierr = PetscCitationsRegister(FECitation,&FEcite);CHKERRQ(ierr); *fem = NULL; ierr = PetscFEInitializePackage();CHKERRQ(ierr); ierr = PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);CHKERRQ(ierr); f->basisSpace = NULL; f->dualSpace = NULL; f->numComponents = 1; f->subspaces = NULL; f->invV = NULL; f->T = NULL; f->Tf = NULL; f->Tc = NULL; ierr = PetscArrayzero(&f->quadrature, 1);CHKERRQ(ierr); ierr = PetscArrayzero(&f->faceQuadrature, 1);CHKERRQ(ierr); f->blockSize = 0; f->numBlocks = 1; f->batchSize = 0; f->numBatches = 1; *fem = f; PetscFunctionReturn(0); } /*@ PetscFEGetSpatialDimension - Returns the spatial dimension of the element Not collective Input Parameter: . fem - The PetscFE object Output Parameter: . dim - The spatial dimension Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim) { DM dm; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(dim, 2); ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); ierr = DMGetDimension(dm, dim);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ PetscFESetNumComponents - Sets the number of components in the element Not collective Input Parameters: + fem - The PetscFE object - comp - The number of field components Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp) { PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); fem->numComponents = comp; PetscFunctionReturn(0); } /*@ PetscFEGetNumComponents - Returns the number of components in the element Not collective Input Parameter: . fem - The PetscFE object Output Parameter: . comp - The number of field components Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp) { PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(comp, 2); *comp = fem->numComponents; PetscFunctionReturn(0); } /*@ PetscFESetTileSizes - Sets the tile sizes for evaluation Not collective Input Parameters: + fem - The PetscFE object . blockSize - The number of elements in a block . numBlocks - The number of blocks in a batch . batchSize - The number of elements in a batch - numBatches - The number of batches in a chunk Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches) { PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); fem->blockSize = blockSize; fem->numBlocks = numBlocks; fem->batchSize = batchSize; fem->numBatches = numBatches; PetscFunctionReturn(0); } /*@ PetscFEGetTileSizes - Returns the tile sizes for evaluation Not collective Input Parameter: . fem - The PetscFE object Output Parameters: + blockSize - The number of elements in a block . numBlocks - The number of blocks in a batch . batchSize - The number of elements in a batch - numBatches - The number of batches in a chunk Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches) { PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); if (blockSize) PetscValidPointer(blockSize, 2); if (numBlocks) PetscValidPointer(numBlocks, 3); if (batchSize) PetscValidPointer(batchSize, 4); if (numBatches) PetscValidPointer(numBatches, 5); if (blockSize) *blockSize = fem->blockSize; if (numBlocks) *numBlocks = fem->numBlocks; if (batchSize) *batchSize = fem->batchSize; if (numBatches) *numBatches = fem->numBatches; PetscFunctionReturn(0); } /*@ PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution Not collective Input Parameter: . fem - The PetscFE object Output Parameter: . sp - The PetscSpace object Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp) { PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(sp, 2); *sp = fem->basisSpace; PetscFunctionReturn(0); } /*@ PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution Not collective Input Parameters: + fem - The PetscFE object - sp - The PetscSpace object Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 2); ierr = PetscSpaceDestroy(&fem->basisSpace);CHKERRQ(ierr); fem->basisSpace = sp; ierr = PetscObjectReference((PetscObject) fem->basisSpace);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product Not collective Input Parameter: . fem - The PetscFE object Output Parameter: . sp - The PetscDualSpace object Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp) { PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(sp, 2); *sp = fem->dualSpace; PetscFunctionReturn(0); } /*@ PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product Not collective Input Parameters: + fem - The PetscFE object - sp - The PetscDualSpace object Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 2); ierr = PetscDualSpaceDestroy(&fem->dualSpace);CHKERRQ(ierr); fem->dualSpace = sp; ierr = PetscObjectReference((PetscObject) fem->dualSpace);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ PetscFEGetQuadrature - Returns the PetscQuadrature used to calculate inner products Not collective Input Parameter: . fem - The PetscFE object Output Parameter: . q - The PetscQuadrature object Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q) { PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(q, 2); *q = fem->quadrature; PetscFunctionReturn(0); } /*@ PetscFESetQuadrature - Sets the PetscQuadrature used to calculate inner products Not collective Input Parameters: + fem - The PetscFE object - q - The PetscQuadrature object Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q) { PetscInt Nc, qNc; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); if (q == fem->quadrature) PetscFunctionReturn(0); ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); ierr = PetscQuadratureGetNumComponents(q, &qNc);CHKERRQ(ierr); 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); ierr = PetscTabulationDestroy(&fem->T);CHKERRQ(ierr); ierr = PetscTabulationDestroy(&fem->Tc);CHKERRQ(ierr); ierr = PetscObjectReference((PetscObject) q);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&fem->quadrature);CHKERRQ(ierr); fem->quadrature = q; PetscFunctionReturn(0); } /*@ PetscFEGetFaceQuadrature - Returns the PetscQuadrature used to calculate inner products on faces Not collective Input Parameter: . fem - The PetscFE object Output Parameter: . q - The PetscQuadrature object Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFEGetFaceQuadrature(PetscFE fem, PetscQuadrature *q) { PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(q, 2); *q = fem->faceQuadrature; PetscFunctionReturn(0); } /*@ PetscFESetFaceQuadrature - Sets the PetscQuadrature used to calculate inner products on faces Not collective Input Parameters: + fem - The PetscFE object - q - The PetscQuadrature object Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFESetFaceQuadrature(PetscFE fem, PetscQuadrature q) { PetscInt Nc, qNc; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); ierr = PetscQuadratureGetNumComponents(q, &qNc);CHKERRQ(ierr); 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); ierr = PetscTabulationDestroy(&fem->Tf);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&fem->faceQuadrature);CHKERRQ(ierr); fem->faceQuadrature = q; ierr = PetscObjectReference((PetscObject) q);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ PetscFECopyQuadrature - Copy both volumetric and surface quadrature Not collective Input Parameters: + sfe - The PetscFE source for the quadratures - tfe - The PetscFE target for the quadratures Level: intermediate .seealso: PetscFECreate(), PetscFESetQuadrature(), PetscFESetFaceQuadrature() @*/ PetscErrorCode PetscFECopyQuadrature(PetscFE sfe, PetscFE tfe) { PetscQuadrature q; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(sfe, PETSCFE_CLASSID, 1); PetscValidHeaderSpecific(tfe, PETSCFE_CLASSID, 2); ierr = PetscFEGetQuadrature(sfe, &q);CHKERRQ(ierr); ierr = PetscFESetQuadrature(tfe, q);CHKERRQ(ierr); ierr = PetscFEGetFaceQuadrature(sfe, &q);CHKERRQ(ierr); ierr = PetscFESetFaceQuadrature(tfe, q);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFEGetNumDof - Returns the number of dofs (dual basis vectors) associated to mesh points on the reference cell of a given dimension Not collective Input Parameter: . fem - The PetscFE object Output Parameter: . numDof - Array with the number of dofs per dimension Level: intermediate .seealso: PetscFECreate() @*/ PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(numDof, 2); ierr = PetscDualSpaceGetNumDof(fem->dualSpace, numDof);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFEGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points on the reference cell Not collective Input Parameter: . fem - The PetscFE object Output Parameter: . T - The basis function values and derivatives at quadrature points Note: $T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c$ 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 $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 Level: intermediate .seealso: PetscFECreateTabulation(), PetscTabulationDestroy() @*/ PetscErrorCode PetscFEGetCellTabulation(PetscFE fem, PetscTabulation *T) { PetscInt npoints; const PetscReal *points; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(T, 2); ierr = PetscQuadratureGetData(fem->quadrature, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr); if (!fem->T) {ierr = PetscFECreateTabulation(fem, 1, npoints, points, 1, &fem->T);CHKERRQ(ierr);} *T = fem->T; PetscFunctionReturn(0); } /*@C PetscFEGetFaceTabulation - Returns the tabulation of the basis functions at the face quadrature points for each face of the reference cell Not collective Input Parameter: . fem - The PetscFE object Output Parameters: . Tf - The basis function values and derviatives at face quadrature points Note:$ 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 $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$ 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 Level: intermediate .seealso: PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy() @*/ PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscTabulation *Tf) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(Tf, 2); if (!fem->Tf) { const PetscReal xi0[3] = {-1., -1., -1.}; PetscReal v0[3], J[9], detJ; PetscQuadrature fq; PetscDualSpace sp; DM dm; const PetscInt *faces; PetscInt dim, numFaces, f, npoints, q; const PetscReal *points; PetscReal *facePoints; ierr = PetscFEGetDualSpace(fem, &sp);CHKERRQ(ierr); ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); ierr = DMPlexGetConeSize(dm, 0, &numFaces);CHKERRQ(ierr); ierr = DMPlexGetCone(dm, 0, &faces);CHKERRQ(ierr); ierr = PetscFEGetFaceQuadrature(fem, &fq);CHKERRQ(ierr); if (fq) { ierr = PetscQuadratureGetData(fq, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr); ierr = PetscMalloc1(numFaces*npoints*dim, &facePoints);CHKERRQ(ierr); for (f = 0; f < numFaces; ++f) { ierr = DMPlexComputeCellGeometryFEM(dm, faces[f], NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); for (q = 0; q < npoints; ++q) CoordinatesRefToReal(dim, dim-1, xi0, v0, J, &points[q*(dim-1)], &facePoints[(f*npoints+q)*dim]); } ierr = PetscFECreateTabulation(fem, numFaces, npoints, facePoints, 1, &fem->Tf);CHKERRQ(ierr); ierr = PetscFree(facePoints);CHKERRQ(ierr); } } *Tf = fem->Tf; PetscFunctionReturn(0); } /*@C PetscFEGetFaceCentroidTabulation - Returns the tabulation of the basis functions at the face centroid points Not collective Input Parameter: . fem - The PetscFE object Output Parameters: . Tc - The basis function values at face centroid points Note: $T->T[0] = Bf[(f*pdim + i)*Nc + c] is the value at point f for basis function i and component c Level: intermediate .seealso: PetscFEGetFaceTabulation(), PetscFEGetCellTabulation(), PetscFECreateTabulation(), PetscTabulationDestroy() @*/ PetscErrorCode PetscFEGetFaceCentroidTabulation(PetscFE fem, PetscTabulation *Tc) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(Tc, 2); if (!fem->Tc) { PetscDualSpace sp; DM dm; const PetscInt *cone; PetscReal *centroids; PetscInt dim, numFaces, f; ierr = PetscFEGetDualSpace(fem, &sp);CHKERRQ(ierr); ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); ierr = DMPlexGetConeSize(dm, 0, &numFaces);CHKERRQ(ierr); ierr = DMPlexGetCone(dm, 0, &cone);CHKERRQ(ierr); ierr = PetscMalloc1(numFaces*dim, ¢roids);CHKERRQ(ierr); for (f = 0; f < numFaces; ++f) {ierr = DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);CHKERRQ(ierr);} ierr = PetscFECreateTabulation(fem, 1, numFaces, centroids, 0, &fem->Tc);CHKERRQ(ierr); ierr = PetscFree(centroids);CHKERRQ(ierr); } *Tc = fem->Tc; PetscFunctionReturn(0); } /*@C PetscFECreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided. Not collective Input Parameters: + fem - The PetscFE object . nrepl - The number of replicas . npoints - The number of tabulation points in a replica . points - The tabulation point coordinates - K - The number of derivatives calculated Output Parameter: . T - The basis function values and derivatives at tabulation points Note:$ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c $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$ 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 Level: intermediate .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy() @*/ PetscErrorCode PetscFECreateTabulation(PetscFE fem, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T) { DM dm; PetscDualSpace Q; PetscInt Nb; /* Dimension of FE space P */ PetscInt Nc; /* Field components */ PetscInt cdim; /* Reference coordinate dimension */ PetscInt k; PetscErrorCode ierr; PetscFunctionBegin; if (!npoints || !fem->dualSpace || K < 0) { *T = NULL; PetscFunctionReturn(0); } PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(points, 4); PetscValidPointer(T, 6); ierr = PetscFEGetDualSpace(fem, &Q);CHKERRQ(ierr); ierr = PetscDualSpaceGetDM(Q, &dm);CHKERRQ(ierr); ierr = DMGetDimension(dm, &cdim);CHKERRQ(ierr); ierr = PetscDualSpaceGetDimension(Q, &Nb);CHKERRQ(ierr); ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); ierr = PetscMalloc1(1, T);CHKERRQ(ierr); (*T)->K = !cdim ? 0 : K; (*T)->Nr = nrepl; (*T)->Np = npoints; (*T)->Nb = Nb; (*T)->Nc = Nc; (*T)->cdim = cdim; ierr = PetscMalloc1((*T)->K+1, &(*T)->T);CHKERRQ(ierr); for (k = 0; k <= (*T)->K; ++k) { ierr = PetscMalloc1(nrepl*npoints*Nb*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);CHKERRQ(ierr); } ierr = (*fem->ops->createtabulation)(fem, nrepl*npoints, points, K, *T);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFEComputeTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided. Not collective Input Parameters: + fem - The PetscFE object . npoints - The number of tabulation points . points - The tabulation point coordinates . K - The number of derivatives calculated - T - An existing tabulation object with enough allocated space Output Parameter: . T - The basis function values and derivatives at tabulation points Note: $T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c$ 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 $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 Level: intermediate .seealso: PetscFEGetCellTabulation(), PetscTabulationDestroy() @*/ PetscErrorCode PetscFEComputeTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T) { PetscErrorCode ierr; PetscFunctionBeginHot; if (!npoints || !fem->dualSpace || K < 0) PetscFunctionReturn(0); PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(points, 3); PetscValidPointer(T, 5); if (PetscDefined(USE_DEBUG)) { DM dm; PetscDualSpace Q; PetscInt Nb; /* Dimension of FE space P */ PetscInt Nc; /* Field components */ PetscInt cdim; /* Reference coordinate dimension */ ierr = PetscFEGetDualSpace(fem, &Q);CHKERRQ(ierr); ierr = PetscDualSpaceGetDM(Q, &dm);CHKERRQ(ierr); ierr = DMGetDimension(dm, &cdim);CHKERRQ(ierr); ierr = PetscDualSpaceGetDimension(Q, &Nb);CHKERRQ(ierr); ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 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); if (T->Nb != Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nb %D must match requested Nb %D", T->Nb, Nb); if (T->Nc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation Nc %D must match requested Nc %D", T->Nc, Nc); if (T->cdim != cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Tabulation cdim %D must match requested cdim %D", T->cdim, cdim); } T->Nr = 1; T->Np = npoints; ierr = (*fem->ops->createtabulation)(fem, npoints, points, K, T);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscTabulationDestroy - Frees memory from the associated tabulation. Not collective Input Parameter: . T - The tabulation Level: intermediate .seealso: PetscFECreateTabulation(), PetscFEGetCellTabulation() @*/ PetscErrorCode PetscTabulationDestroy(PetscTabulation *T) { PetscInt k; PetscErrorCode ierr; PetscFunctionBegin; PetscValidPointer(T, 1); if (!T || !(*T)) PetscFunctionReturn(0); for (k = 0; k <= (*T)->K; ++k) {ierr = PetscFree((*T)->T[k]);CHKERRQ(ierr);} ierr = PetscFree((*T)->T);CHKERRQ(ierr); ierr = PetscFree(*T);CHKERRQ(ierr); *T = NULL; PetscFunctionReturn(0); } PETSC_EXTERN PetscErrorCode PetscFECreatePointTrace(PetscFE fe, PetscInt refPoint, PetscFE *trFE) { PetscSpace bsp, bsubsp; PetscDualSpace dsp, dsubsp; PetscInt dim, depth, numComp, i, j, coneSize, order; PetscFEType type; DM dm; DMLabel label; PetscReal *xi, *v, *J, detJ; const char *name; PetscQuadrature origin, fullQuad, subQuad; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fe,PETSCFE_CLASSID,1); PetscValidPointer(trFE,3); ierr = PetscFEGetBasisSpace(fe,&bsp);CHKERRQ(ierr); ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); ierr = PetscDualSpaceGetDM(dsp,&dm);CHKERRQ(ierr); ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); ierr = DMPlexGetDepthLabel(dm,&label);CHKERRQ(ierr); ierr = DMLabelGetValue(label,refPoint,&depth);CHKERRQ(ierr); ierr = PetscCalloc1(depth,&xi);CHKERRQ(ierr); ierr = PetscMalloc1(dim,&v);CHKERRQ(ierr); ierr = PetscMalloc1(dim*dim,&J);CHKERRQ(ierr); for (i = 0; i < depth; i++) xi[i] = 0.; ierr = PetscQuadratureCreate(PETSC_COMM_SELF,&origin);CHKERRQ(ierr); ierr = PetscQuadratureSetData(origin,depth,0,1,xi,NULL);CHKERRQ(ierr); ierr = DMPlexComputeCellGeometryFEM(dm,refPoint,origin,v,J,NULL,&detJ);CHKERRQ(ierr); /* CellGeometryFEM computes the expanded Jacobian, we want the true jacobian */ for (i = 1; i < dim; i++) { for (j = 0; j < depth; j++) { J[i * depth + j] = J[i * dim + j]; } } ierr = PetscQuadratureDestroy(&origin);CHKERRQ(ierr); ierr = PetscDualSpaceGetPointSubspace(dsp,refPoint,&dsubsp);CHKERRQ(ierr); ierr = PetscSpaceCreateSubspace(bsp,dsubsp,v,J,NULL,NULL,PETSC_OWN_POINTER,&bsubsp);CHKERRQ(ierr); ierr = PetscSpaceSetUp(bsubsp);CHKERRQ(ierr); ierr = PetscFECreate(PetscObjectComm((PetscObject)fe),trFE);CHKERRQ(ierr); ierr = PetscFEGetType(fe,&type);CHKERRQ(ierr); ierr = PetscFESetType(*trFE,type);CHKERRQ(ierr); ierr = PetscFEGetNumComponents(fe,&numComp);CHKERRQ(ierr); ierr = PetscFESetNumComponents(*trFE,numComp);CHKERRQ(ierr); ierr = PetscFESetBasisSpace(*trFE,bsubsp);CHKERRQ(ierr); ierr = PetscFESetDualSpace(*trFE,dsubsp);CHKERRQ(ierr); ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); if (name) {ierr = PetscFESetName(*trFE, name);CHKERRQ(ierr);} ierr = PetscFEGetQuadrature(fe,&fullQuad);CHKERRQ(ierr); ierr = PetscQuadratureGetOrder(fullQuad,&order);CHKERRQ(ierr); ierr = DMPlexGetConeSize(dm,refPoint,&coneSize);CHKERRQ(ierr); if (coneSize == 2 * depth) { ierr = PetscDTGaussTensorQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);CHKERRQ(ierr); } else { ierr = PetscDTStroudConicalQuadrature(depth,1,(order + 1)/2,-1.,1.,&subQuad);CHKERRQ(ierr); } ierr = PetscFESetQuadrature(*trFE,subQuad);CHKERRQ(ierr); ierr = PetscFESetUp(*trFE);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&subQuad);CHKERRQ(ierr); ierr = PetscSpaceDestroy(&bsubsp);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode PetscFECreateHeightTrace(PetscFE fe, PetscInt height, PetscFE *trFE) { PetscInt hStart, hEnd; PetscDualSpace dsp; DM dm; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fe,PETSCFE_CLASSID,1); PetscValidPointer(trFE,3); *trFE = NULL; ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); ierr = PetscDualSpaceGetDM(dsp,&dm);CHKERRQ(ierr); ierr = DMPlexGetHeightStratum(dm,height,&hStart,&hEnd);CHKERRQ(ierr); if (hEnd <= hStart) PetscFunctionReturn(0); ierr = PetscFECreatePointTrace(fe,hStart,trFE);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ PetscFEGetDimension - Get the dimension of the finite element space on a cell Not collective Input Parameter: . fe - The PetscFE Output Parameter: . dim - The dimension Level: intermediate .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension() @*/ PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); PetscValidPointer(dim, 2); if (fem->ops->getdimension) {ierr = (*fem->ops->getdimension)(fem, dim);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@C PetscFEPushforward - Map the reference element function to real space Input Parameters: + fe - The PetscFE . fegeom - The cell geometry . Nv - The number of function values - vals - The function values Output Parameter: . vals - The transformed function values Level: advanced Note: This just forwards the call onto PetscDualSpacePushforward(). Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. .seealso: PetscDualSpacePushforward() @*/ PetscErrorCode PetscFEPushforward(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) { PetscErrorCode ierr; PetscFunctionBeginHot; ierr = PetscDualSpacePushforward(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFEPushforwardGradient - Map the reference element function gradient to real space Input Parameters: + fe - The PetscFE . fegeom - The cell geometry . Nv - The number of function gradient values - vals - The function gradient values Output Parameter: . vals - The transformed function gradient values Level: advanced Note: This just forwards the call onto PetscDualSpacePushforwardGradient(). Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension. .seealso: PetscFEPushforward(), PetscDualSpacePushforwardGradient(), PetscDualSpacePushforward() @*/ PetscErrorCode PetscFEPushforwardGradient(PetscFE fe, PetscFEGeom *fegeom, PetscInt Nv, PetscScalar vals[]) { PetscErrorCode ierr; PetscFunctionBeginHot; ierr = PetscDualSpacePushforwardGradient(fe->dualSpace, fegeom, Nv, fe->numComponents, vals);CHKERRQ(ierr); PetscFunctionReturn(0); } /* Purpose: Compute element vector for chunk of elements Input: Sizes: Ne: number of elements Nf: number of fields PetscFE dim: spatial dimension Nb: number of basis functions Nc: number of field components PetscQuadrature Nq: number of quadrature points Geometry: PetscFEGeom[Ne] possibly *Nq PetscReal v0s[dim] PetscReal n[dim] PetscReal jacobians[dim*dim] PetscReal jacobianInverses[dim*dim] PetscReal jacobianDeterminants FEM: PetscFE PetscQuadrature PetscReal quadPoints[Nq*dim] PetscReal quadWeights[Nq] PetscReal basis[Nq*Nb*Nc] PetscReal basisDer[Nq*Nb*Nc*dim] PetscScalar coefficients[Ne*Nb*Nc] PetscScalar elemVec[Ne*Nb*Nc] Problem: PetscInt f: the active field f0, f1 Work Space: PetscFE PetscScalar f0[Nq*dim]; PetscScalar f1[Nq*dim*dim]; PetscScalar u[Nc]; PetscScalar gradU[Nc*dim]; PetscReal x[dim]; PetscScalar realSpaceDer[dim]; Purpose: Compute element vector for N_cb batches of elements Input: Sizes: N_cb: Number of serial cell batches Geometry: PetscReal v0s[Ne*dim] PetscReal jacobians[Ne*dim*dim] possibly *Nq PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq PetscReal jacobianDeterminants[Ne] possibly *Nq FEM: static PetscReal quadPoints[Nq*dim] static PetscReal quadWeights[Nq] static PetscReal basis[Nq*Nb*Nc] static PetscReal basisDer[Nq*Nb*Nc*dim] PetscScalar coefficients[Ne*Nb*Nc] PetscScalar elemVec[Ne*Nb*Nc] ex62.c: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[], const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]), void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[]) ex52.c: 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) 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) ex52_integrateElement.cu __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec) PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], PetscLogEvent event, PetscInt debug, PetscInt pde_op) ex52_integrateElementOpenCL.c: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[], PetscLogEvent event, PetscInt debug, PetscInt pde_op) __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec) */ /*@C PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration Not collective Input Parameters: + prob - The PetscDS specifying the discretizations and continuum functions . field - The field being integrated . Ne - The number of elements in the chunk . cgeom - The cell geometry for each cell in the chunk . coefficients - The array of FEM basis coefficients for the elements . probAux - The PetscDS specifying the auxiliary discretizations - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements Output Parameter: . integral - the integral for this field Level: intermediate .seealso: PetscFEIntegrateResidual() @*/ PetscErrorCode PetscFEIntegrate(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) { PetscFE fe; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); if (fe->ops->integrate) {ierr = (*fe->ops->integrate)(prob, field, Ne, cgeom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@C PetscFEIntegrateBd - Produce the integral for the given field for a chunk of elements by quadrature integration Not collective Input Parameters: + prob - The PetscDS specifying the discretizations and continuum functions . field - The field being integrated . obj_func - The function to be integrated . Ne - The number of elements in the chunk . fgeom - The face geometry for each face in the chunk . coefficients - The array of FEM basis coefficients for the elements . probAux - The PetscDS specifying the auxiliary discretizations - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements Output Parameter: . integral - the integral for this field Level: intermediate .seealso: PetscFEIntegrateResidual() @*/ PetscErrorCode PetscFEIntegrateBd(PetscDS prob, PetscInt field, void (*obj_func)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), PetscInt Ne, PetscFEGeom *geom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) { PetscFE fe; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); if (fe->ops->integratebd) {ierr = (*fe->ops->integratebd)(prob, field, obj_func, Ne, geom, coefficients, probAux, coefficientsAux, integral);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@C PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration Not collective Input Parameters: + prob - The PetscDS specifying the discretizations and continuum functions . field - The field being integrated . Ne - The number of elements in the chunk . cgeom - The cell geometry for each cell in the chunk . coefficients - The array of FEM basis coefficients for the elements . coefficients_t - The array of FEM basis time derivative coefficients for the elements . probAux - The PetscDS specifying the auxiliary discretizations . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements - t - The time Output Parameter: . elemVec - the element residual vectors from each element Note:$ Loop over batch of elements (e): $Loop over quadrature points (q):$ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q $Call f_0 and f_1$ Loop over element vector entries (f,fc --> i): $elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u) Level: intermediate .seealso: PetscFEIntegrateResidual() @*/ PetscErrorCode PetscFEIntegrateResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) { PetscFE fe; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); if (fe->ops->integrateresidual) {ierr = (*fe->ops->integrateresidual)(prob, field, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@C PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary Not collective Input Parameters: + prob - The PetscDS specifying the discretizations and continuum functions . field - The field being integrated . Ne - The number of elements in the chunk . fgeom - The face geometry for each cell in the chunk . coefficients - The array of FEM basis coefficients for the elements . coefficients_t - The array of FEM basis time derivative coefficients for the elements . probAux - The PetscDS specifying the auxiliary discretizations . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements - t - The time Output Parameter: . elemVec - the element residual vectors from each element Level: intermediate .seealso: PetscFEIntegrateResidual() @*/ PetscErrorCode PetscFEIntegrateBdResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) { PetscFE fe; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); if (fe->ops->integratebdresidual) {ierr = (*fe->ops->integratebdresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@C PetscFEIntegrateHybridResidual - Produce the element residual vector for a chunk of hybrid element faces by quadrature integration Not collective Input Parameters: + prob - The PetscDS specifying the discretizations and continuum functions . field - The field being integrated . Ne - The number of elements in the chunk . fgeom - The face geometry for each cell in the chunk . coefficients - The array of FEM basis coefficients for the elements . coefficients_t - The array of FEM basis time derivative coefficients for the elements . probAux - The PetscDS specifying the auxiliary discretizations . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements - t - The time Output Parameter . elemVec - the element residual vectors from each element Level: developer .seealso: PetscFEIntegrateResidual() @*/ PetscErrorCode PetscFEIntegrateHybridResidual(PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) { PetscFE fe; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); if (fe->ops->integratehybridresidual) {ierr = (*fe->ops->integratehybridresidual)(prob, field, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, elemVec);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@C PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration Not collective Input Parameters: + prob - The PetscDS specifying the discretizations and continuum functions . jtype - The type of matrix pointwise functions that should be used . fieldI - The test field being integrated . fieldJ - The basis field being integrated . Ne - The number of elements in the chunk . cgeom - The cell geometry for each cell in the chunk . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point . coefficients_t - The array of FEM basis time derivative coefficients for the elements . probAux - The PetscDS specifying the auxiliary discretizations . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements . t - The time - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) Output Parameter: . elemMat - the element matrices for the Jacobian from each element Note:$ Loop over batch of elements (e): $Loop over element matrix entries (f,fc,g,gc --> i,j):$ Loop over quadrature points (q): $Make u_q and gradU_q (loops over fields,Nb,Ncomp)$ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) $+ \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)$ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) $+ \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) Level: intermediate .seealso: PetscFEIntegrateResidual() @*/ PetscErrorCode PetscFEIntegrateJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) { PetscFE fe; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); if (fe->ops->integratejacobian) {ierr = (*fe->ops->integratejacobian)(prob, jtype, fieldI, fieldJ, Ne, cgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@C PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration Not collective Input Parameters: + prob - The PetscDS specifying the discretizations and continuum functions . fieldI - The test field being integrated . fieldJ - The basis field being integrated . Ne - The number of elements in the chunk . fgeom - The face geometry for each cell in the chunk . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point . coefficients_t - The array of FEM basis time derivative coefficients for the elements . probAux - The PetscDS specifying the auxiliary discretizations . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements . t - The time - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) Output Parameter: . elemMat - the element matrices for the Jacobian from each element Note:$ Loop over batch of elements (e): $Loop over element matrix entries (f,fc,g,gc --> i,j):$ Loop over quadrature points (q): $Make u_q and gradU_q (loops over fields,Nb,Ncomp)$ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) $+ \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)$ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) $+ \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) Level: intermediate .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() @*/ PetscErrorCode PetscFEIntegrateBdJacobian(PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) { PetscFE fe; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); if (fe->ops->integratebdjacobian) {ierr = (*fe->ops->integratebdjacobian)(prob, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@C PetscFEIntegrateHybridJacobian - Produce the boundary element Jacobian for a chunk of hybrid elements by quadrature integration Not collective Input Parameters: . prob - The PetscDS specifying the discretizations and continuum functions . jtype - The type of matrix pointwise functions that should be used . fieldI - The test field being integrated . fieldJ - The basis field being integrated . Ne - The number of elements in the chunk . fgeom - The face geometry for each cell in the chunk . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point . coefficients_t - The array of FEM basis time derivative coefficients for the elements . probAux - The PetscDS specifying the auxiliary discretizations . coefficientsAux - The array of FEM auxiliary basis coefficients for the elements . t - The time - u_tShift - A multiplier for the dF/du_t term (as opposed to the dF/du term) Output Parameter . elemMat - the element matrices for the Jacobian from each element Note:$ Loop over batch of elements (e): $Loop over element matrix entries (f,fc,g,gc --> i,j):$ Loop over quadrature points (q): $Make u_q and gradU_q (loops over fields,Nb,Ncomp)$ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q) $+ \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)$ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q) \$ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) Level: developer .seealso: PetscFEIntegrateJacobian(), PetscFEIntegrateResidual() @*/ PetscErrorCode PetscFEIntegrateHybridJacobian(PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) { PetscFE fe; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1); ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); if (fe->ops->integratehybridjacobian) {ierr = (*fe->ops->integratehybridjacobian)(prob, jtype, fieldI, fieldJ, Ne, fgeom, coefficients, coefficients_t, probAux, coefficientsAux, t, u_tshift, elemMat);CHKERRQ(ierr);} PetscFunctionReturn(0); } /*@ PetscFEGetHeightSubspace - Get the subspace of this space for a mesh point of a given height Input Parameters: + fe - The finite element space - height - The height of the Plex point Output Parameter: . subfe - The subspace of this FE space Note: For example, if we want the subspace of this space for a face, we would choose height = 1. Level: advanced .seealso: PetscFECreateDefault() @*/ PetscErrorCode PetscFEGetHeightSubspace(PetscFE fe, PetscInt height, PetscFE *subfe) { PetscSpace P, subP; PetscDualSpace Q, subQ; PetscQuadrature subq; PetscFEType fetype; PetscInt dim, Nc; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1); PetscValidPointer(subfe, 3); if (height == 0) { *subfe = fe; PetscFunctionReturn(0); } ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); ierr = PetscFEGetFaceQuadrature(fe, &subq);CHKERRQ(ierr); ierr = PetscDualSpaceGetDimension(Q, &dim);CHKERRQ(ierr); 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);} if (!fe->subspaces) {ierr = PetscCalloc1(dim, &fe->subspaces);CHKERRQ(ierr);} if (height <= dim) { if (!fe->subspaces[height-1]) { PetscFE sub = NULL; const char *name; ierr = PetscSpaceGetHeightSubspace(P, height, &subP);CHKERRQ(ierr); ierr = PetscDualSpaceGetHeightSubspace(Q, height, &subQ);CHKERRQ(ierr); if (subQ) { ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), &sub);CHKERRQ(ierr); ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) sub, name);CHKERRQ(ierr); ierr = PetscFEGetType(fe, &fetype);CHKERRQ(ierr); ierr = PetscFESetType(sub, fetype);CHKERRQ(ierr); ierr = PetscFESetBasisSpace(sub, subP);CHKERRQ(ierr); ierr = PetscFESetDualSpace(sub, subQ);CHKERRQ(ierr); ierr = PetscFESetNumComponents(sub, Nc);CHKERRQ(ierr); ierr = PetscFESetUp(sub);CHKERRQ(ierr); ierr = PetscFESetQuadrature(sub, subq);CHKERRQ(ierr); } fe->subspaces[height-1] = sub; } *subfe = fe->subspaces[height-1]; } else { *subfe = NULL; } PetscFunctionReturn(0); } /*@ PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more sparsity). It is also used to create an interpolation between regularly refined meshes. Collective on fem Input Parameter: . fe - The initial PetscFE Output Parameter: . feRef - The refined PetscFE Level: advanced .seealso: PetscFEType, PetscFECreate(), PetscFESetType() @*/ PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef) { PetscSpace P, Pref; PetscDualSpace Q, Qref; DM K, Kref; PetscQuadrature q, qref; const PetscReal *v0, *jac; PetscInt numComp, numSubelements; PetscInt cStart, cEnd, c; PetscDualSpace *cellSpaces; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr); /* Create space */ ierr = PetscObjectReference((PetscObject) P);CHKERRQ(ierr); Pref = P; /* Create dual space */ ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr); ierr = PetscDualSpaceSetType(Qref, PETSCDUALSPACEREFINED);CHKERRQ(ierr); ierr = DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);CHKERRQ(ierr); ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr); ierr = DMPlexGetHeightStratum(Kref, 0, &cStart, &cEnd);CHKERRQ(ierr); ierr = PetscMalloc1(cEnd - cStart, &cellSpaces);CHKERRQ(ierr); /* TODO: fix for non-uniform refinement */ for (c = 0; c < cEnd - cStart; c++) cellSpaces[c] = Q; ierr = PetscDualSpaceRefinedSetCellSpaces(Qref, cellSpaces);CHKERRQ(ierr); ierr = PetscFree(cellSpaces);CHKERRQ(ierr); ierr = DMDestroy(&Kref);CHKERRQ(ierr); ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr); /* Create element */ ierr = PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);CHKERRQ(ierr); ierr = PetscFESetType(*feRef, PETSCFECOMPOSITE);CHKERRQ(ierr); ierr = PetscFESetBasisSpace(*feRef, Pref);CHKERRQ(ierr); ierr = PetscFESetDualSpace(*feRef, Qref);CHKERRQ(ierr); ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); ierr = PetscFESetNumComponents(*feRef, numComp);CHKERRQ(ierr); ierr = PetscFESetUp(*feRef);CHKERRQ(ierr); ierr = PetscSpaceDestroy(&Pref);CHKERRQ(ierr); ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr); /* Create quadrature */ ierr = PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);CHKERRQ(ierr); ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr); ierr = PetscFESetQuadrature(*feRef, qref);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFECreateDefault - Create a PetscFE for basic FEM computation Collective Input Parameters: + comm - The MPI comm . dim - The spatial dimension . Nc - The number of components . isSimplex - Flag for simplex reference cell, otherwise its a tensor product . prefix - The options prefix, or NULL - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree Output Parameter: . fem - The PetscFE object Note: Each object is SetFromOption() during creation, so that the object may be customized from the command line. Level: beginner .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() @*/ PetscErrorCode PetscFECreateDefault(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem) { PetscQuadrature q, fq; DM K; PetscSpace P; PetscDualSpace Q; PetscInt order, quadPointsPerEdge; PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; PetscErrorCode ierr; PetscFunctionBegin; /* Create space */ ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr); ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr); ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr); ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr); /* Create dual space */ ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr); ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); ierr = DMDestroy(&K);CHKERRQ(ierr); ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr); ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr); ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); /* Create element */ ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr); ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr); ierr = PetscFESetUp(*fem);CHKERRQ(ierr); ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); /* Create quadrature (with specified order if given) */ qorder = qorder >= 0 ? qorder : order; ierr = PetscObjectOptionsBegin((PetscObject)*fem);CHKERRQ(ierr); ierr = PetscOptionsBoundedInt("-petscfe_default_quadrature_order","Quadrature order is one less than quadrature points per edge","PetscFECreateDefault",qorder,&qorder,NULL,0);CHKERRQ(ierr); ierr = PetscOptionsEnd();CHKERRQ(ierr); quadPointsPerEdge = PetscMax(qorder + 1,1); if (isSimplex) { ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); } else { ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); } ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ PetscFECreateLagrange - Create a PetscFE for the basic Lagrange space of degree k Collective Input Parameters: + comm - The MPI comm . dim - The spatial dimension . Nc - The number of components . isSimplex - Flag for simplex reference cell, otherwise its a tensor product . k - The degree k of the space - qorder - The quadrature order or PETSC_DETERMINE to use PetscSpace polynomial degree Output Parameter: . fem - The PetscFE object Level: beginner Notes: 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. .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() @*/ PetscErrorCode PetscFECreateLagrange(MPI_Comm comm, PetscInt dim, PetscInt Nc, PetscBool isSimplex, PetscInt k, PetscInt qorder, PetscFE *fem) { PetscQuadrature q, fq; DM K; PetscSpace P; PetscDualSpace Q; PetscInt quadPointsPerEdge; PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; char name[64]; PetscErrorCode ierr; PetscFunctionBegin; /* Create space */ ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr); ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); /* Create dual space */ ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); ierr = DMDestroy(&K);CHKERRQ(ierr); ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr); ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); /* Create element */ ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); ierr = PetscSNPrintf(name, 64, "P%d", (int) k);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) *fem, name);CHKERRQ(ierr); ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr); ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); ierr = PetscFESetUp(*fem);CHKERRQ(ierr); ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); /* Create quadrature (with specified order if given) */ qorder = qorder >= 0 ? qorder : k; quadPointsPerEdge = PetscMax(qorder + 1,1); if (isSimplex) { ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); } else { ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); } ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C PetscFESetName - Names the FE and its subobjects Not collective Input Parameters: + fe - The PetscFE - name - The name Level: intermediate .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate() @*/ PetscErrorCode PetscFESetName(PetscFE fe, const char name[]) { PetscSpace P; PetscDualSpace Q; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFEGetBasisSpace(fe, &P);CHKERRQ(ierr); ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) P, name);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) Q, name);CHKERRQ(ierr); PetscFunctionReturn(0); } 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[]) { PetscInt dOffset = 0, fOffset = 0, f; PetscErrorCode ierr; for (f = 0; f < Nf; ++f) { PetscFE fe; const PetscInt cdim = T[f]->cdim; const PetscInt Nq = T[f]->Np; const PetscInt Nbf = T[f]->Nb; const PetscInt Ncf = T[f]->Nc; const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; PetscInt b, c, d; ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; for (b = 0; b < Nbf; ++b) { for (c = 0; c < Ncf; ++c) { const PetscInt cidx = b*Ncf+c; u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; } } ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); if (u_t) { for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; for (b = 0; b < Nbf; ++b) { for (c = 0; c < Ncf; ++c) { const PetscInt cidx = b*Ncf+c; u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; } } ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); } fOffset += Ncf; dOffset += Nbf; } return 0; } 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[]) { PetscInt dOffset = 0, fOffset = 0, g; PetscErrorCode ierr; for (g = 0; g < 2*Nf-1; ++g) { if (!T[g/2]) continue; { PetscFE fe; const PetscInt f = g/2; const PetscInt cdim = T[f]->cdim; const PetscInt Nq = T[f]->Np; const PetscInt Nbf = T[f]->Nb; const PetscInt Ncf = T[f]->Nc; const PetscReal *Bq = &T[f]->T[0][(r*Nq+q)*Nbf*Ncf]; const PetscReal *Dq = &T[f]->T[1][(r*Nq+q)*Nbf*Ncf*cdim]; PetscInt b, c, d; fe = (PetscFE) ds->disc[f]; for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0; for (d = 0; d < cdim*Ncf; ++d) u_x[fOffset*cdim+d] = 0.0; for (b = 0; b < Nbf; ++b) { for (c = 0; c < Ncf; ++c) { const PetscInt cidx = b*Ncf+c; u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b]; for (d = 0; d < cdim; ++d) u_x[(fOffset+c)*cdim+d] += Dq[cidx*cdim+d]*coefficients[dOffset+b]; } } ierr = PetscFEPushforward(fe, fegeom, 1, &u[fOffset]);CHKERRQ(ierr); ierr = PetscFEPushforwardGradient(fe, fegeom, 1, &u_x[fOffset*cdim]);CHKERRQ(ierr); if (u_t) { for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0; for (b = 0; b < Nbf; ++b) { for (c = 0; c < Ncf; ++c) { const PetscInt cidx = b*Ncf+c; u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b]; } } ierr = PetscFEPushforward(fe, fegeom, 1, &u_t[fOffset]);CHKERRQ(ierr); } fOffset += Ncf; dOffset += Nbf; } } return 0; } PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[]) { PetscFE fe; PetscTabulation Tc; PetscInt b, c; PetscErrorCode ierr; if (!prob) return 0; ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); ierr = PetscFEGetFaceCentroidTabulation(fe, &Tc);CHKERRQ(ierr); { const PetscReal *faceBasis = Tc->T[0]; const PetscInt Nb = Tc->Nb; const PetscInt Nc = Tc->Nc; for (c = 0; c < Nc; ++c) {u[c] = 0.0;} for (b = 0; b < Nb; ++b) { for (c = 0; c < Nc; ++c) { const PetscInt cidx = b*Nc+c; u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx]; } } } return 0; } PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) { const PetscInt dE = T->cdim; /* fegeom->dimEmbed */ const PetscInt Nq = T->Np; const PetscInt Nb = T->Nb; const PetscInt Nc = T->Nc; const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; PetscInt q, b, c, d; PetscErrorCode ierr; for (b = 0; b < Nb; ++b) elemVec[b] = 0.0; for (q = 0; q < Nq; ++q) { for (b = 0; b < Nb; ++b) { for (c = 0; c < Nc; ++c) { const PetscInt bcidx = b*Nc+c; tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; } } ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); for (b = 0; b < Nb; ++b) { for (c = 0; c < Nc; ++c) { const PetscInt bcidx = b*Nc+c; const PetscInt qcidx = q*Nc+c; elemVec[b] += tmpBasis[bcidx]*f0[qcidx]; for (d = 0; d < dE; ++d) elemVec[b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; } } } return(0); } PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE fe, PetscTabulation T, PetscInt r, PetscScalar tmpBasis[], PetscScalar tmpBasisDer[], PetscFEGeom *fegeom, PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[]) { const PetscInt dE = T->cdim; const PetscInt Nq = T->Np; const PetscInt Nb = T->Nb; const PetscInt Nc = T->Nc; const PetscReal *basis = &T->T[0][r*Nq*Nb*Nc]; const PetscReal *basisDer = &T->T[1][r*Nq*Nb*Nc*dE]; PetscInt q, b, c, d, s; PetscErrorCode ierr; for (b = 0; b < Nb*2; ++b) elemVec[b] = 0.0; for (q = 0; q < Nq; ++q) { for (b = 0; b < Nb; ++b) { for (c = 0; c < Nc; ++c) { const PetscInt bcidx = b*Nc+c; tmpBasis[bcidx] = basis[q*Nb*Nc+bcidx]; for (d = 0; d < dE; ++d) tmpBasisDer[bcidx*dE+d] = basisDer[q*Nb*Nc*dE+bcidx*dE+d]; } } ierr = PetscFEPushforward(fe, fegeom, Nb, tmpBasis);CHKERRQ(ierr); ierr = PetscFEPushforwardGradient(fe, fegeom, Nb, tmpBasisDer);CHKERRQ(ierr); for (s = 0; s < 2; ++s) { for (b = 0; b < Nb; ++b) { for (c = 0; c < Nc; ++c) { const PetscInt bcidx = b*Nc+c; const PetscInt qcidx = (q*2+s)*Nc+c; elemVec[Nb*s+b] += tmpBasis[bcidx]*f0[qcidx]; for (d = 0; d < dE; ++d) elemVec[Nb*s+b] += tmpBasisDer[bcidx*dE+d]*f1[qcidx*dE+d]; } } } } return(0); } 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[]) { const PetscInt dE = TI->cdim; const PetscInt NqI = TI->Np; const PetscInt NbI = TI->Nb; const PetscInt NcI = TI->Nc; const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; const PetscInt NqJ = TJ->Np; const PetscInt NbJ = TJ->Nb; const PetscInt NcJ = TJ->Nc; const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; PetscInt f, fc, g, gc, df, dg; PetscErrorCode ierr; for (f = 0; f < NbI; ++f) { for (fc = 0; fc < NcI; ++fc) { const PetscInt fidx = f*NcI+fc; /* Test function basis index */ tmpBasisI[fidx] = basisI[fidx]; for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; } } ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); for (g = 0; g < NbJ; ++g) { for (gc = 0; gc < NcJ; ++gc) { const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ tmpBasisJ[gidx] = basisJ[gidx]; for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; } } ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); for (f = 0; f < NbI; ++f) { for (fc = 0; fc < NcI; ++fc) { const PetscInt fidx = f*NcI+fc; /* Test function basis index */ const PetscInt i = offsetI+f; /* Element matrix row */ for (g = 0; g < NbJ; ++g) { for (gc = 0; gc < NcJ; ++gc) { const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ const PetscInt j = offsetJ+g; /* Element matrix column */ const PetscInt fOff = eOffset+i*totDim+j; elemMat[fOff] += tmpBasisI[fidx]*g0[fc*NcJ+gc]*tmpBasisJ[gidx]; for (df = 0; df < dE; ++df) { elemMat[fOff] += tmpBasisI[fidx]*g1[(fc*NcJ+gc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(fc*NcJ+gc)*dE+df]*tmpBasisJ[gidx]; for (dg = 0; dg < dE; ++dg) { elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((fc*NcJ+gc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; } } } } } } return(0); } 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[]) { const PetscInt dE = TI->cdim; const PetscInt NqI = TI->Np; const PetscInt NbI = TI->Nb; const PetscInt NcI = TI->Nc; const PetscReal *basisI = &TI->T[0][(r*NqI+q)*NbI*NcI]; const PetscReal *basisDerI = &TI->T[1][(r*NqI+q)*NbI*NcI*dE]; const PetscInt NqJ = TJ->Np; const PetscInt NbJ = TJ->Nb; const PetscInt NcJ = TJ->Nc; const PetscReal *basisJ = &TJ->T[0][(r*NqJ+q)*NbJ*NcJ]; const PetscReal *basisDerJ = &TJ->T[1][(r*NqJ+q)*NbJ*NcJ*dE]; const PetscInt Ns = isHybridI ? 1 : 2; const PetscInt Nt = isHybridJ ? 1 : 2; PetscInt f, fc, g, gc, df, dg, s, t; PetscErrorCode ierr; for (f = 0; f < NbI; ++f) { for (fc = 0; fc < NcI; ++fc) { const PetscInt fidx = f*NcI+fc; /* Test function basis index */ tmpBasisI[fidx] = basisI[fidx]; for (df = 0; df < dE; ++df) tmpBasisDerI[fidx*dE+df] = basisDerI[fidx*dE+df]; } } ierr = PetscFEPushforward(feI, fegeom, NbI, tmpBasisI);CHKERRQ(ierr); ierr = PetscFEPushforwardGradient(feI, fegeom, NbI, tmpBasisDerI);CHKERRQ(ierr); for (g = 0; g < NbJ; ++g) { for (gc = 0; gc < NcJ; ++gc) { const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ tmpBasisJ[gidx] = basisJ[gidx]; for (dg = 0; dg < dE; ++dg) tmpBasisDerJ[gidx*dE+dg] = basisDerJ[gidx*dE+dg]; } } ierr = PetscFEPushforward(feJ, fegeom, NbJ, tmpBasisJ);CHKERRQ(ierr); ierr = PetscFEPushforwardGradient(feJ, fegeom, NbJ, tmpBasisDerJ);CHKERRQ(ierr); for (s = 0; s < Ns; ++s) { for (f = 0; f < NbI; ++f) { for (fc = 0; fc < NcI; ++fc) { const PetscInt sc = NcI*s+fc; /* components from each side of the surface */ const PetscInt fidx = f*NcI+fc; /* Test function basis index */ const PetscInt i = offsetI+NbI*s+f; /* Element matrix row */ for (t = 0; t < Nt; ++t) { for (g = 0; g < NbJ; ++g) { for (gc = 0; gc < NcJ; ++gc) { const PetscInt tc = NcJ*t+gc; /* components from each side of the surface */ const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ const PetscInt j = offsetJ+NbJ*t+g; /* Element matrix column */ const PetscInt fOff = eOffset+i*totDim+j; elemMat[fOff] += tmpBasisI[fidx]*g0[sc*NcJ*Nt+tc]*tmpBasisJ[gidx]; for (df = 0; df < dE; ++df) { elemMat[fOff] += tmpBasisI[fidx]*g1[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisDerJ[gidx*dE+df]; elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g2[(sc*NcJ*Nt+tc)*dE+df]*tmpBasisJ[gidx]; for (dg = 0; dg < dE; ++dg) { elemMat[fOff] += tmpBasisDerI[fidx*dE+df]*g3[((sc*NcJ*Nt+tc)*dE+df)*dE+dg]*tmpBasisDerJ[gidx*dE+dg]; } } } } } } } } return(0); } PetscErrorCode PetscFECreateCellGeometry(PetscFE fe, PetscQuadrature quad, PetscFEGeom *cgeom) { PetscDualSpace dsp; DM dm; PetscQuadrature quadDef; PetscInt dim, cdim, Nq; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); ierr = PetscDualSpaceGetDM(dsp, &dm);CHKERRQ(ierr); ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); ierr = PetscFEGetQuadrature(fe, &quadDef);CHKERRQ(ierr); quad = quad ? quad : quadDef; ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); ierr = PetscMalloc1(Nq*cdim, &cgeom->v);CHKERRQ(ierr); ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->J);CHKERRQ(ierr); ierr = PetscMalloc1(Nq*cdim*cdim, &cgeom->invJ);CHKERRQ(ierr); ierr = PetscMalloc1(Nq, &cgeom->detJ);CHKERRQ(ierr); cgeom->dim = dim; cgeom->dimEmbed = cdim; cgeom->numCells = 1; cgeom->numPoints = Nq; ierr = DMPlexComputeCellGeometryFEM(dm, 0, quad, cgeom->v, cgeom->J, cgeom->invJ, cgeom->detJ);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode PetscFEDestroyCellGeometry(PetscFE fe, PetscFEGeom *cgeom) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFree(cgeom->v);CHKERRQ(ierr); ierr = PetscFree(cgeom->J);CHKERRQ(ierr); ierr = PetscFree(cgeom->invJ);CHKERRQ(ierr); ierr = PetscFree(cgeom->detJ);CHKERRQ(ierr); PetscFunctionReturn(0); }