Gaussian process regression techniques have been used in fluid mechanics for the reconstruction of flow fields from a reduction-of-dimension perspective. A main ingredient in this setting is the construction of adapted covariance functions, or kernels, to obtain such estimates. In this paper, we present a general method for constraining a prescribed Gaussian process on an arbitrary compact set. The kernel of the pre-defined process must be at least continuous and may include other information about the studied phenomenon. This general boundary-constraining framework can be implemented with high flexibility for a broad range of engineering applications. From this, we derive physics-informed kernels for simulating two-dimensional velocity fields of an incompressible (divergence-free) flow around aerodynamic profiles. These kernels allow to define Gaussian process priors satisfying the incompressibility condition and the prescribed boundary conditions along the profile in a continuous manner. We describe an adapted numerical method for the boundary-constraining procedure parameterized by a measure on the compact set. The relevance of the methodology and performances are illustrated by numerical simulations of flows around a cylinder and a NACA 0412 airfoil profile, for which no observation at the boundary is needed at all.