 Solution

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.