Combining multiple multivariate integrals into a single function Let $p$ be a multivariate polynomial (with property that for any $x\in[0, 1]^n$, $p(x) \in [0,1]$). If I want to find a $g$ function such that $\frac{\partial g}{\partial x_1} = p(x)$ then clearly $g=\int p(x) dx_1$. 
Q1: Suppose, instead I have two direction vectors $u$ and $v$, and I want $(v \cdot \nabla g(x)) = p(x)$ and $(u \cdot \nabla g(x)) = p(x)$, then what should g be? I can do this for one directional derivative, but not sure how to combine for multiple directional derivatives. 
Q2: Can I go one step further, where given $q:[0,1]^n \rightarrow [0, 1]^n$, define function $g$ such that $((q(x) - x)\cdot \nabla g(x))=p(x)$? Here each $q_i$ is a polynomial again.  
 A: If you have the equation
$$
u \cdot \nabla g(x) = p(x),
$$
this can be written (as you said) as
$$
D_u g(x) = p(x),
$$
where $D_u$ is the directional derivative. This suggests that we go on that direction and see what happens. Let $x_1 = u_1 \eta$, $x_2 = u_2 \eta$ and $x_3 = u_3 \eta$, then,
$$
\frac{d g}{d\eta} = u \cdot \nabla g = p(\eta)
$$
wich is, sort of, what you require. Of course, there is more to it, as what you'r trying to solve is an uncoupled system of PDEs (one for $u$ and one for $v$), but you're lacking a second function (i.e. $g_1$ and $g_2$), so, there are going to be very strong compatibility conditions, which will condition the existence of the function $g$.
For the PDE
$$
(q(x) - x) \cdot \nabla g(x) = p(x),
$$
if $g(x)$ is a solution, then the vector $\begin{pmatrix}\nabla g\\-1\end{pmatrix}$ is normal to the surface generated by function $g$. This means that
$$
\begin{pmatrix}q(x) - x \\ p(x)\end{pmatrix} \cdot \begin{pmatrix}\nabla g\\-1\end{pmatrix} = 0.
$$
In other words, the vector $\begin{pmatrix}q(x) - x \\ p(x)\end{pmatrix}$ is tangent to it. So, we want to find integral surfaces in $\mathbb{R}^{n+1}$ such that their tangents are $\begin{pmatrix}q(x) - x \\ p(x)\end{pmatrix}$ and pass trough an initial surface 
$$
\begin{pmatrix}x_1(\xi_1,\ldots,\xi_{n-1})\\\vdots\\x_n(\xi_1,\ldots,\xi_{n-1})\\f(\xi_1,\ldots,\xi_{n-1})\end{pmatrix}.$$
This means that
$$
\begin{align}
\frac{d x_1}{d \eta} &= q_1(x_1,\ldots,x_n) - x_1, \\
&\qquad\qquad\vdots\\
\frac{d x_n}{d \eta} &= q_n(x_1,\ldots,x_n) - x_n, \\
\frac{d g}{d \eta} &= p(x_1,\ldots,x_n),
\end{align}
$$
where
$$
\begin{align}
x_1\big|_{\eta = 0} &= x_1(\xi_1,\ldots,\xi_{n-1}),\\
&\,\,\vdots\\
x_n\big|_{\eta = 0} &= x_n(\xi_1,\ldots,\xi_{n-1}),\\
g\big|_{\eta = 0} &= f(\xi_1,\ldots,\xi_{n-1}).
\end{align}
$$
In order to have a solution, you need:


*

*That the system of ODEs is solvable (existence).

*That the transformation $(x_1,\ldots,x_n) \to (\xi_1,\ldots,\xi_{n-1},\eta)$ is invertible (this condition ensures that you actually have a surface [i.e. its normal is well defined] and that you can go from the parametric surface representation of the solution to the function $g = g(x_1,\ldots,x_n)$).


Is pretty straight forward to prove that if 1. and 2. are satisfied, then $g$ is the solution (modulo uniqueness) of the PDE.
You can check this for $q(x) - x = u$.
