Subject: [Fenics] The order of test and trial functions

From: Anders Logg <[email-masked]>
Date: Thu, 29 Apr 2010 12:03:03 +0200

As you are all aware, FEniCS uses a somewhat unorthodox convention for
the order of test and trial function, which is to place the test
function first:

  a(v, u)

At some point in history, u was chosen as a name for the unknown and v
was chosen as a name for the test function. At another point in
history, i was chosen as the name of the row index in a matrix and j
was chosen as a column index. So people tend to write

  a(u, v)
  A_ij

This is somewhat unfortunate, since i (rows, equations) corresponds to
the test function v and j corresponds to the trial function u.

So we must define the matrix A by

  A_ij = a(phi_j, phi_i)

Notice the swapping of i and j here. This might not be much of a
problem but it's certainly a problem when we consider multilinear
forms and their corresponding tensors (as we do frequently in FEniCS):

  A_i = a(..., phi_i3, phi_i2, phi_i1)

For this reason, Rob and I decided to change the order of u and v and
write

  a(v, u)

when we wrote the first FFC paper back in 2005. Rob is not so much to
blame for this, but I've been pushing this quite hard since then with
the result that v is before u in the current draft version of the
FEniCS book, as well as in all demos.

The question is now, should we stick by this? The reason for doing so
are:

  1. It the logic thing: v <--> i and u <--> j

  2. UFL, UFC and the DOLFIN assembler are all based on this

The reasons for not doing so are:

  1. a(u, v) looks better (alphabetical ordering)

  2. It's what most people are used to

If we choose this second option, we need to define the order of
arguments in a form to start from the right in all text that relates
to finite element assembly (in the FEniCS book) and in all code that
deals with forms (UFL, UFC, DOLFIN). Note that this means we don't
need to change any code at all, but the text would need changing.
In particular, the UFC function create_finite_element(0) would still
return the finite element for the test space so that the assembler
would assemble A and not A transposed.

So for a bilinear form

  a(u, v)

we would need to say that v is the first argument (numbered 0) and u
would be the second argument (numbered 1).

And for a semilinear form (see FEniCS book chapter on forms), we would
need to write

  F(v_r, ..., v_3, v_2, v_1; w_1, w_2, w_3, ..., w_n)

where {v_i} are arguments (basis functions) and {w_i} are
coefficients.

This is a bit awkward (and the reason for changing the order) in the
first place, but this would only appear deep down in the more
technical parts of the text and so most users would not see it. They
will only see the swapping of i and j when defining the matrix for
Poisson.

I hope we can have a long and good discussion on this so we settle
this once and for all. All comments/votes are welcome.

--
Anders

---

From: Hans Petter Langtangen <[email-masked]>
Date: Thu, 29 Apr 2010 12:52:25 +0200

I am probably the one to blame for introducing the confusion of two
notations in a FEniCS context. In the tutorial (to appear first in the
FEniCS book), I initially adopted the a(v,u) notation, simply because
that was the notation used in the book's intro chapter to FEM and
because I immediately realized the logical arguments that Anders
states for this choice.

However, the tutorial was used in a FEM course by Doug Arnold, and he
and his students complained about a(v,u). Most books write a(u,v), and
it was confusing to see that even the very first FEniCS example used a
different notation.  There will be many newcomers to FEniCS that are
also newcomers to FEM.  I think these users should be given the most
widespread and traditional notation so that the first experience with
FEniCS maximizes good vibrations. This argument was strong enough to
make me shift from a(v,u) to a(u,v).

The chapters in the FEniCS book where the a(v,u) notation is
preferable can use this notation, with a little argument the first
time the notation appears. Obviously, consistency throughout a book is
always desirable, but in my opinion this argument is weaker than
pleasing newcomers to FEniCS and minimizing potential confusion during
the first 30 min trial of the software. Remember that many (most?)
people will stop there and not have time to become heavy users, but
they will then have a very good impression of the software - this is
extremely important for widespread use and recognition in the long run.

Hans Petter

> _______________________________________________
> Mailing list: https://launchpad.net/~fenics
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~fenics
> More help   : https://help.launchpad.net/ListHelp

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 29 Apr 2010 14:19:44 +0100

