I was thinking if we have Dirac distribution as boundary condition, then what will be regularity of solution. Problem is following, $$ \left\{\begin{matrix} \nabla_x(\gamma(x)\nabla_x u(x,y))=0 & in~\Omega\\ u(x,y)=\delta_y(x) & on ~\partial \Omega \end{matrix}\right. $$ where $x,y\in \mathbb R^n , \gamma(x)>c>0 $. $\Omega$ is with smooth boundary.
Any help or reference is greatly appreciated.
For a regular enough boundary function $\psi$ the solution of the first BVP can be written as $$ u(x)=\int_{\partial \Omega}\frac{\partial G(x,y)}{\partial \bar n_y}\gamma(y)\psi(y)\,dy, $$ where $G(x,y)$ is the Green function of the corresponding first BVP and $\bar n_y$ is a unit normal vector to $\partial\Omega$. It follows from the Green's third identity written out for your operator.
Substituting $\psi(x)=\delta_y(x)$ gives $u(x)=\gamma(y)\frac{\partial G(x,y)}{\partial n_y}$. In some cases this function is also called the Poisson kernel. Estimates for the Green function are known (for smooth enough coefficients and bounded domains they are the same as for a fundamental solution) and it follows that $$ |u(x,y)|\le C|x-y|^{1-n} \quad x\in \Omega,\ y\in \partial \Omega. $$