  In this paper, we demonstrate that many of the computational tools for
univariate orthogonal polynomials have analogues for a family of bivariate
orthogonal polynomials on the triangle, including Clenshaw's algorithm and
sparse differentiation operators. This allows us to derive a practical spectral
method for solving linear partial differential equations on triangles with
sparse discretizations. We can thereby rapidly solve partial differential
equations using polynomials with degrees in the thousands, resulting in sparse
discretizations with as many as several million degrees of freedom.
