This paper is devoted to the study of non-conforming h-p spectral element methods for three dimensional elliptic problems on non-smooth domains using parallel computers. To overcome the singularities which arise in the neighbourhoods of vertices, edges and vertex-edges we use local systems of coordinates together with geometric meshes which become finer near corners and edges. We provide estimates for the error and the computational complexity of our method and present results of numerical simulations that have been performed to validate the theoretical estimates obtained in [22]. © 2014 Elsevier Inc. All rights reserved.