How to calculate radial coordinate angle with uflacs?

Question

It seems like uflacs doesn't support the atan_2 function, as I get an error for the code below. Do you know another way to calculate an angle coordinate that will not crash if x= 0?

from dolfin import *
parameters["form_compiler"]["representation"] = "uflacs"

mesh = UnitSquareMesh(2,2)
VV = VectorFunctionSpace(mesh, "CG", 1)
V = VV.sub(0).collapse()

r2 = 1.0/sqrt(2)
e = interpolate(Constant([r2,r2]), VV)

phi = project(atan_2(e[1], e[0]), V)
plot(phi, interactive = True)

Comments

Hi, I am not getting any errors running your code with FEniCS stack compiled on 04/05/2016. What is your setup?

---

This has already been fixed in the development version and will be part of the next release.

Answer 1

For those who are running the current release here is a workaround (thanks Øyvind Evju) and a few tests.

from dolfin import *
parameters["form_compiler"]["representation"] = "uflacs"

def atan_2(y,x):
    return conditional(ne(y, 0), 2*atan((sqrt(x**2+y**2)-x)/y), conditional(gt(x,0), 0, conditional(lt(x,0), pi, 1.0e20)))

def check_value(x,y, atanval, VV, V):
    e = interpolate(Constant([x, y]), VV)
    theta = project(atan_2(e[1], e[0]), V)
    assert near(theta.vector().array()[0], atanval, 1.0e-13)

def test_atan2():
   mesh = UnitIntervalMesh(1)
   VV = VectorFunctionSpace(mesh, "CG", 1, dim = 2)
   V = VV.sub(0).collapse()

   r2 = 1.0/sqrt(2)
   check_value(r2, 0, 0.0, VV, V)
   check_value(r2, r2, pi/4.0, VV, V)
   check_value(0, r2, pi/2.0, VV, V)
   check_value(-r2, r2, 3*pi/4.0, VV, V)
   check_value(-r2, 0, pi, VV, V)
   check_value(-r2, -r2, -3*pi/4.0, VV, V)
   check_value(0, -r2, -pi/2.0, VV, V)
   check_value(r2, -r2, -pi/4.0, VV, V)

   test_atan2()