ffc team mailing list archive

[Alessio Quaglino <alessio@xxxxxx>]

Thank you for the reply.
So you mean that such limitation in present only for triangles, while for tets you can take gradients of n-vectors? Then I don't understand, also because the workaround I used in the code I posted treats grad(phi) and grad(b) as 3x2 matrices and converts them to vectors, and I don't get any "index out of range" error.

Alessio 


________________________________________
Da: Garth N. Wells [gnw20@xxxxxxxxx]
Inviato: venerdì 28 novembre 2008 9.41
A: Alessio Quaglino
Cc: ffc-dev@xxxxxxxxxx
Oggetto: Re: [FFC-dev] bug with triangles?

Alessio Quaglino wrote:
> I have the following code:
>
>   q = 2
>   EL1 = VectorElement("Lagrange", "triangle", q, 3)
>   EL2  = VectorElement("Lagrange", "triangle", q - 1, 3)
>   Mh = EL1 + EL2
>   (psi, phi) = TestFunctions(Mh)
>   (b, w) = TrialFunctions(Mh)
>   k = dot( grad(phi),grad(b) )
>
> which gives me the error:
>
>   k = dot( grad(phi),grad(b) )
>   File "/Library/Python/2.5/site-packages/ffc/compiler/language/operators.py", line 102, in dot
>   form = form + v[i][j]*w[i][j]
>   IndexError: list index out of range
>
> I don't think I hven't understood what I supposed to give to ffc, since if I define the following I get the desired behavior:
>
>   def vect(a):
>     return [ a[i][j] for j in range(2) for i in range(3) ]
>   gphi = vect( grad(phi) )
>   gb = vect( grad(b) )
>   k = dot( gphi,gb )
>
> Is there a problem in handling (n x d) matrices with 2D meshes?

The problem is that FFC doesn't know how to take the gradient of a
vector function of length 3 on a 2D element (triangle). The below works

   q = 2
   EL1 = VectorElement("Lagrange", "triangle", q, 2)
   EL2  = VectorElement("Lagrange", "triangle", q - 1, 2)
   Mh = EL1 + EL2
   (psi, phi) = TestFunctions(Mh)
   (b, w) = TrialFunctions(Mh)
   k = dot( grad(phi),grad(b) )

The '2' at the end of the element definition is redundant in this case
since FFC will default to 2 for 2D elements.

Garth

> Alessio
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/ffc-dev