You could use
Cubic Hermite Splines to interpolate the values at the vertices. The interpolated points will be linear in the 1st order (connect along the edges), but not in the 2nd order since the surface through each triangle's vertices is underconstrained.
Wikipedia gives a 2d case. For the 3d case you need to compute the weighted sum of 10 control points. The weight is given by the Bernstein function, generalized for a barycentric coordinate zyx:
Babc = (6 * z^a * x^b * y^c) / (a!*b!*c!).
for x+y+z=1 and a+b+c = 3.
So for instance at point A(z=1) so B300 = 6* 1^3 * 0^0 * 0^0 / 6 = 1. So trivially at the control point of vertex A, the value is 1* the value of A.
Generate 6 additional control points, by using the normal at the vertices to compute two slopes along each edge. The normal can be computed from all neighbours of every vertex. These control points (for instance A+dA/3, where dA is the slope at A along the edge to B) are multiplied by Bernstein function B210 and so on, analogous to the 2d case on wiki.
Finally you need a middle control point multiplied by B111. This is the parameter that is underconstrained. It's like a tetrahedron added to the spline, the apex of which is the control point. If it's the average of A,B,C then the tetrahedron will be flat. This control point won't influence the value at the edges, but it will influence the slope, so if you want the result to be 'smooth' you'll have to figure out this control point relative to the neighbouring triangles.