I just fixed up the quadrature rules for pyramids, which have been
wrong since I first coded them up years ago!! If you've ever done
anything with pyramids in the library, it was probably at the very
least inaccurate and (most likely) totally wrong :)
Anyway, I actually wanted to start a thread about (the lack of)
pyramid refinement in the library. I've been working with some viz
folks here at TACC recently, and after telling them "yeah we can do
pyramids" I actually went back and checked... Anyway, there appears
to be at least a couple issues to discuss:
1.) A possible refinement pattern for the pyramid involves 10
children: 4 sub-pyramids on the base, one in the apex, and one pyramid
upside-down below the apex pyramid. Then, there are 4 tets filling in
the gaps. (This is a little hard to visualize, but I can send you an
image if you want to see what I mean.) One bad thing about this
refinement pattern is that under repeated refinements, you end up with
a little "ribbon" of pyramids running through the original element.
So, we probably need to discuss other possibilities, but this one
would at least get us started.
2.) *Any* refinement pattern for the pyramid will likely involve a
splitting where some children have a different type than the parent.
Since this would be the only element with such an inhomogeneous
splitting, I'd like to discuss the best way to do this in the present
context of the embedding matrix. Since embedding_matrix() is a
virtual function, it can be redefined in Pyramid sub-classes to do
whatever we want, e.g. we could place -1 entries in the actual tables
and just ensure, through embedding_matrix(), that those entries are
never accessed. The real issue seems to be determining where in the
library, if anywhere, we've assumed that all the children have the
same geometric type. I'd be glad to hear your comments there.
3.) The pyramid basis functions are a little weird (rational), and,
for example, their evaluation at the apex is a little touchy. I need
to do more testing on these, including checking the accuracy of
integration for e.g. a mass matrix, and making sure the element
Jacobian is well-defined everywhere. Anyone know anything about the
accuracy of quadrature for functions which are ratios of polynomials?
I asked Roy about some sort of macroelement splitting to get away
from the rational basis functions, but I don't think that idea can
really go anywhere. If you have read anything about other types of
pyramid finite elements besides the one we've got, please let me know.
Finally, to be more useful in incompressible Navier-Stokes
calculations, were are gonna need a quadratic pyramid at some point.
I have a paper on how to define arbitrary-order Lagrange pyramids with
rational basis functions, so it's on my long term development plan.
Your comments on any of the above are welcomed.