Legendre transformation: non-convex/non-convave functions
As long as long-range interactions (e.g., gravity) are negligible, the internal energy $U$ of a homogeneous system is an extensive quantity. In such a case, the second law of thermodynamics in fact guarantees that $U$ is a convex function of other extensive quantities it depends on, e.g., $S$, $V$, and $N$.
An important observation is that given a thermodynamic system, a tiny part of it should be in equilibrium with the rest. Suppose that we take this tiny part to be negligibly small compared to the entire system, and simply call it "the system." The rest of the entire system acts as a reservoir with some fixed temperature $T_{0}$, pressure $P_{0}$, and chemical potential $\mu_{0}$, so let's call this "the reservoir." When the system and the reservoir are in thermodynamic equilibrium, the second law of thermodynamics demands that the total entropy be maximized. Notice that the change in the total entropy is given by \begin{equation} \Delta S_{\mathrm{total}} = \Delta S + \frac{Q_{\mathrm{reservoir}}}{T_{0}} = \Delta S - \frac{\Delta U + P_{0}\Delta V - \mu_{0}\Delta N}{T_{0}} = -\frac{1}{T_{0}} \Delta \Omega, \end{equation} where $\Omega = U-T_{0}S + P_{0}V-\mu_{0}N$. Hence, maximizing $S_{\mathrm{total}}$ is equivalent to minimizing $\Omega$. Here, $\Omega$ has three, as opposed to four, independent variables because there is a relation between $U$, $S$, $V$, and $N$. Taking the independent parameters to be $S$, $V$, and $N$, we have \begin{equation} \Omega(S,V,N) = U(S,V,N) - T_{0}S + P_{0} V - \mu_{0} N. \end{equation}
A consistency condition is that for an arbitrary $T_{0}$, $P_{0}$, and $\mu_{0}$, there exists a unique minimum of $\Omega(S,V,N)$ with respect to its independent variables, viz., the system should be able to equilibrate with the reservoir. It then follows that the function $U(S,V,N)$ is convex.
So far, $S$, $V$, $N$, and $U(S,V,N)$ are properties of a small system that is a part of a much larger body. However, for a homogeneous system, the internal energy of the entire body must have the same form as any small part of it. Therefore, what I have shown equally holds for an arbitrary homogeneous system.
[Taking the independent variables of $\Omega$ to be $U$, $V$, and $N$ would lead to the concavity of the function $S(U,V,N)$.]
The procedure $f(x)\rightarrow g(u)=xu-f(x)$ obtained by inverting $u=f'(x)$ is the "wrong" definition of the Legendre transform and that's what gets you in trouble when the function is not convex or twice differentiable. The correct definition in the one given by Werner Fenchel, namely, $$ g(u)=\sup_{x}\ (\ xu-f(x)\ )\ . $$ It is always convex and well defined even if $f$ is not convex and it stays the same if one replaces $f$ by its convex envelope. You can learn all about it in Yvan's book (I don't know why he didn't mention it in his comment above), pp. 480-488.
If you want more thorough references on convexity in statistical mechanics, the canonical ones are:
"Convexity in the Theory of Lattice Gases" by Robert Israel. The introduction by Arthur Wightman is a real beauty and in particular it discusses the argument due to Gibbs about how the convexity of $U$ as a function of $S$ essentially follows from the Second Law of Thermodynamics (see also the answer by higgsss).
"Convexity: An Analytic Viewpoint" by Barry Simon.