UFL condition in if statement

Question

Hi All,

thanks in advance for your help.

I have problems understanding how UFL conditions should be handled correctly.
Here is an example

import dolfin as dlfn
ScalarSpace = dlfn.FunctionSpace(mesh, 'DG', 0)
Yinit_value = 1.0
Yinit = dlfn.project(Yinit_value, ScalarSpace)
Y = dlfn.project(1.0/wf*(heavisideFunction(sigmaH)*K/2.0*eps_solve[k,k]**2
                +mu*eps_solve_dev[i,j]*eps_solve_dev[j,i]), ScalarSpace)

Y is computed from a solution obtained solving a weak form. Results are as follows

Y.vector().array()
Out[80]: 
array([ 3.59723591,  3.59723591,  3.59723591,  3.59723591,  3.59723591,
        3.59723591])
Yinit.vector().array()
Out[81]: array([ 1.,  1.,  1.,  1.,  1.,  1.])

I now want to do something like

if dlfn.conditional(dlfn.gt(Y,Yinit), 1, 0)==1:
    print "True"
else:
    print "False"

Unfortunately this gives "False" although Y > Yinit...

I also have a workaround

diff=Y.vector().array()-Yinit.vector().array()
diff[diff<0]=0 # set all values in diff which are <0 to zero
if diff.any()==True:
        print "True"

My question is now, how does this work with UFL conditionals?
(without using arrays)
Also how would pointwise conditional operations on tensors of rank 2 work?
(is it necessary to formulate such for each component in an loop?)

Thanks, Philipp

Answer 1

Hi, dlfn.conditional(dlfn.gt(Y,Yinit), 1, 0) does not evaluate to 1 or 0 (check it's type). This is why your if test returns False. Your conditional is something that when compiled with FFC will act in a form as lamba x: 1 if Y(x) > Yinit(x) else 0. Note, that comparison works iff Y, Yinit are scalar valued. Conditional branches can be tensor valued expressions (the rank of the two must be identical).

from dolfin import *

mesh = UnitSquareMesh(100, 10)

x = SpatialCoordinate(mesh)
f = conditional(gt(x[0], 0.5), x[0]+x[1], Constant(0))
# In the generated code f will behave as (x, y) -> 1 if x+y > 0.5 else 0
V = FunctionSpace(mesh, 'DG', 1)
fh = project(f, V)

plot(fh, interactive=True)

# Make up tensor
s = outer(grad(f), grad(f))

# Conditionals can only compare scalar valued expression
A = Constant(((0, 0), (0, 0)))
try: # Fail comparing two matrices
    conditional(eq(s, A), Constant(1), Constant(2))
except dolfin.UFLException as e: 
    pass
# Compare S, A vie (S-A):(S-A)
g = conditional(le(abs(inner(s-A, s-A)), Constant(0.1)), Constant((0, 1)), Constant((0, 0)))

W = VectorFunctionSpace(mesh, 'DG', 0)
gh = project(g, W)
plot(gh, interactive=True)