I have not implemented these elements before, but I like their reduced cardinality compared to (e.g.) a tensor product of lagrange interpolants, which is very "overcomplete" (especially for orders>2)

I think they are not complete. My test program takes a quad and builds the 14x14 mass (projection?) matrix consisting the 8 serendipity functions + the 6 functions that comprise the second order polynomial basis (1 x y xx xy yy). Matlab reports that matrix has a rank of 9 (not 8), indicating there is a function inside the polynomial basis that the serendipity functions fail to capture. The thing is, the 9'th eigenvalue is like 10^5 down from the others, the serendipity basis "almost" captures all the polynomials, just barely doesn't. When I deliberately make the quad rectangular, the rank is reported as 8, the 9th eigenvalue vanishes to machine precision. That 9th eigenvalue grows as the quad is skewed more and more.

My suspicion has grown further upon reexamination of my reference material more closely (Huebner's "The Finite Element Method for Engineers", pp 167-171) - the serendipity functions are in the section entitled "Rectangular Elements" (e.g. orthogonal edges), and there's no mention of them in the "Arbitrary Quadrilateral Elements".

I have been googling and regoogling but thought somebody here might know off the top of their head, being more familiar with the space (seems like its more common in structures - I'm more an electromagnetics guy).

I'm not a fan of serendipity spaces. The number of degrees of freedom saved is not very high (especially in 3D) and their effect is easily removed in the first step of a solver, so they don't actually cost much. For example, a sparse direct solver essentially does static condensation of the internal nodes in its first step. The overall cost does not depend on those degrees of freedom because it does not change the minimal vertex separators. Those rows of the matrix also have fewer degrees of freedom than the interface nodes. Also, leaving the full tensor product enables performance optimizations in local computations.

Unfortunately these are going to be used in a boundary element method (no condensation... right?) and I was trying to pinch every penny. Thanks nonetheless for the very informative response!
–
rchilton1980Jan 31 '13 at 23:17