This paper presents a novel method for solving partial differential equations on three-dimensional CAD geometries by means of immersed isogeometric discretizations that do not require quadrature schemes. It relies on a newly developed technique for the evaluation of polynomial integrals over spline boundary representations that is exclusively based on analytical computations. First, through a consistent polynomial approximation step, the finite element operators of the Galerkin method are transformed into integrals involving only polynomial integrands. Then, by successive applications of the divergence theorem, those integrals over B-Reps are transformed into the first surface and then line integrals with polynomials integrands. Eventually, these line integrals are evaluated analytically with machine precision accuracy. The performance of the proposed method is demonstrated by means of numerical experiments in the context of 2D and 3D elliptic problems, retrieving optimal error convergence order in all cases. Finally, the methodology is illustrated for 3D CAD models with an industrial level of complexity.