We describe a unified computational framework for polynomial fitting method based on weighted least squares approximations. This framework is first motivated by computing normals and curvatures on discrete surface mesh, then it is extended to construct the continuous global surface and calculate high order surface integration. Finally, the idea of generalized finite difference method is extracted by inverting the Vandermonde matrix from local polynomial fitting, it will be used to solve various geometric geometric partial differential equations and an elastic membrane problem. Different applications require different techniques and impose different challenges including simplicity, accuracy, continuity, robustness and efficiency. Our research focus on high order accuracy, but also give considerations to other requirements. For normal and curvature calculations, surface reconstruction and integration, we demonstrate order of accuracy up to 6 . For numerical solution of geometric PDEs and the elastic membrane problem, we introduce an accurate spatial discretization over triangular surface mesh and our semi-implicit schemes can achieve at least quadratic convergence rate while being much more accurate and stable than using explicit schemes.