Singular right-hand side

variational formulation
1

Hello,
I want to solve a PDE with a rhs f that is singular at the origin (see example below). Increasing the number of quadrature points does not help for coarse grids. Is there a way to apply some composite quadrature rule? Moreover, can a increase the number of quadrature points locally (more quadrature points for triangles close to the origin)? Best regards, Johannes

from fenics import *

p = 16./15
exp_str = ‘’’-(pow(2,-1 + p)(-(ppow(x[0],2)) + pow(x[0],4) + (-p + 10*(-1 + p)pow(x[0],2) + (8 - 7p)pow(x[0],4) + (-1 + p)pow(x[0],6))
pow(x[1],2) + (1 + (8 - 7p)pow(x[0],2) + (-6 + 4p)*pow(x[0],4))*pow(x[1],4) + (-1 + p)pow(x[0],2)pow(x[1],6))
pow(pow(x[1],2) + pow(x[0],2)(1 + (-4 + pow(x[0],2))*pow(x[1],2) + pow(x[1],4)) + pow(10,-14),(-4 + p)/2.))’’’
f = Expression(exp_str,p = p,degree = 5)

mesh= UnitSquareMesh(1,1)

print(assemble(f2*dx(mesh,metadata={‘quadrature_degree’: 20})))
print(assemble(f2*dx(mesh,metadata={‘quadrature_degree’: 1})))

for ref in range(0,9):
mesh = refine(mesh)
print(assemble(f**2*dx(mesh,metadata={‘quadrature_degree’: 0})))

2

Expressions get interpolated to FE-spaces, in your case a CG space of degree 5. So the expression no longer has a singularity. Try writing f as a ufl function and use it directly in the integrals, e.g.

from dolfin import *

mesh = UnitSquareMesh(10,10)

x = SpatialCoordinate(mesh)
f = x[0]**(-2) + x[1]**(-2)

print(assemble(f**2 * dx))

Quadrature degree will automatically be estimated.

You can try locally refining the mesh near the singularity to improve quadrature locally.