Here is my question
Given a point $P$ in space, and given a piece of malleable material of constant density, how should you shape and place the material in order to create the largest possible gravitational field at $P$?
Label the points on the surface by their distance $r$ from $P$, and by the angle $θ$ that the line of this distance subtends with the $x$-axis. Then a small mass $dm$ on the surface provides an $x$-component of the gravitational field equal to
$$F_{x} = \cos θ \frac {G\; dm}{r^2} $$
Any help is appreciated!
Before we start, let us rephrase the question to a way I understand:
The problem has a lot of symmetry, it is invariant under translation of $P$ and rotation of directional axis. WOLOG, we only need to consider the case $P$ is the origin and the coordinate axis is chosen so that the "largest" gravitation force is pointing towards the $z$ direction. i.e.
$$\text{maximize}\;\frac{F_z}{Gm\rho} = \int_{\Omega} \frac{z}{(x^2+y^2+z^2)^{3/2}} dx dy dz \quad\text{ subject to }\quad \verb/Volume/(\Omega) = V$$ Let $1_\Omega(x,y,z)$ be the indicator function for $\Omega$, it will be a solution to following problem
$$\text{maximize}\;\int f(x,y,z)\frac{z}{(x^2+y^2+z^2)^{3/2}}dxdydz \quad\text{ subject to }\quad \begin{cases}\int f(x,y,z) dxdydz = V\\ f(x,y,z) = 0 \text{ or } 1 \end{cases} $$
By a variant of Rearrangement inequality for integral, the LHS is maximized when $f$ has the form
$$f(x,y,z) = \begin{cases} 1, & \frac{z}{r^3} > \frac{1}{a^2}\\ 0, & \frac{z}{r^3} < \frac{1}{a^2} \end{cases} $$ for some suitable chosen $a$. In terms of polar coordinate $$[0,\infty) \times [0,\pi] \times [ 0,2\pi ] \ni (r,\theta,\phi) \quad\mapsto\quad (x,y,z) = (r\sin\theta\cos\phi,r\sin\theta\sin\phi,r\cos\theta) \in \mathbb{R}^3 $$ $\Omega$ is the region such that $$\frac{z}{(x^2+y^2+z^2)^{3/2}} = \frac{\cos\theta}{r^2} \ge \frac{1}{a^2}\quad\iff\quad r \le a\sqrt{\cos\theta}, \theta \in [0,\frac{\pi}{2}]$$ It is a solid of revolution which can be generated by rotating the dipole curve $(x^2+z^2)^3 = z^2$ in $xz$-plane with respect to $z$-axis and then keep the component bounded by the resulting surface above the $xy$-plane.
For any $a$, the volume of $\Omega$ equals to
$$\begin{align} V = \int_{\Omega} dxdydz &= \int_0^{2\pi}\int_0^{\pi/2}\int_0^{a\sqrt{\cos\theta}} r^2 \sin\theta dr d\theta d\phi\\ &= \frac{2\pi a^3}{3} \int_0^{\pi/2}(\cos\theta)^{3/2} \sin\theta d\theta = \frac{2\pi a^3}{3} \int_0^1 u^{3/2} du = \frac{4\pi a^3}{15} \end{align} $$ This allow us to fix $a$ as $\displaystyle\left(\frac{15 V}{4\pi}\right)^{1/3} = \left(\frac{15M}{4\pi\rho}\right)^{1/3}$. The maximum value of $|\vec{F}|$ satisfies $$\begin{align} \frac{|\vec{F}|_{max}}{Gm\rho} = \int_{\Omega} \frac{\cos\theta}{r^2} dxdydz &= \int_0^{2\pi}\int_0^{\pi/2}\int_0^{a\sqrt{\cos\theta}} \cos\theta \sin\theta dr d\theta d\phi\\ &= 2\pi a \int_0^{\pi/2}(\cos\theta)^{3/2}\sin\theta d\theta = 2\pi a \int_0^1 u^{3/2} du = \frac{4\pi a}{5}\\ &= \frac{3V}{a^2} \end{align} $$ This leads to $\displaystyle\;|\vec{F}|_{max} = \frac{3GMm}{a^2} = 3GMm\left(\frac{4\pi\rho}{15M}\right)^{2/3}$.