  We develop a new $C^{\,0}$-continuous Petrov-Galerkin spectral element method
for one-dimensional fractional elliptic problems of the form
${}_{0}{\mathcal{D}}_{x}^{\alpha} u(x) - \lambda u(x) = f(x)$, $\alpha \in
(1,2]$, subject to homogeneous boundary conditions. We employ the standard
(modal) spectral element bases and the Jacobi poly-fractonomials as the test
functions [1]. We formulate a new procedure for assembling the global linear
system from elemental (local) mass and stiffness matrices. The Petrov-Galerkin
formulation requires performing elemental (local) construction of mass and
stiffness matrices in the standard domain only once. Moreover, we efficiently
obtain the non-local (history) stiffness matrices, in which the non-locality is
presented analytically for uniform grids. We also investigate two distinct
choices of basis/test functions: i) local basis/test functions, and ii) local
basis with global test functions. We show that the former choice leads to a
better-conditioned system and accuracy. We consider smooth and singular
solutions, where the singularity can occur at boundary points as well as in the
interior domain. We also construct two non-uniform grids over the whole
computational domain in order to capture singular solutions. Finally, we
perform a systematic numerical study of non-local effects via full and partial
history fading in order to further enhance the efficiency of the scheme.
