  We design and numerically validate a recovery based linear finite element
method for solving the biharmonic equation. The main idea is to replace the
gradient operator $\nabla$ on linear finite element space by $G(\nabla)$ in the
weak formulation of the biharmonic equation, where $G$ is the recovery operator
which recovers the piecewise constant function into the linear finite element
space. By operator $G$, Laplace operator $\Delta$ is replaced by $\nabla\cdot
G(\nabla)$. Furthermore the boundary condition on normal derivative $\nabla
u\cdot \pmb{n}$ is treated by the boundary penalty method. The explicit matrix
expression of the proposed method is also introduced. Numerical examples on
uniform and adaptive meshes are presented to illustrate the correctness and
effectiveness of the proposed method.
