Subject: Re: [Ufl] [Dolfin] [Question #152702]: indices in an expression?

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 13 Apr 2011 19:24:33 +0100

The permutation symbol should be added to UFL.

Garth

On 13/04/11 19:20, Martin Sandve Alnæs wrote:
> Python does not know anything about implicit index summation, so in your
> code e() is only called once, and reaches the else condition and thus
> returns 0.0.
> 
> Similar to as_vector, you can use the as_tensor function, since you need
> three dimensions.
> 
> Martin
> 
> On 13 April 2011 19:04, B. Emek Abali
> <[email-masked]
> <mailto:[email-masked]>> wrote:
> 
>     New question #152702 on DOLFIN:
>     https://answers.launchpad.net/dolfin/+question/152702
> 
>     May be a non-needed trial for some, but I try to define an
>     expression (permutation symbol)
> 
>     def e(i,j,k):
>            if (i,j,k) == (1.0,2.0,3.0): return 1.0
>            elif (i,j,k) == (2.0,3.0,1.0): return 1.0
>            elif (i,j,k) == (3.0,1.0,2.0): return 1.0
>            elif (i,j,k) == (3.0,2.0,1.0): return -1.0
>            elif (i,j,k) == (2.0,1.0,3.0): return -1.0
>            elif (i,j,k) == (1.0,3.0,2.0): return -1.0
>            else: return 0.0
> 
>     and use (with summation over indices, it is a curl operator
>     contracted with w) in the form like
> 
>     f=v[j]*w[i]*e(i,j,l)*v[l]*dx
> 
>     which raises:
>     ufl.log.UFLException: Trying to integrate expression of rank 0 with
>     free indices
> 
>     Do I have to use as_vector(...) type definition to sum with indices
>     or is there a nice possibility to let it sum over the arguments of
>     e(...) expression?
> 
>     --
>     You received this question notification because you are a member of
>     DOLFIN Team, which is an answer contact for DOLFIN.
> 
>     _______________________________________________
>     Mailing list: https://launchpad.net/~dolfin
>     Post to     : [email-masked]
>     <mailto:[email-masked]>
>     Unsubscribe : https://launchpad.net/~dolfin
>     More help   : https://help.launchpad.net/ListHelp
> 
> 
> 
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Wed, 13 Apr 2011 20:31:22 +0200

That should be fairly easy to do, can be handled similarly to Identity.

I = Identity(cell.d)
f = I[i,j] * ...

e = Permutation(cell.d)
f = e[i,j,l] * v[j]*w[i]*v[l]*dx

Other suggestions for naming?

Martin

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 13 Apr 2011 19:38:19 +0100

> That should be fairly easy to do, can be handled similarly to Identity.

Yes - I had a look in the code this afternoon because I thought that it
may have already been implemented somewhere close to Identity.

> Other suggestions for naming?

Maybe

   e = PermutationSymbol()

? 'Permutation' sounds a bit generic.

I don't think that it needs the cell dimension since it's only defined
for d = 3.

Garth

> Martin
> 
> On 13 April 2011 20:24, Garth N. Wells <[email-masked]
> <mailto:[email-masked]>> wrote:
> 
>     <mailto:[email-masked]>
>     > <mailto:[email-masked]
>     <mailto:[email-masked]>>> wrote:
>     >
>     >     <mailto:[email-masked]
>     <mailto:[email-masked]>>

---

From: Douglas N Arnold <[email-masked]>
Date: Wed, 13 Apr 2011 13:57:33 -0500

Why only three dimensions?  That seems arbitrary. The mathematical 
definition is of course

  epsilon(i_1,...,i_n) =  1, -1, or 0

according to whether (i_1,...,i_n) is an even permutation
of (1,...,n), an odd permutation of (1,...,n), or neither.
It is defined in all dimensions n\ge 1, and is very useful
for all n>1.  The 2-d case is often used in describing finite
element methods for elasticity with weakly imposed symmetry.

The wikipedia article on this is "Levi-Civita symbol",
and the synonyms for the name are "permutation symbol", "antisymmetric
symbol", or "alternating symbol".  I would think

  e = PermutationSymbol()   or  e = PermutationSign()

would be relativitely intuitive.

On 04/13/2011 01:38 PM, Garth N. Wells wrote:
>> On 13 April 2011 20:24, Garth N. Wells<[email-masked]
> Mailing list: https://launchpad.net/~ufl
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ufl

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Thu, 14 Apr 2011 10:10:22 +0200

It's working, from the test:

Note that e here will be a 3x3x3 tensor which must be indexed to get eps_ijk
        e = PermutationSymbol(3)
        self.assertEqual(e.shape(), (3, 3, 3))

Evaluation with int indices produce correct values for eps_ijk:
        for i in range(3):
            for j in range(3):
                for k in range(3):
                    value = (j-i)*(k-i)*(k-j)/2
                    self.assertEqual(e[i, j, k], value)

If e is indexed with free indices, we get an Indexed object as usual in ufl
        i, j, k = indices(3)
        self.assertIsInstance(e[i,j,k], Indexed)

and evaluating the product eps_ijk eps_ijk with implicit summation gives the
right sum
        x = (0,0,0)
        self.assertEqual((e[i,j,k] * e[i,j,k])(x), 6)

I haven't tested with n != 3, but the formula Garth copied should handle
that.

Martin

On 14 April 2011 09:19, Mikael Mortensen <[email-masked]> wrote:

> The source of truth suggests Levi-Civita, alternating or permutation
> symbol: http://en.wikipedia.org/wiki/Levi-Civita_symbol
>
> Mikael
>
>
> On 13 April 2011 20:31, Martin Sandve Alnæs <[email-masked]> wrote:
>

---

From: Mikael Mortensen <[email-masked]>
Date: Thu, 14 Apr 2011 09:19:03 +0200

The source of truth suggests Levi-Civita, alternating or permutation symbol:
http://en.wikipedia.org/wiki/Levi-Civita_symbol

Mikael