I'm trying to produce a set of permutations of 3 axis Cartesian coordinates which are limited by a specific geometric constraint:
Perhaps the simplest way to visualize this is with a Sudoku-like game where the player must fill an n square grid (where n>1) with all integers up to n so that when we add up the integers on each row or on each column the result is n!. (the reason I call these sets 'cubes' is explained better, with images below.)
eg: for a 2 cube:
0 1
1 0
And
1 0
0 1
and one of the n = 3 cubes:
1 2 0
2 0 1
0 1 2
an n = 4 cube:
3 2 1 0
2 3 0 1
1 0 3 2
0 1 2 3
or
3 1 2 0
1 0 3 2
2 3 0 1
0 2 1 3
or
3 1 0 2
0 3 2 1
2 0 1 3
1 2 3 0
A more visual explanation using vectors (x,y,z) as the centers of unit cubes:
lets say that in the above Sudoku-like examples, the coordinates of the grids' cells define the x and y components and the integers within each cell, the z component of a set or array of vectors.
let's start with a unit cube and label it with coordinates (x,y,z).Next we 'cube' the unit cube (0,0,0) so that we end up with 8 cubes labeled: (0,0,0); (1,0,0); (0,1,0); (1,1,0); (0,0,1); (1,0,1); (0,1,1); (1,1,1).
Now from we remove enough cubes so that each coordinate axis contains only one unit cube.
eg: in this case there are only two possibilities: (0,0,0); (1,1,0); (0,1,1);(1,0,1) and a rotation of Pi/2 radians about any axis, say the Z axis (1,0,1); (0,1,0); (0,0,1); (1,1,1).

In the case of a 3 cube we also get more than one possibility:(1,0,0); (2,1,0); (0,2,0); (0,0,1); (1,1,1); (2,2,1); (2,0,2); (0,1,2); (1,2,2); (2,0,2) and like before a PI/2 radians rotation around any axis:
A 4 cubed version can be constructed by simply replacing each cube of the two cube version with a scaled down version of the whole:

Likewise, a 6 cubed version can be constructed by replacing each cube of the two cubed version with a scaled down 3 cube.

Here and in the previous case I use a fractional 'scale down' for brevity however, any scaling will introduce fractional coordinates which will unnecessarily complicate things. What should actually happen then is after scaling down, the whole should be scaled up the same amount to restore coordinates to integer values.
To sum up: I tried using Python 'itertools' library, but I haven't been successful yet. The sets of coordinates that I get all have more than n^2 vectors. My coding skills are rather rudimentary and limited to mostly trial and error. I deduced the examples above on paper using the Sudoku-like approach, as well as a parametric 3d modeling program to manipulate cubes in 3D space. I initially suspected that what I'm after is an already defined kind of permutation but again, I've not been able to find it.
If it isn't obvious above, each set of cubes I show look exactly the same when viewed along each Cartesian axis and there's never a cube in front of any other.