In this paper, we propose a progressive iterative approximation (PIA) algorithm for B-spline surface approximation with normal constraint. On the one hand, the discrete data points of the tangent vector, normal vector, curvature and other geometric characteristics are fully applied to the approximation problem of discrete data points, using the two directions of the tangent vector to construct normal constraint can avoid unnecessary fluctuations, and obtain better approximation effect. On the other hand, the number of selecting feature points is less than the number of data points, so the PIA algorithm can be used for the approximation of the mass of discrete data points. The steps in the process of each iteration of the algorithm are independent, which is easy to be applied to the parallel computing, which greatly improve the computational efficiency. Some examples are given to show the validity of the algorithm.