# Newtonian Potential

**URL:** <https://fenicsproject.discourse.group/t/newtonian-potential/9176>\
**Category:** Site Feedback\
**Created:** [August 31, 2022, 2:06pm UTC](https://fenicsproject.discourse.group/t/newtonian-potential/9176 "2022-08-31T14:06:00Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![JanBohn](https://avatars.discourse-cdn.com/v4/letter/j/848f3c/32.png) [@JanBohn](https://fenicsproject.discourse.group/u/JanBohn)\
**Post date:** [August 31, 2022, 2:06pm UTC](https://fenicsproject.discourse.group/t/newtonian-potential/9176/1 "2022-08-31T14:06:00Z")

</div>

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:

```auto
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

---

<div class="post-metadata">

**Author:** ![PeculiarCarrot](https://avatars.discourse-cdn.com/v4/letter/p/4af34b/32.png) [@PeculiarCarrot](https://fenicsproject.discourse.group/u/PeculiarCarrot)\
**Post date:** [September 2, 2022, 4:32pm UTC](https://fenicsproject.discourse.group/t/newtonian-potential/9176/2 "2022-09-02T16:32:13Z")

</div>

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