It is well known that using the conventional non-Hermitian but ${\cal PT}-$symmetric Bose-Hubbard Hamiltonian with real spectrum one can realize the Bose-Einstein condensation (BEC) process in an exceptional-point limit of order $N$. Such an exactly solvable simulation of the BEC-type phase transition is, unfortunately, incomplete because the standard version of the model only offers an extreme form of the limit characterized by a minimal geometric multiplicity $K=1$. In our paper we describe a rescaled and partitioned direct-sum modification of the linear version of the Bose-Hubbard model which remains exactly solvable while admitting any value of $K\geq 1$. It offers a complete menu of benchmark models numbered by a specific combinatorial scheme. In this manner, an exhaustive classification of the general BEC patterns with any geometric multiplicity is obtained and realized in terms of an exactly solvable generalized Bose-Hubbard model.