Calculating in closed form $\int _0^1\int _0^1\frac{1}{1+x y (x+y)} \ dx \ dy$ Integrating with respect to a variable and then to the other one, things look pretty
complicated, but I'm sure you have ideas that might simplify the job to do here.
This time we're talking about 
$$\int _0^1\int _0^1\frac{1}{1+x y (x+y)} \ dx \ dy$$

The bounty moment: after 2 years and 8 months from the releasing moment of the question, it's time for a 300 points bounty for finding the simplest closed-form of the integral!

Supplementary question: Here is an extension of the question for those with a need for more challenging questions.
Calculate
$$\int _0^1\int _0^1\cdots\int _0^1\frac{1}{1+x_1 x_2\cdots x_n (x_1+x_2+\cdots +x_n)} \ \textrm{d}x_1 \ \textrm{d}x_2\cdots \textrm{d}x_n, \ n\ge 2.$$
Last but not least, special greetings to Cleo!
 A: By symmetry (i.e. by exploiting the fact that our integral is twice the integral over the sub-region $0\leq y\leq x\leq 1$) we just have to compute:
$$ I=2 \int_{0}^{1}\int_{0}^{1}\frac{x}{1+x^3 y(1+y)}\,dy\,dx = 2 \int_{0}^{1}\int_{0}^{2}\frac{x}{(1+x^3y)\sqrt{1+4y}}\,dy\,dx$$
Integrating with respect to $x$ first,
$$\begin{eqnarray*} I &=& 2\int_{0}^{2}\frac{2 \sqrt{3}\arctan\left(\frac{\sqrt{3}}{1-2y^{1/3}}\right)-2\log\left(1+y^{1/3}\right)+\log\left(1-y^{1/3}+y^{2/3}\right)}{6 y^{2/3}\sqrt{1+4y}}\,dy\\&=&\int_{0}^{2^{1/3}}\frac{2\sqrt{3}\arctan\left(\frac{\sqrt{3}}{1-2y}\right)-3\log(1+y)+\log(1+y^3)}{\sqrt{1+4y^3}}\,dy\end{eqnarray*}$$
but the resulting integrals in just one variable do not look so appealing.
Am I missing some crucial simplification that follows from replacing $y$ with a Jacobi elliptic function (maybe $\text{dn}$) or with the Weierstrass elliptic function $\wp(z)$ (corresponding to $g_2=0,g_3=-1$) then exploiting some weird/mystical product formulas?
A: $\def\l{\ell}$I have found (many) integral forms no nicer than those in the other answers.
Here I give a series expansion.
We have
\begin{align*}
\frac{1}{1+xy(x+y)}
 &= \frac{1}{2\left(1+\frac{xy(x+y)-1}{2}\right)} \\
    &= \sum_{k=0}^\infty \frac{(-1)^k}{2^{k+1}} [xy(x+y)-1)]^k \\
    &= \sum_{k=0}^\infty \frac{(-1)^k}{2^{k+1}}
    \sum_{\l=0}^k {k\choose \l} [xy(x+y)]^\l (-1)^{k-\l} \\
    &= \sum_{k=0}^\infty \frac{(-1)^k}{2^{k+1}}
    \sum_{\l=0}^k (-1)^{k-\l}{k\choose \l}  x^\l y^\l
    \sum_{m=0}^\l {\l\choose m} x^{\l-m}y^m \\
    &= \sum_{k=0}^\infty \sum_{\l=0}^k \sum_{m=0}^\l
    \frac{(-1)^\l}{2^{k+1}}{k\choose \l}{\l\choose m}
    x^{2\l-m}y^{\l+m}.
\end{align*}
Thus,
\begin{align}\tag{1}
\int_0^1\int_0^1 \frac{1}{1+xy(x+y)} dx\, dy
 &= \sum_{k=0}^\infty \sum_{\l=0}^k \sum_{m=0}^\l
    \frac{(-1)^\l}{2^{k+1}}{k\choose \l}{\l\choose m}
    \frac{1}{(2\l-m+1)(\l+m+1)}.
