Motivated by the magneto hydrodynamic (MHD) simulation for Tokamaks with Isogeometric analysis, we present a new type of splines defined over a rectangular mesh with arbitrary topology, which are piecewise polynomial functions of bidegree (d,d) and C^r parameter continuity. In particular, We compute their dimension and exhibit basis functions called Hermite bases for bicubic spline spaces. We investigate their potential applications for solving partial differential equations (PDEs) over a complex physical domain in the framework of Isogeometric analysis. In particular, we analyze the property of approximation of these spline spaces for the L2-norm. Despite the fact that the basis functions are singular at extraordinary vertices, we show that the optimal approximation order and numerical convergence rates are reached by setting a proper parameterization.