  The finite-size renormalization-group approach for isotropic O$(n)$-symmetric
systems introduced previously [V. Dohm, Phys. Rev. Lett. {\bf 110}, 107207
(2013)] is extended to weakly anisotropic O$(n)$-symmetric systems. Our theory
is formulated within the $\varphi^4$ model with lattice anisotropy in a
$d$-dimensional block geometry with periodic boundary conditions. It describes
the crossover from low- to high-temperature fluctuations including
Goldstone-dominated and critical fluctuations for $1\leq n \leq \infty$ in
$2<d<4$ dimensions. An exact representation is derived for the large-distance
behavior of the bulk correlation function of anisotropic systems in terms of
the principal correlation lengths and an anisotropy matrix ${\bf \bar A}$. This
includes the long-ranged correlations with an anisotropic algebraic decay at
low temperatures due to the Goldstone modes for $n>1$. We calculate the
finite-size scaling functions of the excess free energy and thermodynamic
Casimir force. Exact results are derived in the large-$n$ limit. Applications
are given for $L_\parallel^{d-1} \times L$ slab geometries with a finite aspect
ratio $\rho=L/L_\parallel$ as well as for the film limit $\rho \to 0$. For
weakly anisotropic systems two-scale-factor universality is replaced by
multiparameter universality. This implies a substantial reduction of the
predictive power of bulk and finite-size theory for anisotropic systems as
compared to isotropic systems. The validity of multiparameter universality is
confirmed analytically for a nontrivial example of the $d=2,n=1$ universality
class. Anisotropy-dependent minima of the Casimir force scaling function are
found below $T_c$. Both the sign and magnitude of the Casimir amplitude in the
Goldstone and critical regimes are affected by the lattice anisotropy.
Quantitative predictions are made that can be tested by Monte Carlo
simulations.
