ufl team mailing list archive

[Martin Sandve Alnæs <martinal@xxxxxxxxx>]

On Tue, Jun 16, 2009 at 7:12 PM, Kristian
Oelgaard<k.b.oelgaard@xxxxxxxxxx> wrote:
> Quoting Martin Sandve Alnæs <martinal@xxxxxxxxx>:
>
>> On Tue, Jun 16, 2009 at 9:35 AM, Kristian
>> Oelgaard<k.b.oelgaard@xxxxxxxxxx> wrote:
>> > Quoting Martin Sandve Alnæs <martinal@xxxxxxxxx>:
>> >
>> >> 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.
>> >
>> > 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'm not into this notation, do you have any paper defining it?
>>
>> > 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
>>
>> 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.
>
>
> 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