While the triangular element is simple and can be used to mesh any two-dimensional geometry, its accuracy is limited by the assumption of constant strain. To improve upon this, we introduce the 4-node rectangular element. This element allows for a more complex strain distribution, leading to more accurate results.

1. Displacement Field
For this element, the displacement field (e.g., for the u-displacement) is approximated using a four-term polynomial. A crucial addition is the xy term, which allows for a non-constant strain distribution.
This can be written in matrix form:
By enforcing this equation at each of the four nodes, we can solve for the coefficients a in terms of the nodal displacements u. This leads to the familiar relationship
2. Shape Functions
Instead of a brute-force matrix inversion, the shape functions for a rectangular element can be constructed elegantly through the product of one-dimensional linear interpolation functions. Consider a rectangle with dimensions a and b. We can define simple linear functions in each direction:
- In the x-direction:
- In the y-direction:
The two-dimensional shape functions are then formed by taking products of these 1D functions. For a node i, the shape function is the product of the 1D functions that are equal to 1 at that node.

3. Strain-Displacement Matrix and Stiffness
Since the shape functions now contain x and y terms, their derivatives are no longer constant. For example, for
The B matrix will now contain terms that are functions of x and y. This means that the strain, given by
The element stiffness matrix is calculated using the standard formula:
Because the B matrix is a function of x and y, the integrand is no longer a constant and the integration must be performed explicitly.