fenics team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jun 01, 2010 at 12:29:32AM +0200, Douglas N Arnold wrote:
> 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).
> >
>
> 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.

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

Attachment: signature.asc
Description: Digital signature