Evaluating $\int_0^1\frac{x^{2/3}(1-x)^{-1/3}}{1-x+x^2}dx$
Update
I have now finally found a way to take my hypergeometric solution all the way to its final elementary form.
Let $$I = \int_0^1 \frac{x^{2/3}}{\sqrt[3]{1 - x} (1 - x +x^2)} \, dx. \tag1$$ Enforcing a substitution of $x \mapsto 1 - x$ leads to $$I = \int_0^1 \frac{(1 - x)^{2/3}}{\sqrt[3]{x} (1 - x + x^2)} \, dx. \tag2$$ Adding (1) to (2) produces $$I = \frac{1}{2} \int_0^1 \frac{dx}{\sqrt[3]{x - x^2} (1 - x + x^2)}.$$
Expanding the second term appearing in the denominator in terms of a geometric series, we have \begin{align} I &= \frac{1}{2} \int_0^1 \frac{dx}{(x - x^2)^{1/3} [1 - (x - x^2)]}\\ &= \frac{1}{2} \int_0^1 \frac{1}{(x - x^2)^{1/3}} \sum_{n = 0}^\infty (x - x^2)^n \, dx\\ &= \frac{1}{2} \sum_{n = 0}^\infty \int_0^1 x^{n - 1/3} (1 - x)^{n - 1/3} \, dx\\ &= \frac{1}{2} \sum_{n = 0}^\infty \operatorname{B} \left (n + \frac{2}{3}, n + \frac{2}{3} \right ), \tag3 \end{align} where $\operatorname{B}(x,y)$ is the Beta function. Making use of the result $$\operatorname{B} (x,x) = \frac{\sqrt{\pi} 2^{1 - 2x} \Gamma (x)}{\Gamma \left (x + \frac{1}{2} \right )},$$ the sum in (3) can be written as \begin{align} I &= \frac{\sqrt{\pi}}{2 \sqrt[3]{2}} \sum_{n = 0}^\infty \frac{\Gamma \left (n + \frac{2}{3} \right )}{\Gamma \left (n + \frac{7}{6} \right ) 4^n}\\ &= \frac{\sqrt{\pi}}{2 \sqrt[3]{2}} \cdot \frac{\Gamma (\frac{2}{3})}{\Gamma (\frac{7}{6})} \sum_{n = 0}^\infty \frac{\left (\frac{2}{3} \right )_n (1)_n}{\left (\frac{7}{6} \right )_n 4^n n!}\\ &= \frac{\sqrt{\pi}}{2 \sqrt[3]{2}} \cdot \frac{\Gamma (\frac{2}{3})}{\Gamma (\frac{7}{6})}\ _2F_1 \left (\frac{2}{3}, 1; \frac{7}{6}; \frac{1}{4} \right ),\tag4 \end{align} where $_2F_1 (a,b;c;x)$ is the Gauss hypergeometric function.
To reduce the hypergeometric function that appears in (4) into elementary form, we proceed as follows.
Firstly, since $_2F_1 (a,b;c;x) =\ _2F_1 (b,a;c;x)$ on applying the second of Pfaff's transformations, namely $$_2F_1 (a,b;c;x) = (1 - x)^{-a}\ _2F_1 \left (a,c-b;c;\frac{x}{x - 1} \right ),$$ to the hypergeometric function, we have $$_2F_1 \left (1, \frac{2}{3}; \frac{7}{6}; \frac{1}{4} \right ) = \frac{4}{3}\ _2F_1 \left (1, \frac{1}{2}; \frac{7}{6}; -\frac{1}{3} \right ).\tag5$$ Next, applying Euler's transformation, namely $$_2F_1 (a,b;c;x) = (1 - x)^{c - a - b}\ _2F_1 (c-a,c-b;c;x),$$ we have $$_2F_1 \left (1, \frac{2}{3}; \frac{7}{6}; \frac{1}{4} \right ) = \frac{4^{2/3}}{3^{2/3}}\ _2F_1 \left (\frac{1}{6}, \frac{2}{3}; \frac{7}{6}; -\frac{1}{3} \right ). \tag6$$
Finally, from DLMF: 15.4.31 we see that $$_2F_1 \left (a, \frac{1}{2} + a; \frac{3}{2}-2a; -\frac{1}{3} \right ) = \left (\frac{8}{9} \right )^{-2a} \frac{\Gamma (\frac{4}{3}) \Gamma (\frac{3}{2} - 2a)}{\Gamma (\frac{3}{2}) \Gamma (\frac{4}{3} - 2a)}.$$ Setting $a = 1/6$ leads to $$_2F_1 \left (\frac{1}{6}, \frac{2}{3}; \frac{7}{6}; -\frac{1}{3} \right ) = \frac{\sqrt[3]{9}}{\sqrt{\pi}} \Gamma \left (\frac{4}{3} \right ) \Gamma \left (\frac{7}{6} \right ).$$ Thus (6) becomes $$_2F_1 \left (1, \frac{2}{3}; \frac{7}{6}; \frac{1}{4} \right ) = \frac{2^{4/3}}{\sqrt{\pi}} \Gamma \left (\frac{4}{3} \right ) \Gamma \left (\frac{7}{6} \right ).$$ On substituting this result into (4), one has \begin{align} I &= \Gamma \left (\frac{2}{3} \right ) \Gamma \left (\frac{4}{3} \right )\\ &= \frac{1}{3} \Gamma \left (\frac{2}{3} \right ) \Gamma \left (\frac{1}{3} \right )\\ &= \frac{1}{3} \Gamma \left (1 - \frac{1}{3} \right ) \Gamma \left (\frac{4}{3} \right )\\ &= \frac{\pi}{3 \sin (\frac{\pi}{3})}\\ &= \frac{2\pi}{3 \sqrt{3}}, \end{align} as required.
The solution heavily exploits symmetry of the integrand.
Let $$I = \int_0^1\frac{x^{2/3}(1-x)^{-1/3}}{1-x+x^2} dx $$ Replace $x$ by $1-x$ and sum up gives $$\tag{1} 2I = \int_0^1 \frac{x^{2/3}(1-x)^{-1/3} + (1-x)^{2/3}x^{-1/3}}{1-x+x^2} dx = \int_0^1 \frac{x^{-1/3}(1-x)^{-1/3}}{1-x+x^2} dx$$
Let $\ln_1$ be complex logarithm with branch cut at positive real axis, while $\ln_2$ be the one whose cut is at negative real axis. Denote $$f(z) = \frac{2}{3}\ln_1(x) - \frac{1}{3}\ln_2 (1-x)$$ Then $f(z)$ is discontinuous along the positive axis, but have different jump in $\arg$ across intervals $[0,1]$ and $[1,\infty)$.
Now integrate $g(z) = e^{f(z)}/(1-z+z^2)$ using keyhole contour. Let $\gamma_1$ be path slightly above $[0,1]$, $\gamma_4$ below. $\gamma_2$ be path slightly above $[1,\infty)$, $\gamma_3$ below. It is easily checked that $$\int_{\gamma 1} g(z) dz = I \qquad \qquad \int_{\gamma 4} g(z) dz = I e^{4\pi i/3}$$ $$\int_{\gamma 2} g(z) dz = e^{\pi i/3} \underbrace{\int_1^\infty \frac{x^{2/3}(x-1)^{-1/3}}{1-x+x^2} dx}_J\qquad \int_{\gamma 3} g(z) dz = e^{\pi i} J$$
If we perform $x\mapsto 1/x$ on $J$, we get $\int_0^1 x^{-1/3}(1-x)^{-1/3}/(1-x+x^2)dx$, thus $J = 2I$ by $(1)$.
Therefore $$I(1-e^{4\pi i/3}) + 2I(e^{\pi i / 3} - e^{\pi i}) = 2\pi i\times \text{Sum of residues of } g(z) \text{ at } e^{\pm 2\pi i /3}$$ From which I believe you can work out the value of $I$.
There is a completely elementary way to solve this. In the end, I do not see how to find the elementary primitive in a simple and intuitive way (if others do, then please edit the answer accordingly), but Rubi helped me. For this reason I post this but make it cw. I also would like to thank @JanG who pointed me to the question and who actually was the one doing the first changes of variables.
Set $$ I=\int_0^1\frac{x^{2/3}(1-x)^{-1/3}}{1-x(1-x)}\,dx. $$ By doing $x\mapsto 1-x$ and adding, others have found that $$ I= \frac{1}{2}\int_0^1\frac{x^{-1/3}(1-x)^{-1/3}}{1-x(1-x)}\,dx. $$ Next, let $x=(1+y)/2$. Then the integral becomes $$ I=2^{2/3}\int_{-1}^1\frac{1}{(1-y^2)^{1/3}(3+y^2)}\,dy= 2^{5/3}\int_{0}^1\frac{1}{(1-y^2)^{1/3}(3+y^2)}\,dy. $$ This is very similar to (and just a $y=\cos t$ away from) the integral @Zacky observes in a comment to the question.
This can be put into Rubi, and surprisinlgy the result is elementary, $$ I=\biggl[\frac{1}{\sqrt{3}}\arctan\Bigl(\frac{\sqrt{3}}{y}\Bigr) +\frac{1}{\sqrt{3}}\arctan\Bigl(\frac{\sqrt{3}\bigl(1-(2-2y^2)^{1/3}\bigr)}{y}\Bigr)-\frac{1}{3}\text{artanh}\,y+\text{artanh}\,\Bigl(\frac{y}{1+(2-2y^2)^{1/3}}\Bigr)\biggr]_0^1 $$ Inserting the boundarys, the upper one gives (using a limit) $2\pi/(3\sqrt{3})$ and the lower one gives $0$. Hence $$ I=\frac{2\pi}{3\sqrt{3}}. $$