Subject: [UFL-dev] Extending UFL

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 15 Jun 2009 14:33:14 +0200

Hi all,
thanks for all the nice comments and feedback at the workshop.
There was in particular some interest in how to extend UFL
with experimental features without having to interfere with the
main development branch. I'd like to compile a list of the kind
of features that are most interesting to people, so we can figure
out what to make extensible. Then I'll see what I can do.

Martin

---

From: Shawn Walker <[email-masked]>
Date: Mon, 15 Jun 2009 09:40:15 -0400 (EDT)

One experimental feature would be to add higher order 'curved' elements 
(e.g. quadratic triangles).  I know this is sort of in place, but it has 
not been implemented.  Maybe this belongs more with FFC.

- Shawn

On Mon, 15 Jun 2009, Martin Sandve Alnæs wrote:

> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev
>

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 15 Jun 2009 16:49:20 +0200

I'm only talking about UFL here, not FEniCS in general.
UFL can also be used by non-FEniCS software, which
may need certain extensions that are too specialized
to be added to FEniCS in general. This was requested
at the FEniCS workshop.

Martin

On Mon, Jun 15, 2009 at 3:40 PM, Shawn Walker<[email-masked]> wrote:
> (e.g. quadratic triangles).  I know this is sort of in place, but it has not
> been implemented.  Maybe this belongs more with FFC.
>
> - Shawn
>

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 16 Jun 2009 09:35:16 +0200

> out what to make extensible. Then I'll see what I can do.

Extending UFL with lifting operators.

Recall that for the Poisson equation we have:

a = \int (grad(v) + R(v))(grad(u) + R(u)) d\Omega
    + \sum_e \int (\eta r^e(v) r^e(u) d\Omega

with

R(u) = \sum_e r^e(u)
and
\int_{E = E^{+} + E^{-}} w r^e(u) d\Omega = - \int_e {w}[u] d\Gamma

where the last equation is the lifting operation.

I think we need to introduce two new functions in order to distinguish terms
like grad(v)R(u) and the stabilisation term r^e(v)r^e(u).
As mentioned at the workshop, we probably also need a new integral type.

So (some of) the terms in the Poisson equation could be represented as:

element  = FiniteElement("Discontinuous Lagrange", triangle, 1)
l_space = VectorElement("Discontinuous Lagrange", triangle, 0)

v = TestFunction(element)
u = TrialFunction(element)
R = LiftingFunction(l_space)
r = LiftingOperator(l_space)

a = inner(grad(v), R(u))*dE + inner(r(v), r(u))*dE

For this to work, we probably have to introduce a new type of integral in UFC
and modify the assembly in DOLFIN.

Kristian

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 16 Jun 2009 09:03:39 +0100

> and modify the assembly in DOLFIN.

I think that this can be cast in a more generic framework of static 
condensation. The question then is how UFL should let a form compiler 
know that a dof should be eliminated at some local level (typically at 
the level of a cell or the union of two cells). This may also require an 
extension of UFC (to return/accept local and global parts of the tensor).

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 16 Jun 2009 10:30:16 +0200

> extension of UFC (to return/accept local and global parts of the tensor).

OK, so you want to define e.g., v and r(v) on a mixed element and then eliminate
r(v)?
I guess, in addition to where, we also need to define how we eliminate the dofs
i.e.

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Tue, 16 Jun 2009 19:04:18 +0200

> where the last equation is the lifting operation.

I'm not into this notation, do you have any paper defining it?

> a = inner(grad(v), R(u))*dE + inner(r(v), r(u))*dE

I thought a bit about it, and I can make this notation work in UFL.
I'll represent the result of R(u) as an expression node, suggestion
for the name of this type? Same for r(u) I guess.
But it would help to understand roughly what r(u) does,
in particular to update some algorithms.

> and modify the assembly in DOLFIN.

Of course.

Martin

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Thu, 18 Jun 2009 18:02:38 +0200

On Tue, Jun 16, 2009 at 7:12 PM, Kristian
Oelgaard<[email-masked]> wrote:
> Note that Garth suggested implementing this in a more general way by
> dof-elimination. Have you thought about how that can be done? Do you even think
> this is possible?
>
> Kristian

My suggested solution (first considering the lifting operator only)
is something like this. Feel free to suggest better names if
something is misleading!

class LiftingOperatorResult(Operator):
    def __init__(self, operator, operand):
        Operator.__init__(self)
        self._operator = operator
        self._operand = operand

    def operands(self):
        return (self._operator, self._operand)
    # TODO: implement some functions from the Expr interface here

class LiftingOperator(Terminal):
    def __init__(self, element):
        Terminal.__init__(self)
        self._element = element

    def __call__(self, *args, **kwargs):
        if len(args) == 1 and len(kwargs) == 0:
            a, = args
            if isinstance(a, Expr):
                return LiftingOperatorResult(self, a)
        return Terminal.__call__(self, *args, **kwargs)
    # TODO: implement some functions from the Expr interface here

cell = triangle

l_space = VectorElement("DG", cell, 0)
r = LiftingOperator(l_space)

u_space = VectorElement("DG", cell, 1)
u = Function(u_space)

ru = r(u)

In addition, a LiftingFunction / LiftingFunctionResult
can be defined, and the L...Operator and L...Function
can share a baseclass and some common code.

The classes for R and r must be Terminal subclasses
because they depend on non-Expr data (l_space).
When applied to an Expr such as in r(u), the result
is an expression node of an Operator type, since
it depends on other Expr instances.

Since r and R are affine operators, most algorithms will be simple to update.

Derivatives are easy, we can just define d[r(u)]/dv = r(du/dv).

Is it ok to simply disallow free indices in the argument to r and R?

The value shape of r(u) == the value shape of jump(u, n).

I didn't do much here that is particular to the
lifting operators, but the types carry a lot of
information about what should be done.
I don't know enough about static condensation
to define a generalization as Garth suggested.

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 19 Jun 2009 10:32:33 +0200

> The value shape of r(u) == the value shape of jump(u, n).

Looks good. What about the integral/measure? Did you reserve a new token for the
macro domain integral.

> to define a generalization as Garth suggested.

Me neither, I'd settle for the above as a first try.

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Fri, 19 Jun 2009 12:17:18 +0200

> Looks good. What about the integral/measure? Did you reserve a new token for the

I implemented a way to add new measure types
(simply identified by a string in UFL) last week.
What should we call it? "macro cell integral"?
Is "dE" generally used for this?
I guess this leads to adding ufc::macro_cell_integral as well?

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 19 Jun 2009 13:20:59 +0200

Quoting Martin Sandve Alnæs <[email-masked]>:

> >> >> > I think we need to introduce two new functions in order to
> distinguish
> >> >> terms
> >> > dof-elimination. Have you thought about how that can be done? Do you
> even
> >> think
> >> > this is possible?
> >> >
> >> > Kristian
> >>
> >>

I don't know, but it looks OK to me.

> I guess this leads to adding ufc::macro_cell_integral as well?

Yes, but we don't need that until we start generating code for these operators.
So we can add it once we know how the interface should look like.

Kristian