This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Question on quadrature method for singular integrands

+4 votes

Dear All

For my current research project the potential has Coulomb singularities of the
type 1/r near the nuclei. To make these integrals converge in principle I will make
use of a coordinate transformation from x,y,z to sqrt(x),sqrt(y) and sqrt(z)
that results in a custom integration rule as discussed in the
previous questions.
However, as I understand the standard method of evaluating
matrix elements the potentials or other coefficient functions are
assumed to continuous on the whole element and the accuracy
can be controlled only via specifying a degree parameter for the expression.
As discussed in the first paragraph above this won't work for me
Is there any other possibility in fenics at this point.
This is indeed very important for me.
In two dimensions the volume element 2 pi rho d rho dz
makes the integrals still converge , however in three dimensions
this is not the case.
Thanks for all you help
I really do appreciate it !

regards

Moritz

asked Jan 30, 2014 by moritzbraun FEniCS User (1,390 points)

Gauss quadrature can handle $1/r$ as long as there is no quadrature point at $r=0$ which is often the case. Nevertheless Gauss quadrature can't be exact for non-polynomials and the degree of $p/q$ for $p,q$ polynomials is estimated by UFL as $deg(p)+deg(q)$ which is heuristic. So you may need to increase quadrature degree manually.

commented Jan 30, 2014 by Jan Blechta FEniCS Expert (51,420 points)

Dear Jan

I am still very confused.
It would help if you could show me with an example
code, how to instruct UFL and FFC to evaluate the potential at the grid points instead of
interpolating it.

regards

Moritz

commented Jan 30, 2014 by moritzbraun FEniCS User (1,390 points)

I'm confused also - what do you mean by potential and grid points.

Using Gauss quadrature, integrand is never evaluated at vertices and facets. So if you have singularity at vertex (which is the case for 1/r term and vertex at [0, 0, 0]), integrand is not evaluated there and you don't obtain infinities in your tensors.

commented Jan 30, 2014 by Jan Blechta FEniCS Expert (51,420 points)

Dear Jan

I am still confused as well.

I hope that the following code snippet helps us to make progress:

import numpy as np
from dolfin import *
N=2
mesh=UnitCubeMesh(1,1,1)
V=FunctionSpace(mesh,"Lagrange",2)
coords=mesh.coordinates()
print coords
ListOfPoints=[]
class Pot(Expression):
def eval(self,values,x):
r=sqrt(x[0]x[0]+x[1]x[1]+x[2]x[2])
ListOfPoints.append(list(x))
values[0]=-2.0/(r+1e-6)
pot=Pot()
u=TrialFunction(V)
v=TestFunction(V)
P=upotvdx
PMAT=assemble(P)
print "NPt=",len(ListOfPoints)
ListOfPoints=np.array(ListOfPoints)
print ListOfPoints
VMAT=PMAT.array()
N=VMAT.shape[0]
VMAT.shape=N*N
print max(VMAT)

Sorry that the code is not properly indented ; I do not know how to fix it!
My output for this program is

