A framework for the construction of stable spectral methods on arbitrary domains with unstructured grids is presented. Although most of the developments are of a general nature, an emphasis is placed on schemes for the solution of partial differential equations defined on the tetrahedron. In the first part the question of well-behaved multivariate polynomial interpolation on the tetrahedron is addressed, and it is shown how to extend the electrostatic analogy of the Jacobi polynomials to problems beyond the line. This allows for the identification of nodal sets suitable for polynomial interpolation within the tetrahedron and, subsequently, for the formulation of stable spectral schemes on such unstructured nodal sets. The second part of this work is devoted to a discussion of weakly imposed boundary conditions, and energy-stable schemes are formulated for a wide class of problems, exemplified by advection problems, advection-diffusion problems, and linear symmetric hyperbolic systems. Finally, in the third part, issues related to computational efficiency and implementation of the schemes are discussed. The spectral accuracy of the approximation is confirmed through an example, and factorization methods for the efficient computation of derivatives on the general nodal sets within the d-simplex are developed, ensuring that the proposed schemes are competitive with tensor-product-based methods. In this last part we also show that the advective operator results in an O (n(-2)) restriction on the time-step, similar to that of spectral collocation methods employing a tensor-product-based approximation. The performance of the proposed scheme is illustrated by solving a wave problem on a triangulated domain, confirming the expected accuracy and stability.