While the finite element isogeometric analysis (FE IGA) has received most of the attention of the scientific community, only few results are available for boundary element methods (BEM) and boundary integral equations (BIE) in general. In 2D, there are some works on potential flows, such as [36,44], and on elastostatic analysis [40]. Three dimensional appli cations are of great interest in the maritime community, and there are some works on NURBS based panel methods to study, for example, marine propellers [25] or the wavemaking resistance problem [10]. Three dimensional IGA BEM linear elasticity has recently been studied in [20] while some work on shape optimization has been presented in [29]. Boundary element isogeometric analysis (IGA BEM) is very attractive for the solution of homogeneous elliptic PDEs, both in the interior and in the exterior cases, since it requires the solution of integral equations only on the boundary of the enclos ing (or enclosed) domain, which is typically the only information that is generated by standard CAGD tools. In contrast, FE IGA requires the extension of the computational domain inside (or outside) the enclosing CAGD surface, an outstand ing problem in IGA.While the dimensional reduction of the domain diminishes significantly both the number of degrees of freedom and the computational cost associated with the handling of the geometry, there are two main drawbacks which require special attention, and are common to all boundary integral discretizations: (i) the system matrices that are produced with a bound ary integral method are full, and (ii) the integrands contain ð1=rÞ singularities or worse (in three dimensions) around the source point, arising from the fundamental solutions in formulating the boundary integrals. Both issues have perhaps con tributed in keeping the IGA community away from boundary integral formulations.In the literature, the first issue is usually tackled by acceleration techniques such as the fast multipole method [39] (see [31] for a survey of this technique applied to BEM), which allow one to solve the system without assembling the full matri ces, and at the same time reduce the cost of the matrix vector multiplication from Oðn 2 Þ to OðnÞ (this technique has already been applied to Laplace problems in two dimensional IGA BEM in [44]). The second issue is more delicate and its solution depends highly on the degree of the discrete functional spaces and of the geometry description, remaining one of the most difficult aspects of the implementation of BEM. Different approaches are possible, (see, for example, [23]) but high order methods and curved geometries still remain challenging.Popular approaches to tackle singular integration involve a local change of variables for the quadrature rules, following the principle that the Jacobian of the transformation can be used to eliminate locally the kernel singularities. Examples of these methods are given by Lachat and Watson [28] and by Telles [45]. Both methods are elegant and produce accurate re sults, but t...