Generating a divergence-free basis I'm trying to derive an algorithm that would generate a set of divergence-free vectors, which shall be used as basis vectors later on. Using a simple example, a 2D second-order divergence-free basis would look like this:
$$
\begin{bmatrix}
1 \\
0
\end{bmatrix}, 
\begin{bmatrix}
0 \\
1
\end{bmatrix},
\begin{bmatrix}
x \\
-y
\end{bmatrix},
\begin{bmatrix}
y \\
0
\end{bmatrix},
\begin{bmatrix}
0 \\
x
\end{bmatrix},
\begin{bmatrix}
x^2 \\
-2xy
\end{bmatrix},
\begin{bmatrix}
xy \\
-\frac{1}{2} y^2
\end{bmatrix},
\begin{bmatrix}
y^2 \\
0
\end{bmatrix},
\begin{bmatrix}
0 \\
x^2
\end{bmatrix}.
$$
So far, I have manually implemented these calculations, but it will be quite tedious to work it out for 3D cases and for higher order. As of now, my pseudocode for 2D case looks like this:
switch (order):
  order=1:
     // manually calculate vectors for first order
  order=2:
     // manually calculate vectors for second order
  [...]

I would imagine things to be nastier in 3D... Currently, I am trying both to derive the whole basis manually and to come up with an algorithm to do it. Any suggestions?
 A: For dimension $d$, the first $d-1$ entries of your vector $X$ can be arbitrary, and then you want
$$ X_d = \int \left(-\sum_{j=1}^{d-1} \dfrac{\partial X_j}{\partial x_j}\right) \; dx_d + c(x_1, \ldots, x_{d-1})$$
where $c$ is an arbitrary function of $x_1, \ldots, x_{d-1}$.
For a basis of the vector space of divergence-free $d$-dimensional vectors whose entries are polynomials in $x_1, \ldots, x_d$ of total degree $\le p$, consider the cases where one of $X_1, \ldots, X_{d-1}$ is a monomial of total degree $\le p$, the others are $0$, and $X_d$ is the appropriate monomial, or all of $X_1, \ldots, X_{d-1}$ are $0$ and $X_d$ is a monomial that doesn't depend on $x_d$.
For example, the case $d=3$ and order $2$ is
$$\eqalign{&\left[ \matrix{ 1 \cr 0 \cr 0 } \right] \left[ \matrix{ z \cr 0 \cr 0 } \right] \left[ \matrix{ z^2 \cr 0 \cr 0 } \right] \left[ \matrix{ y \cr 0 \cr 0 } \right] \left[ \matrix{ yz \cr 0 \cr 0 } \right] \left[ \matrix{ y^2 \cr 0 \cr 0 } \right] \left[ \matrix{ x \cr 0 \cr -z } \right] \left[ \matrix{ xz \cr 0 \cr -z^2/2 } \right]\cr & \left[ \matrix{ xy \cr 0 \cr -yz } \right] \left[ \matrix{ x^2 \cr 0 \cr -2xz } \right] \left[ \matrix{ 0 \cr 1 \cr 0 } \right] \left[ \matrix{ 0 \cr z \cr 0 } \right] \left[ \matrix{ 0 \cr z^2 \cr 0 } \right] \left[ \matrix{ 0 \cr y \cr -z } \right] \left[ \matrix{ 0 \cr yz \cr -z^2/2 } \right]\cr & \left[ \matrix{ 0 \cr y^2 \cr -2yz } \right] \left[ \matrix{ 0 \cr x \cr 0 } \right] \left[ \matrix{ 0 \cr xz \cr 0 } \right] \left[ \matrix{ 0 \cr xy \cr -xz } \right] \left[ \matrix{ 0 \cr x^2 \cr 0 } \right] \left[ \matrix{ 0 \cr 0 \cr 1 } \right] \left[ \matrix{ 0 \cr 0 \cr x } \right] \left[ \matrix{ 0 \cr 0 \cr x^2 } \right]}
  $$
