  In this work we present an adaptive boundary element method for computing the
electromagnetic response of wave interactions in hyperbolic metamaterials. One
unique feature of hyperbolic metamaterial is the strongly directional wave in
its propagating cone, which induces sharp transition for the solution of the
integral equation across the cone boundary when wave starts to decay or grow
exponentially. In order to avoid a global refined mesh over the whole boundary,
we employ a two-level a posteriori error estimator and an adaptive mesh
refinement procedure to resolve the singularity locally for the solution of the
integral equation. Such an adaptive procedure allows for the reduction of the
degree of freedom significantly for the integral equation solver while
achieving desired accuracy for the solution. In addition, to resolve the fast
transition of the fundamental solution and its derivatives accurately across
the propagation cone boundary, adaptive numerical quadrature rules are applied
to evaluate the integrals for the stiff matrices. Finally, in order to
formulate the integral equations over the boundary, we also derive the limits
of layer potentials and their derivatives in the hyperbolic media when the
target points approach the boundary.
