We already know the relationship between Green's function and solution to elliptic partial differential equation, i.e $$u(y)=\int\_{\partial \Omega}u\frac{\partial G}{\partial n} ds+\int\_\Omega G\Delta u dx $$ where $n$ is the unit outward normal , $G$ is Green's function on $\Omega$, and $p=\frac{\partial G}{\partial n}$ is called the Poisson kernel. 


So now I'm wondering the case if we have a partial differential operator $$Lu(t,x)=\partial\_tu(t,x)+\sum\_{i,j}^da\_{i,j}(t,x)\frac{\partial^2u(t,x)}{\partial x\_i\partial x\_j}$$ 
I know that in this case we call $G$ the fundamental solution if $$LG(s,y;t,x)=\delta(s-t)\delta(x-y)$$ 
Then we can find that $$u(t,x)=\int G(s,y;t,x)Lu(s,y)dsdy$$ without boundary conditions, if my understanding is correct. 


Are there the similar relationship between fundamental solution and Poisson kernel (if this notion is correct) on the boundary? In one paper, $$L=\frac{\partial}{\partial t}+\operatorname{div}\big(A(x,t)\nabla\_x\big)$$ where $A$ is a $d\times d$ real symmetric matrix, uniformly elliptic on $\Omega=D\times(0,\infty)$. They have the corresponding Poisson kernel $$p(x,t;y,s)=\frac{\partial G}{\partial N(y,s)}(x,t;y,s)$$ where $N(y,s)=A(y,s)n(y)$ with $n(y)$ is the unit inner normal to $\partial D$ at $y$. It confused me for several days. I really want to figure out how it was calculated. 



Green's functions/fundamental solution to a non-constant coefficients pde
First, unless I misinterpret what you wrote, your integral representation for $u(t,x)$ is not correct. At the very least, it's missing boundary terms. Without boundary terms the formula obviously fails for any $u(t,x)\ne 0$ with $Lu(t,x)=0$. Whatever the right boundary term is, that will give you the $L$-analog of the Poisson kernel.


Here's how to find it. First, you need the formal adjoint $L^\*$ of the operator $L$. Everybody knows how to find the formal adjoint, by using integration by parts. But what that really implies is that there exists a pair of bidifferential operators $V[-,-]$ and $W[-,-]$, linear in each argument such that
$$
 v L[u] - L^\*[v] u = \partial\_t V[v,u] + \operatorname{div}\_x W[v,u] .
$$
Basically, $V$ and $W$ collect all the boundary terms that you get by moving all the derivatives from $u$ to $v$ using integration by parts. If that's confusing, see this recent [answer](https://mathoverflow.net/a/289070) for some explanation. In your case,
\begin{align\*}
 L&=\frac{\partial}{\partial t}+\operatorname{div}\_x\big(A(x,t)\nabla\_x\big) , \\
 L^\*&=-\frac{\partial}{\partial t}+\operatorname{div}\_x\big(A(x,t)\nabla\_x\big) , \\
 V[v,u] &= v u , \\
 W[v,u] &= v A(t,x) \nabla\_x u - u A(t,x) \nabla\_x v .
\end{align\*}


Next, you actually need $G(s,y;t,x)$ to satisfy
$$
 L\_{s,y}^\*[G](s,y;t,x) = \delta(s-t) \delta(y-x) .
$$
Then, for $(t,x) \in \Omega = [0,T] \times D$, and assuming that $G$ and $u$ are sufficiently regular for all the differentiations to be well defined (perhaps in a weak sense),
\begin{align}
 u(t,x)
 &= \int\_\Omega \delta(s-t) \delta(y-x) u(s,y) \, ds dy \\
 &= \int\_\Omega L\_{s,y}^\*[G](s,y; t,x) u(s,y) \, ds dy \\
 &= \int\_\Omega \left( \partial\_s V[u,v] + \operatorname{div}\_y W[u,v] \right) \, ds dy
 + \int\_\Omega G(s,y; t,x) L\_{s,y}[u](s,y) \, ds dy \\
 &= \left. \int\_D V[G,u] \right|\_{s=0}^T \, dy
 + \int\_{\partial D}\left( \int\_0^T n(y)\cdot W[G,u] \, ds \right) dy\_{\|}
 + \int\_\Omega G(s,y; t,x) L\_{s,y}[u](s,y) \, ds dy ,
\end{align}
where $n(y)$ is the outer normal vector to the boundary $\partial D$ and $dy\_{\|}$ is the surface area element of $\partial D$.


So, in your case, the final integral representation of $u(t,x)$ inside $\Omega$ is
\begin{multline\*}
 u(t,x) = \left. \int\_D G(s,y; t,x) u(s,y) \right|\_{s=0}^T \, dy \\
 + \int\_{\partial D} \int\_0^T \left( G(s,y; t,x) \, n(y)\cdot A(s,y)\nabla\_y u(s,y) - u(s,y)\, n(y)\cdot A(s,y) \nabla\_y G(s,y; t,x) \right) \, ds dy\_{\|} \\
 + \int\_\Omega G(s,y; t,x) L[u](s,y) \, ds dy .
\end{multline\*}
Depending on what properties you require of $u$ and $G$ in the interior and on the boundary of $\Omega$, some of the terms on the right hand side of that identity will vanish. For instance, if $L[u] = 0$, then you are left with a boundary integral representation of the form
$$
 u(t,x) = \int\_{\partial\Omega} p(s,y; t,x) \, (ds dy)\_{\|} ,
$$
where again the $\|$ subscript denotes the appropriate surface area element on $\partial \Omega$. I guess $p(s,y; t,x)$ would be your $L$-analog of a Poisson kernel, which you can read off from the preceding identity.


