We present a framework for registering and analyzing functional neuroimaging data constrained to the cortical surface of the brain. We assume as input a set of labeled data points that lie on a set of parameterized topologically spherical surfaces that represent the cortical surfaces of multiple subjects. To perform analysis across subjects, we first co-register the coordinates from each surface to a cortical atlas using labeled sulcal maps as constraints. The registration minimizes a thin plate spline energy function on the deforming surface using covariant derivatives to solve the associated PDEs in the intrinsic geometry of the individual surface. The resulting warps are used to bring the functional data for multiple subjects into a common surface atlas coordinate system. We then present a novel method for performing statistical analysis of points on this atlas surface. We use the Green's function of the heat equation on the surface to model probability distributions and thus demonstrate the use of PDEs for statistical analysis in Riemannian manifolds. We describe methods for estimating the mean and variance of a set of points, such that the mean also lies in the manifold. We demonstrate the utility of this framework in the development of a maximum likelihood classifier for parcellation of somatosensory cortex in the atlas based on current dipole fits to MEG data, simulated to represent a somatotopic mapping of S1 sensory areas in multiple subjects.