  We analyze theoretically and experimentally vortex configurations in
mesoscopic superconducting squares. Our theoretical approach is based on the
analytical solution of the London equation using Green's-function method. The
potential-energy landscape found for each vortex configuration is then used in
Langevin-type molecular-dynamics simulations to obtain stable vortex
configurations. Metastable states and transitions between them and the ground
state are analyzed. We present our results of the first direct visualization of
vortex patterns in micrometer-sized Nb squares, using the Bitter decoration
technique. We show that the filling rules for vortices in squares with
increasing applied magnetic field can be formulated, although in a different
manner than in disks, in terms of formation of vortex "shells".
