  Photonic graphene, a photonic crystal with honeycomb structures, has been
intensively studied in both theoretical and applied fields. Similar to graphene
which admits Dirac Fermions and topological edge states, photonic graphene
supports novel and subtle propagating modes (edge modes) of electromagnetic
waves. These modes have wide applications in many optical systems. In this
paper, we propose a novel gradient recovery method based on Bloch theory for
the computation of topological edge modes in the honeycomb structure. Compared
to standard finite element methods, this method provides higher order accuracy
with the help of gradient recovery technique. This high order accuracy is
highly desired for constructing the propagating electromagnetic modes in
applications. We analyze the accuracy and prove the superconvergence of this
method. Numerical examples are presented to show the efficiency by computing
the edge mode for the $\mathcal{P}$-symmetry and $\mathcal{C}$-symmetry
breaking cases in honeycomb structures.
