This is not an answer as such, and noch much progress is made. The argument is one I know almost nothing about, but I found the problem intriguing, and I hope someone can expand/rubbish this idea. Apologies if this is not ok.
Let $f \colon [0,1]^d \to \mathbb{R}$ be a continuous function. Approximate $|f|$ with the multivariate Bernstein polynomials (writing $f_{k_1\cdots k_d} = f\bigl(\tfrac{k_1}{n}, \ldots \tfrac{k_d}{n}\bigr)$ for short),
$$B_{|f|,n}(x) = \sum_{k_1 \cdots k_d} |f_{k_1\cdots k_d}| \prod_{i=1}^{d} \binom{n}{k_i} \, x_i^{k_i} (1-x_i)^{n-k_i} \;.$$
Let $S_{lm} = E[X^m (1-X)^{l-m}]$. The expectation of the approximant is then
$$E[B_{|f|,n}(X)] = \sum_{k_1 \cdots k_d} |f_{k_1\cdots k_d}| \prod_{i=1}^{d} \binom{n}{k_i} \, S_{nk_i} \;,$$
which should be maximal. The value of $S_{nk}$ is a linear combination on the first $n$ moments $M$ of the resulting distribution, and by a result of Hausdorff, there is indeed a unique distribution for a given set of moments iff all $S_{lm}$ (ie. not just for $n, k$) are positive. Expanding
$$S_{lm} = \sum_{r=0}^{l} \binom{l-m}{r-m} (-1)^{r-m} \, M_r \;,$$
the expectation in terms of the moments is
$$E[B_{|f|,n}(X)] = \sum_{k_1 \cdots k_d} \sum_{r_1 \cdots r_d} |f_{k_1\cdots k_d}| \prod_{i=1}^{d} \binom{n}{r_i} \binom{r_i}{k_i} (-1)^{r_i-k_i} \, M_{r_i} \;.$$
In the general case, the problem at this point seems to be the constrained multilinear optimisation
$$ \begin{aligned} \max_{M_1 \cdots M_n} \;&\; \sum_{k_1 \cdots k_d} \sum_{r_1 \cdots r_d} |f_{k_1\cdots k_d}| \prod_{i=1}^{d} \binom{n}{r_i} \binom{r_i}{k_i} (-1)^{r_i-k_i} \, M_{r_i} \;, \\ \text{s.t.} \;&\; \sum_{r=0}^{l} \binom{l-m}{r-m} (-1)^{r-m} \, M_r \ge 0 \;, \\ &\; M_r \in [0,1] \;, \quad M_0 = 1 \;, \end{aligned} $$
which appears generally ill-behaved from a cursory Google search (but I have no idea, really). Nevertheless, it is a recipe for calculations: For the baby problem, I find the moments (with $n = 12$) 0.5, 0.40042, 0.35063, 0.323633, 0.308033, 0.298489, 0.292331, 0.288158, 0.2852, 0.283021, 0.28136, 0.280058 using a quick Mathematica test run, and the numbers are relatively robust against changes in $n$ and match the distributions suggestions from the comments. It should be possible to fit $n$ point masses to create the distribution.
Perhaps a little more progress could be made in special cases, even though I am leaning very far outside my window now. Expanding the expectation further almost factors the expression,
$$ E[B_{|f|,n}(X)] = \sum_{{k_1 \cdots k_d}} \frac{|f_{k_1 \cdots k_d}| \, (-1)^{k_1 + \ldots + k_d}}{k_1! \cdots k_d!} \sum_{{r_1 \cdots r_d}} \left\{\prod_{i=1}^{d} \frac{n! \, (-1)^{r_i}}{(n-r_i)! \, (r_i-k_i)!} \, M_{r_i}\right\} \;. $$
If $f$ is symmetric, and the symbol $f_{k_1 \cdots k_d}$ is invariant under exchange of the indices, and it feels like something should give. Assuming that the remaining term could somehow be factored, the summation could be rewritten in multinomial terms,
$$ \begin{aligned} E[B_{|f|,n}(X)] &= \sum_{{k_1 \cdots k_d}} \frac{|f_{k_1 \cdots k_d}| \, (-1)^{k_1 + \ldots + k_d}}{k_1! \cdots k_d!} \sum_{|p| = d} \binom{d}{p_0 \,\cdots\, p_n} \prod_{i=0}^{n} \bigl(a_{i,k?} \, M_i\bigr)^{p_i} \\ &= \sum_{{k_1 \cdots k_d}} \frac{|f_{k_1 \cdots k_d}| \, (-1)^{k_1 + \ldots + k_d}}{k_1! \cdots k_d!} \left(a_{0,k_?} \, M_0 + \ldots + a_{n, k?} \, M_n\right)^d \;, \end{aligned} $$ which looks like it might lead somewhere nice.
Update: Fitting a polynomial distribution
It's quite natural to try and approximate the recovered distribution $g$ itself as a Bernstein polynomial of order $n$,
$$B_{g,n}(x) = \sum_{k=0}^{n} \binom{n}{k} \, g_k \, x^k (1-x)^{n-k} \;,$$
with the function values $g_k = g\bigl(\frac{k}{n}\bigr)$ the unknown parameters of the fit. The approximate moments are
$$ \begin{aligned} \tilde{M}_r &= \sum_{k=0}^{n} \binom{n}{k} \, g_k \int_0^1 \! x^{k+r} (1-x)^{n-k} \, dx \\ &= \sum_{k=0}^{n} \frac{n! \, (k+r)!}{k! \, (n+r+1)!} \, g_k \;, \end{aligned} $$
This is readily inverted to give the function values $g_k$. Here is a sample plot for $n = 10$, where you can see that Runge's phenomenon is quite pronounced for such low values of $n$.
Update 2: It actually does not matter if $f$ is symmetric or not.
One can permute the indices $k_1 \cdots k_d$ without changing the product in the expectation, so defining the symmetric version of $|f_{k_1 \cdots k_d}|$,
$$|f^S_{k_1 \cdots k_d}| = \frac{1}{d!} \sum_{\sigma \in S_d} |f_{k_{\sigma(1)} \cdots k_{\sigma(d)}}| \;,$$
the expectation can be rewritten as
$$ E[B_{|f|,n}(X)] = \sum_{k_1 \cdots k_d} \sum_{r_1 \cdots r_d} |f^S_{k_1 \cdots k_d}| \prod_{i=1}^{d} \binom{n}{r_i} \binom{r_i}{k_i} (-1)^{r_i-k_i} \, M_{r_i} \;. $$