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.

Use the usual bounded Lipschitz domain assumptions so that the Sobolev trace is defined. Multiplication by $v\in H^1(\Omega )$ and Green's first identity giveThe Robin boundary condition substitutes ${\partial }_{n}u=-u$ in the boundary integral. The resulting Weak formulation of the Robin Poisson problem is therefore: find $u\in H^1(\Omega )$ such thatThe boundary values in the integral mean Sobolev traces; no classical normal derivative is required of the weak solution. For uniqueness, subtract two solutions and test by their difference $w$. This yields ${\int }_{\Omega }|\nabla w{|}^2+{\int }_{\partial \Omega }{w}^2=0$. Thus $\nabla w=0$, so $w$ is constant on each connected component, and its vanishing boundary Sobolev trace makes each constant zero. The weak solution is unique.

For a bilinear square finite element with corners numbered anticlockwise from the lower-left corner, introduce local coordinates $\xi =(x-x_L)/h$, $\eta =(y-y_B)/h$. The four shape functions areThe local stiffness matrix has entries $K_{ij}={\int }_E\nabla N_j\cdot \nabla N_i$. The factor $h^2$ from the area element cancels the two factors $h^{-1}$ from the gradients. Symmetry leaves three types: diagonal, adjacent corners, and opposite corners. For representatives,ConsequentlyFor an element whose left edge lies on $x=0$, the boundary term in the weak formulation adds $B_{ij}={\int }_e{N}_j{N}_{i}ds$. Only $N_1=1-\eta$ and $N_4=\eta$ are nonzero on that edge. Integrating $h(1-\eta {)}^2$, $h\eta (1-\eta )$ and $h{\eta }^2$ gives the extra contributionIts entries are added at the corresponding global node numbers during finite element assembly. Additional boundary edges contribute their own analogous edge matrices. The Robin boundary condition is represented by these integrals, so boundary nodes remain unknowns.

A shape function is a local basis function associated with a degree of freedom of a finite element. For nodal degrees of freedom it is one at its own node and zero at the other nodes, allowing local interpolation from the nodal values. For a bilinear square finite element, the four products of linear coordinate functions provide the nodal shape functions.