Double integral over a triangular domain

Setup.
$$\iint_D (x^2+y^2) \, dx \, dy = \int_0^1 \left( \int_0^{1-x} (x^2+y^2) \, dy \right) dx.$$

Integration in $y$.
$$\int_0^{1-x} (x^2+y^2) \, dy = x^2(1-x) + \frac{(1-x)^3}{3}.$$

Integration in $x$.
$$\int_0^1 \left( x^2(1-x) + \frac{(1-x)^3}{3} \right) dx = \left[ \frac{x^3}{3} - \frac{x^4}{4} - \frac{(1-x)^4}{12} \right]_0^1.$$
At $x=1$: $\frac13 - \frac14 - 0 = \frac{1}{12}$. At $x=0$: $0 - 0 - \frac{1}{12} = -\frac{1}{12}$. Hence
$$\iint_D (x^2+y^2) \, dx \, dy = \frac{1}{12} - \left(-\frac{1}{12}\right) = \frac16.$$

Remark. This is consistent with the area of the triangle ($1/2$) and the symmetry of the domain.