This paper presents the higher order extension of Particle Discretization Scheme (PDS) and its implementation in FEM framework (PDS-FEM) to solve boundary value problems of linear elastic solids, including brittle cracks. Higher order PDS defines an approximation fd(x) of a function f(x), defined over domain Ω, as the union of local polynomial approximation of f(x) over each Voronoi tessellation elements of Ω. The support of the local polynomial bases being confined to the domain of each Voronoi element, fd(x) consists of discontinuities along each Voronoi boundaries. Considering local polynomial approximations over elements of Delaunay tessellation, PDS define bounded derivatives for this discontinuous fd(x). Utilizing the inherent discontinuities in fd(x), PDS-FEM proposes a numerically efficient treatment for modeling cracks. This novel use of local polynomial approximations in FEM is verified with a set of linear elastic problems, including mode-I crack tip stress field.