  In this paper we present a method to accurately compute the energy of the
magnetostatic interaction between linearly (or uniformly, as a special case)
magnetized polyhedrons. The method has applications in finite element
micromagnetics, or more generally in computing the magnetostatic interaction
when the magnetization is represented using the finite element method (FEM).
  The magnetostatic energy is described by a six-fold integral that is singular
when the interaction regions overlap, making direct numerical evaluation
problematic. To resolve the singularity, we evaluate four of the six iterated
integrals analytically resulting in a 2d integral over the surface of a
polyhedron, which is nonsingular and can be integrated numerically. This
provides a more accurate and efficient way of computing the magnetostatic
energy integral compared to existing approaches.
  The method was developed to facilitate the evaluation of the demagnetizing
interaction between neighouring elements in finite-element micromagnetics and
provides a possibility to compute the demagnetizing field using efficient fast
multipole or tree code algorithms.
