Dealing with removable singularities in cylindrical coordinates

variational formulation
1

When working in cylindrical coordinates, one often has to deal with removable singularities (e.g. in linear elasticity): for instance, the divergence of a vector field u contains the term u_r/r. Mathematically, this is not an issue since u_r must vanish at r=0.

However, numerically this can be a problem, since the fraction of very small numbers is computed. Is there a canonical way to deal with such terms in fenics?

In the code in this example the singularities are simply ignored. Is this safe to do? Is fenics dealing with this issue automatically?

I am using fenics version 2018.1.0.

2

Hi,

indeed I did not discuss this issue in the example. Here are some thoughts on that:

  1. The most singular term is epsilon is in 1/r, the elastic energy being in 1/r^2 with an integration measure r dr so that in 1D, integrals over 0 to some distance of dr/r have to be computed.

  2. For a fixed mesh size, such integrals are well defined except on the first element close to zero. FEniCS uses quadrature points to compute such integrals with a precomputed or user-defined quadrature degree. Now you assemble for instance u_/x[0]*dx with u_ being a TestFunction, the only dof which is not well defined with respect to the quadrature degree is the first one located at r=0 all the other converge to a given finite value.
    If you impose a DirichletBC on the symmetry axis, this dof will be suppressed so that there are no problems any more.

  3. Convergence with respect to the mesh size is then ensured by the mathematical properties of the elasticity problem.

1 Like
3

Thank you for your helpful answer!

I understand that (for reasonable triangulations and realistic discretization parameters h) evaluations of the function (r,z) \mapsto u_r(r,z)/r are not an issue for quadrature points NOT lying on the symmetry axis r=0.

However, when evaluating ON the symmetry axis, the term u_r(r,z)/r|_{r=0} must be replaced by \partial_r u_r(0,z) since u_r(0,z)=0. So, is fenics doing this on its own?

4

Indeed, it seems that fenics can handle this. In the following very simple example, the correct result is returned:

import dolfin as df
import numpy as np

mesh = df.UnitSquareMesh(1, 1)
V = df.FunctionSpace(mesh, ā€˜Pā€™, 1)
T = df.Function(V)
T.vector()[:] = np.array([0, 0, 1, 1])

r = df.Expression(ā€˜x[0]ā€™, degree = 1)
print(df.assemble(T/r * df.dx))

Nevertheless, it would be interesting to know whatā€™s happening behind the scenes. It seems to me that fenics is not evaluating the function at r=0!?

5

The integrand is evaluated only at quadrature points, which, for typical Gaussian quadrature rules, are in the interior of each element. If you instead use a (non-default) quadrature rule with points at vertices, you will have problems:

# Default quadrature works with singularity:
print(df.assemble(T/r * df.dx))

# Vertex rule works without singularity:
print(df.assemble(1.0 * df.dx(domain=mesh),
                  form_compiler_parameters={'quadrature_rule': 'vertex'}))

# Vertex rule gives nan with singularity:
print(df.assemble(T/r * df.dx,
                  form_compiler_parameters={'quadrature_rule': 'vertex'}))
5 Likes
6

Thank you, this finally answers my question!

7

The singularity seems to be not properly processed in FEniCS.

Let us try the integral : Integral[ 1/Sqrt[x x+y y], {x,0,1},{y,0,1}] = 1.76275
Below is the dolfin code:

local_mesh = UnitSquareMesh(10,10)

singular_f1 = Expression("1/sqrt(x[0]*x[0]+x[1]*x[1])", degree = 8)
print( assemble( singular_f1*dx(local_mesh)) )

singular_f2 = Expression("sqrt(x[0]*x[0]+x[1]*x[1])", degree = 8)
print( assemble( Constant(1.0)/singular_f2*dx(local_mesh)) )

The first way to do integral gives: nan
The second way gives raw result: 1.75835412.

Is there any idea to the integrally correctly.

Update: The problem comes to me when I want to define an extended FEM with a base function of singularity. For example, for the Poisson equation over an L-shaped domain. The domain has re-entry corner at (0,0).

V = {CG FEM} + { phi }
phi = r^{2/3} sin(2/3 \theta) * CuttingFunction
8

One trick that might be helpful is to transform away the singularity by making a change of variable:
\hat{u}_r = u_r / \sqrt{r}.
This way, the transformed Lagrangian will not have any singular terms and can be integrated using the usual quadrature methods.

9

@Xuefeng_LIU
You can try the following:

from dolfin import *

mesh = UnitSquareMesh(10,10)
x = SpatialCoordinate(mesh)

print(assemble(1/sqrt(x[0]*x[0]+x[1]*x[1]) * dx, form_compiler_parameters={'quadrature_degree': 100}))

Of course playing with the quadrature degree and mesh size parameter changes the accuracy.

10

Many thanks. I have found the another approach, which gives the same value as your proposed one.
It seems that the integral is performed with the same quadrature using the value of 1/r over inertial point. But, it is obvious that if the integral is done under polar coordinate, then a better result can be obtained.

I am wondering that is it possible to let FEniCS perform the integral under polar coordinate?

from dolfin import *
mesh = UnitSquareMesh(10,10)
x = SpatialCoordinate(mesh)
print("Integral by spatial coordinate: ", assemble(1/sqrt(x[0]*x[0]+x[1]*x[1]) * dx, form_compiler_parameters={ā€˜quadrature_degreeā€™: 20}))

Q=FiniteElement(ā€œQuadratureā€,triangle,degree=20,quad_scheme=ā€˜defaultā€™)
sf = Expression("1/sqrt(x[0]*x[0]+x[1]x[1])", element=Q)
print("Integral by Quadrature element: ", assemble( sfdx(mesh,degree=20)) )

11

hi, I am trying to solve laplace equation in cylindrical coordinate with theta symmetry rāˆ‚āˆ‚r(rāˆ‚Vāˆ‚r)+āˆ‚2Vāˆ‚z2 =0. The discontinuity at r =0 in this equation how to remove in defining variational problem ? I am using rectangular mesh (0,1) ,(4,0). Vr =0 at r =0

12

Is this okay to do? It will change the stress/strain calculation, or?

One trick that might be helpful is to transform away the singularity by making a change of variable:
^ur=ur/āˆšr.
This way, the transformed Lagrangian will not have any singular terms and can be integrated using the usual quadrature methods.

13

No, if you transform back using
u_r = \hat{u}_r \sqrt{r}.

You will need to apply this to all derivatives, i.e. when finding strain components. Example:
\frac{\partial u_r}{\partial r} = \sqrt{r}\frac{\partial \hat{u}_r}{\partial r}+\frac{1}{2}r^{-1/2} \hat{u}_r.