@HansLundmark has already given away the answer. To prove it we divide the region $R=\{(x,y)|0\le x\le1,0\le y\le1\}$ into $R_1=\left\{(x,y)|\frac12-r<x<\frac12+r,\frac12-r<y<\frac12+r\right\}$ and $R_2=R\setminus R_1$. Then
$$\begin{align}I=\int\int_R\left[xy(1-x)(1-y)\right]^nf(x,y)d^2A&=\int\int_R\left[xy(1-x)(1-y)\right]^nf\left(\frac12,\frac12\right)d^2A\\
&\quad+\int\int_{R_1}\left[xy(1-x)(1-y)\right]^n\left[f(x,y)-f\left(\frac12,\frac12\right)\right]d^2A\\
&\quad+\int\int_{R_2}\left[xy(1-x)(1-y)\right]^n\left[f(x,y)-f\left(\frac12,\frac12\right)\right]d^2A\\
&=I_0+I_1+I_2\end{align}$$
Now,
$$I_0=f\left(\frac12,\frac12\right)\left[B(n+1,n+1)\right]^2=\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2f\left(\frac12,\frac12\right)$$
Since $f(x,y)$ is continuous at $(x,y)=\left(1/2,1/2\right)$ we know that for any $\epsilon_1>0$ there is a $\delta_1(\epsilon_1)>0$ such that $\left|f(x,y)-f\left(\frac12,\frac12\right)\right|<\epsilon_1$ whenever
$$\sqrt{\left(x-\frac12\right)^2+\left(y-\frac12\right)^2}<\delta_1$$
So we choose
$$0<r<\min\left(\frac{\delta_1\left(\frac{\epsilon}2\right)}{\sqrt2},\frac12\right)$$
So that
$$\begin{align}\left|I_1\right|&\le\int\int_{R_1}\left[xy(1-x)(1-y)\right]^n\left|f(x,y)-f\left(\frac12,\frac12\right)\right|d^2A\\
&<\int\int_{R_1}\left[xy(1-x)(1-y)\right]^n\left(\frac{\epsilon}2\right)d^2A\\
&<\int\int_{R}\left[xy(1-x)(1-y)\right]^n\left(\frac{\epsilon}2\right)d^2A\\
&=\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2\left(\frac{\epsilon}2\right)\end{align}$$
If $n\ge2$ then $\frac83n=2n+\frac23n>2n+1$ so
$$\begin{align}\frac{(2n+1)!}{n\cdot2^{2n}\left(n!\right)^2}&=\frac{(2n+1)}n\prod_{k=1}^n\frac{(2k-1)(2k)}{(2k)(2k)}=\frac{(2n+1)}n\left(\frac12\right)\left(\frac34\right)\prod_{k=3}^n\frac{(2k-1)}{(2k)}\\
&=\frac{(2n+1)}{\frac83n}\prod_{k=3}^n\frac{(2k-1)}{(2k)}\le1\end{align}$$
Also the continuity of $f$ over a closed set $R$ implies that it is bounded on that set so there is some $M$ such that $\left|f(x,y)\right|\le M$ for all $(x,y)\in R$, so by the triangle inequality
$$\left|f(x,y)-f\left(\frac12,\frac12\right)\right|\le2M$$
For $(x,y)\in R\supset R_2$. Also $xy(1-x)(1-y)$ attains its maximum absolute value in $R_2$ at $\left(\frac12\pm r,\frac12\right)$ and at $\left(\frac12,\frac12\pm r\right)$ of $\frac1{16}\left(1-4r^2\right)$. Then
$$\begin{align}\left|I_2\right|&\le\int\int_{R_2}\left[xy(1-x)(1-y)\right]^n\left|f(x,y)-f\left(\frac12,\frac12\right)\right|d^2A\\
&\le\int\int_{R_2}\frac1{2^{4n}}\left(1-4r^2\right)^n(2M)d^2A\\
&=\int\int_{R_2}\left[\frac{(2n+1)!}{n\cdot2^{2n}\left(n!\right)^2}\right]^2\left[\frac{n\cdot\left(n!\right)^2}{(2n+1)!}\right]^2\left(1-4r^2\right)^n(2M)d^2A\\
&\le\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2n^2\left(1-4r^2\right)^n(2M)\int\int_{R_2}d^2A\\
&\le\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2n^2\left(1-4r^2\right)^n(2M)\int\int_{R}d^2A\\
&=\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2(2M)n^2\left(1-4r^2\right)^n\end{align}$$
We may establish via L'Hopital's Rule that for $0<r<\frac12$
$$\lim_{n\rightarrow\infty}g(n,r)=\lim_{n\rightarrow\infty}n^2\left(1-4r^2\right)^n=0$$
Which means that for any $\epsilon_2>0$ there is an $N_2\left(\epsilon_2,r\right)$ such that $g(n,r)<\epsilon_2$ whenever $n>N_2\left(\epsilon_2,r\right)$ so we choose
$$n>\max\left(N_2\left(\frac{\epsilon}{4M},r\right),2\right)$$
Which yields the desired
$$\left|I_2\right|\le\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2\left(\frac{\epsilon}2\right)$$
Given such an $n$ we can now put together
$$\begin{align}\left|\left[\frac{(2n+1)!}{\left(n!\right)^2}\right]^2I-f\left(\frac12,\frac12\right)\right|&=\left|\left[\frac{(2n+1)!}{\left(n!\right)^2}\right]^2\left(I_0+I_1+I_2\right)-f\left(\frac12,\frac12\right)\right|\\
&\le\left|\left[\frac{(2n+1)!}{\left(n!\right)^2}\right]^2\left(\left|I_1\right|+\left|I_2\right|\right)\right|\\
&<\left|\frac{\epsilon}2+\frac{\epsilon}2\right|=\epsilon\end{align}$$
Let's see... we divided by $4M$ earlier, but if $4M=0$ then the limit is trivially $f\left(1/2,1/2\right)=0$. So this limit gives us another realization of the Dirac $\delta$ function.