[[ 0. 0. 0.]
[ 1. 0. 0.]
[ 0. 1. 0.]
[ 1. 1. 0.]
[ 0. 0. 1.]
[ 1. 0. 1.]
[ 0. 1. 1.]
[ 1. 1. 1.]]
NPt= 60
[[ 0. 0. 0. ]
[ 1. 0. 0. ]
[ 1. 1. 0. ]
[ 1. 1. 1. ]
[ 1. 1. 0.5]
[ 1. 0.5 0.5]
[ 1. 0.5 0. ]
[ 0.5 0.5 0.5]
[ 0.5 0.5 0. ]
[ 0.5 0. 0. ]
[ 0. 0. 0. ]
[ 1. 0. 0. ]
[ 1. 0. 1. ]
[ 1. 1. 1. ]
[ 1. 0.5 1. ]
[ 1. 0.5 0.5]
[ 1. 0. 0.5]
[ 0.5 0.5 0.5]
[ 0.5 0. 0.5]
[ 0.5 0. 0. ]
[ 0. 0. 0. ]
[ 0. 0. 1. ]
[ 1. 0. 1. ]
[ 1. 1. 1. ]
[ 1. 0.5 1. ]
[ 0.5 0.5 1. ]
[ 0.5 0. 1. ]
[ 0.5 0.5 0.5]
[ 0.5 0. 0.5]
[ 0. 0. 0.5]
[ 0. 0. 0. ]
[ 0. 1. 0. ]
[ 1. 1. 0. ]
[ 1. 1. 1. ]
[ 1. 1. 0.5]
[ 0.5 1. 0.5]
[ 0.5 1. 0. ]
[ 0.5 0.5 0.5]
[ 0.5 0.5 0. ]
[ 0. 0.5 0. ]
[ 0. 0. 0. ]
[ 0. 0. 1. ]
[ 0. 1. 1. ]
[ 1. 1. 1. ]
[ 0.5 1. 1. ]
[ 0.5 0.5 1. ]
[ 0. 0.5 1. ]
[ 0.5 0.5 0.5]
[ 0. 0.5 0.5]
[ 0. 0. 0.5]
[ 0. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 1. 1. ]
[ 1. 1. 1. ]
[ 0.5 1. 1. ]
[ 0.5 1. 0.5]
[ 0. 1. 0.5]
[ 0.5 0.5 0.5]
[ 0. 0.5 0.5]
[ 0. 0.5 0. ]]
6349.00815294

This clearly shows, that the Expression pot is evaluated at [0,0,0]

in contract to what you said yesterday

I will appreciate any hints on how to change the settings in FEniCS
to prevent that!

regards

Moritz

commented Jan 31, 2014 by moritzbraun FEniCS User (1,390 points)

Discussing a code is definitely a good idea. Nevertheless you must format the code properly. I'm not able to use it this way either because of indentation and asterisks. Read http://fenicsproject.org/qa/50/does-one-write-code-with-syntax-highlighting-fenics-q%26a-forum

commented Jan 31, 2014 by Jan Blechta FEniCS Expert (51,420 points)

Dear Jan
Dear Colleagues

Thanks for the hint
I did not know where to find the guideslines!

regards

M Braun

commented Jan 31, 2014 by moritzbraun FEniCS User (1,390 points)

Dear Jan
Dear Friends

Here we go with the proper formatting.

import numpy as np
from dolfin import *
N=2
mesh=UnitCubeMesh(1,1,1)
V=FunctionSpace(mesh,"Lagrange",2)
coords=mesh.coordinates()
print coords
ListOfPoints=[]
class Pot(Expression):
    def eval(self,values,x):
        r=sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])
        ListOfPoints.append(list(x))
        values[0]=-2.0/(r+1e-6)
pot=Pot()
u=TrialFunction(V)
v=TestFunction(V)
P=u*pot*v*dx
PMAT=assemble(P)
print "NPT=",len(ListOfPoints)
ListOfPoints=np.array(ListOfPoints)
print ListOfPoints

The output for this code is:

