We derive the exact boundary, in any direction, of functional (semi-infinite) multivariate inequality constraints on the parameter space of a regression model. A parameter space subject to multiple constraints is also possible. The calculated boundaries inform the use of an MCMC hit-and-run algorithm used to determine the posterior distribution of the constrained regression models. The method is applied to a forensic morphometric analysis of human skulls where monotonicity is required in more than one dimension which, to our knowledge, alternate methods are yet to achieve. Methods are compared to those currently available with one dimensional constraints and the results discussed.