  In our previous work [Phys. Rev. B 101,245134 (2020)], we constructed a
series of higher-dimensional exactly solvable fracton ordered models with
spatially extensive excitations whose mobility and deformability are
restricted. In this paper, we proceed further to propose a combinatorial method
to calculate their ground state degeneracy (GSD). By decomposing a state into a
collection of subsystem ground state sectors, we find that a ground state can
be reconstructed from a consistent collection of subsystem data. As a result,
we find that GSD formulas (after taking logarithm) have very diverse
polynomial-dependency on the sizes. Inspired from known results about the
relation between GSD of X-cube models and the first Betti numbers with
$\mathbb{Z}_2$ coefficients of leaves, we expect that the coefficients in the
polynomials of these higher-dimensional models may reflect some mathematical
properties in higher dimensions.
