mlarson:
You made the statement: "Wouldn't you miss what is happening in the center of that element since p-elements would only compute along the edges?" In general this statement is not true--in p-element codes, you can compute engineering data such as stresses, etc. along the edges and in the internal regions, and there is no such intrinsic limitation in the p-element method. Not being familiar with Mechanica, I can only guess that you can also compute internal engr. data with Mechanica, but perhaps you just can't figure out how to get the post processor to do it for you (it happens to the best of us!

). I know internal stresses can be computed with StressCheck, another p-element code, and I computed internal stresses myself in two research quality p-element codes I wrote.
One point about the elements themselves--in isoparametric h-elements, the nodes and the shape functions are tied directly together, so that the shape functions are defined as value of 1.0 at the assigned node, and 0.0 as the shape function passes through all the other nodes. In p-elements, only the corner nodes have shape functions assigned directly to them; the edges do not have nodes, but there are shape functions associated to the edges. In addition, there are shape functions associated with the center region of the elements. If you restrict the polynomial level 'p' to 2 (that is, quadratic), then in principle you can have the same number of degrees of freedom (DOFs) on the p-elements as you do in the h-elements (for 2D quadratic elements, that's 16 or 18 DOFs per each element, no matter the method used, p-element or h-element). However, in the p-element method, you can keep increasing polynomial level of the shape functions to arbitrarily large value--normally you are restricted only by the size and speed of your computer. So why choose one over the other? If you care anything at all about the reliability and accuracy of the Finite element solution, then pick a p-element code, unless the p-element code does not model the physics of your problem well (which is a valid criterion for selecting any particular piece of engineering software).
One last point about integration points. The integration points are the Gauss points for the integration of all the stiffness integrals you need to compute over each of the elements; you are merely approximating the exact integration of these integrals with a finite sum of pairs of weight fcns multiplied the integral argument evaluated at the Gauss points. In principle then there is no limitation on the number of Gauss points you could use (or, for that matter, the polynomial order, p--typically the maximum 'p' is 8 or 9, but I've seen up to p=14!). It just so happens the optimal number of Gauss points is related to the polynomial order of the integral arguments, but that there is no reason other than personal preference of the code developers for restricting the number of Gauss points to say 9 along the edges.