\end{align}
Similarly,
\begin{align*}
\int_0^1 \cdots \int_0^1 & \frac{1}{1+x_1\cdots x_n(x_1+\ldots +x_n)} dx_1\cdots dx_n \\
 &= \sum_{k=0}^\infty \sum_{\l=0}^k
    \sum_{{m_1,\ldots,m_n}\atop{m_1+\ldots +m_n=\l}}
    (-1)^\l\frac{(n-1)^{k-\l}}{n^{k+1}}
    {k\choose \l} \frac{\l!}{m_1!\cdots m_n!}
    \frac{1}{(\l+m_1+1)\cdots (\l+m_n+1)}.
\end{align*}
It is possible to write the sum in equation (1) in terms of a single (infinite) sum over special functions, but it is not particularly illuminating.
I am curious to see if there is some nice closed form at all! 
A: Another single integral form,
$$ 2 \int_0^1 \frac{\arctan{\big(\sqrt{y(4-y^3)}/(y^2+2)\big)}}{\sqrt{y(4-y^3)}}\,dy. $$
It doesn't look simpler than the first 1D integral given, however, so I'll skip the proof.
A: Starting in the polar coordinates,
\begin{align}
&J=\int\limits_0^1\int\limits_0^1 \dfrac{dx\,dy}{1+xy(x+y)} = 2\int\limits_0^{\dfrac\pi4}\int\limits_0^\dfrac1{\cos\varphi}\dfrac{\rho\,d\rho\,d\varphi}{1+\rho^3\cos\varphi\sin\varphi(\cos\varphi+\sin\varphi)} = 2\int\limits_0^{\pi/4}F(\rho,\varphi)\,d\varphi,\\
\end{align}
and substitution
\begin{align}
&\rho=\dfrac1{v \cos\varphi}\tag1
\end{align}
gives
\begin{align}
&F(\rho,\varphi) = \int\limits_1^\infty\dfrac1 {1+\dfrac1{v^3\cos^3\varphi}\cos\varphi\sin\varphi(\cos\varphi+\sin\varphi)}\,\dfrac{dv}{v^3\cos^2\varphi} =\\ &\dfrac1{\cos^2\varphi}\int\limits_1^\infty\dfrac{dv}{\tan\varphi(\tan\varphi+1)+ v^3},\quad\text{or}\\
&J=2\int\limits_1^\infty G(v,\varphi)\, dv,\quad\text{where}\\
&G(v,\varphi) = \int\limits_0^{\pi/4}\dfrac{dv}{\tan\varphi(\tan\varphi+1)+v^3}\cdot\dfrac{d\varphi}{\cos^2\varphi}=\int\limits_0^1\dfrac{dt}{t^2+t+v^3} = \\
&\dfrac1{\sqrt{v^3-\dfrac14}}\left.\arctan\dfrac{t+\dfrac12}{\sqrt{v^3-\dfrac14}}\right|_{t=0}^1 = \dfrac1{\sqrt{v^3-\dfrac14}}\left(\arctan\dfrac{\dfrac32}{\sqrt{v^3-\dfrac14}}-\arctan\dfrac{\dfrac12}{\sqrt{v^3-\dfrac14}}\right),\tag2\\[4pt]
&G(v,\varphi) = \dfrac1{\sqrt{v^3-\dfrac14}}\arctan\dfrac{\dfrac1{\sqrt{v^3-\dfrac14}}}{1+\dfrac3{4\left(v^3-\dfrac14\right)}} = \dfrac2{\sqrt{4v^3-1}}\arctan\dfrac{2\sqrt{4v^3-1}}{4v^3+2}.\tag3
\end{align}
Therefore,
\begin{align}
&J=\int\limits_1^\infty \dfrac4{\sqrt{4v^3-1}}\arctan\dfrac{2\sqrt{4v^3-1}}{4v^3+2}\, dv\approx 0.798965\tag4
\end{align}
(see also Wolfram Alpha).
Substitution
$$u=\sqrt{4v^3-1}$$
leads to
\begin{align}
&v=\dfrac{\sqrt[3]2}2\sqrt[3]{u^2+1},\quad dv=\dfrac{2\sqrt[3]2}3\dfrac1{\sqrt[3]{(u^2+1)^2}}u\,du,\\
&J=\dfrac{4\sqrt[3]2}3\int\limits_\sqrt3^\infty\arctan\dfrac {2u}{3+u^2}\dfrac{du}{\sqrt[3]{(u^2+1)^2}}\tag5\\
&J=\dfrac{4\sqrt[3]2}3\int\limits_\sqrt3^\infty\arctan\dfrac {u-\dfrac13u}{1+u\cdot\dfrac13u}\dfrac{du}{\sqrt[3]{(u^2+1)^2}},\\
&J=\dfrac{4\sqrt[3]2}3\int\limits_\sqrt3^\infty\left(\arctan u - \arctan\dfrac u3\right)\dfrac{du}{\sqrt[3]{(u^2+1)^2}} = J_1-J_2, \tag6\\
\end{align}
(see also Wolfram Alpha), wherein
$$J_1=\dfrac{4\sqrt[3]2}3\int\limits_{\pi/3}^{\pi/2}\dfrac{z\,dz}{\sqrt[3]{\cos^2z}}\\$$
has the closed form.
But the further attempts to get a closed form weren't successful. 
A: Why do you think that there is a closed form?  Most integrals do not have one, and Maple cannot find one.  The numerical value is 0.79896482380785081628946784922318984550713669761340, which is not recognised by the Inverse Symbolic Calculator (https://isc.carma.newcastle.edu.au/advancedCalc).  Given these facts, you should assume that there is no closed form unless you have a very good reason to think otherwise.
A: Another 1-D integral form is the following.
Making the substitution
$$
\left\{ \matrix{  x = s - u \hfill \cr y = s + u \hfill \cr}  \right.
  \quad \left\{ \matrix{  s = \left( {y + x} \right)/2 \hfill \cr   u = \left( {y - x} \right)/2 \hfill \cr}  \right.
$$
the Jacobian of which is $2$, we can write
$$
\eqalign{
  & I = \int_{x = 0}^1 {\int_{y = 0}^1 {{1 \over {1 + xy\left( {x + y} \right)}}dxdy} }
  = 2\left( {\int_{x = 0}^1 {\int_{y = 0}^1 {{1 \over {1 + 2\left( {s^{\,2}  - u^{\,2} } \right)s}}dsdu} } } \right) =   \cr 
  &  = 4\left( {\int_{s = 0}^{1/2} {\int_{u = 0}^s {{1 \over {1 + 2\left( {s^{\,2}  - u^{\,2} } \right)s}}dsdu} }
  + \int_{s = 1/2}^1 {\int_{u = 0}^{1 - s} {{1 \over {1 + 2\left( {s^{\,2}  - u^{\,2} } \right)s}}dsdu} } } \right) =   \cr 
  &  = 4\left( {\int_{s = 0}^{1/2} {\int_{u = 0}^s {{1 \over {1 + 2\left( {s^{\,2}  - u^{\,2} } \right)s}}dsdu} }
  + \int_{s = 0}^{1/2} {\int_{u = 0}^s {{1 \over {1 + 2\left( {\left( {1 - s} \right)^{\,2}  - u^{\,2} } \right)\left( {1 - s} \right)}}dsdu} } } \right) =   \cr 
  &  = 4\int_{s = 0}^{1/2} {\int_{u = 0}^s {\left( {g(s,u) + g(1 - s,u)} \right)dsdu} }  \cr} 
$$
where
$$
\eqalign{
  & g(s,u) = {1 \over {1 + 2\left( {s^{\,2}  - u^{\,2} } \right)s}} = {1 \over {1 + 2s^{\,3}  - 2su^{\,2} }} =   \cr 
  &  = {1 \over {\left( {1 + 2s^{\,3} } \right)\left( {1 - \left( {\left( {\sqrt {{{2s} \over {1 + 2s^{\,3} }}} } \right)u} \right)^{\,2} } \right)}} \cr} 
$$
Putting
$$
A(s) = 1 + 2s^{\,3} \quad a\left( s \right) = \sqrt {{{2s} \over {1 + 2s^{\,3} }}} 
$$
it is easily verified that
$$
0 \le u \le s \le 1/2\quad  \Rightarrow \quad \left\{ \matrix{
  0 \le a\left( s \right) \le a\left( {1 - s} \right) < 1 \hfill \cr 
  0 \le u\,a\left( s \right) \le u\,a\left( {1 - s} \right) < 1/2 \hfill \cr 
  0 \le s\,a\left( s \right) \le s\,a\left( {1 - s} \right) < 1/2 \hfill \cr}  \right.
$$
so we can make the substitution
$$
au = \sin t\quad du = \,{1 \over a}\cos tdt
$$
to get
$$
\eqalign{
  & \int_{u = 0}^s {g(s,u)du}  = {1 \over {a\left( s \right)A(s)}}\int_{t = 0}^{\arcsin (as)} {{1 \over {\cos t}}dt}
  = {1 \over {2a\left( s \right)A(s)}}\ln \left( {{{1 + sa(s)} \over {1 - sa(s)}}} \right)  \cr 
  & \int_{u = 0}^s {g(1 - s,u)du}  = {1 \over {2a\left( {1 - s} \right)A(1 - s)}}\ln \left( {{{1 + sa(1 - s)} \over {1 - sa(1 - s)}}} \right) \cr} 
$$
that is:
$$ \bbox[lightyellow] {  
I = 4\int_{s = 0}^{1/2} {\left( {{1 \over {2a\left( s \right)A(s)}}\ln \left( {{{1 + sa(s)} \over {1 - sa(s)}}} \right)
 + {1 \over {2a\left( {1 - s} \right)A(1 - s)}}\ln \left( {{{1 + sa(1 - s)} \over {1 - sa(1 - s)}}} \right)} \right)ds} 
}$$
which can be checked to give the correct result.
In spite of its look, the integrand is a "well behaved" function , and since either
$s\,a(s)$ and $s\,a(1-s)$ are bounded as above, it might be possible to develop it
into a series with coefficients expressable in a closed form.
The first  summand in fact can be expressed as
$$
\eqalign{
  & {1 \over {2a\left( s \right)A(s)}}\ln \left( {{{1 + sa(s)} \over {1 - sa(s)}}} \right) =   \cr 
  &  = s\sum\limits_{0\, \le \,n} {{{\left( {2s^{\,3} } \right)^{\,n} } \over {2n + 1}}{1 \over {\left( {1 + 2s^{\,3} } \right)^{\,n + 1} }}}  =   \cr 
  &  = s\sum\limits_{0\, \le \,n} {{{\left( {2s^{\,3} } \right)^{\,n} } \over {2n + 1}}\sum\limits_{0\, \le \,m} {\left( { - 1} \right)^{\,m} \left( \matrix{
  n + m \cr 
  n \cr}  \right)\left( {2s^{\,3} } \right)^{\,m} } }  =   \cr 
  &  = s\sum\limits_{0\, \le \,n} {\left( {\sum\limits_{0\, \le \,k\,\left( { \le \,n} \right)} {{{\left( { - 1} \right)^{\,k} } \over {2k + 1}}\left( \matrix{
  n \cr 
  k \cr}  \right)} } \right)\left( { - 2s^{\,3} } \right)^{\,n} }  =   \cr 
  &  = \sum\limits_{0\, \le \,n} {\left( {{{\left( { - 1} \right)^{\,n} 2^{\,3n} \left( {n!} \right)^{\,2} } \over {\left( {2n + 1} \right)!}}} \right)s^{\,3n + 1} }  \cr} 
$$
which is easily integrable and converging fast enough.
I am looking for a suitable form for the second as well.
