Constraint
$$((b^3 - abc) + ab(a + b) - c(a^2 + b^2))((abc - b^3) - bc(b + c) + a(b^2 + c^2)) \le 0,$$
or
$$(b^3-a^3)(b^3 - c^3)\ge 0,$$
is equivalent to
$$b^2-(a+c)b+ac\ge 0\tag1.$$
The goal function
$$f(a,b,c)=c^2 - ca + a^2 - \frac{(b + 1)^2(b - 2)}{c + a}$$
under the conditions
$$a^3 + c^3 = m\tag2$$
can be presented in the form of
$$f(a,b,c)=\dfrac1m(c^2-ca+a^2)(m-(b^3-3b-2)).\tag3$$
$\color{brown}{\textbf{At the edges}}.$
$$\max f(0,b,c) = \max f(a,b,0)
= \begin{cases}
\dfrac{m+4}{\sqrt[3]m},\quad\text{if}\quad m\in(0,1]\\
\dfrac{3\sqrt[3]m+2}{\sqrt[3]m},\quad\text{if}\quad m\in\left(1,3\sqrt3\right]\\
\dfrac{m+2}{\sqrt[3]m},\quad\text{if}\quad m\in\left(3\sqrt3,\infty\right).
\end{cases}\tag4$$

The greatest value of
$$g(a,c)=c^2-ca+a^2$$
under the condition $(2)$ can be achieved at the edges of the area or in the inner stationary points.
At the edges,
$$g_e = g\left(0,\sqrt[3]{m}\right) = g\left(0,\sqrt[3]{m}\right) = \sqrt[3]{m^2}.\tag5$$
The inner stationary points
of $g(a,c)$ under the condition $(2),$ in accordance with Lagrange multupliers method, corresponds with the stationary points of
$$G(a,c,\lambda) = a^2-ac+c^2+\lambda(m-a^3-c^3),$$
which can be found from the system $G'_a = G'_c = G'(\lambda) = 0,$ or
$$
\begin{cases}
2a-c-3\lambda a^2=0\\
-a+2c-3\lambda c^2 = 0\\
a^3+c^3=m
\end{cases}\ \Rightarrow
\begin{cases}
(2a-c)c^2=(-a+2c)a^2\\
a^3+c^3=m
\end{cases}\ \Rightarrow
\begin{cases}
(a-c)(a^2-ac+c^2)=0\\
a^3+c^3=m,
\end{cases}
$$
with the solution
$$g_s=g\left(\sqrt[3]{\dfrac m2},\sqrt[3]{\dfrac m2}\right)=\sqrt[3]{\dfrac {m^2}4}.\tag6$$
The global maximum is
$$\max f(a,0,c)=\dfrac{m+2}{\sqrt[3]m}\max(g_e,g_s) = \dfrac{m+2}{\sqrt m}\ \text{ at }\
\dbinom ac\in\left\{\dbinom{\sqrt[3]m}0,\dbinom0{\sqrt[3]m}\right\}.\tag6$$
$\color{brown}{\textbf{The inner stationary points}}$
of $f(a,b,c)$ under the conditions $(1),(2),$ in accordance with Lagrange multupliers method, corresponds with the stationary points of
$$F(a,b,c,\lambda,\mu) = (c^2-ca+a^2)(m-a^3+3a+2)+\lambda(m-a^3-c^3)+\mu(b^2-(a+c)b+ac),$$
which can be found from the system
$$F'_a = F'_b = F'_c = F'(\lambda) = F'_{\mu} = 0,$$
and can be solved similarly.
If the inner maxima are less than on the edges then $(4)$ defines the the global maxima.