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

From: Douglas N Arnold <[email-masked]>
Date: Tue, 01 Jun 2010 00:29:32 +0200

On Wed, May 05, 2010 at 10:42:32AM +0100, Anders Logg wrote:
 >
> On Wed, May 05, 2010 at 09:05:19AM +0100, Garth N. Wells wrote:
>>
>>
>> On 04/05/10 20:32, Anders Logg wrote:
>> >On Mon, May 03, 2010 at 11:08:13AM +0200, Kent Andre wrote:
>> >>On Mon, 2010-05-03 at 09:36 +0100, Garth N. Wells wrote:
>> >>>
>> >>>On 03/05/10 08:48, Kent Andre wrote:
>> >>>>
>> >>>>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.
>> >>>>
>> >>>
>> >>>They do need to worry about it on the DOLFIN side with the creation of a
>> >>>Form:
>> >>>
>> >>>     Poissin::Bilinearform a(V0, V1);
>> >>>
>> >>>Garth
>> >>>
>> >>
>> >>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 :)
>> >
>> >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?
>> >
>>
>> (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
>
> 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

I am coming to this discussion rather late, but I cast
a vote for a(u,v) where u is the trial function and v is
the test function.  It does not make much difference
to a sophisticated user, but I think to students it does.
Almost all of the books and papers on the mathematical
side of finite elements use the following notation:

   Galerkin's method:
   Find u in Vh such that  a(u,v) = F(v)  for all v in Vh.

   Stiffness matrix:
   A_ij = a(phi_j,phi_i)

   Load vector: b_i = F(phi_i)

Since a big selling point of FEniCS is that you enter
a bilinear form just as you would write it in a paper
or read it in a book, I find this a fairly strong argument
in favor of changing.  Whether it is worth the effort
and the side-effects I cannot judge.

   -- Doug

---

From: Anders Logg <[email-masked]>
Date: Tue, 1 Jun 2010 08:23:11 +0200

> and the side-effects I cannot judge.

I see the point. I often find myself in the position where I need to
convert others (like my students) to a(v, u). They will not be happy
if we change back to a(u, v) but that's my problem (and theirs). :-)

But what about the transpose?

  A_ij = a(phi_j, phi_i)

Why is there a mismatch? Is there any argument one could make which
one is more "correct", either the standard bilinear form notation
a(u, v) or the standard matrix notation A_ij?

--
Anders

---

From: Anders Logg <[email-masked]>
Date: Tue, 1 Jun 2010 16:14:44 +0200

Here's a summary of opinions (in order of appearance):

  Anders: ambivalent
  Garth: neutral, maybe a(u, v)
  Hans Petter: a(u, v)
  Kent: a(u, v)
  Marie: a(v, u)
  Patrick: a(v, u)
  Doug: a(u, v)

More thoughts/opinions?

--
Anders

---

From: Anders Logg <[email-masked]>
Date: Tue, 1 Jun 2010 18:01:47 +0200

ok, so one more for a(u, v). It's not surprising or you would need to
rewrite that book. ;-)

Current results:

  Anders: ambivalent
  Ridg: a(u, v)

--
Anders

On Tue, Jun 01, 2010 at 07:40:26AM -0700, Ridgway Scott wrote:
> I am formally ambivalent in the sense that it is just notation and thus
> there cannot be any mathematical argument one way or the other.
> But I support Doug's analysis: a(u,v) is the way it is commonly used
> even though it is backwards with respect to the standard matrix
> numbering. So I have to say: why change? So a(u,v) for me.
>
> Ridg
>
> On Jun 1, 2010, at 7:14 AM, Anders Logg wrote:
>

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Jun 2010 11:25:46 +0200

Perhaps we don't have an "overwhelming majority" but the indication
seems pretty clear. I've added this as a blueprint for 0.9.9:

https://blueprints.launchpad.net/dolfin/+spec/order-of-test-and-trial-functions

--
Anders

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