Subject: [UFL-dev] Index renumbering

From: Anders Logg <[email-masked]>
Date: Sat, 21 Feb 2009 21:24:49 +0100

Is index renumbering implemented in UFL? I'd like to renumber all
indices in a form such that the first index starts at zero. I'd also
like to be able to have separate index types/flags for argument basis
functions and other indices.

For Poisson, FFC renumbers from 

  (dXa10/dx0)(dXa11/dx0) | ((d/dXa10)vi-2)*((d/dXa11)vi-1)*dX(0) +
  (dXa12/dx1)(dXa13/dx1) | ((d/dXa12)vi-2)*((d/dXa13)vi-1)*dX(0)
  
to
  
  (dXa0/dx0)(dXa1/dx0) | ((d/dXa0)vi0)*((d/dXa1)vi1)*dX(0) +
  (dXa0/dx1)(dXa1/dx1) | ((d/dXa0)vi0)*((d/dXa1)vi1)*dX(0)

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Sun, 22 Feb 2009 10:35:06 +0100

No, index renumbering isn't implemented.
And basis function indices are nonexisting,
I didn't think we needed them but maybe
you do for the tensor representation.
Lets discuss it next week.

Martin

On Sat, Feb 21, 2009 at 9:24 PM, Anders Logg <[email-masked]> wrote:
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkmgYxEACgkQTuwUCDsYZdEMCACfbxiYdEQrXpr9Ex5axYMqQ/Ja
> PbAAnjPaJkzuhXyEwBxKnDi1peZrKx+Y
> =OyGI
> -----END PGP SIGNATURE-----
>
> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev
>
>

---

From: Anders Logg <[email-masked]>
Date: Sun, 22 Feb 2009 11:15:55 +0100

I thought about it some more and I think I might be able to generate
the indices on the fly as part of the monomial extraction so it might
not be necessary.

So if it's not needed for quadrature it can be part of the monomial
extraction. Otherwise, it seems pretty simple to make a transformer
for index renumbering.

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Sun, 22 Feb 2009 13:11:33 +0100

We need some kind of canonical numbering of indices (and form
arguments) for repr to give a unique ID for the cache.

Martin

On Sun, Feb 22, 2009 at 11:15 AM, Anders Logg <[email-masked]> wrote:
> iEYEARECAAYFAkmhJdsACgkQTuwUCDsYZdEo0gCfXNGRlgpfNBABWgAVLeMexyrW
> ZJkAn2f85OXCiyBjVnNl+CftHHG3irJH
> =/4QY

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 23 Feb 2009 11:11:25 +0100

> extraction. Otherwise, it seems pretty simple to make a transformer

I don't think index renumbering is needed for quadrature. A thing that is
needed, however, is the component number of functions/basisfunctions but is this
alway represented by index?

Kristian

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 23 Feb 2009 11:23:30 +0100

> alway represented by index?

Sorry, 'index' should be Indexed in the above.

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 23 Feb 2009 11:30:10 +0100

> Sorry, 'index' should be Indexed in the above.

This question is too inaccurate for me to know what to say... But I'll try.

A (basis) function never _has_ an index, but expressions
may refer to a component of a function using Indexed.
However, you can also have e.g. sums of Functions
on vector elements, and Indexed outside that.

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 23 Feb 2009 11:44:56 +0100

Quoting Martin Sandve Alnæs <[email-masked]>:

> >> needed, however, is the component number of functions/basisfunctions but
> is
> >> this

Does this means that for

a = V[0].dx(1)*dx

one will have to return something like [V[0].dx(1), V[1].dx(1)] from
BasisFunction and then let Indexed pick the correct value? This looks very
inefficient and it makes it more complicated to know which variable names are
going to be used in the code.

BTW index renumbering breaks for

V[0]*V[1]*dx

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 23 Feb 2009 12:01:50 +0100

> going to be used in the code.

I'm not sure, but if I understand you right, you're talking about using
the transformer framework. As I've suggested before, state can be
kept on a stack, updated by Indexed and SpatialDerivative, and
then basis_function can simply return the current component (0)
and derivative (1).

But as I've said before, I intend to use linearized computational
graphs as an intermediate step for quadrature code generation.
I'm soon getting to updating the graph tools since differentiation
seems to work ok now. (I plan to get back to improving differentiation
later with better algorithms (reverse mode AD), but want to get
something up and running first!)

> BTW index renumbering breaks for

Fixed. (not tested)

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 23 Feb 2009 12:37:16 +0100

> then basis_function can simply return the current component (0)

Yes, and I get all that. I just wanted to be sure that Indexed always referred
to components of basisfunctions. But I see now that Indexed is also used for the
components of ComponentTensor but in this case a second mapping of index values
is needed.

E.g.,

a = inner(grad(V), grad(U))*dx

IndexSum                         
(                                
  IndexSum                     
  (                            
    Product                  
    (                        
      Indexed              
      (                    
        ComponentTensor  
        (                
          Indexed      
          (            
            SpatialDerivative
            (                
              BasisFunction(VectorElement('Lagrange',1, 2), -2)
              MultiIndex((Index(16),))
            )
            MultiIndex((Index(17),))
          )
          MultiIndex((Index(16), Index(17)))
        )
        MultiIndex((Index(18), Index(19)))
      )
[snip]

Here the first Indexed will set the values of index18 and index19, which
ComponentTensor should use to map the values of index16 and index17 such that
BasisFunction can generate the correct component V[index17].dx(index16)

correct?
 
> something up and running first!)

OK, that's understandable.

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 23 Feb 2009 12:59:33 +0100

> components of ComponentTensor but in this case a second mapping of index values

I think you've got it right, but just to be sure: It's a bit more
general than this.

Indexed represents components of <any tensor-valued expression>,
which can also be for example a Sum of other tensor-valued expressions,
a ComponentTensor, or a ListTensor.

ComponentTensor represents a tensor-valued expression
with components given by a scalar expression with free indices
that are "bound" to tensor axes.

In a sense ComponentTensor and Indexed are opposite operations,
Indexed mapping from tensor valued to scalar valued with free indices,
and ComponentTensor mapping from scalar valued with free indices to
tensor valued.

I plan to document this in more detail with mathematical notation...

> BasisFunction can generate the correct component V[index17].dx(index16)

Seems right.

Having Indexed(ComponentTensor(expr, ii), jj) might seem like
additional garbage in the expression tree, since this is basically
represents a relabeling of the free indices ii in expr to jj, but this is
a result of keeping subexpressions unmodified as long as possible.

---

From: Anders Logg <[email-masked]>
Date: Mon, 23 Feb 2009 17:42:44 +0100

> and ComponentTensor mapping from scalar valued with free indices to

So why then are both needed in the current example? I would imagine
Indexed applied to ComponentTensor would cancel out so that we can get

 IndexSum
 (
  IndexSum
  (
    Product
    (
      Indexed
        (
        SpatialDerivative
        (

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 23 Feb 2009 18:15:21 +0100

On Mon, Feb 23, 2009 at 5:42 PM, Anders Logg <[email-masked]> wrote:
> - Show quoted text -

It's generally a good idea to not modify
input expressions to any operators.
Doing so may have side effects in other
algorithms, and can be done at a later
stage anyway. This is also related to the
index renumbering I mentioned eariler
today, so I guess I'll get back to this later.

Martin