Subject: Re: [FFC-dev] Suggestions for changes [...]

From: <[email-masked]>
Date: Tue, 24 May 2005 14:08:34 -1100

>A Function in FFC is always a member of a finite element 
space, and
>each Function declared in FFC is always connected to a 
finite element
>Function in DOLFIN.
>
>The FFC Function is the restriction of a DOLFIN Function to 
a local
>element.

ok.  I was unaware of this assumption.  Sundance and PETSc 
both allow functions that are not in the finite element 
space.  But FFC only allows things that are strictly 
multilinear, so it requires the mapping into the space.

>Of course, one could change the syntax from
>
>  f = Function(element)
>  L = v*f*dx
>
>to
>
>  f  = Function()
>  pi = Projection(element)
>  L  = v*pi(f)*dx
>
>but that would only make things more complicated to write 
(but perhaps
>it would be more clear what happens).
>

It is more to type, but more clear.  It opens two nice 
avenues:
i.) Adding quadrature-based modules that will *not* do the 
projection explicitly
ii.) Allowing multiple kinds of projections (L2, H1, nodal 
interpolation)

>What I meant was that in some cases, we want to project a 
Function
>locally from one finite element space to another. Take as 
an example
>the stabilization term of Navier-Stokes with linear 
elements on tets:
Good point, I hadn't thought about this.  Another example 
would be any kind of multigrid or multilevel method.

>
>And D(i) would look very much like a gradient:
>
>  D(i)(u)*D(i)(v)
>
I'n not sure what you mean -- this is the dot product of 
gradients of scalar functions.

>I need to think about this. At some point, we should add 
grad,
>div, rot anyway.
I would add these in terms of what you've already done by 
simple definitions, e.g.

grad(u) --> u.dx(i)

But maybe it's two lines -- one if u is a scalar, and 
another if it's a vector

grad(u) --> u[i].dx(j) if u is a vector
grad(u) --> u.dx[i] if u is a scalar.

Rot/curl is a little trickier, as you need to know whether 
you're in 2d or 3d.

Rob

---

From: Anders Logg <[email-masked]>
Date: Mon, 30 May 2005 23:08:50 -0500

> ii.) Allowing multiple kinds of projections (L2, H1, nodal 

I need to think this through. At some point we probably need to add an
option for quadrature, but there may be other ways to handle it, for
example by declaring dx differently. We currently have

dx = Integral("interior")
...
a = u.dx(i)*v.dx(i)*dx

The declaration of dx is normally not visible (it's added by the
simple FFC parser that preprocesses a given .form file and adds the
necessary declarations on top of the file).

For quadrature, one could  do

dx = Integral("interior quadrature")
...
a = u.dx(i)*v.dx(i)*dx

Different terms can be integrated differently (quadrature or tensor)
and combined in the same form.

> Rot/curl is a little trickier, as you need to know whether 

Good idea. I also need to add a scalar product (could overload the ','
operator) that translates

    (grad(u), grad(v))

into

    u.dx(i)*v.dx(i)

Note that the index needs to be the same in both u.dx(i) and v.dx(i)
which means that grad(u) can probably not be mapped to u.dx(i), but
rather to a function that takes the index as an argument (and ',' will
give the same index to both u.dx() and v.dx()).

/Anders

---

From: Robert C. Kirby <[email-masked]>
Date: Wed, 1 Jun 2005 09:20:13 -0500

> option for quadrature, but there may be other ways to handle it, for

This is also possible to allow further experimentation.

> and combined in the same form.

This is a good idea -- I meant by projection the L2 projection
(Pi u,v) = (u,v) or the (local) H1 projection
(grad( Pi u ) , grad( v ) ) = ( grad( u ) , grad( v ) )

as opposed to nodal interpolation.

Of course, opening up the possibility of varying quadrature is an 
orthogonal dimension to these.

> operator) that translates

I wouldn't overload the comma operator since that is used in duality 
pairings and FFC works on integrands / measure.  I don't think you even 
can overload , in Python.

> give the same index to both u.dx() and v.dx()).

You can introduce a class "grad" whose constructor takes a function, 
and overload member __mul__ to take another "grad" instance and return 
u.dx(i) * v.dx(i) for some appropriate index i.  See above comment on 
",".

Rob

---

From: Anders Logg <[email-masked]>
Date: Wed, 1 Jun 2005 13:07:06 -0500

> can overload , in Python.

You're probably right. (But you can in C++...)

> u.dx(i) * v.dx(i) for some appropriate index i.  See above comment on 

Good idea, but I would probably create a class Gradient that overloads
'*' to create a new Index i and return u.dx(i) * v.dx(i) for that
Index. In addition, I will make a function grad() that returns a
Gradient.

/Anders