2d-element constraint can have out-of-bounds indices

[IgnaceBleukx]

The current implementation of the Element constraint has a bug:

import cpmpy as cp

iv = cp.intvar(0,10, shape=(3,3))
a,b = cp.intvar(0,10, shape=2)

m = cp.Model()
m += iv[a,b] == 3

this model can have the solution a = 0, b=4, iv[1,1] = 3 which is not what the user intended.
I think the index variables should be constrained to only take value within the bounds of the repsective axis.

[Dimosts]

Again, we face the partial or total function dilemma.

If we want it as partial function, we can just check the indices when we convert it to 1-D, and throw an exception if they are not within the bounds.

If we want it as total function, we have to return constraints forcing the variables to be within the bounds along with the element constraint. This would also work in reified context.

[Wout4]

Since the multidimensional case gets mapped to the one-dimensional, this will also be 'fixed' by #393 where we disallow partial functions

[IgnaceBleukx]

Re-iterating on this issue. I think the correct way to do things is to make a new global constaint "MultiDimElement" or somtehing, and ensure the decomposition is used when a solver does not support multi-dim element.
That allows to safen the expression too such as we do in #515 .