Calculate in closed form $\int_0^1 \int_0^1 \frac{dx\,dy}{1-xy(1-x)(1-y)}$

Let $x=\sin^2(t)$ and $y=\sin^2(w)$, we obtain \begin{align} I & = \int_0^{\pi/2}\int_0^{\pi/2} \dfrac{4\sin(t)\cos(t)\sin(w)\cos(w)dtdw}{1-\sin^2(t)\cos^2(t)\sin^2(w)\cos^2(w)}\\ & = 4\sum_{k=0}^{\infty}\int_0^{\pi/2}\int_0^{\pi/2} \sin^{(2k+1)}(t) \cos^{(2k+1)}(t)\sin^{(2k+1)}(w) \cos^{(2k+1)}(w)dtdw\\ & = 4 \sum_{k=0}^{\infty} \left(\int_0^{\pi/2}\sin^{(2k+1)}(t) \cos^{(2k+1)}(t)dt\right)^2\\ & = 4 \sum_{k=0}^{\infty} \dfrac1{4^{2k+1}}\left(\int_0^{\pi/2}\sin^{(2k+1)}(2t)dt\right)^2 = \sum_{k=0}^{\infty} \dfrac1{4^{2k}}\left(\int_0^{\pi}\sin^{(2k+1)}(t)\dfrac{dt}2\right)^2\\ & = \sum_{k=0}^{\infty} \dfrac1{4^{2k}}\left(\int_0^{\pi/2}\sin^{(2k+1)}(t)dt\right)^2 \end{align} Recall that $$\int_0^{\pi/2}\sin^{(2k+1)}(t)dt=\dfrac{4^k}{2k+1}\dfrac1{\dbinom{2k}k}$$ This gives us $$I = \sum_{k=0}^{\infty} \dfrac1{(2k+1)^2} \dfrac1{\dbinom{2k}k^2} = _3F_2\left(1,1,1;\frac32,\frac32;\frac1{16}\right)$$

An incomplete answer.

Tackling the inner integral with elementary methods, \begin{align} \int^1_0\frac{{\rm d}x}{1-ax+ax^2} &=\int^1_0\frac{{\rm d}x}{a\left(x-\frac{1}{2}\right)^2+1-\frac{a}{4}}\\&=\int^\frac{1}{2}_{-\frac{1}{2}}\frac{{\rm d}x}{ax^2+1-\frac{a}{4}}\\ &=\frac{1}{\sqrt{a\left(1-\frac{a}{4}\right)}}\int^\frac{\sqrt{a}}{2\sqrt{(1-\frac{a}{4})}}_{-\frac{\sqrt{a}}{2\sqrt{(1-\frac{a}{4})}}}\frac{{\rm d}x}{1+x^2}\\&=\frac{4}{\sqrt{a(4-a)}}\arctan\left(\sqrt{\frac{a}{4-a}}\right) \end{align} Thus \begin{align} \iint\limits_{[0,1]^2}\frac{{\rm d}A}{1-xy(1-x)(1-y)} &=4\int^1_0\arctan\left(\sqrt{\frac{x(1-x)}{x^2-x+4}}\right)\frac{{\rm d}x}{\sqrt{x(1-x)(x^2-x+4)}}\\ &=16\int^1_0\arctan\left(\sqrt{\frac{1-x^2}{15+x^2}}\right)\frac{{\rm d}x}{\sqrt{(1-x^2)(15+x^2)}}\\ &=\frac{16}{\sqrt{15}}\int^1_0\operatorname{F}\left(\arcsin{x}\right)\frac{x\ {\rm d}x}{\sqrt{(1-x^2)(15+x^2)}}\\ &=\frac{8}{15}\int^{\operatorname{K}}_{-\operatorname{K}}u\operatorname{sn}(u)\ {\rm d}u \end{align} where $\operatorname{sn}(u)$ is a Jacobi elliptic function with modulus $k=15^{-1/2}i$.

Let $R$ be the parallelogram with vertices $-\operatorname{K}-2i\operatorname{K}'$, $\operatorname{K}-2i\operatorname{K}'$, $\operatorname{K}+2i\operatorname{K}'$, $-\operatorname{K}+2i\operatorname{K}'$. Parameterising along the contour, \begin{align} \int_{\partial R}z^{\color{red}{2}}\operatorname{sn}(z)\ {\rm d}z &=-8i\operatorname{K}\left(\frac{4}{\sqrt{15}}\right)\int^{\operatorname{K}}_{-\operatorname{K}}u\operatorname{sn}(u)\ {\rm d}u+2i\int^{2\operatorname{K}'}_{-2\operatorname{K}'}(\operatorname{K}^2-t^2)\operatorname{cd}(it)\ {\rm d}t\\ &=2\pi i\sum_{\pm}\operatorname*{Res}_{z=\pm i\operatorname{K}'}z^2\operatorname{sn}(z)=-4\pi\sqrt{15}\operatorname{K}\left(\frac{4}{\sqrt{15}}\right)^2 \end{align} $\operatorname{cd}(t)$ has an anti-derivative, namely $$2i\operatorname{K}^2\int^{2\operatorname{K}'}_{-2\operatorname{K}'}\operatorname{cd}(it)\ {\rm d}t=-4i\sqrt{15}\operatorname{K}^2\left.\ln\left(\operatorname{nd}(t)+\frac{i}{\sqrt{15}}\operatorname{sd}(t)\right)\right|^{2i\operatorname{K}'}_0$$ Otherwise, this seems to be as far as I can go for now.

