Note that I didn't really look to closely at your picture, as it was quite messy, so if this is what you did then you got it right.
Notice that $S$ is the graph of the function
$$f:\Omega\to\mathbb{R},\quad (x,y)\mapsto x^2+y,$$
where
$$\Omega=\{(x,y)\in\mathbb{R}^2 : 0\leq x\leq y\land x+y\leq 1\}.$$
Notice also that
$$\nabla f(x,y)=(2x,1),$$
from which we have that
$$\lVert \nabla f(x,y)\rVert^2=4x^2+1.$$
We then have that
$$I=\iint_S\frac{xy}{\sqrt{1+2x^2}}~\mathrm{d}S=\iint_\Omega \frac{xy}{\sqrt{1+2x^2}}\sqrt{1+4x^2+1}~\mathrm{d}x~\mathrm{d}y=\iint_\Omega xy\frac{\sqrt{2}\sqrt{1+2x^2}}{\sqrt{1+2x^2}}~\mathrm{d}x~\mathrm{d}y=\sqrt{2}\iint_\Omega xy~\mathrm{d}x~\mathrm{d}y.$$
Now let's turn this into an iterated integral. Let's take $x$ as our outer integration variable, and so we can let $y$ depend on $x$ (the other way can also be done, but it gets more messy, and will lead to having to split it into two integrals). Notice that we can rewrite
$$x+y\leq 1$$
as
$$y\leq 1-x.$$
Combining this with the inequality
$$0\leq x\leq y$$
gives us then that $\Omega$ can be expressed as
$$\Omega=\left\{(x,y)\in\mathbb{R}^2 : 0\leq x\leq \frac{1}{2} \land x\leq y\leq 1-x\right\}$$
($x\leq\frac{1}{2}$ as $x> x-1$ otherwise). This finally means that
$$I=\sqrt{2}\int_0^\frac{1}{2}\int_x^{1-x}xy~\mathrm{d}y~\mathrm{d}x.$$
At this point it's just a standard computation so I'll leave the details to you.