  A new class of high-order maximum principle preserving numerical methods is
proposed for solving parabolic equations, with application to the semilinear
Allen--Cahn equation. The proposed method consists of a $k$th-order multistep
exponential integrator in time, and a lumped mass finite element method in
space with piecewise $r$th-order polynomials and Gauss--Lobatto quadrature. At
every time level, the extra values violating the maximum principle are
eliminated at the finite element nodal points by a cut-off operation. The
remaining values at the nodal points satisfy the maximum principle and are
proved to be convergent with an error bound of $O(\tau^k+h^r)$. The accuracy
can be made arbitrarily high-order by choosing large $k$ and $r$. Extensive
numerical results are provided to illustrate the accuracy of the proposed
method and the effectiveness in capturing the pattern of phase-field problems.
