The Earth's tectonic plates are strong, viscoelastic shells which make up the outermost part of a thermally convecting, predominantly viscous layer. In order to build a more realistic simulation of the planet's evolution, the complete viscoelastic convection system ought to be included. A particle-in-cell finite element method is demonstrated which can simulate very large deformation viscoelasticity. This is applied to a plate-deformation problem. Numerical accuracy is demonstrated relative to analytic benchmarks, and the characteristics of the method are discussed.