A general result on the convex function

Let $g(t):=\max(f(t),f(1-t))$. Then the function $g$ is convex on $[0,1]$, symmetric about $1/2$, and hence nonincreasing on $[0,1/2]$. So, letting $h(s):=g(s/2)$ for $s\in[0,1]$, we rewrite the inequality in question as $$\int_0^ 1 h\ge\frac23\,h\Big(\frac23\Big)+\frac13\,h\Big(\frac13\Big) \tag{1}\label{1},$$ where the function $h$ is convex and nonincreasing on $[0,1]$.

So, letting $h_z(s):=\max(0,z-s)$ for some $z\in(0,1]$ and all $s\in[0,1]$, we see that the function $h$ is a nonnegative mixture of a constant and functions $h_z$ -- see the detail at the end of this answer.

So, it is enough to verify \eqref{1} when $h$ is either a constant or of the form $h_z$ for some $z\in(0,1]$. This verification is elementary and straightforward.

Let me comment on the nice comment by Paata Ivanisvili. The first part of the GPT-5 Pro solution is essentially the same as the above reduction of the inequality in question to inequality \eqref{1} for functions $h$ that are convex and nonincreasing on $[0,1]$.

The second parts of the two solutions are quite different, though. That of the GPT-5 Pro solution is rather ad hoc and creative.

On the other hand, the second part above requires no thinking or creativity. Moreover (or in other words), it is much more generally applicable. Indeed, with the same approach, it will be easy to create and then elementary and straightforward to verify any inequality of the form $$\int_0^ 1 h\ge\sum_{i\in I}p_i\, h(b_i), \tag{2}\label{2}$$ where again the function $h$ is convex and nonincreasing on $[0,1]$ and for each $i$ in an at most countable set $I$ we have the following: $A_i$ is a measurable subset of $[0,1]$, $p_i:=\int_{A_i}dx>0$, $b_i\in[a_i,1]$, $a_i:=\int_{A_i}x\,dx/p_i$ is the barycenter of $A_i$, and $(A_i)_{i\in I}$ is a partition of the interval $[0,1]$.

Inequality \eqref{2} holds because for each $i\in I$ we have $\frac1{p_i}\,\int_{A_i} h\ge h(a_i)\ge h(b_i)$, by Jensen's inequality for the convex function $h$ and because $h$ is nonincreasing.

If we only know the $p_i$'s and the $b_i$'s (but not the $A_i$'s), then, to prove inequality \eqref{2}, we can just note that this inequality trivially holds when $h$ is a constant and then verify \eqref{2} for $h$ of the form $h_z$ for some $z\in(0,1]$, which is elementary and straightforward.

Inequality \eqref{1} is a special case of \eqref{2}, with $I=\{1,2\}$, $A_1=(1/3,1]$, $A_2=[0,1/3]$, $a_1=2/3$, $a_2=1/6$, $b_1=2/3\ge a_1$, $b_2=1/3\ge a_2$. In particular, this provides another, and this time immediate, proof of \eqref{1}.

Detail: Our function $h$ is convex and nonincreasing on $[0,1]$. So, letting $h'$ denote the left derivative of the convex function $h$ on $(0,1]$, for any $u\in(0,1]$ we have $$h'(1)-h'(u)=\int_{[u,1)}\mu(dz),$$ where $\mu$ is the (nonnegative) Lebesgue--Stieltjes measure corresponding to the nondecreasing left-continuous function $h'$. So, for any $s\in(0,1]$ $$h(s)=h(1)-\int_s^1 du\,h'(u) \\ =h(1)-h'(1)(1-s)+\int_s^1 du\,\int_{[u,1)}\mu(dz) \\ =h(1)-h'(1)(1-s)+\int_{[s,1)}\mu(dz)\int_{(s,z]} du \\ =h(1)-h'(1)h_1(s)+\int_{[s,1)}\mu(dz)h_z(s) \\ =h(1)+|h'(1)|h_1(s)+\int_{(0,1)}\mu(dz)h_z(s);$$ here it was taken into account that $h$ is nonincreasing and hence $h'\le0$.

So, indeed, $h$ is a nonnegative mixture of the constant $h(1)$ and functions $h_z$ for $z\in(0,1]$ -- that is, for some nonnegative measure $\nu$ and all $s\in(0,1]$ we have $$h(s)=h(1)+\int_{(0,1]}\nu(dz)h_z(s).$$