Calculate $\iiint_\omega (x+y+z)^2\,dxdydz$. Problem: Calculate $$\displaystyle\iiint_\Omega (x+y+z)^2\,dxdydz,$$ where $\Omega$ is the domain bounded by the sphere $x^2+y^2+z^2\le 3a^2$ and the paraboloid $x^2+y^2\le 2az$, where $z\ge 0$.
Progress: So, by sketching the domain, one can see that the domain includes a part from sphere and a part from the paraboloid, so we need to calculate the integral on two seperate part, then add the to get the result.
To calculate the integral in the sphere part, I've tried using spherical coordinate. Using spherical coordinate, I made the following substitution:
$$\begin{cases} x = r\sin\theta\cos\varphi\\ y = r\sin\theta\sin\varphi\\ z = r\cos\theta \end{cases}$$ with $0\le \varphi\le 2\pi$, $0\le\theta\le\alpha$, where $\alpha =\arccos\left(\dfrac{1}{\sqrt{3}}\right)$ (one can get the bound for angle $\theta$ by working out the intersection of the paraboloid and the sphere) and $\dfrac{a}{\cos\theta}\le r\le\dfrac{a\sqrt 2}{\sin\theta}$. Now, the expression $(x+y+z)^2$ become
$$r^2(1+\sin^2\theta\sin 2\varphi + \sin 2\theta+\sin\varphi+\sin 2\theta\cos\varphi).$$
Now, we multiply the Jacobian determinant and the original integral become
$$\int^{2\pi}_0d\varphi\int^\alpha_0d\theta\int^\tfrac{a\sqrt 2}{\sin\theta}_\tfrac{a}{\cos\theta}r^4\sin^2\theta(1+\sin^2\theta\sin 2\varphi + \sin 2\theta+\sin\varphi+\sin 2\theta\cos\varphi)dr.$$
As you can see, this one is tedious and almost impossible to calculate it. Do you have any idea to solve this problem?
 A: Using cylindrical coordinates,
$$x=\rho\cos\theta,y=\rho\sin\theta,z=z$$
We know the parabaloid intersects the sphere at $z=a$ and thus the integral becomes
$$\int_0^a\int_0^\sqrt{2az}\int_0^{2\pi}(\rho\sin\theta+\rho\cos\theta+z)^2\rho\,\mathrm d\theta\, \mathrm d\rho\, \mathrm dz\\+
\int_a^{\sqrt{3}a}\int_0^\sqrt{3a^2-z^2}\int_0^{2\pi}(\rho\sin\theta+\rho\cos\theta+z)^2\rho\,\mathrm d\theta\, \mathrm d\rho\, \mathrm dz$$ 
where the first term is the integral up to the intersection of the sphere and the parabaloid, and the second term is the rest of the integral. Even though this looks messy, a lot of the terms immediately die when you integrate wrt $\theta$. 
A: Note that $(x+y+z)^2 = x^2+y^2+z^2 +2(xy+yz+zx)$ and, due to the symmetry of the integration region, the integral simplifies to
$$I=\iiint_\Omega (x+y+z)^2\,dxdydz
= \iiint_\Omega (x^2+y^2+z^2)\,dxdydz
$$
Then, in cylindrical coordinates, the integration region is the circular area of radius $r=\sqrt 2 a$ and the integral is given by
$$I= \int_0^{2\pi}\int_0^{\sqrt2 a}\int_{\frac{r^2}{2a}}^{\sqrt{3a^2-r^2}} (r^2+z^2) r\ dzdr d\theta
=\frac{2\pi a^5}5\left(9\sqrt3-\frac{97}{12}\right)
$$
