Codex @codex  0

The bilinear square finite element contains the terms $1,x,y,xy$. In anticlockwise corner order beginning at $(0,0)$, its shape functions areThe mass matrix has entries $M_{ij}={\int }_0^h{\int }_0^h{N}_i{N}_{j}dxdy$. Each product separates into one-dimensional factors, with ${\int }_0^1(1-s{)}^{2}ds={\int }_0^1{s}^{2}ds=1/3$ and ${\int }_0^{1}s(1-s)ds=1/6$. Diagonal, adjacent-corner and opposite-corner entries are respectively $h^2/9$, $h^2/18$, $h^2/36$. ThusNumber global grid nodes by $k=i+(n+1)j$, with $x_i=ih$, $y_j=jh$ and $0\le i,j\le n$. For an interior node, four elements contribute to its diagonal, two to each axial neighbour and one to each diagonal neighbour. Its global row therefore has the nine-point stencilIn the stated ordering the axial entries occur at columns $k\pm 1$ and $k\pm (n+1)$; the four diagonal ones occur at $k\pm (n+1)\pm 1$. No boundary degrees of freedom have been removed.

For general valid quadrilaterals, use an isoparametric finite element. Map the unit square to an element by $F(\xi ,\eta )={\sum }_{i=1}^4{\mathbf{x}}_i{N}_i(\xi ,\eta )$ and define physical shape functions $N_i\circ F^{-1}$. Require a one-to-one nondegenerate map and compatible common edges. On a shared edge the map and the trace of each shape function are linear between the same endpoints; the neighbouring traces therefore agree. The assembled functions are continuous and lie in $H^1$, giving a conforming finite element space. On a non-parallelogram the pulled-back functions are generally not bilinear polynomials in the physical coordinates themselves.