We fully generalize a previously-developed computational geometry tool [1] to perform large-scale simulations of arbitrary two-dimensional faceted surfaces $z = h(x,y)$. Our method uses a three-component facet/edge/junction storage model, which by naturally mirroring the intrinsic surface structure allows both rapid simulation and easy extraction of geometrical statistics. The bulk of this paper is a comprehensive treatment of topological events, which are detected and performed explicitly. In addition, we also give a careful analysis of the subtle pitfalls associated with time-stepping schemes for systems with topological changes. The method is demonstrated using a simple facet dynamics on surfaces with three different symmetries. Appendices detail the reconnection of holes left by facet removal and a strategy for dealing with the inherent kinematic non-uniqueness displayed by several topological events. [1] S.A. Norris and S.J. Watson, Acta Mat. 55 (2007) p. 6444