This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

UFL conditional with Coefficient

+2 votes

Hi fenics users!
I am trying to define a form in which one of the integrals is on a part of the domain defined by $ E = (x : a > \epsilon ) $ and $ a $ is a = Coefficient(Q) where Q is some scalar Finite Element space.

I thought to integrate on the whole domain the characteristic function of $ E $, using the conditional function provided by UFL like this

Pen = conditional(gt(a,eps),1,0)*dx

but ffc complains that this is an invalid expression, in particular the gt(a,eps).

What's my best option? How should I do this?
Thanks in advance!
Massimiliano

closed with the note: Solved
asked Nov 10, 2014 by Massimiliano Leoni FEniCS User (1,760 points)
closed Nov 12, 2014 by Massimiliano Leoni

I found out additional infos: the line

Pen = 1.0/gammaB*conditional(gt(ac,1-tol),1,0)*dx

compiles well, but

Pen = 1.0/gammaB*(1 - ac)**2*conditional(gt(ac,1-tol),1,0)*dx

doesn't, and ffc complains about the conditional, which is the same as in the first line.

For some reason adding the (1 - ac)**2 makes it freak out.
The error is ufl.log.UFLException: Found Argument in Conditional(GT(Argument(FiniteElement('Lagrange', Domain(Cell('triangle', 2), label=None, data=None), 2, quad_scheme=None), 1, None), Sum(IntValue(1, (), (), {}), Product(IntValue(-1, (), (), {}), Constant(Domain(Cell('triangle', 2), label=None, data=None), 13)))), IntValue(1, (), (), {}), Zero((), (), {})), this is an invalid expression.

Any idea?

Hi, can you post more of the ufl file? This compiles with FFC 1.4.0+

element = FiniteElement('Lagrange', triangle, 2)

c = Coefficient(element)
alpha = Constant(triangle)

form = 1/alpha*(1-c)**2*conditional(gt(c, 0.2), 1, 0)*dx
forms = [form]

Wow, this is really crazy. My complete UFL file includes the same commands you listed, and if I delete all but those, it works O.O

In particular, the following snippet alone works,

cell = triangle
Q = FiniteElement("Lagrange", cell, 2)

ac = Coefficient(Q)
gammaB = Constant(cell)
tol = Constant(cell)
Pen = 1.0/gammaB*(1 - ac)**2*conditional(gt(ac,1-tol),1,0)*dx

forms = [Pen]

whereas the whole script doesn't :S

cell = triangle

V = VectorElement("Lagrange", cell, 2)
Q = FiniteElement("Lagrange", cell, 2)

u = TrialFunction(V)
v = TestFunction(V)
a = TrialFunction(Q)
b = TestFunction(Q)

f = Coefficient(V)
l = Constant(cell)
w1 = Constant(cell)

uc = Coefficient(V)
ac = Coefficient(Q)

E = (1 - ac)**2
w = w1 * ac

gammaB = Constant(cell)
tol = Constant(cell)
Pen = 1.0/gammaB*(1 - ac)**2*conditional(gt(ac,1-tol),1,0)*dx

# dimensionless version
P = (0.5*inner(E*sym(grad(uc)),sym(grad(uc))) + w +
     0.5*w1*l**2*inner(grad(ac),grad(ac)) - inner(f,uc))*dx \
    + Pen
Pu = derivative(P, uc, v)
Pa = derivative(P, ac, b)

Pu = replace(Pu, {uc:u})
Pa = replace(Pa, {ac:a})

au = lhs(Pu)
Fu = rhs(Pu)

aa = lhs(Pa)
Fa = rhs(Pa)


# nonlinear version
Panl = derivative(P, ac, b)         ## ~ E_alpha
PanlJ = derivative(Panl, ac, a)     ## ~ E_alpha_alpha


element = [V, Q]
forms = [au, aa, Fu, Fa, Panl, PanlJ]

EDIT: if I use another acc = Coefficient(Q) in the definition of Pen, it works fine. Could it be some mess with the replacing?? That's the only thing that affects ac.
Also, can automatic differentiation handle this?

Hi,

the problem is that the call

Pa = replace(Pa, {ac:a})

puts a TrialFunction into the conditional, i.e. it becomes conditional(gt(a,1-tol),1,0) and that is illegal. If you define the Pen with acc, which does not get substituted in the call to replace things will work which is what you observed.

Yes, that's what I was thinking! Thank you for your help!

...