  This paper is concerned with finite element approximations of $W^{2,p}$
strong solutions of second-order linear elliptic partial differential equations
(PDEs) in non-divergence form with continuous coefficients. A nonstandard
(primal) finite element method, which uses finite-dimensional subspaces
consisting globally continuous piecewise polynomial functions, is proposed and
analyzed. The main novelty of the finite element method is to introduce an
interior penalty term, which penalizes the jump of the flux across the interior
element edges/faces, to augment a nonsymmetric piecewise defined and
PDE-induced bilinear form. Existence, uniqueness and error estimate in a
discrete $W^{2,p}$ energy norm are proved for the proposed finite element
method. This is achieved by establishing a discrete Calderon-Zygmund-type
estimate and mimicking strong solution PDE techniques at the discrete level.
Numerical experiments are provided to test the performance of proposed finite
element method and to validate the convergence theory.
