1/x near boundary, NaNs where there shouldn't be NaNs

Question

I have a PDE in cylindrical coordinates on the unit square domain that contains the term 1 / x[0] on the right-hand side. One might assume this is problematic given the fact that the domain includes x[0]=0; surely there will be NaNs. And indeed there are for all points at the boundary. Luckily, those are ironed out by Dirichlet conditions later on.

Unfortunately, FEniCS also produces NaNs for points that are only part of an element that touches the boundary. Mathematically, its shape function is linear in x[0] so the singularity cancels out nicely. FEniCS however has no way of knowing that, encountering a 0 / 0 and spitting out a NaN.

Here's an MWE:

from dolfin import (
    parameters, UnitSquareMesh, Expression, as_backend_type,
    FunctionSpace, dx, assemble, TestFunction
    )
parameters['linear_algebra_backend'] = 'Eigen'

mesh = UnitSquareMesh(1, 1)

V = FunctionSpace(mesh, 'Lagrange', 1)
v = TestFunction(V)

one_over_x = Expression('1.0 / x[0]', degree=5)
f = one_over_x * v * dx

b = assemble(f)

b = as_backend_type(b)
for k in range(4):
    print(b[k])

All four elements are NaNs while mathematically, only the two x=0-boundary elements should be.

Curiously, when replacing the one_over_x term by

one_over_x = 1.0 / Expression('x[0]', degree=1)

_none_ of the elements are NaNs. No idea why.

The actual RHS in the application comes as an C-expression from a sympy computation, so I guess I won't be able to use that workaround.

Any other ideas? (Perhaps forcing a quadrature that doesn't include boundary points?)

Answer 1

This happens because the Expression is interpolated to continuous Lagrange elements instead of being directly evaluated in the quadrature points. You can avoid this by specifying a quadrature element. For example,

degree = 5
Q = FiniteElement("Quadrature", triangle, 
                  degree = degree, quad_scheme = "default")
one_over_x = Expression('1.0 / x[0]', element = Q)
f = one_over_x * v * dx(degree = degree)

Sometimes it is possible to formulate the coefficient as a function of the SpatialCoordinate instead of using Expression, for example

x = SpatialCoordinate(mesh)
f = (1 / x[0]) * v * dx(degree = degree)

Comments

Ah cool, I didn't about this. Thanks! I do get a bunch of errors like

Exception: Points must be equal to coordinates of quadrature points

now though. Will I have to carry degree=degree through the entire code base?

---

Also, isn't the degree setting wrong in the dx? It should be a degree higher since the (linear) v joined the game, shouldn't it?

---

When you specify the element for an Expression, this element will be used instead of the default Lagrange element family when the Expression is interpolated. The "dofs" for the quadrature element are evaluation in the quadrature points, but this element does not have a defined basis. This means that a Function with this element type cannot be evaluated outside of the quadrature points.

The quadrature rule defining the quadrature element therefore needs to be the same as the quadrature rule used to assemble the form. To ensure that this is the case, you should specify the degree of the measure dx or pass it as a form compiler parameter when assembling.