PS:
For example, like BDM element in a rectangle, one adds two vector basis ( -x**{k+1}, (k+1)*x**k*y ) and ( (k+1)*x*y**k, -y**{k+1} ) to the standard Pk*Pk.
Now, we want to add ( -x**{k+1}, (k+1)*x**k*y ) and ( (k+1)*x*y**k, -y**{k+1} ) to Qk*Qk element-by-element in a quadrilateral mesh. How can we finish it ? I know it is a tricky work.