ufl team mailing list archive

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

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

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