We describe a local Lagrange interpolation method using cubic (i.e. non-tensor product) C 1 splines on cube partitions with five tetrahedra in each cube. We show, by applying a complex proof, that the interpolation method is local, stable, has optimal approximation order and linear complexity. Since no numerical results on trivariate cubic C 1 spline interpolation are known from the literature, the steps of the algorithm, which are different from those of the known methods, are focused on its implementation. In this way, we are able to describe the first implementation of a trivariate C 1 spline interpolation method, run numerical tests and visualize the corresponding isosurfaces. These tests with up to 5.5 × 10 11 data confirm the efficiency of the algorithm.