Subject: [UFL-dev] Restriction syntax

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 14 May 2009 07:15:56 +0100

Is the correct syntax for the restriction of the gradient (or any other 
operator) of a function to one side a facet

     grad(u)('+')

or

     grad(u('+'))

?

Garth

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Thu, 14 May 2009 09:31:24 +0200

Is there a mathematical difference?

Martin

On Thu, May 14, 2009 at 8:15 AM, Garth N. Wells <[email-masked]> wrote:
> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev
>

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 14 May 2009 08:40:41 +0100

> Is there a mathematical difference?

Splitting hairs, I would say yes and that grad(u('+')) is correct and 
grad(u)('+') is not because grad(u) is not defined in a classical sense. 
I prefer grad(u('+')), but FFC used to use grad(u)('+'). My preference 
is however not strong.

My immediate question is will the generated code differ for grad(u('+')) 
and grad(u)('+')?

Garth

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Thu, 14 May 2009 09:42:36 +0200

On Thu, May 14, 2009 at 9:40 AM, Garth N. Wells <[email-masked]> wrote:
> grad(u)('+') is not because grad(u) is not defined in a classical sense. I
> prefer grad(u('+')), but FFC used to use grad(u)('+'). My preference is
> however not strong.
>
> My immediate question is will the generated code differ for grad(u('+')) and
> grad(u)('+')?

SFC doesn't have interior facet integrals,
but there should be no problem doing this
whatever way the right way is.

Martin

---

From: Anders Logg <[email-masked]>
Date: Thu, 14 May 2009 09:43:26 +0200

> My immediate question is will the generated code differ for grad(u('+')) 

No, since the restriction will be propagated to the leaves... :-)

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Thu, 14 May 2009 09:55:10 +0200

> No, since the restriction will be propagated to the leaves... :-)

The leaves are terminal expressions, that is literal constants,
geometric quantities, functions and basis functions.
"Propagating to the leaves" means modifying for example

  (grad(f) + grad(g))('+')
to
  (grad(f('+')) + grad(g('+')))

such that the restrictions apply directly to the terminal expressions.
That will simplify the work of the form compilers, and allow
for example taking the jump of anything.

Basically, I think the canonical representation after some algorithms should be
  u('+').dx(i)
and with indexing either
  u('+').dx(i)[j]
or
  u('+')[j].dx(i)

If restrictions, (indexing) and derivatives are applied to a
(non-scalar) function (terminal).
Currently, both
  u.dx(i)[j]
and
  u[j].dx(i)
can occur, and restrictions are not handled.

Martin

---

From: Anders Logg <[email-masked]>
Date: Thu, 14 May 2009 10:35:34 +0200

> Basically, I think the canonical representation after some algorithms should be

I prefer this one:

  u('+')[j].dx(i)

Incidentally, the monomial extraction in FFC represents basis
functions or functions as

  class MonomialFactor:

      def __init__(self, arg=None):
          if isinstance(arg, MonomialFactor):
              self.function = arg.function
              self.components = arg.components
              self.derivatives = arg.derivatives
              self.restriction = arg.restriction
          ...

-- 
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 14 May 2009 16:59:06 +0100

> "Propagating to the leaves" means modifying for example

So the above works now? It would be super useful for me.

Garth

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Thu, 14 May 2009 23:09:16 +0200

> So the above works now? It would be super useful for me.

I don't know what FFC does.

But I implemented the propagation right now, only tested
for a simple case but it's a simple algorithm:

# test_restrictions.py:
from ufl import *

e = FiniteElement("DG", triangle, 1)
f = Function(e)
g = Function(e)
a = (grad(f) + grad(g))('-')
print a

from ufl.algorithms import propagate_restrictions

b = propagate_restrictions(a)
print b

martinal@martinal-znote:~/dev/fenics/ufl/sandbox/algorithmtests$
python test_restrictions.py
((grad(w_0)) + (grad(w_1)))('-')
(grad((w_0)('-'))) + (grad((w_1)('-')))

FFC just needs to apply propagate_restrictions to the integrand of
interior facet integrals.

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 15 May 2009 09:56:39 +0200

> interior facet integrals.

I just tried it and it appears to be working. However, I got the error:

*** FFC: Not expecting twice restricted expression.
/home/oelgaard/fenics/clone/ffc/scripts/ffc:131: DeprecationWarning:
BaseException.message has been deprecated as of Python 2.6
  error(exception.message)
*** FFC: Not expecting twice restricted expression.
*** FFC: To get more information about this error, rerun FFC with --debug.

for the following form (from Garth)

sigma_a = -b*u - gamma*b*(uhat - u)
sigma_d = kappa*grad(u) + (alpha/h)*kappa*(uhat -u)*n('+')
sigma_i = dot(sigma_a + sigma_d, n)

sigma_plus = sigma_i('+')

a = v('+')*sigma_plus*dS

The reason is that sigma_d already contains n('+'). What should happen in a case
like this? For complex forms it is sometimes easier to declare temporary
expressions and in these cases one might accidentally (or intentionally?)
restrict a function multiple times. Should a restriction just be overwritten or
what would you expect Garth for this form?

Kristian

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 15 May 2009 08:58:50 +0100

> what would you expect Garth for this form?

I would like to get an error. n('+') shouldn't be in sigma_d (my mistake).

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 15 May 2009 10:02:43 +0200

> I would like to get an error. n('+') shouldn't be in sigma_d (my mistake).

OK, fine by me. I have updated FFC to use propagate_restrictions() it should be
working now.

Kristian