On 29/04/10 11:03, Anders Logg wrote:
>    1. It the logic thing: v<-->  i and u<-->  j
>

I don't feel particularly strongly either way, but have a slight leaning 
to go with a(u, v). Both a(u, v) and a(v, u) are commonly used, and 
which is more common depends on the field.

Garth

---

From: Anders Logg <[email-masked]>
Date: Thu, 29 Apr 2010 15:41:34 +0200

On Thu, Apr 29, 2010 at 02:19:44PM +0100, Garth N. Wells wrote:
> I don't feel particularly strongly either way, but have a slight
> leaning to go with a(u, v). Both a(u, v) and a(v, u) are commonly
> used, and which is more common depends on the field.

Where is a(v, u) used? I only hear comments that a(u, v) is used
*everywhere* else and we need to change.

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 29 Apr 2010 14:53:16 +0100

> *everywhere* else and we need to change.

Here are two 'a(v, u)' books from a quick scan of my book shelves

   - Hughes, TJR, The Finite Element Method
   - Donea and Huerta, Finite Element Methods for Flow Problems

Engineering numerical analysis and mechanics book often use a(v, u), I 
assume because when considering the weight function as a 'variation', 
it's common to use the order

    \delta u u

Garth

---

From: Anders Logg <[email-masked]>
Date: Thu, 29 Apr 2010 16:04:51 +0200

On Thu, Apr 29, 2010 at 02:53:16PM +0100, Garth N. Wells wrote:
> Engineering numerical analysis and mechanics book often use a(v, u),
> I assume because when considering the weight function as a
> 'variation', it's common to use the order
>
>    \delta u u
>
> Garth

ok, so they did the right thing... :-)

Another thought is this: u should come before v since u is more
important. Same goes for matrices, columns before rows since a matrix
is really a list of column vectors (the range of the matrix is
spanned by its columns), so we should really be changing to A_ji to
match a(u, v)...

--
Anders

---

From: Anders Logg <[email-masked]>
Date: Mon, 3 May 2010 09:23:41 +0200

There doesn't seem to be any overwhelming pressure to change the
notation from a(v, u) to a(u, v). More opinions?

In summary, the options are:

1. Keep a(v, u) so that A_ij = a(phi_i, phi_j)

2. Change to a(u, v) and do some kind of swapping of indices

3. Change to a(u, v) and change to A_ij = a(phi_i, phi_j)

Option 3 seems unreasonable (although one could argue this is where
the problem lies, columns should be numbered first).

--
Anders

---

From: Kent Andre <[email-masked]>
Date: Mon, 03 May 2010 09:48:10 +0200

I am in favour of a(u,v), simply because it is the way I am used to 
read and write it. Also, very few users need to care about this due to
the high-level syntax in UFL which implicitly takes care of the
ordering.

Kent

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 03 May 2010 09:36:32 +0100

> the high-level syntax in UFL which implicitly takes care of the

They do need to worry about it on the DOLFIN side with the creation of a 
Form:

    Poissin::Bilinearform a(V0, V1);

Garth

---

From: Kent Andre <[email-masked]>
Date: Mon, 03 May 2010 11:08:13 +0200

>     Poissin::Bilinearform a(V0, V1);

Yes, of course, sorry, in C++ this is important,  

Anyway, I think both alternatives have their advantages. 
But, as Godel put it, you must choose between consistency and
conformity :)  

Kent

---

From: Anders Logg <[email-masked]>
Date: Tue, 4 May 2010 21:32:39 +0200

> But, as Godel put it, you must choose between consistency and

I've been thinking some more, and I'm getting more convinced that the
problem is really the ordering of rows and columns in matrices, which
should be (column, row) rather than (row, column) because a matrix is
really a set of columns that span the range of A, so A_i should be
column i and A_ij should be element j of that column.

So then a(u, v) would be consistent with *that* notation but it would
be a pain for us to have an unorthodox matrix index notation.

Which leaves us with the following options:

1. Keep a(v, u) as we have now and not change anything.

