We propose an efficient and flexible scheme to fairly interpolate or approximate the vertices of a given triangular mesh. Instead of generating a piecewise polynomial representation, our output will be a refined mesh with vertices lying densely on a surface with minimum bending energy. To obtain those, we generalize the finite differences technique to parametric meshes. The use of local parameterizations (charts) makes it possible to cast the minimization of non-linear geometric functionals into solving a sparse linear system. Efficient multi-grid solvers can be applied which leads to fast algorithms that generate surfaces of high quality.