Models for the geometrical contact of wheel and rail are a basic component of multibody system (MBS) models for wheel-rail systems. Approximations are used to get sufficiently differentiable contact conditions that can be evaluated efficiently. We discuss an approximation that is essentially based on polynomial 2D-tensorproduct splines and minimizes a functional that combines a weighted least squares approximation of a rigid contact model with a smoothing term. An efficient algorithm to compute the tensorproduct spline is developed, the parallelization on a cluster of workstations is discussed. We report on results of both the sequential and the parallel algorithm and give simulation results for a rigid wheelset on a straight track.