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[0]-r0,2)+pow(x[1]-r1,2)+pow(x[2]-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[2]
s3=np.sin(2*np.pi*x)
c3=np.cos(2*np.pi*x)
value[0] = s3
value[1] = c3
value[2] = 0.0
def value_shape(self):
return (3,)
maprfkt = interpolate(MyExpressionM(degree=3), V3)
res= assemble(maprfkt[0]*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.
Thanks in advance,
Jan