suggested improvement for pyramidal elements

Development
1

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.

Please, let me know what are your thoughts on this, and if you have any more questions about this proposal.

In any case, thank you for your attention

1 Like
2

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
3

(Edited for clarity)

4

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

5

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.

Thank you for your attention

6

Greetings Dr. westphal,

I havent finished the code on VTK yet: all my testing has been done with Octave’s symbolical features, and so far i’ve only added/modified things related to the Biquadratic pyramid (PYR14). Also, due to compilation issues with the Git version, I’m using an older VTK 9.0.1 extracted from the .tar files, so it’s a bit out-of date, though most of my changes are quite self contained. I couldnt find any Gaussian quadrature functions on VTK, could you tell me if they exist and where would they be so i can provide all necessary changes? What I can provide at this time is an octave file that symbolically computes all sets of (N, dN) at the reference pyramid K, the unit cube Q and the reference cube Qr (-1,1), while also providing comparative integration tests. Please let me know your opinion on this.

7

Please open an issue. I have not idea how to answer the rest of your questions sadly.

8

I will, thank you very much. I have been quite busy with this, so i’ve had no time to open the issue regarding compilation. Regarding the quadratures, I’ll look into the unity tests and examples to see if i find something related to this. About the code sharing, as soon as i have the VTK files up and running properly (before today i hope :smiley: ) I’ll submit them through the proper channel, but for now i’d be happy to share my algebraic results so that you can verify them.