  In this article, quadrature rules for the efficient computation of the
stiffness matrix for the fractional Laplacian in three dimensions are
presented. These rules are based on the Duffy transformation, which is a common
tool for singularity removal. Here, this transformation is adapted to the needs
of the fractional Laplacian in three dimensions. The integrals resulting from
this Duffy transformation are regular integrals over less-dimensional domains.
We present bounds on the number of Gauss points to guarantee error estimates
which are of the same order of magnitude as the finite element error. The
methods presented in this article can easily be adapted to other singular
double integrals in three dimensions with algebraic singularities.