From hypergeometric function to elliptic integrals

The hypergeometric function $_3F_2$ from your link (or user17762's answer) may be rewritten (using following integral representation and $_2F_1\left(1,1;\frac32; t\right)=\frac d{dt}\left[\left(\arcsin\sqrt{t}\,\right)^2\right]$) as : \begin{align} \tag{1}I&=\;_3F_2\left(1,1,1;\frac32,\frac32;\frac1{16}\right)=\frac 12\int_0^1\frac{_2F_1\left(1,1;\frac32;\frac x{16}\right)}{\sqrt{1-x}}dx\\ &=2\int_0^1\frac{\arcsin\frac{\sqrt{x}}4}{\sqrt{x\,(1-x)\left(1-\frac x{16}\right)}}dx\\ \tag{2}&=16\int_0^{1/4}\frac{\arcsin u}{\sqrt{(1-16 u^2)\left(1-u^2\right)}}\,du\\ &\tag{3}=16\int_0^{\arcsin(1/4)} \frac {y}{\sqrt{1-4^2\sin(y)^2)}}\;dy\\ &=\left.16\;y\;F\left(y, 4\right)\right|_0^{\arcsin(1/4)}-16\int_0^{\arcsin(1/4)}\;F\left(y, 4\right)\;dy\\ \tag{4}&=4\;\arcsin(1/4)\;K\left(\frac 1{4}\right)-16\int_0^{\arcsin(1/4)}\;F\left(y, 4\right)\;dy\\ \end{align}

With $\;\displaystyle F(\phi,k):=\int_0^{\phi} \frac {dy}{\sqrt{1-k^2\sin(y)^2)}},\;K(k):=F\left(\frac{\pi}2,k\right)\;$ the incomplete and complete elliptic integral of the first kind $F$ and $K$ (one should take care of replacing $F(\phi,k)$ by EllipticF[$\phi,k^2$] and $K(k)$ by EllipticK[$k^2$] while using Alpha/Mathematica).

'$I$' may be expressed with elliptic integrals but this doesn't seem to help for a closed form... Observing that $\;\displaystyle K\left(\frac 1{4}\right)=\frac{\pi}2\sum_{k=0}^\infty\frac{\binom{2k}{k}^2}{2^{8\,k}}\;$ let's start it all again using the

Squared central binomial series

We want $\;\displaystyle I=\sum_{k=0}^\infty \frac 1{(2k+1)^2\binom{2k}{k}^2}\;$ but Nected Batir proved in 2004 following more general formula (for $n\in\mathbb{N},\; n\ge3$) :

\begin{align} \tag{5}I_n(x)&:=\sum_{k=0}^\infty \frac{x^{2k+1}}{(2k+1)^n\,\binom{2k}{k}^2}\\ &=\frac{4\,(-1)^{n-3}}{(n-3)!}\int_0^{\pi/2}\int_0^{\arcsin(x\cos(y)/4)}\frac t{\sin(t)}\;\log^{n-3}\left(\frac{4\sin(t)}{x\cos(y)}\right)\;dt\,dy\\ \end{align} For $\,n=3\,$ this is simply $\;\displaystyle I_3(x)=\int_0^{\pi/2}\int_0^{\arcsin(x\cos(y)/4)}\frac t{\sin(t)}\;dt\,dy\;\;$ and after derivation \begin{align} I_2(x)&=x\,I_3(x)'=4\,x\int_0^{\pi/2} \frac {\arcsin(x\cos(y)/4)}{x\cos(y)/4} \frac {\partial}{\partial x}\left[\arcsin(x\cos(y)/4)\right]\;dy\\ \tag{6}I_2(x)&=4\int_0^{\pi/2} \frac {\arcsin(x\cos(y)/4)}{\sqrt{1-(x\cos(y)/4)^2}} \;dy\\ \end{align} Unfortunately this is merely $(2)$ in the case $x=1$ (the integral for $x=4$ is divergent and $I_3(4)=8\pi G-14\zeta(3)\;$ but this won't help here...) so nothing complete here either...

The journey may be of interest anyway...