  A method for analytical computation of the double surface integrals for all
layer potential kernels associated with the Laplace Green's function, in the
Galerkin boundary element method (BEM) in $\mathbb{R}^3$ using piecewise
constant flat elements is presented. The method uses recursive dimensionality
reduction from 4D ($\mathbb{R}^2\times\mathbb{R}^2$) based on Gauss' divergence
theorem. Computable analytical expressions for all cases of relative location
of the source and receiver triangles are covered for the single and double
layer potentials and their gradients with analytical treatment of the singular
cases are presented. A trick that enables reduction of the case of gradient of
the single layer to the same integrals as for the single layer is introduced
using symmetry properties. The method was confirmed using analytical benchmark
cases, comparisons with error-controlled computations of regular
multidimensional integrals, and a convergence study for singular cases.
