  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.
