"Mollification" of a convex function with a finite set of points unchanged $f:\mathbb R\to\mathbb R$ is a convex continuous function. We have a finite or a countable set of triples: $\{(x_n,f(x_n),D_n)\}_{n\in N}$, where $D_n$ is the slope of a tangent line $L_n$ at $x_n$ (if at a point $f$ is not differentiable, then multiple lines can be tangents; $L_n$ is just one of those lines).
Assuming that, for any $n,m,k$, the intersection of $L_n$ and $L_m$ cannot be the point $(x_k, f(x_k))$, then we want to prove that there exists a smooth function $g$ such that $g(x_n)=f(x_n)$ and $g'(x_n)=D_n$ for any $n$.

The original problem that I am trying to solve involves multi-dimensional manifolds, but I think it is easy to generalize the 2-dimensional case.
By the mollification theorem, a smooth function approximating $f$ must exist, but can it contain a set of points that corresponds precisely to the points on the graph of $f$?
 A: No, as long as the number of triples you are allowing is infinite. Pick any strictly convex function that is smooth anywhere but at zero (say $f(x) = |x|+|x|^2$). And take as your set of $x_n$ to be the set $\{1/m, m \in \mathbb Z \setminus \{0\}\}$.
Since there's points near zero where the slope is set to be (essentially) $1$ and points where it is set to be (essentially) $-1$ you cannot find a smooth extension.
A: Making quantitative the assumption "the intersection of $L_n$ and $L_m$ cannot be the point $(x_k, f(x_k))$", it is possible to construct such a function $g$ (as pointed out by Jaume, the nonquantitative assumption is not sufficient).
Let us consider the problem in $\mathbb R^n$.
Given a family of indices $I$, let $(x_i)_{i\in I}\subset \mathbb R^n$ be a family of points and let $(L_i)_{i\in I}$ be a family of affine functions. We assume that there is a convex function $f:\mathbb R^n\to\mathbb R$ such that $f(x_i) = L_i(x_i)$ and $f\ge L_i$.
Quantitative assumption: There is an $\varepsilon > 0$ such that for any $i\not=j$ we have $L_i(x)\ge L_j(x)$ for any $x\in B(x_i, \varepsilon)$ (when $x = x_i$ this follows from the convexity of $f$) (this assumption is equivalent to the original one if $I$ is finite).
Let $h=\sup_{i\in I} L_i$. The function $h$ is convex and $h\le f < \infty$.
Let $\rho\in C^{\infty}_c(B(0,\varepsilon))$ be a convolution kernel ($\rho\ge 0$ and $\int\rho=1$).
Define $g=h\star\rho$. It is not hard to check that $g$ is smooth, $g(x_i)=L_i(x_i)$ and $g\ge L_i$ for each $i\in I$.
Relaxing the assumption: It is possible, with an appropriate partition of unity argument, to prove the result even if the value of $\varepsilon$ is allowed to depend on $i$ (so we would have $\varepsilon_i$), provided that it remains locally bounded away from $0$ (i.e., $\inf_{|x_i|<R}\varepsilon_i > 0$).
