Spline interpolation for vector-valued data in 3D space

I have output from a 3D linear elasticity finite element simulation which uses linear tetrahedral elements, such that the displacement is continuous over the nodes but the gradient is not ($C_0$ continuous). I wish to "upgrade" the displacement data so that it can be reported as a $C_1$ continuous quantity which I can directly differentiate to produce the gradient, which would yield a consistency that I find appealing. Instead of directly interpolating gradient data over the nodes (or using a higher order basis function, which of course would make this exercise moot), I wish to interpolate the primary variable, displacement, as a quadratic polynomial with splines.

For 1D quadratic splines, each element or segment takes the form:
$$f_e(x)=a_e+b_ex+c_ex^2$$
$$nabla f_e(x)=b_e+2 c_ex$$
To form the system of equations, the function value and gradient are matched at end points. And since this is 1D, each end point connects to at most one other segment.
$$f_e(x_i)=f_{e+1}(x_i)$$ and $$nabla f_e(x_i)=nabla f_{e+1}(x_i)$$
Such is the predictable nature of this arrangement that you are guaranteed to have $3(n-1)$ variables and at least $n+2(n-1)=3n-2$ equations. This can be supplemented with an additional assumption or solved in a linear least squares fashion. Awesome.

For the 3D case, it looks similar:
$$f_e(x)=a_e+x^Tb_e+x^Tc_ex$$ and $$nabla f_e(x)=b_e+2c_ex$$where $x$ is a $3$ x $1$ vector, $a_e$ is a scalar, $b_e$ is a $3$ x $1$ vector, and $c_e$ is a $3$ x $3$ matrix. Because $x^Tc_ex$ is a quadratic form, only the symmetric part of $c_e$ is of any consequence, so it can be safely assumed $c_e$ is a symmetric matrix with 6 independent variables:

$$f_e(x)=a_e+x^Tb_e+begin{pmatrix}x_0x_0 & x_0x_1 & x_0x_2 & x_1x_1 & x_1x_2 & x_2x_2end{pmatrix}^Tbegin{pmatrix}c_{e00} & c_{e01} & c_{e02} & c_{e11} & c_{e12} & c_{e22}end{pmatrix}$$
$$nabla f_e(x)=b_e+begin{pmatrix}2x_0 & x_1 & x_2 & 0 & 0 & 0 0 & x0 & 0 & 2x_1 & x_2 & 0 0 & 0 & x_0 & 0 & x_1 & 2x_2end{pmatrix}^Tbegin{pmatrix}c_{e00} & c_{e01} & c_{e02} & c_{e11} & c_{e12} & c_{e22}end{pmatrix}$$

The hard part is that each node is not limited to 2 neighboring elements, in some cases it can be upward of 20! In one very modest case (I usually work with much finer grids), I have approximately 2000 tetrahedra and 500 nodes, with 51,000 (!) combinations where I must enforce gradient continuity. For example: if two tetrahedra (T1, T2) share 3 nodes (A, B, C), I must enforce gradient continuity for at A, B and C using the spline variables for T1 and T2.

This yields a sparse matrix ($10^{-4}$ sparsity) with $4t+3*51,000 = 161,000$ equations and $10t=20,000$ variables, and a $161,000$ x $3$ right-hand side (each component of displacement is independently interpolated). This is obviously way more complicated than calculating the original finite element system, which seems counterintuitive to me. Nevertheless, it does work…slowly. In the case above, I solved the original finite element system in 500 milliseconds and the interpolation took 30 seconds. I didn’t even have the memory needed (it blew up to 120 gigs or something silly like that) to interpolate my other grids which takes a whopping 4 seconds to solve per time step. I am using SuiteSparse’s SPQR solver, which produces a least squares solution using a QR factorization of the giant sparse system.

Questions:

1. What do you call this interpolation procedure? Splines? Super-splines? B-Splines? I can’t make heads or tails of what I read in literature on anything beyond 1D splines.

2. Am I "double counting" some conditions and making the system overly complicated? I thought I could ignore a few conditions, such as pairs of tetrahedra that touch at one point only, but ignoring even that destroys the interpolation. I am dumbfounded by how much more complicated this interpolation procedure is than the original finite element calculations.

3. Any ideas on how to simplify this system? I am confident that the algorithms I am using are top notch — I’m just giving them too much work.

Examples plotting pressure (trace of the gradient) without interpolation:

and with the procedure above: