  Gaussian processes (GPs) provide a powerful framework for extrapolation,
interpolation, and noise removal in regression and classification. This paper
considers constraining GPs to arbitrarily-shaped domains with boundary
conditions. We solve a Fourier-like generalised harmonic feature representation
of the GP prior in the domain of interest, which both constrains the GP and
attains a low-rank representation that is used for speeding up inference. The
method scales as $\mathcal{O}(nm^2)$ in prediction and $\mathcal{O}(m^3)$ in
hyperparameter learning for regression, where $n$ is the number of data points
and $m$ the number of features. Furthermore, we make use of the variational
approach to allow the method to deal with non-Gaussian likelihoods. The
experiments cover both simulated and empirical data in which the boundary
conditions allow for inclusion of additional physical information.
