UFL conditional with Coefficient

Question

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

Comments

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!