[[ 0. 0. 0.]
[ 1. 0. 0.]
[ 0. 1. 0.]
[ 1. 1. 0.]
[ 0. 0. 1.]
[ 1. 0. 1.]
[ 0. 1. 1.]
[ 1. 1. 1.]]
NPT= 60
[[ 0. 0. 0. ]
[ 1. 0. 0. ]
[ 1. 1. 0. ]
[ 1. 1. 1. ]
[ 1. 1. 0.5]
[ 1. 0.5 0.5]
[ 1. 0.5 0. ]
[ 0.5 0.5 0.5]
[ 0.5 0.5 0. ]
[ 0.5 0. 0. ]
[ 0. 0. 0. ]
[ 1. 0. 0. ]
[ 1. 0. 1. ]
[ 1. 1. 1. ]
[ 1. 0.5 1. ]
[ 1. 0.5 0.5]
[ 1. 0. 0.5]
[ 0.5 0.5 0.5]
[ 0.5 0. 0.5]
[ 0.5 0. 0. ]
[ 0. 0. 0. ]
[ 0. 0. 1. ]
[ 1. 0. 1. ]
[ 1. 1. 1. ]
[ 1. 0.5 1. ]
[ 0.5 0.5 1. ]
[ 0.5 0. 1. ]
[ 0.5 0.5 0.5]
[ 0.5 0. 0.5]
[ 0. 0. 0.5]
[ 0. 0. 0. ]
[ 0. 1. 0. ]
[ 1. 1. 0. ]
[ 1. 1. 1. ]
[ 1. 1. 0.5]
[ 0.5 1. 0.5]
[ 0.5 1. 0. ]
[ 0.5 0.5 0.5]
[ 0.5 0.5 0. ]
[ 0. 0.5 0. ]
[ 0. 0. 0. ]
[ 0. 0. 1. ]
[ 0. 1. 1. ]
[ 1. 1. 1. ]
[ 0.5 1. 1. ]
[ 0.5 0.5 1. ]
[ 0. 0.5 1. ]
[ 0.5 0.5 0.5]
[ 0. 0.5 0.5]
[ 0. 0. 0.5]
[ 0. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 1. 1. ]
[ 1. 1. 1. ]
[ 0.5 1. 1. ]
[ 0.5 1. 0.5]
[ 0. 1. 0.5]
[ 0.5 0.5 0.5]
[ 0. 0.5 0.5]
[ 0. 0.5 0. ]]

From the above it is evident, that the code
evaluates the Expression pot only on the
nodes of the form functions and
not on the gauss points!
I have been looking for a solution to this problem since June last year!

Thanks for all your help and patience!

regards

Moritz

commented Jan 31, 2014 by moritzbraun FEniCS User (1,390 points)

1 Answer

+2 votes

To your example
Note that there are kwargs element, degree for Expression constructor. By default, CG element of degree 1 is used to interpolate expression in UFC/DOLFIN evaluation chain to obtain finite element representation of expression in form. This representation is then really evaluated at Gauss quadrature points. To do different interpolation you can do

elem = FiniteElement('DG', tetrahedron, 0)
pot = Pot(element=elem)

so that expression is regarded as cell-wise constant. Or you can do

elem = FiniteElement('Q', tetrahedron, 4)
pot = Pot(element=elem)

but you must be able to determine a degree of guadrature element needed which is, I guess, here 4 as u*v term in the form is of polynomial degree 4. But you should check this precisely.

Other way
You may be able to define your singular potentials directly with UFL language avoiding usage of Expression. For example

x = V.cell().x
# Spherical coordinates
r     = x[0]
theta = x[1]
phi   = x[2]

# Singular potential
Phi = 1./r

# Now use Phi in some forms.
# But note that polynomial degree of Phi estimated to 1 by UFL
# so you may want to increase quadrature degree manually.

But this approach may not be of any use for you if you have many singular potentials seeded through whole domain. It is of use when formulating equations in cylindrical/spherical coordinates and 1/r term in Jacobian rises. The term is never evaluated in singularity (with properly aligned mesh).

answered Jan 31, 2014 by Jan Blechta FEniCS Expert (51,420 points)
edited Jan 31, 2014 by Jan Blechta

Maybe I misinterpreted a degree of quadrature element. Thinking about it, it maybe works the way that by picking up degree 4 quad element, you just tell to UFL that whole term u*v*pot or any other term with pot will be integrated by degree 4 quadrature rule which might be insufficient here (i.e. yield the same result as DG0 element).

commented Jan 31, 2014 by Jan Blechta FEniCS Expert (51,420 points)

Dear Jan

Thanks a lot again!

I have now figured out that I need
to define to FunctionSpaces
V of type Lagrange for the finite element functions
and
V1 of type "Quadrature" for the Expression
both with the same deree
If I do it like that the finite element calculation
works as intended in 3D for the hydrogen atom with the Coulomb potential
The next step for me will be to add custom integration for
those elements that included the point (0,0,0)
by extending the python script in the ffc sources accordingly.

regards

Moritz

commented Jan 31, 2014 by moritzbraun FEniCS User (1,390 points)
...