  Elliptic partial differential equations arise in many fields of science and
engineering such as steady state distribution of heat, fluid dynamics,
structural/mechanical engineering, aerospace engineering and seismology etc. In
three dimensions it is well known that the solutions of elliptic boundary value
problems have singular behavior near the corners and edges of the domain. The
singularities which arise are known as vertex, edge, and vertex-edge
singularities. We propose a nonconforming h-p spectral element method to solve
three dimensional elliptic boundary value problems on non-smooth domains to
exponential accuracy. To overcome the singularities which arise in the
neighbourhoods of the vertices, vertex-edges and edges we use local systems of
coordinates. These local coordinates are modified versions of spherical and
cylindrical coordinate systems in their respective neighbourhoods. Away from
these neighbourhoods standard Cartesian coordinates are used. In each of these
neighbourhoods we use a geometrical mesh which becomes finer near the corners
and edges. We then derive differentiability estimates in these new set of
variables and a stability estimate on which our method is based for a
non-conforming h-p spectral element method. The Sobolev spaces in vertex-edge
and edge neighbourhoods are anisotropic and become singular at the corners and
edges.
  The method is essentially a least-squares collocation} method and a solution
can be obtained using Preconditioned Conjugate Gradient Method (PCGM). To solve
the minimization problem we need to solve the normal equations for the
least-squares problem. The residuals in the normal equations can be obtained
without computing and storing mass and stiffness matrices. Computational
results for a number of model problems confirm the theoretical estimates
obtained for the error and computational complexity.
