The present book develops the mathematical and numerical analysis of linear, elliptic and parabolic partial differential equations (PDEs) with coefficients whose logarithms are modelled as Gaussian random fields (GRFs), in polygonal and polyhedral physical domains. Both, forward and Bayesian inverse PDE problems subject to GRF priors are considered. Adopting a pathwise, affine-parametric representation of the GRFs, turns the random PDEs into equivalent, countably-parametric, deterministic PDEs, with nonuniform ellipticity constants. A detailed sparsity analysis of Wiener-Hermite polynomial chaos expansions of the corresponding parametric PDE solution families by analytic continuation into the complex domain  is developed, in corner- and edge-weighted function spaces on the physical domain. The presented Algorithms and results are relevant for the mathematical analysis of many approximation methods for PDEs with GRF inputs, such as model order reduction, neural network and tensor-formatted surrogates of parametric solution families. They are expected to impact computational uncertainty quantification subject to GRF models of uncertainty in PDEs, and are of interest for researchers and graduate students in both, applied and computational mathematics, as well as in computational science and engineering.