### How to compute the integral of function with singularity?

155

views

0

I tried to solve a problem where one of the function is L2 function but has a singularity at the origin. When I tried to use assemble to find the integration, the result is nan. Why does that happen? I thought all the integradtion was done by using quadrature points and weight, which should not cause this problem. Is there any way to solve it?

Here is a simple example：

Here is a simple example：

```
from __future__ import print_function
from fenics import *
from dolfin import *
import numpy as np
from numpy import *
N=4*2**0
intDegree = 2
mesh = RectangleMesh(Point(-1., -1.), Point(1., 1.), N, N)
V = FunctionSpace(mesh, 'CG', intDegree)
n = FacetNormal(mesh)
f = Expression('pow(x[0]*x[0]+x[1]*x[1],-1/8.)',degree=intDegree)
fp=project(f,V)
A=assemble(fp*dx)
print(A)
```

Community: FEniCS Project

### 1 Answer

4

The

When you try to assemble

Note, that there always will be a degree of freedom (for CG spaces) at point (0, 0) for your mesh, that is why you get nan=1/0.0. In your case you can bypass the

By the virtue of the process, you do not need any function space. Function

But be careful. You might introduce huge error when integrating such singularity via quadrature rule. Very nice approach is to subtract asymptotically same function which you can then integrate by hand, see https://www.johndcook.com/blog/2012/02/21/care-and-treatment-of-singularities/

`Expression`

works in a bit less intuitive way. When you try to assemble

`Expression(..., degree=3)`

it evaluates the C++ code at degrees of freedom of CG3 space! and then interpolates locally to get values at quadrature points which are used for assembly.Note, that there always will be a degree of freedom (for CG spaces) at point (0, 0) for your mesh, that is why you get nan=1/0.0. In your case you can bypass the

`Expression`

class and prepare purely UFL expression, e.g.```
from __future__ import print_function
import dolfin as d
import numpy as np
from numpy import *
N=4*2**0
intDegree = 10
mesh = d.RectangleMesh(d.Point(-1., -1.), d.Point(1., 1.), N, N)
# Get symbolic UFL representation of spatial coordinate
x = d.SpatialCoordinate(mesh)
# Prepare the form based solely on this symbolic repre,
# no CPP code
f0 = (x[0]*x[0] + x[1]*x[1])**(-1/8.)
A = d.assemble(f0*d.dx(degree=intDegree))
print(A)
```

By the virtue of the process, you do not need any function space. Function

`(x[0]*x[0] + x[1]*x[1])**(-1/8.)`

is evaluated at every quadrature point of order `intDegree`

and so is integral computed.But be careful. You might introduce huge error when integrating such singularity via quadrature rule. Very nice approach is to subtract asymptotically same function which you can then integrate by hand, see https://www.johndcook.com/blog/2012/02/21/care-and-treatment-of-singularities/

Please login to add an answer/comment or follow this question.