Question on Derivative and Conversion on Chebyshev^2

I am somewhat puzzled about what exactly is happening when computing derivatives and conversions on 2d Chebyshev product spaces. Some enlightenment would be greatly appreciated. This is also related to #725.

My question circles around the fact that Chebyshev(0..1)^2 is defined on Padua points, whereas derivatives and projections thereof, living in Ultraspherical(1, 0..1) ⊗ Chebyshev(0..1), Chebyshev(0..1) ⊗ Ultraspherical(1, 0..1) or Ultraspherical(1, 0..1)^2, are defined on tensor product meshes of Chebyshev points.
I am wondering how exactly this works and how this is dealt with internally. I found that there are special methods for points(), transform() and itransform(), dispatching on TensorSpace{<:Tuple{<:Chebyshev{<:ChebyshevInterval},<:Chebyshev{<:ChebyshevInterval}}}, but I could not find similar methods for Derivative() or Conversion(), taking care of the special nature of the Padua points.

For context: I have a problem (namely computing Poincaré integral invariants), where

• I construct a Chebyshev product space,

C2 = Chebyshev(0..1)^2

• extract the Padua points and evaluate some scalar function on those points,

p = points(C2)
v = f.(p)

• transform those to coefficients and create a fun to work with, e.g., to compute derivatives,

plan = plan_transform(C2, v)
c = plan * v
f = Fun(C2, c)
# ...

In the end I need to extract all values, derivatives, etc., on a common grid, which is why everything that does not already live there is projected onto Ultraspherical(1, 0..1)^2.

I implemented this using both, ApproxFun and OrthogonalPolynomialsQuasi, where the former uses Padua points on the Chebyshev()^2 space and the latter uses a Chebyshev tensor grid. Not surprisingly, the ApproxFun implementation gives better results for a comparable number of points. However, even when using the full grid corresponding to the Padua grid, I cannot get comparable results in the OrthogonalPolynomialsQuasi implementation (accuracy is much worse). Trying to understand the reason for this and debugging my code lead to the above questions.

This also raises the question if using Padua points is possible together with OrthogonalPolynomialsQuasi (although this might not be the right place to ask this question)?

Start by reading the existing points(), transform(), and itransform() dispatches for the Chebyshev tensor space, then trace the Derivative() and Conversion() paths mentioned in the issue. Compare the documented behavior with the C2, plan_transform(), and Fun examples; done would be a clear explanation of the Padua-point handling and whether any documentation or implementation change is needed.