How to compute the integral of function with singularity?


155
views
0
5 months ago by
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:
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
5 months ago by
The 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.

Similar posts:
Search »