  The cosmic web consists of a complex configuration of voids, walls,
filaments, and clusters, which formed under the gravitational collapse of
Gaussian fluctuations. Understanding under what conditions these different
structures emerge from simple initial conditions, and how different
cosmological models influence their evolution, is central to the study of the
large-scale structure. Here, we present a general formalism for setting up
initial random density and velocity fields satisfying non-linear constraints
for specialized N-body simulations. These allow us to link the non-linear
conditions on the eigenvalue and eigenvector fields of the deformation tensor,
as specified by caustic skeleton theory, to the current-day cosmic web. By
extending constrained Gaussian random field theory, and the corresponding
Hoffman-Ribak algorithm, to non-linear constraints, we probe the statistical
properties of the progenitors of the walls, filaments, and clusters of the
cosmic web. Applied to cosmological N-body simulations, the proposed techniques
pave the way towards a systematic investigation of the evolution of the
progenitors of the present-day walls, filaments, and clusters, and the embedded
galaxies, putting flesh on the bones of the caustic skeleton. The developed
nonlinear constrained random field theory is valid for generic cosmological
conditions. For ease of visualization, the case study presented here probes the
two-dimensional caustic skeleton.
