Summary: Motivated by models for multiway comparison data, we consider the problem of estimating a coordinatewise isotonic function on the domain \([0,1]^d\) from noisy observations collected on a uniform lattice, but where the design points have been permuted along each dimension. While the univariate and bivariate versions of this problem have received significant attention, our focus is on the multivariate case \(d\ge 3\). We study both the minimax risk of estimation (in empirical \(L_2\) loss) and the fundamental limits of adaptation (quantified by the adaptivity index) to a family of piecewise constant functions. We provide a computationally efficient Mirsky partition estimator that is minimax optimal while also achieving the smallest adaptivity index possible for polynomial time procedures. Thus, from a worst-case perspective and in sharp contrast to the bivariate case, the latent permutations in the model do not introduce significant computational difficulties over and above vanilla isotonic regression. On the other hand, the fundamental limits of adaptation are significantly different with and without unknown permutations: Assuming a hardness conjecture from average-case complexity theory, a statistical-computational gap manifests in the former case. In a complementary direction, we show that natural modifications of existing estimators fail to satisfy at least one of the desiderata of optimal worst-case statistical performance, computational efficiency and fast adaptation. Along the way to showing our results, we improve adaptation results in the special case \(d=2\) and establish some properties of estimators for vanilla isotonic regression, both of which may be of independent interest.