We present the implementation of a fast estimator for the full dark matter bispectrum of a three-dimensional particle distribution relying on a separable modal expansion of the bispectrum. The computational cost of accurate bispectrum estimation is negligible relative to simulation evolution, so the isotropic bispectrum can be used as a standard diagnostic whenever the power spectrum is evaluated. As an application we measure the evolution of gravitational and primordial dark matter bispectra in $N$-body simulations with Gaussian and non-Gaussian initial conditions of the local, equilateral, orthogonal and flattened shape. The results are compared to theoretical models using a 3D visualisation, 3D shape correlations and the cumulative bispectrum signal-to-noise, all of which can be evaluated extremely quickly. Our measured bispectra are determined by $mathcal{O}(50)$ coefficients, which can be used as fitting formulae in the nonlinear regime and for non-Gaussian initial conditions. In the nonlinear regime with $k<2h,mathrm{Mpc}^{-1}$, we find an excellent correlation between the measured dark matter bispectrum and a simple model based on a `constant bispectrum plus a (nonlinear) tree-level gravitational bispectrum. In the same range for non-Gaussian simulations, we find an excellent correlation between the measured additional bispectrum and a constant model plus a (nonlinear) tree-level primordial bispectrum. We demonstrate that the constant contribution to the non-Gaussian bispectrum can be understood as a time-shift of the constant mode in the gravitational bispectrum, which is motivated by the one-halo model. The final amplitude of this extra non-Gaussian constant contribution is directly related to the initial amplitude of the constant mode in the primordial bispectrum. We also comment on the effects of regular grid and glass initial conditions on the bispectrum.