Surface meshes are widely used by many numerical methods for solving partial differential equations. They not only represent computational grids for various discretization methods, but also are numerical objects in themselves. The accuracy of numerical methods, especially high-order methods, are highly dependent on the geometrical accuracy of the mesh as well as on that of differential or integral quantities defined over them. The situation is further complicated if the surface mesh does not have an underlying CAD model representation or is evolving. As a result one generally has access only to a discrete surface for numerical computations. We describe a computational framework based on local polynomial fittings using a weighted least squares approach. Using this framework, we can obtain a piecewise C0 continuous high-order approximations over discrete surface, which subsequently can be used in various mesh-based numerical computation. A least-squares based approach over interpolation allows greater flexibility and stability without any loss of accuracy. This vastly increases the applicability of our methods to a number of applications. We demonstrate the application of our framework to solve problems related to computation of differential quantities, surface integrals, surface optimization, surface adaptation, numerical PDE’s, etc.