How do I proceed towards the evaluation of this surface integral? 
Evaluate $\oint (xy\, dx + xy^2 \,dy)$ by Stokes' Theorem where the
  surface $S$ is a square with the vertices at $(0,1)$ , $(1,0)$ ,
  $(0,-1)$ and $(-1,0)$.


My Attempt: The square is composed of $4$ lines namely -


*

*$y=-x+1$

*$y=x-1$

*$y=-x-1$

*$y=x+1$


From Stokes' Theorem, I get from the given line integral, $$\iint_S (y^2-x) \,dx\, dy$$
Now since the square is symmetrically divided into $4$ parts by the $X$ and $Y$ axes, we can consider only the triangle in the first quadrant and  we get $$\iint_S (y^2-x) \,dx\, dy=4\cdot \int_{x=0}^{x=1}\int_{y=1}^{y=1-x} (y^2-x) \,dx\, dy \tag1$$
Evaluating this, I get the answer as $\boxed{\color{red}{4\cdot \frac{1}{12}= \frac{1}{3}}}$
However, a friend of mine told me that although the square is symmetric, the vector field given by $(xy\hat i+xy^2\hat j)$ is not symmetrically divided into $4$ parts by the $X$ and $Y$ axes. So I cannot take into account only one part and multiply by $4$. He said that I should evaluate the integral for all $4$ parts (or triangles).
I did that and as answer got $-\frac{1}{3}$, which may have the negative sign since I traversed in clockwise direction. But direction should not matter in surface integrals. I am also in doubt of my working process in this latter attempt. In this attempt, I used the last mentioned integral $(1)$ and changing its limits from the relations of the $4$ straight lines, evaluated the surface integrals.
Finally, I traversed the path in a clockwise manner and evaluated the original line integral as $-\frac{1}{3}$.

My question is: Was my original effort via Stokes' Theorem a correct approach? Or was my friend correct on how I should use the Stokes' Theorem?
 A: Denote
$$D_1=\{(x,y)|0\leq y\leq 1,y-1\leq x\leq 1-y,\},$$
$$D_2=\{(x,y)|-1\leq y\leq 0,-y-1\leq x\leq 1+y,\},$$
$$D=D_1+D_2,$$
by Green's Theorem,
$$\oint_L xydx+xy^2dy=\iint_D(y^2-x)d\sigma=\iint_Dy^2d\sigma-\iint_Dxd\sigma.$$
For $D$ is axisymmetric about $x$ and $y$, we got
$$\iint_Dxd\sigma=0,\quad \iint_Dy^2d\sigma=2\iint_{D_1}y^2d\sigma,$$
Then
$$\oint_L xydx+xy^2dy=2\iint_{D_1}y^2d\sigma=2\int_0^1\int_{y-1}^{1-y}y^2dxdy\\
=2\int_0^1y^2(2-2y)dy=2\left(\frac23-\frac24\right)=\frac13.$$
[hint]: If $D$ is axisymmetric about $x$ (or $y=0$), and $f(x,-y)=-f(x,y)$, then 
$$\iint_Df(x,y)d\sigma=0.$$
In fact, let $D=D_1+D_2$, where $D_1, D_2$ are the parts of $D$ with $y\geq 0$ and $y\leq 0$ respectively. Assume we can denote
$$D_1=\{(x,y)|a\leq x\leq b,0\leq y\leq y_1(x)\},$$
$$D_2=\{(x,y)|a\leq x\leq b,-y_1(x)\leq y\leq 0\},$$
then
$$\iint_{D_1}f(x,y)d\sigma=\int_a^b\int_0^{y_1(x)}f(x,y)dydx,$$
let $t=-y$, then
$$\iint_{D_1}f(x,y)d\sigma=\int_a^b\int_0^{y_1(x)}f(x,y)dydx\\
=\int_a^b\int_0^{-y_1(x)}f(x,-t)d(-t)dx\\
=-\int_a^b\int_0^{-y_1(x)}f(x,-t)dtdx\\
=\int_a^b\int_0^{-y_1(x)}f(x,t)dtdx\\
=-\int_a^b\int_{-y_1(x)}^0f(x,t)dtdx\\
=-\iint_{D_2}f(x,y)d\sigma,$$
hence we get
$$\iint_{D}f(x,y)d\sigma=\iint_{D_1}f(x,y)d\sigma+\iint_{D_2}f(x,y)d\sigma=0$$
