ufl team mailing list archive

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

On Fri, Jun 19, 2009 at 10:32 AM, Kristian
Oelgaard<k.b.oelgaard@xxxxxxxxxx> wrote:
> Quoting 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).
>
> Looks good. What about the integral/measure? Did you reserve a new token for the
> macro domain integral.

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


>> 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.
>
> Me neither, I'd settle for the above as a first try.
>
> Kristian
>
>> Martin
>> _______________________________________________
>> UFL-dev mailing list
>> UFL-dev@xxxxxxxxxx
>> http://fenics.org/mailman/listinfo/ufl-dev
>>
>
>
> _______________________________________________
> UFL-dev mailing list
> UFL-dev@xxxxxxxxxx
> http://fenics.org/mailman/listinfo/ufl-dev
>