Newtonian Potential

Site Feedback
1

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

2

Hi Jan,

My thoughts on this: either you remove a small sphere around the origin from your domain so that your calculation won’t evaluate it there, or you add a min function to cap the expression at a certain value.

Best,
Carrot