ufl team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxxx>]

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 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
>