The analysis of the joint cumulative distribution function (CDF) with bivariate event time data is a challenging problem both theoretically and numerically. This paper develops a tensor spline-based sieve maximum likelihood estimation method to estimate the joint CDF with bivariate current status data. The I-splines are used to approximate the joint CDF in order to simplify the numerical computation of a constrained maximum likelihood estimation problem. The generalized gradient projection algorithm is used to compute the constrained optimization problem. Based on the properties of B-spline basis functions it is shown that the proposed tensor spline-based nonparametric sieve maximum likelihood estimator is consistent with a rate of convergence potentially better than $n^{1/3}$ under some mild regularity conditions. The simulation studies with moderate sample sizes are carried out to demonstrate that the finite sample performance of the proposed estimator is generally satisfactory.