Implement caching of the element stiffness matrices

Since all element geometries are identical, the element stiffness
matrices only depend on the phases at the quadrature nodes.
We use those to access a cache.  This leads to tremendous time savings.