def _get_basis_product_indices(basis, dimensions, W):
    """
    Returns the indicies for the product between the weights for each
    interpolant for basis functions.
    """
    BPInd = None
    if dimensions == 1:
        return None
    elif dimensions == 2:
        if basis[0][0] == 'T':
            return None
        elif len(basis) == 2:
            if (basis[0][0] == 'L' and basis[1][0] == 'L') or \
               (basis[0][0] == 'L' and basis[1][0] == 'H') or \
               (basis[0][0] == 'H' and basis[1][0] == 'L'):
                BPInd = []
                for ind1 in range(W[1].shape[1]):
                    for ind0 in range(W[0].shape[1]):
                        BPInd.append([ind0, ind1])
            elif basis == ['H3', 'H3']:
                BPInd = [[0, 0], [1, 0], [0, 1], [1, 1],
                      [2, 0], [3, 0], [2, 1], [3, 1],
                      [0, 2], [1, 2], [0, 3], [1, 3],
                      [2, 2], [3, 2], [2, 3], [3, 3]]
            else:
                raise ValueError('Basis combination not supported')
    elif dimensions == 3:
        if len(basis) == 3:
            if (basis[0][0] == 'L' and basis[1][0] == 'L' and basis[2][0] == 'L'):
                BPInd = []
                for ind2 in range(W[2].shape[1]):
                    for ind1 in range(W[1].shape[1]):
                        for ind0 in range(W[0].shape[1]):
                            BPInd.append([ind0, ind1, ind2])
            elif basis == ['H3', 'H3', 'H3']:
                BPInd = [
                    [0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0],
                    [0, 0, 1], [1, 0, 1], [0, 1, 1], [1, 1, 1],
                    [2, 0, 0], [3, 0, 0], [2, 1, 0], [3, 1, 0],
                    [2, 0, 1], [3, 0, 1], [2, 1, 1], [3, 1, 1],
                    [0, 2, 0], [1, 2, 0], [0, 3, 0], [1, 3, 0],
                    [0, 2, 1], [1, 2, 1], [0, 3, 1], [1, 3, 1],
                    [2, 2, 0], [3, 2, 0], [2, 3, 0], [3, 3, 0],
                    [2, 2, 1], [3, 2, 1], [2, 3, 1], [3, 3, 1],
                    
                    [0, 0, 2], [1, 0, 2], [0, 1, 2], [1, 1, 2],
                    [0, 0, 3], [1, 0, 3], [0, 1, 3], [1, 1, 3],
                    [2, 0, 2], [3, 0, 2], [2, 1, 2], [3, 1, 2],
                    [2, 0, 3], [3, 0, 3], [2, 1, 3], [3, 1, 3],
                    [0, 2, 2], [1, 2, 2], [0, 3, 2], [1, 3, 2],
                    [0, 2, 3], [1, 2, 3], [0, 3, 3], [1, 3, 3],
                    [2, 2, 2], [3, 2, 2], [2, 3, 2], [3, 3, 2],
                    [2, 2, 3], [3, 2, 3], [2, 3, 3], [3, 3, 3]]
            else:
                raise ValueError('Basis combination not supported')
        else:
            raise ValueError('Basis combination not supported')
    else:
        raise ValueError('%d dimensions not supported' % (len(basis)))
        
    return BPInd