Shape functions do not determine the shape of the element. Rather, shape functions are determined (at least partially) by the shape of the element.
The set of shape functions are ultimately determined by the degrees of freedom needed to cover interpolation of a given polynomial order for all points within the element as well as on the boundaries of the element.
For instance, consider a polynomial space of order 1 (linear function space). In 1d, a linear function space is spanned by 2 functions (e.g. 1 and x, but lagrange functions are more useful). Thus there are 2 degrees of freedom, corresponding to the coefficients assigned to the two functions that span this space.
In 1D, all “elements” have the same shape: a straight line. This makes it very easy to make the interpolation cover the entire element. If we choose the degrees of freedom to correspond to the nodal values at the boundary points (i.e. using lagrange functions), then interpolation between these two points covers the entire element.
But in 2d, things change a lot. A linear space is now spanned by 3 functions (i.e. 1, x, and y). So, a space spanned by 3 functions requires 3 degrees of freedom. On a 2D triangle, we need to choose these degrees of freedom in such a way that interpolation across the entire element, particularly across all the boundaries of the element. If we chose the degrees of freedom to be nodal values at the midpoints between the corners instead of corner points, this would not work. Interpolation between the degrees of freedom would only cover the middle of the element, not the corners. One would have to choose the degrees of freedom at the 3 corner points of the element. This establishes linear interpolation within the element, as well as the boundaries of the element.
But a triangle is not the only possible element shape in 2D. Now contrast the triangular element with a quadrilateral element. Technically speaking, a linear function space is still spanned by 3 functions. But no matter which 3 corners of the element you choose, it is impossible to establish linear interpolation to cover the remaining corner. To overcome this, we modify the function space to include an additional basis function while still preserving the polynomial degree order. For quadrilaterals, we use what is called “bilinear interpolation”. It is a space spanned by four functions (e.g. 1, x, y, and x*y) whose degrees of freedom can be chosen at the corners of the element so that interpolation can span all points within the element and along the boundaries of the element. Note that in either spatial dimension (x or y), the polynomial order is still linear.
We can certainly choose a higher order polynomial space to do interpolation, which would require more functions and therefore more degrees of freedom. The degrees of freedom do not necessarily have to all be corresponding to point values on the boundaries. But the number of shape functions depends on the spatial dimension size, and the polynomial functions used for interpolation. These polynomial functions are chosen in such a way as to cover the interpolation order across the entire shape of the element.
Interpolation on other shapes (e.g. pentagons, hexagons, etc.) and in 3 dimensions (e.g. tetrahedra, hexahedra, etc) is possible but requires similar considerations: the polynomial interpolation functions (and degrees of freedom) need to be chosen to preserve the polynomial order interpolation within the element and along its boundaries.
So what exactly are shape functions? They are the set of functions who, together, accomplish interpolation of a given order across the entire shape of the element (interior and boundary points).
As for the differences among “approximation” functions, “interpolation” functions and “shape” functions, your confusion is well merited. I have seen each of these terms used synonymously. Though the math remains the same, different authors emphasize different labels for different contexts. I believe most (not necessarily all) refer to the “shape” functions as the basis functions for a reference element as opposed to the basis function on the global coordinate system.
