In the comments, Tom Copeland asked if, six years later, I had any further insights into the diagrammatics I proposed near the end of my question. So I figured I'd mention what I know, which is yet another (marginally) diagrammatic description of BCH multiplication and in particular the map $mathcal U mathfrak g to mathcal S mathfrak g$. What I'll describe is some part of http://dx.doi.org/10.1090/pspum/088 (also available at http://arxiv.org/abs/1307.5812). Some pictures are available in the last chapter of my thesis, but note that that chapter has a subtle error (the details of which I haven't yet sussed out) which I have corrected in the linked paper; the error does not affect the case of (duals of) Lie algebras, but does affect Poisson structures whose Taylor expansions include quadratic or higher terms.
I will describe some homological algebra, and then unpack into diagrams. Let $mathfrak G = mathfrak g otimes Omega_{mathrm{cpt}}(mathbb R)[1]$ denote the cochain complex of $mathfrak g$-valued compactly supported de Rham forms on $mathbb R$, shifted so that its cohomology is concentrated in degree $0$. Then $mathrm H^0(mathfrak G) = mathfrak g$, and the projection to homology is given by integrating de Rham forms — there is a cohomological-degree-$0$ map $int : mathfrak G to mathfrak g$. Note also that $Omega_{mathrm{cpt}}(mathbb R)$ has a non-unital graded-commutative multiplication ($wedge$), through which $mathfrak G$ picks up a graded-symmetric (!) cohomological-degree-$(+1)$ map $delta : mathfrak G otimes mathfrak G to mathfrak G$ by $delta(xotimes alpha,yotimes beta) = [x,y] otimes (alpha wedge beta)$ (up to a sign that depends on a choice of conventions about how to handle elements of shifted complexes). Consider the (graded-commutative) symmetric algebra $mathcal S mathfrak G$. I will denote its usual differential (which is the extension of de Rham-form differentiation as a derivation) by $partial_{mathrm{dR}}$.
Integration of de Rham forms extends to an algebra homomorphism $int : (mathcal Smathfrak G, partial_{mathrm{dR}}) to mathcal S mathfrak g$.
The map $delta$ has a canonical extension to a second-order differential operator on $mathcal S mathfrak G$ which I will also call $delta$. (A second order differential operator on a symmetric algebra is unique determined by its values on constant, linear, and quadratic terms. We declare that $delta$ vanishes on constants and linears, and set it to be the original $delta$ on quadratics.) Because $wedge$ plays well with $partial_{mathrm{dR}}$, $delta$ and $partial_{mathrm{dR}}$ graded-commute. The Jacobi identity implies $delta^2 = 0$. So we have a new differential $partial_q = partial_{mathrm{dR}} + delta$ on $mathcal S mathfrak G$.
Pick any $alpha in Omega_{mathrm{cpt}}^1(mathbb R)$ with $intalpha = 1$. (Henceforth I will call such one-forms bumps.) The map $x mapsto x otimes alpha$ from $mathfrak g to mathfrak G$ extends to an algebra homomorphism $mathcal S mathfrak g to mathcal S mathfrak G$ splitting $int$. There is a unique contracting homotopy $eta_alpha : Omega_{mathrm{cpt}}(mathbb R) to Omega_{mathrm{cpt}}(mathbb R)$ (of cohomological degree $-1$) such that the graded commutator $[partial_{mathrm{dR}},eta_{alpha}] = mathrm{id} - (alpha otimes) circ int$; it vanishes on $Omega_{mathrm{cpt}}^0$ and on $Omega_{mathrm{cpt}}^1$ satisfies $eta_alpha(beta) = partial_{mathrm{dR}}^{-1}(beta - (int beta)alpha)$; note that $int(beta - (int beta)alpha) = 0$, so the one-form $beta - (int beta)alpha$ has a unique antiderivative among compactly-supported functions.
We can extend $eta_alpha$ to $mathcal Smathfrak G$ in many ways, and the choice can be proven not to matter. To make a choice, declare that on the constants $mathcal S^0mathfrak G$ we have $eta_alpha = 0$, and on $mathcal S^nmathfrak G$ we have
$$ eta_alpha(beta_1 odot dots odot beta_n) = frac1n sum_i beta_1 odot dots odot eta_alpha(beta_i) odot dots odot beta_n, $$
where $odot$ denotes the symmetric multiplication in $mathcal S$. Note that this is not the extension of $eta_alpha$ as a derivation.
Since $eta_alpha$ always attaches 1-forms and $delta$ involves wedge multiplication, $eta_alpha : mathcal S mathfrak g to (mathcal S mathfrak G,partial_q)$ is a chain map. The integration map $int : mathcal S mathfrak G to mathcal S mathfrak g$ splitting this map is not a chain map from $(mathcal S mathfrak G,partial_q)$. But with the above choices we can choose a different splitting, namely $int circ (mathrm{id} - delta eta_alpha)^{-1}$. (I have a 50% chance of getting that minus sign wrong.) I will leave checking that this is a chain map splitting $eta_alpha$ to you. Note that $(mathrm{id} - delta eta_alpha)^{-1} = sum_N (delta eta_alpha)^N$ converges on $mathcal S mathfrak G$, since $delta eta_alpha$ drops polynomial degree by $1$.
Pick bumps $alpha_1, dots,alpha_n$ such that the support of $alpha_i$ is in $[i-1,i]$, and pick one final bump $alpha$ arbitrarily. One can prove that the map $mathcal U mathfrak g to mathcal S mathfrak g$ is given on monomials $x_1 dots x_n$ (with multiplication in $mathcal U$) by
$$ mathcal U mathfrak g ni x_1 dots x_n mapsto int circ (mathrm{id} - delta eta_alpha)^{-1} bigl( (alpha_1 otimes x_1) odot dots odot (alpha_n otimes x_n) bigr) in mathcal S mathfrak g$$
(or I might be off by a sign somewhere). In general, similar formulas describe the entire product on $mathcal Smathfrak g$ given by transporting the one from $mathcal U mathfrak g$ along the symmetrization isomorphism.
Let me now unpack this formula, or rather give the answer after some unpacking. (Proving that this is a valid unpacking is straightforward: you need to track the numerical factors coming from $eta_alpha$, understand how to apply a second-order differential operator to a monomial, and also include a brief "degree reasons" argument to get $eta_alpha$ and $delta$ to apply always to the same things at the same time.)
Define an $n$-leaf binary heighted forest, abbreviated forest, to be set of binary rooted trees whose leaves are put in bijection with the set ${1,dots,n}$ and whose nodes are totally ordered (I mean: totally order the set of all nodes) such that in a given tree, and path from root to leaf is increasing for the total ordering. Arbitrarily choose for each node which of its two branches is left and which right (the choice will cancel out).
Given a forest and the list $x_1,dots,x_n$ of elements in $mathfrak g$, there is an obvious element of $mathcal S mathfrak g$ given by putting the $x_i$s at the leaves and reading the forest as instructions of who to bracket with whom (then multiply the "root" outputs).
Now I will describe, for each forest, how to compute a number. Consider the map $Omega_{mathrm{cpt}}^1(mathbb R) otimes Omega_{mathrm{cpt}}^1(mathbb R) to Omega_{mathrm{cpt}}^1(mathbb R)$ given by $beta_1 otimes beta_2 mapsto beta_1 wedge eta_alpha (beta_2) - beta_2 wedge eta_alpha(beta_1) $. Place this map at each vertex, and $alpha_i$ at the $i$th leaf, and let the forest tell you how to apply this map to end up with a bunch of $1$-forms at the roots. Then integrate all these $1$-forms to get numbers, and multiply those numbers together.
Finally, suppose there are $k$ roots (and hence $(n-k)$ nodes). Then multiply by $frac1 n frac1{n-1} dots frac1{n-k}$.
Note that neither the number nor the element of $mathcal S mathfrak g$ determined by a forest depends on the height ordering. But I now want you to sum over all forests with total node-ordering of the product of these two numbers. That sum computes the map $mathcal U mathfrak g to mathcal S mathfrak g$ above. (If I had a good way to count the number of total orderings of the nodes for a given un-heighted forest, I would have used it.)
This is similar to, but not the same as, Kontsevich's star product. In particular, note that my forests have no wheels, whereas Kontsevich does not describe the symmetrization map $mathcal S mathfrak g cong mathcal U mathfrak g$, but rather this map twisted by some traces in the adjoint representation.