Subject: [FIAT-dev] Evaluation of functionals

From: Marie Rognes <[email-masked]>
Date: Mon, 10 Dec 2007 15:00:24 +0100

We would like to evaluate the degrees of freedom on a given function in 
ffc. In order to do that with relative easy, it would be very handy if 
we could extend the current information about the dual bases in FIAT.

We were thinking that if we could represent functionals in terms of 
three attributes, 'points', 'components' and 'weights', then it would be 
very easy to generate code for the evaluation of these in ffc.

Some examples for the motivation for the choice of representation

    Example 1: L = Point evaluation at component k:

       weights = [1]
       components = [k]
       points = [p_0]

    Then given a function f, we want

       L(f) = 1*f[k](p_0)

    Example 2: L = Normal component at diagonal edge:

       weights = [1, 1]
       components = [0, 1]
       points  = [p]

    Then

       L(f) = [1, 1]*(f[0](p), f[1](p)) = 1*f[0](p) + 1*f[1](p)

    Example 3: L = Integral over interior of component k:

       weights, points = make_quadrature(...)
       components = [k]

       L(f) = [w_0, ... w_m] (f[k](p_0), .... ,f[k](p_m))
            = sum_j w_j f[k](p_j)

Maybe the finite elements/dual bases in FIAT could supply this 
information in addition to the 'entity_ids'?

I was thinking something along the lines of letting the Functional class 
have the additional attributes, and initializing these when constructing 
the dual bases.

If there are some attributes that don't apply (such as weights and 
component for scalar Lagrange elements), these could of course be left 
empty.

How does this sound?

-- 
Marie E. Rognes
Ph.D Fellow, 
Centre of Mathematics for Applications, 
University of Oslo
http://folk.uio.no/meg

---

From: "Robert Kirby" <[email-masked]>
Date: Mon, 10 Dec 2007 09:08:33 -0600

By "given function", do you mean an arbitrary input function (something that
supports __call__) or something in the finite element space?

Functionals know how to apply themselves to the finite-dimensional space
members via __call__ (it's just a dot product).

There is some preliminary support for the types of functionals in the
"get_type()" method.  It would need to be expanded to include some of this
other information.

On Dec 10, 2007 8:00 AM, Marie Rognes <[email-masked]> wrote:

> _______________________________________________
> FIAT-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/fiat-dev
>

---

From: Marie Rognes <[email-masked]>
Date: Mon, 10 Dec 2007 16:44:11 +0100

Robert Kirby wrote:
> By "given function", do you mean an arbitrary input function 
> (something that supports __call__) or something in the finite element 
> space?
>

Arbitrary input functions, or more(?) precisely, a function that we know 
how to evaluate at a point but not much more.

> Functionals know how to apply themselves to the finite-dimensional 
> space members via __call__ (it's just a dot product).
>

Yes, but we typically want to use this in order to interpolate functions 
onto the element spaces...

> "get_type()" method.  It would need to be expanded to include some of 
> this other information.

Yes, that is a starting point. If you think the idea sounds sensible, 
then I can take a shot at making it work?

---

From: Anders Logg <[email-masked]>
Date: Mon, 10 Dec 2007 18:42:13 +0100

> very easy to generate code for the evaluation of these in ffc.

Just a comment. We need this because we generate code for the UFC
function evaluate_dof which is called at run-time to interpolate
user-defined functions (like coefficients, right-hand sides) to the
finite element space.

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Mon, 10 Dec 2007 20:51:47 +0100

> supports __call__) or something in the finite element space?

A function that can be evaluated at any point. For example, someone
has specified a right-hand side by "return sin(x)".

-- 
Anders

> There is some preliminary support for the types of functionals in the "get_type
> ()" method.  It would need to be expanded to include some of this other
> information.

---

From: "Robert Kirby" <[email-masked]>
Date: Mon, 10 Dec 2007 18:47:27 -0600

1.) This is ill-posed, in general.  You can't integrate a generic function
just by using it as a black-box.  This handles Lagrange interpolation, but
not Raviart-Thomas, which has integral moments against some functions as
some of the nodes.
2.) Why don't you do L2 projection instead?  Generally, you get a mass
matrix that is so much cheaper to solve than anything else you're doing that
it doesn't cost very much.  Everything is already in place for that.

