This paper is concerned with the derivation of a generic sampling technique for a class of non-Gaussian vector-valued random fields. Such an issue typically arises in uncertainty quantification for complex systems, where the input coefficients associated with the elliptic operators must be identified by solving statistical inverse problems. Specifically, we consider the case of non-Gaussian random fields with values in some arbitrary bounded or semibounded subsets of ℝn. The approach involves two main features. The first is the construction of a family of random fields converging, at a user-controlled rate, toward the target random field. Each of these auxialiary random fields can be subsequently simulated by solving a family of Itô stochastic differential equations. The second ingredient is the definition of an adaptive discretization algorithm. The latter allows refining the integration step on-the-fly and prevents the scheme from diverging. The proposed strategy is finally exemplified on three examples, each serving as a benchmark, either for the adaptivity procedure or for the convergence of the diffusions.