The Padua points $\mathrm{Pad}_{n} \subset [-1,1]^{2}$ are a unisolvent pointset with optimal growth of Lebesgue constant, described in detail here. With some work they can be used to generate a multivariate Chebyshev series representation of the function they attempt to interpolate.
I would like to use the Padua points to generate a multivariate Chebyshev series of a function singular along the line $x=y$, but this is impossible because two of the Padua points are consecutive vertices of the square.
Is there a way to generate a multivariate Chebyshev polynomial interpolant of a function $f\colon [-1,1]^2\to \mathbb{R}$ which is singular on the line $x=y$, perhaps using a modification of the Padua points?
Update: The goal is to write a piece of code which projects as many 2D functions defined on rectangular regions as possible onto Chebyshev polynomials and allow for interpolation and recovery of the coefficients for perhaps other uses. This is similar to the goals of chebfun in 2D.