Strange behavior of project vs interpolate for C++/CPP vs UFL Expression (azimuth, radial or polar angle)

1

Hi guys,

For a cylindrical geometry, I want to get the azimuth (i.e. radial or polar) angle \phi evaluated at all mesh coordinates.

I tried four approaches with interpolate/project and cpp/ufl expressions. However, if I use the ufl.atan_2 and try to project onto a simple Lagrange finite element function space, the angle is not correct at the \phi = \pi discontinuity (see screenshot). However, if a construct the angle with an cpp expression, it works using both interpolate and project.

So, the main question is, why the ufl expression behaves different than the cpp expression when using the project method (case 1 vs 3 in the MWE)?

13
Bildschirmfoto 2019-07-07 um 15.18.13.png2044×1452 285 KB

MWE:

from dolfin import *
import matplotlib.pyplot as plt
import ufl

mesh = RectangleMesh(Point(-1.0, -1.0), Point(1.0, 1.0), 10, 10)

el = FiniteElement("Lagrange", mesh.ufl_cell(), degree=1)
V = FunctionSpace(mesh, el)
x = SpatialCoordinate(mesh)

# Approach 1: Fails with project
phi = project(ufl.atan_2(x[1], x[0]), V)
# phi = interpolate(ufl.atan_2(x[1], x[0]), V) # fails, no eval() method..
plt.subplot(221)
p1 = plot(phi, title="phi1")
plt.colorbar(p1)
file = File("phi1.pvd")
file.write(phi)

# Approach 2: Works with interpolate and UserExpression (has eval method)
class RadialAngle(UserExpression):
    def eval(self, values, x):
        values[0] = ufl.atan_2(x[1], x[0])
phi2 = interpolate(RadialAngle(), V)
plt.subplot(222)
p2 = plot(phi2, title="phi2")
plt.colorbar(p2)
file2 = File("phi2.pvd")
file2.write(phi2)

# Approach 3: Works with project
phi3 = project(Expression("std::atan2(x[1],x[0])", element=V.ufl_element()), V)
plt.subplot(223)
p3 = plot(phi3, title="phi3")
plt.colorbar(p3)
file3 = File("phi3.pvd")
file3.write(phi3)

# Approach 4: Works with interpolate
phi4 = interpolate(Expression("std::atan2(x[1],x[0])", element=V.ufl_element()), V)
plt.subplot(224)
p4 = plot(phi4, title="phi4")
plt.colorbar(p4)
file4 = File("phi4.pvd")
file4.write(phi4)

plt.rcParams['figure.figsize'] = [15, 10]
plt.show()
3

The C++ Expressions are already interpolated in V (or, more generally, whatever element you pass to the constructor). Projecting them onto V is therefore a do-nothing operation. The UFL expression is being directly evaluated at Gauss points, so you are getting the overshoot/undershoot that you expect from L^2 projection of a discontinuous function onto a continuous space (in contrast to interpolation, which captures the discontinuity in a monotone way).

4

Thank you! Is there an option to interpolate a UFL expression without the need to define a UserExpression, simply like with the C++ expressions?

Something like that:

interpolate(ufl.atan_2(x[1], x[0]), V)
5

No, although, at least with CG1 elements, you can use lumped-mass projection to get monotone approximations of discontinuities (and avoid the linear solve in a consistent-mass projection):

from dolfin import *
mesh = UnitIntervalMesh(100)
V = FunctionSpace(mesh,"Lagrange",1)

def lumpedProject(f,V):
    v = TestFunction(V)
    lhs = assemble(inner(Constant(1.0),v)*dx)
    rhs = assemble(inner(f,v)*dx)
    u = Function(V)
    as_backend_type(u.vector())\
        .vec().pointwiseDivide(as_backend_type(rhs).vec(),
                               as_backend_type(lhs).vec())
    return u

# Compare results:
from matplotlib import pyplot as plt
x = SpatialCoordinate(mesh)
f = 0.2*sin(2.0*pi*x[0]) + conditional(gt(x[0],0.5),1.0,Constant(0.0))
plot(project(f,V))
plot(lumpedProject(f,V))
plt.show()
6

Thank you, that works!