Refactor `expand_dim()` method of the `Grid` class

Currently, when expanding the dimension of instances of the Grid class means expanding the dimension of the multi-index set by appending zero exponent to the new dimensions and followed by one of these two:

• if generating_function is not None, then the generating points can be created on the fly according to the multi-index set.
• if generating_function is None, then the generating points are appended by zero column-wise.

The second outcome is problematic because all zeros as the generating points of the new dimensions are inconsistent with what the array actually means. Each column is one-dimensional interpolation points in that dimension; if there are four points, it means that a cubic one-dimensional polynomial in that dimension can be supported by the grid. If the column values are all zeros then no such polynomials can be supported. Furthermore, appending zeros will violate the condition that the values in each column must be unique to be valid generating points.

Consider the following grid whose multi-index set is:

[
  [0],
  [1],
  [2],
]

with the generating points

array([[ 1.000000e+00],
       [-1.000000e+00],
       [-6.123234e-17]])

After expanding the dimension of the grid to two, the multi-index set becomes

[
  [0, 0],
  [1, 0],
  [2, 0],
]

as expected. The generating points become:

array([[ 1.000000e+00,  0.000000e+00],
       [-1.000000e+00,  0.000000e+00],
       [-6.123234e-17,  0.000000e+00]])

Now if we add an exponent to the grid (currently a valid operation), say, [0, 1], the multi-index set becomes:

array([[0, 0],
       [1, 0],
       [2, 0],
       [0, 1]]),

as expected. But the unisolvent nodes are:

array([[ 1.000000e+00,  0.000000e+00],
       [-1.000000e+00,  0.000000e+00],
       [-6.123234e-17,  0.000000e+00],
       [ 1.000000e+00,  0.000000e+00]])

which are not unisolvent for the two-dimensional polynomial prescribed by the multi-index set because the values of the second dimension are all identical.

I propose removing the possibility of expanding the dimension when there are only generating points defined. This will simplify the code and avoid unwanted outcomes downstream. Besides, without a generating function we have no way of knowing which generating points the new dimensions (post expansion) to use.

This issue is part of the refactoring of the Grid class (Issue #135).