I don't know if this is an easy question for specialists in the field. Consider 
the following interpolation problem : let $\varepsilon >0$, $X$ be a finite
set of real numbers and $g$ be a real-valued function
on $X$. The goal is to find a function $f$, defined on an interval
containing $X$, that coincides with $g$ on $X$, and admits
a continuous second derivative $f''$ bounded by $\varepsilon$ :
$|f''|< \varepsilon$. 


A necessary condition is that, for any $x< y < z$ in $X$,  

$(\*) |\frac{(y-x)g(z)-(z-x)g(y)+(z-y)g(x)}{(z-x)(z-y)(y-x)}| < \frac{\varepsilon}{2}$. 


Is that condition also sufficient ? If $X$ has just three elements then (\*) says
precisely that the Lagrange interpolating polynomial $f$ satisfies $|f''|< \varepsilon$,
so that equivalence holds in this case.



Interpolation by a function whose second derivative is bounded
It isn't sufficient. Suppose that the elements of $X$ are 
$$x\_1 < x\_2 < \cdots < x\_n,$$
and suppose for simplicity that $\epsilon = 1$. Then the values $f(x\_k)$ carry the same information as the integrals
$$\langle p, f'' \rangle = \int\_{x\_1}^{x\_n} p(x) f''(x) dx,$$
where $p(x)$ is a continuous, piecewise-linear function with $p(x\_1) = p(x\_n) = 0$. You can convert between such an integral and linear combinations of values of $f(x)$ by integrating by parts twice. For any such test function $p$, the inequality
$$|\langle p, f'' \rangle| \le \int\_{x\_1}^{x\_n} |p(x)| dx \qquad \qquad$$
holds. Moreover, many of these inequalities are logically independent. More precisely, the extremal choices of $p(x)$ are those for which $p(x\_k)$ alternate in sign for $i \le k \le j$, and $p(x\_k) = 0$ when $k < i$ or $k > j$. The idea to show that such as $p$ is extremal is to make a dual choice of $f''$ that switches between $1$ and $-1$ when $p$ crosses the $x$ axis. (I apologize for skipping some of the logic in the extremality calculation, but I think that this is correct.) The resulting set of values for $f(\vec{x})$ is not a polytope, but a certain convex region with both flat and curved boundary.


For example, consider your statistic,
$$f\_k = \frac{(x\_{k+1}-x\_k)f(x\_{k+1})-(x\_{k+1}-x\_{k-1})f(x\_k)+(x\_k-x\_{k+1})f(x\_{k-1})}{(x\_{k+1}-x\_k)(x\_{k+1}-x\_{k-1})(x\_k-x\_{k-1})},$$ 
for consecutive triples of points, and suppose that $n=4$. Then the pair $(f\_2,f\_3)$ must lie inside the square $[-\frac12,\frac12]^2$; you identified two other pairs of inequalities, but they do not clip the square further. What actually happens is that two of the corners of the square are rounded by parabolas. The supporting lines of the parabolas come from the test function $p$ that is a piecewise-linear interpolation of $(0,-1,a,0)$ with $a > 0$.