On Dec 10, 2007 1:51 PM, Anders Logg <[email-masked]> wrote:

> "get_type

---

From: Anders Logg <[email-masked]>
Date: Tue, 11 Dec 2007 10:25:59 +0100

On Mon, Dec 10, 2007 at 06:47:27PM -0600, Robert Kirby wrote:
> 1.) This is ill-posed, in general.  You can't integrate a generic function just
> by using it as a black-box.  This handles Lagrange interpolation, but not
> Raviart-Thomas, which has integral moments against some functions as some of
> the nodes.

Can't we just include the value of the function we integrate against
in the weights?

> 2.) Why don't you do L2 projection instead?  Generally, you get a mass matrix
> that is so much cheaper to solve than anything else you're doing that it
> doesn't cost very much.  Everything is already in place for that.

L2 projection between which spaces? What we want to do is to
interpolate (evaluate the dofs) locally on each element. We don't
solve a global system.

-- 
Anders

---

From: "Robert Kirby" <[email-masked]>
Date: Tue, 11 Dec 2007 07:41:24 -0600

L2 projection into the finite element space.  On each element, or globally,
solve( pi u - u , w ) = 0 for a function pi u for all of the test functions
in the finite element space.

All you need is point evaluation on the input function u for quadrature, and
the weighted values of u become the rhs of the system.

We could include the function value being integrated against, but it's going
to be inexact.  I see creating interpolants as quite a bit of work with
little payback when the infrastructure to do projection is already in place
and typically is what is required by theorems anyway (e.g. project initial
conditions for parabolic problems into FE space).

On Dec 11, 2007 3:25 AM, Anders Logg <[email-masked]> wrote:

> matrix
> (something
> >     that
> space

---

From: Marie Rognes <[email-masked]>
Date: Tue, 11 Dec 2007 17:16:21 +0100

Robert Kirby wrote:
> L2 projection into the finite element space.  On each element, or 
> globally, solve
> ( pi u - u , w ) = 0 for a function pi u for all of the test functions 
> All you need is point evaluation on the input function u for 
> quadrature, and the weighted values of u become the rhs of the system.
>
> We could include the function value being integrated against, but it's 
> going to be inexact.  I see creating interpolants as quite a bit of 
> work with little payback when the infrastructure to do projection is 
> already in place and typically is what is required by theorems anyway 
> ( e.g. project initial conditions for parabolic problems into FE space).

But say that you want to enforce essential boundary conditions on your 
spaces through replacing the equations in the linear system with the 
values of the degrees of freedom on the boundary.  Then you would need 
to know how to evaluate the degrees of freedom, right...?

---

From: Anders Logg <[email-masked]>
Date: Tue, 11 Dec 2007 21:54:16 +0100

> to know how to evaluate the degrees of freedom, right...?

I think so. For right-hand sides, we could do quadrature as suggested
or we could first evaluate the Lagrange nodes and then do a linear
transform from the Lagrange coefficients to say the BDM coefficients.

But for boundary conditions, we only know the trace of some given
function. Then I don't see how to do either projection (quadrature) or
evaluating Lagrange nodes and then projection. But direct evaluation
of the functionals (like computing the normal component at a point)
would work.

-- 
Anders

---

From: "Robert Kirby" <[email-masked]>
Date: Wed, 12 Dec 2007 13:31:50 -0600

I see -- we need some kind of information for assembly.
"Direct" evaluation of the functionals is not possible except for things
that rely on black-box __call__() -- point evaluation and normal components
are good, integral moments are bad.

Marie, it's fine if you want to take a stab at getting this information in
as you suggested.