# Newtonian Potential

Hello,

I would like to evaluate the Newtonian potential
N(u)(r) = \int_\Omega \frac{1}{|x-r|} u(x) dx
(not over the whole \mathbb{R}^3 but only a bounded domain \Omega).
This is the code I use:

from dolfin import *
import numpy as np

h=5
meshh=UnitCubeMesh(h, h, h)

Pr3 = VectorElement('Lagrange', meshh.ufl_cell(), 1, dim=3);
V3 = FunctionSpace(meshh, Pr3)
DG0 = FunctionSpace(meshh, 'DG', 0)

kernelexp = "pow(pow(x-r0,2)+pow(x-r1,2)+pow(x-r2,2),-0.5)"
kernel = Expression(kernelexp,degree=3,r0=0.0,r1=0.0,r2=0.0)

class MyExpressionM(UserExpression):
def eval(self, value, x):
x=x
s3=np.sin(2*np.pi*x)
c3=np.cos(2*np.pi*x)
value = s3
value = c3
value = 0.0

def value_shape(self):
return (3,)
maprfkt = interpolate(MyExpressionM(degree=3), V3)

res= assemble(maprfkt*kernel*dx)

print(res)


If I evaluate this on a grid node, e.g. r=(0,0,0) the result is nan because the
quadrature rule evaluates in the singularity. Is there any possibility to overcome this problem?
I tried to modify the kernel by \frac{1}{|x-r|+\epsilon} for small \epsilon>0, but then
the result depends heavily on \epsilon.