2. Change to a(u, v) and comment on the need for swapping the indices
when we define the matrix (tensor) in the text (book and manual).

Alternative (2) would imply the following changes:

(i) Redefine the UFC interface so that the test space is the last
space, not the first.

(ii) Swapping of indices inside Assembler.cpp.

Further comments?

--
Anders

---

From: Marie Rognes <[email-masked]>
Date: Tue, 04 May 2010 21:42:46 +0200

> (ii) Swapping of indices inside Assembler.cpp.

Vote in favor of (1).

--
Marie

> --
> Anders
>   
> ------------------------------------------------------------------------
>

---

From: Hans Petter Langtangen <[email-masked]>
Date: Wed, 5 May 2010 08:26:52 +0200

> >(ii) Swapping of indices inside Assembler.cpp.

I also vote in favor of (1) for the source code and technical documentation.

However, I still repeat my arguments for a(u, v) in the FEniCS tutorial,
at least in the beginning. We may well introduce a(v, u) later in the tutorial
and explain the u,v-v,u problem, or postpone that to the introductory chapter
on the finite element method. 

Regardless of what is mostly used in the FEniCS book and technical
documentation, there will probably be a lot of application
contributions over the next years where people write a(u,v).
Therefore I think we should have a small section in the book explaining
the issue of a(u,v) vs a(v,u), so we have a reference for resolving
potential confusion.

Hans Petter

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 05 May 2010 09:05:19 +0100

> (ii) Swapping of indices inside Assembler.cpp.

(1) sounds best/easiest.

In terms of user code, it's not such a big deal since it only becomes 
important when declaring in C++ a bilinear form which uses different 
spaces for the test and trial function.

Garth

---

From: Anders Logg <[email-masked]>
Date: Wed, 5 May 2010 10:42:32 +0200

On Wed, May 05, 2010 at 09:05:19AM +0100, Garth N. Wells wrote:
> In terms of user code, it's not such a big deal since it only
> becomes important when declaring in C++ a bilinear form which uses
> different spaces for the test and trial function.
>
> Garth

But my biggest concern is consistency. There are many pages in the
book that deal with representation of forms, linearization of forms,
assembly of forms; a(v, u), F(v; u) appear in I don't know how many
places. I think it would be very confusing to write a(u, v) and then
change to a(v, u) at some arbitrary point in the text. It would then
not be enough to comment on the issue of the order of v and u in one
place, but we would need to do it in many places.

I'm open for (2), especially if my conclusions of the order of indices
in matrices is correct. The change we would need to make in the text
is to define

  A_i^T = a(\phi_{i_1}, \phi_{i_2}, \phi_{i_3}, ...)

that is, we add a transpose to the tensor A (which we would take to
mean reverse the indices).

We would also need to deal with the transpose operation in the
assembler (easy fix) and redefine the UFC interface (harder).

--
Anders

---

From: Patrick Riesen <[email-masked]>
Date: Wed, 5 May 2010 17:51:52 +0200

>   2. UFL, UFC and the DOLFIN assembler are all based on this

i don't know if my argumentation does help, but i think you do a better 
job with a(v, u) since
all the books on fem and variational/weak formulations i came across 
write something like,

'...multiply by a testfunction (v in our case) and integrat over the 
domain ....'

and what they do is put the testfunction v in front of your term to 
integrate over when forming the linear, bilinear or even multilinear 
forms, so the most stuff i came across looks always something like

\int_G v:u dG  or int_G v:f dG etc.

in the full write out. but the 'short hand' notation, using a(_,_) it is 
reversed, the a(u, v) way. i don't know why they do this, but so it 
appears to me that you (fenics) are consistent and the others are not.
i vote for consistency and logics, thus i like your writing with a(v, 
u). that's the way how to do it and not care about 'what most people are 
used to'.
i always have the most problems figuring out some mathematics in paper 
or book, because in there, the people tend to use 'what most people are 
used to' and just leave some definitions or obvious things (which appear 
totally non-obvious to me, especially if you do such stuff for the first 
time) away and the like. that way it takes ++++time to figure out what 
they actually did.  go for logics. use a(v, u) as you have.
best,
patrick