How to compute the integral of function with singularity?
5 months ago by
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*x+x*x,-1/8.)',degree=intDegree) fp=project(f,V) A=assemble(fp*dx) print(A)
Community: FEniCS Project
5 months ago by
Expressionworks 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
Expressionclass 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*x + x*x)**(-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*x + x*x)**(-1/8.)is evaluated at every quadrature point of order
intDegreeand 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.