# suggested improvement for pyramidal elements

Greetings,

I’ve started the task of implementing the Biquadratic pyramid in VTK (vtkBiQuadraticPyramid with 14 nodes), and noticed that both vtkQuadraticPyramid and vtkPyramid have issues when computing the Jacobian inverse at the apex.

This is due to how the derivatives are implemented: define Nk as the interpolation functions over the reference pyramid, and T as the Bedrosian transform.

The functions over a unit cube can be found by taking `Nk o T` (coordinate transform), which is what you’ve done. For the derivatives however, you’ve taken them over the space of the unit cube itself (say `xq`, `yq`, `zq`), leading to your issues with the Jacobian inverse. I propose instead to apply `dNq = grad(Nq,xk)`, that is, take gradients of the cube shape functions on the space of the pyramid (`xk`, `yk`, `zk`), which is accomplished by first computing `dNk = grad(Nk,xk)` and then `dNq = dNk o T` (apply Bedrosian transform over the derivatives in the reference pyramid).

The implication of this method is that now integrals have to be computed according to the paper of M. Bergot: an integral of fk on the pyramidal space is the same as an integral of `4*fq*((1-zq)^2)` on the unit cube, and this integral over the unit cube itself can easily be passed to the more common Gaussian space `(-1,1)`, or isoparametric cube `Qr (xr,yr,zr).`

The advantage is that now it is possible to properly compute the Jacobian and it’s inverse at the apex without any special “tricks”. I have computed `N` and `dN` for both PYR05 and PYR14, and can provide the symbolic results needed.

Regarding references, I’ve used the paper “Higher-order finite elements for hybrid meshes using new nodal pyramidal elements” by M. Bergot, G. Cohen and M. Durufle. specifically, the quadrature rule can be found on chapter 3, equation 12. I have tested this with both exact and approximate quadratures for a Laplacian operator over a reference pyramid for both PYR05 and PYR14, as well as for the creation of the consistent mass matrix for a deformed affine PYR05 element. In both cases, the 3 quadrature forms i propose give the exact same result while still allowing an exact Jacobian inversion at the apex.

In any case, thank you for your attention

1 Like

Wow, the amount of work and thought that went into this is very impressive, thank you. There is no doubt that some of VTK’s higher order cells can be improved, so I personally am excited by any chance to improve what’s there. The biggest concern I have is how it might compare / affect the past implementation, and what impact it might have on any existing users. Secondary concerns include performance impacts, or the addition of dependencies on any external numerical libraries etc - although these seem not to be the case based on what you wrote. But my sense is that if there are no serious objections from the community we should move forward and integrate your work.

1 Like

(Edited for clarity)

Could you open a MR to share your code ?

See here for the process : https://gitlab.kitware.com/vtk/vtk/-/blob/master/Documentation/dev/git/develop.md

Dear Dr. Schroeder,

Thank you for your kind comment. We’ve started some ground work on high-order FEM implementation at BSC, and using VTK/ParaView to visualize the cells help us quite a lot. As far as I know, as long as we provide my suggested chenges to N, dN and the integration rule (if applicable), previous results should not be affected, since `int(dNq/dxq*detJeq(dNq/dxq) dQ) = int(4*dNq/dxk*((1-zq)**2)*detJeq(dNq/dxk) dQ)`. The major difference would be the ability to properly compute jacobians at the apex, without use of approximations and tolerances, although it must be noted that a Jacobian based on my dNq/dxk is different from Je on reference K itself. I do believe, though i have not yet tested, that this would improve the point finding algorithm as well, since in theory we shouldn’t have issues at the apex anymore (no rationals involved). This proposal is also true for the Linear pyramid as well, meaning changes to vtkPyramid, vtkQuadraticPyramid, and introduction of vtkBiquadraticPyramid, which i’m working on still. regarding code submission, I’m working on VTK 9.0.1 since the latest Git version does not compile on my PC for some reason. Still, the functions have not changed much between versions, so i could send the files i’ve been working on. Performancewise, I believe it should be the same, since we are pre-computing these spatial transforms, meaning no extra computation is needed. It does not require any external numerical library, all changes are to VTK files themselves.