It is well known that, using the conventional non-Hermitian but 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≥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.