To implement any of the aforementioned stabilization forms, one must be able to find face neighbors of given elements and assemble contributions from pairs of elements sharing a face. For higher-order elements, the face-based ghost penalty requires computation of higher-order derivatives. The element-based formulation may then be preferable. All the stabilization terms lead to additional coupling and, as a consequence, fill-in in the stiffness matrix [1]}, [2]}. Systems stabilized by the ghost penalty have similar condition number bounds as boundary-fitted formulations, and therefore larger systems (i.e., finer grids) yield larger condition numbers and require more iterations in iterative solvers. This can be overcome by combining the ghost penalty with a multigrid solver; see the discussion on this in Section REF .
