I don't quite get why samples from prior in Gaussian Processes look like they do. I build covariance matrix, for example using radial basis kernel or a polynomial one, then I draw samples from multivariate normal Gaussian with mean $0$ and covariance that I calculate and the result is smooth, look like a continuous function. I don't get why. My intuition, that is obviously wrong, would tell me that there will be some discontinuities. Similarly with periodic kernel, I can visualize the covariance function using some heatmap, I can see the periodic pattern, but I don't understand why when I plug that into a multivariate Gaussian it yields periodic functions in return. I would be grateful for any help.
The covariance function encodes prior beliefs about the nature of the function. It basically says how similar the output of the Gaussian process should be as a function of the input features. If you have a covariance function (such as the RBF or polynomial covariance functions) that give a high value if the input features are very similar, then that means you will get a smooth function. The output values for similar input vectors will be similar because that is what the covariance function specifies.