We prove convergence of renormalized correlations of primary fields, i. e., spins, disorders, fermions and energy densities, in the scaling limit of the critical Ising model in arbitrary finitely connected domains, with fixed (plus or minus) or free boundary conditions, or mixture thereof. We describe the limits of correlations in terms of solutions of Riemann boundary value problems, and prove their conformal covariance. Moreover, we prove fusion rules, or operator product expansions, which describe asymptotics of the scaling limits of the correlations as some of the points collide together. We give explicit formulae for correlations in the case of simply-connected and doubly-connected domains. Our presentation is self-contained, and the proofs are simplified as compared to the previous work where particular cases are treated.