Abstract. A trivariate Lagrange interpolation method based on C 1 cubic splines is described. The splines are defined over a special refinement of the Freudenthal partition of a cube partition. The interpolating splines are uniquely determined by data values, but no derivatives are needed. The interpolation method is local and stable, provides optimal order approximation, and has linear complexity. §1. Introductionbe a set of points in R 3 . In this paper we are interested in the following problem. Problem 1. Find a tetrahedral partition whose set of vertices includes V, an Ndimensional space S of trivariate splines defined on with a prescribed smoothness r, and a set of additional points {η i } N i=n+1 such that for every choice of the data, there is a unique spline s ∈ S satisfying (1.1)We call P :and S a Lagrange interpolation pair. It is easy to solve this problem using C 0 splines; see Remark 1. In this case no additional interpolation points are needed. However, the situation is much more complicated if we want to use C r splines with r ≥ 1. For r = 1 the problem was first solved in [29] for the special case where the points in V are the vertices of a cube partition. The method uses quintic splines, and provides full approximation power six for sufficiently smooth functions. In [23] we recently described two C 1 methods which solve the problem for arbitrary initial point sets V. The first method is based on quadratic splines, and requires that each tetrahedron be split into 24 subtetrahedra. The second method is based on cubic splines, and requires that each tetrahedron be split into 12 subtetrahedra. Both methods provide approximation order three for sufficiently smooth functions. For the quadratic spline case this is the optimal approximation order, but for the cubic case it is suboptimal.The aim of this paper is to describe a new method based on cubic C 1 splines which yields optimal approximation order four. To achieve this, we restrict ourselves