In this paper we study the cluster formation and evolution in spin crossover systems during the thermal transition in the frame of a mechano-elastic model applied to open boundary hexagonal lattices. The switching processes between the high-spin (HS) and low-spin (LS) state are studied by a method combining a Monte Carlo standard procedure on the spin state and the lattice relaxation. In the present study, we adopt the transition probabilities of the spin state taking into account the energy gap between the two states, the degeneracy ratio and the local pressure determined by the elongations of the closest springs. It is found that clusters of molecules in the same state tend to grow starting from corners, as in available experimental data. Some qualitative differences between the processes of cluster formation for the two hysteresis branches, i.e., HS to LS and LS to HS are pointed out. Moreover, we have studied the dependence of cluster formation on the strength of the elastic interactions, and also on the system size. The size dependence of the ratio between the system size and the maximum cluster length is very weak, which indicates the appearance of macroscopic domains.