  We develop a boundary integral equation solver for elliptic partial
differential equations on complex \threed geometries. Our method is efficient,
high-order accurate and robustly handles complex geometries. A key component is
our singular and near-singular layer potential evaluation scheme, \qbkix: a
simple extrapolation of the solution along a line to the boundary. We present a
series of geometry-processing algorithms required for \qbkix to run efficiently
with accuracy guarantees on arbitrary geometries and an adaptive upsampling
scheme based on a iteration-free heuristic for quadrature error. We validate
the accuracy and performance with a series of numerical tests and compare our
approach to a competing local evaluation method.
