Wavelets for the discretization of boundary integral operators usually have fixed order and are constructed in some parameter space of the surface. Here a new approach is presented, where the order is flexible and no parameterizations are needed. The wavelets are restrictions of piecewise polynomial functions in three variables on the boundary manifold. This construction is especially suited for surfaces with complicated geometries. If the polynomial order is suitably adjusted to the level of the wavelet, then the non-standard form of a large class of boundary integral operators can be truncated to contain only O(N) non-vanishing entries while retaining the asymptotic convergence of the full Galerkin scheme. An algorithm which sets up the basis and the non-standard form in O(N) operations will be discussed. The method is applied to problems from potential theory and Stokes flow and compared with the Fast Multipole Method.