how to evaluate the angle from the x-axis of a gradient

Use the math operators from UFL, not NumPy.

On 11/11/10 12:45, Thomas Fraunholz wrote:
> New question #133595 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/133595
>
> I have a time-dependent Cahn-Hilliard-Type problem and I need to
> evaluate the angle from the x-axis of a gradient in order to define
> my variational problem. The angle of the gradient can be calculated
> by using phi = arctan( dc/dy / dc/dx ). This is needed to model an
> anisotropic N-fold surface tension with the help of cos(N*phi).
>
> A very native approach was to use the following lines of code. It
> doesn't work of course. But I think it helps to understand my
> problem.
>
> -------------------- import numpy as np from dolfin import *
>
> ...
>
> V = FunctionSpace(mesh, "CG", 2) w = TestFunction(V)
>
> c = Function (V) # current solution c0 = Function(V) # old
> solution
>
> ...
>
> # Anisotropic N-fold surface tension: def anisotropic_tension(c0):
> return np.cos(N*np.arctan2(c0.dx(1),c0.dx(0)))
>
> L = L - dt*anisotropic_tension(c0)*(inner(div(grad(w)),
> div(grad(c)))*dx)
>
> ...
>
> --------------------
>
> I think one basic problem is the interaction with numpy and the
> encapsuled evaluation of the derivatives. I'm glad for any hint how
> to proper evaluate the angle and how to use it in order to define L.

Try using 'Atan' (a UFL function) directly, e.g.

M = dt*N*Atan(....)* . . . . .

Garth

>

Ok. When I use

L = L - dt*(inner(div(grad(w)), anisotropic_tension(c0)*div(grad(c)))*dx)

and the simplified expression

def anisotropic_tension(c0):
    return cos(c0.dx(1))

it works. But there is still a problem with the arctan I need. If I use the ufl version, e.g.

def anisotropic_tension(c0):
    return cos(N*atan(c0.dx(1)/c0.dx(0)))

I got a native problem when c0.dx(0) equals zero. Then the solver returns in every step 'not a number'. So I'm looking for a possibility like atan2 in the cmath library or the arctan2 numpy version. A first try with an if/else expression didn't helped either. But I'm not really sure if I am allowed to writefor example

def anisotropic_tension(c0):
    if c0.dx(0) > 0.01:
        return cos(N*atan(c0.dx(1)/c0.dx(0)))
    else:
        return cos(N*atan(c0.dx(1)/0.01))

Thanks in advance for your help,

Thomas

Thanks Kristian B. Ølgaard, that solved my question.