Introducing inequality constraints in Gaussian processes can lead to more realistic uncertainties in learning a great variety of real-world problems. We consider the finite-dimensional Gaussian model from Maatouk and Bay [Math. Geosci., 49 (2017), pp. 557--582] which can satisfy inequality conditions everywhere (either boundedness, monotonicity, or convexity). Our contributions are threefold. First, we extend their approach in order to deal with sets of linear inequalities. Second, we explore different Markov chain Monte Carlo (MCMC) methods to approximate the posterior distribution. Third, we investigate theoretical and numerical properties of a constrained likelihood for covariance parameter estimation. According to experiments on both artificial and real data, our framework together with a Hamiltonian Monte Carlo sampler provides efficient results on both data fitting and uncertainty quantification.