Replacing FacetNormal produces conj

[Jørgen Dokken]

By replacing the second argument of an inner product with a FacetNormal, one obtains the conjugate of the integral. This is not correct.

Minimal example:

from dolfin import *
from ufl import replace
mesh =  UnitSquareMesh.create(10,10, CellType.Type.triangle)
u = Constant((2,0))
v = Constant((1,2))
n = FacetNormal(mesh)
L = inner(u,v)*ds(domain=mesh)
Lun = replace(L, {u: n})
Lvn = replace(L, {v: n})
print(Lun)
print(Lvn)

Produces the following output:

{ (n) : (f_7) } * ds(<Mesh #3>[everywhere], {})
{ conj(((n) : (f_6))) } * ds(<Mesh #3>[everywhere], {})

[Stephan Schmidt]

Also, the behavior of inner and conjugate is completely random and unclear. Have a look at this:

from dolfin import *

a = Constant((1.0,1.0))
b = Constant((1.0,1.0))

print(str(inner(a,b)))
print(str(inner(b,a)))

The output is:

(f_0) : (f_1)
conj(((f_0) : (f_1)))

So right now, I see no reason, why apparently "a" should be the first argument and "b" is the second to avoid a conjugation. In any case, the conjugate should be part of "inner" or the ":", so the form/string representation in either case should not get a conjugate, as it is part of the "inner" or the ":" anyway...

[Jørgen Dokken]

It appears that the conjugate is dependent of the order of declaration of variables. That is fucked up. Minimal example:

from dolfin import *

a = Constant((1.0,2.0))
c = Constant((4.0,5.0))
b = Constant((2.0,3.0))

print(inner(a,b))
print(inner(b,a))
print("...")
print(inner(b,c))
print(inner(c,b))

Output:

(f_0) : (f_2)
conj(((f_0) : (f_2)))
...
conj(((f_1) : (f_2)))
(f_1) : (f_2)

[Stephan Schmidt]

Anyway, there should never be a conjugate because it is built in in "inner", "dot" or ":" and "." respectively already....

[Lawrence Mitchell]

So inner(a, b) gets transformed, later in the day, into sum(a[i]*conj(b[i]), i). This is where the conjugate is inserted.

To allow for common subexpression elimination. Where you might have inner(a, b) and inner(b, a), in non-complex UFL these are canonicalised such that the ordering of operands is consistent.

Check out what happens on master:

In [2]: a = Constant((1.0,2.0))
   ...: c = Constant((4.0,5.0))
   ...: b = Constant((2.0,3.0))
   ...: 
   ...: print(inner(a,b))
   ...: print(inner(b,a))
   ...: print("...")
   ...: print(inner(b,c))
   ...: print(inner(c,b))
   ...: 
   ...: 
(w_0) : (w_2)
(w_0) : (w_2)
...
(w_1) : (w_2)
(w_1) : (w_2)

Notice how inner(a, b) and inner(b, a) have the same string representation. This is because they are the same.

To allow for this to still occur in complex mode, we need to sort the operands. But, this might result in an invalid inner. So, we use the identityinner(a, b) == conj(inner(b, a)).

And now we can sort the operands. This sorting is based on the count of a Coefficient, which of course, varies according to the order in which you defined the coefficients.

In summary, not a bug, think harder.

[Jørgen Dokken]

@wence How do you get this string representation. Using docker with latest dev, i do not get this.

In [4]: from dolfin import *
   ...: 
   ...: a = Constant((1.0,2.0))
   ...: c = Constant((4.0,5.0))
   ...: b = Constant((2.0,3.0))
   ...: 
   ...: print(inner(a,b))
   ...: print(inner(b,a))
   ...: print("...")
   ...: print(inner(b,c))
   ...: print(inner(c,b))
   ...: 
   ...: 
(f_9) : (f_11)
conj(((f_9) : (f_11)))
...
conj(((f_10) : (f_11)))
(f_10) : (f_11)

[Stephan Schmidt]

I did think hard about this, and I can confirm inner(a,b) and inner(b,a) should have the same string representation, but they don't... that's the problem....

How does your example above resolve constants to "w_", but in my fenics, they are resolved to "f_"?

[Jan Blechta]

inner(a, b) and inner(b, a) cannot have the same representation. And this can't be stressed more. As Lawrence said inner is roughly equivalent to sum(a[i]*conj(b[i]), i). This is a decision which has been taken for definition of inner. I believe that time to debate about definition of inner is over. Please follow the #complex channel on Slack: https://fenicsproject.slack.com/messages/C93RPKTNC/

[Lawrence Mitchell]

Why should inner(a, b) and inner(b, a) have the same string representation, surely, they should map to a : b and b : a (i.e., the argument order is changed).

w_ is because I happened to run in firedrake (which labels constant differently), but the UFL issues is the same, I think.

I should note that the UFL I was using is a few weeks out of date (so complex is not merged). Hence no conj anywhere.

[Lawrence Mitchell]

Not a bug, but see https://bitbucket.org/fenics-project/ufl/issues/107 for plausible changes to UFL that would perhaps make things easier for Jørgen and Stephan.