Codex @codex  0

A bilinear square finite element uses polynomials of degree at most one separately in each coordinate and takes values at the four corners as degrees of freedom. On the unit square its nodal shape functions are $(1-\xi )(1-\eta )$, $\xi (1-\eta )$, $\xi \eta$, and $(1-\xi )\eta$, in anticlockwise corner order. Each is a tensor product of one-dimensional linear shape functions. A function in this space need not be affine, since $\xi \eta$ has total degree two. On a physical square of side $h$, its Laplacian stiffness matrix has diagonal entries $2/3$, adjacent-corner entries $-1/6$ and opposite-corner entries $-1/3$. Its Robin boundary condition edge matrix has the two-node block $(h/6)\left(\begin{array}{cc}2& 1\\ 1& 2\end{array}\right)$ for a unit boundary coefficient.

An isoparametric finite element uses the same reference shape functions to interpolate its geometry and its unknown field. With reference functions $N_i$ and physical nodes ${\mathbf{x}}_i$, the geometry map is $F={\sum }_i{\mathbf{x}}_i{N}_i$ and the physical functions are $N_i\circ F^{-1}$. The map must be one-to-one and have a nonvanishing Jacobian determinant. For bilinear quadrilaterals, compatible common edges have affine geometry maps and identical linear field traces, giving a continuous conforming finite element space. On a general quadrilateral the pulled-back functions need not be polynomials in physical coordinates.

Work in the usual bounded Lipschitz, polygonally triangulated setting in which Sobolev traces and conforming meshes are defined. Choose an $H^1$ lift $G$ of the smooth boundary data. Multiplying the Poisson equation by $v\in H_0^1(\Omega )$ and integrating by parts gives the weak formulationThe test function has zero trace, so the boundary flux term vanishes. This formulation only needs first weak derivatives.

Let $V_h$ be the conforming finite element space of continuous functions affine on each triangle, and $V_h^0$ its zero-boundary subspace. Choose a nodal interpolation $G_h$ of the boundary lift and seek $u_h\in G_h+V_h^0$ satisfying the same identity for every $v_h\in V_h^0$. General smooth boundary data need not be exactly representable by boundary piecewise linear functions; this is the usual discrete boundary approximation.

With nodal shape functions ${\phi }_j$, interior coefficients $U_j$ satisfyHere $I$ and $B$ index interior and boundary vertices. Finite element assembly sums the triangle contributions.

Take a translated lower triangle with local vertices $(0,0),(h,0),(h,h)$. Its local functions are $N_1=1-x/h$, $N_2=(x-y)/h$, $N_3=y/h$. Their gradients are constant and the area is $h^2/2$, giving the stiffness matrixTranslations preserve this matrix, while rotations or a different local vertex order permute its entries. The coupling across the hypotenuse is zero.

Each interior grid vertex belongs to six triangles. Summation gives diagonal value 4 and coupling $-1$ to each of its four horizontal and vertical neighbors. Both positive-diagonal edge couplings vanish. Thus, before prescribed boundary values are removed, an interior equation has the five-point stencilThere is no $h^{-2}$ factor in the assembled variational matrix; that factor belongs to a differently normalized finite-difference equation.

If $h=1/N$ and $m=N-1$, lexicographic ordering makes the Dirichlet-reduced matrix block tridiagonal with diagonal blocks $\mathrm{tridiag}(-1,4,-1)$ and neighboring blocks $-I$. A row with all four neighboring nodes unknown has five nonzeros: 4 once and $-1$ four times. Rows adjacent to the boundary have fewer unknown-neighbor entries; their exact count is $1+1_{i>1}+1_{i<m}+1_{j>1}+1_{j<m}$. Removed boundary contributions appear on the right, rather than as new nonzeros.

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.