Subject: [FFC-dev] bug with triangles?

From: Alessio Quaglino <[email-masked]>
Date: Fri, 28 Nov 2008 09:13:25 +0100

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?
Alessio

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 28 Nov 2008 08:41:15 +0000

> 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
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev

---

From: Alessio Quaglino <[email-masked]>
Date: Fri, 28 Nov 2008 10:11:23 +0100

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 [[email-masked]]
Inviato: venerdì 28 novembre 2008 9.41
A: Alessio Quaglino
Cc: [email-masked]
Oggetto: Re: [FFC-dev] bug with triangles?

   q = 2

Garth

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 28 Nov 2008 09:25:52 +0000

Alessio Quaglino wrote:
> 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? 

It's not a limitation and it applies to all spatial dimensions. The 
gradient is an operator that is only defined for  tensors of dimension n 
in a space of dimension n. The gradient of a function with three 
components in two dimensions is not defined.

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.
>

We should add a test that would give an error for your second piece of 
code because the gradient isn't defined for the type of function that 
it's being given.

You can still take derivatives of your function with three components, 
but you'll have to to it using .dx(0), etc or by defining you own 
grad-like operator.

Garth

> Alessio
> 
> 
> ________________________________________ Da: Garth N. Wells
> [[email-masked]] Inviato: venerdì 28 novembre 2008 9.41 A: Alessio
> Quaglino Cc: [email-masked] Oggetto: Re: [FFC-dev] bug with
> triangles?
> 
> Alessio Quaglino wrote:
>> 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) )
>> 
>> 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 )
>> 
> 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 [email-masked] 
> _______________________________________________ FFC-dev mailing list 
> [email-masked] http://www.fenics.org/mailman/listinfo/ffc-dev

---

From: Alessio Quaglino <[email-masked]>
Date: Fri, 28 Nov 2008 10:42:38 +0100

As far as I know the gradient is an opperation on scalar fields giving a vector, while the jacobian is on vector fields giving a tensor. There's no assumption on the dimension of the vector (which you have with the divergence). I then guess page 42 of ffc documentation must be corrected, since there is stated that grad() applied to vectors gives the jacobian matrix.

Alessio

________________________________________
Da: Garth N. Wells [[email-masked]]
Inviato: venerdì 28 novembre 2008 10.25
A: Alessio Quaglino
Cc: [email-masked]
Oggetto: Re: [FFC-dev] R:  bug with triangles?

Garth

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 28 Nov 2008 11:15:10 +0000

Alessio Quaglino wrote:
> As far as I know the gradient is an opperation on scalar fields
> giving a vector, while the jacobian is on vector fields giving a
> tensor. There's no assumption on the dimension of the vector (which
> you have with the divergence). I then guess page 42 of ffc
> documentation must be corrected, since there is stated that grad()
> applied to vectors gives the jacobian matrix.
> 

I had a look at the FFC code and it will indeed compute the Jacobian 
when when you ask for the gradient. The problems arises when you ask FFC 
to compute the inner product. FFC will only compute the inner product 
for scalars, equal length vetors and for square rank-2 tensors. There is 
no check in the FFC code that the tensor is square.

I'll add a test soon which will return a sensible error message.

Garth

> [[email-masked]] Inviato: venerdì 28 novembre 2008 10.25 A: Alessio
> Quaglino Cc: [email-masked] Oggetto: Re: [FFC-dev] R:  bug with
> triangles?
> 
> Alessio Quaglino wrote:
>> 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?
> 
> gradient is an operator that is only defined for  tensors of
> dimension n in a space of dimension n. The gradient of a function
> with three components in two dimensions is not defined.
> 
>> matrices and converts them to vectors, and I don't get any "index
>> out of range" error.
>> 
> 
> We should add a test that would give an error for your second piece
> of code because the gradient isn't defined for the type of function
> that it's being given.
> 
> You can still take derivatives of your function with three
> components, but you'll have to to it using .dx(0), etc or by defining
> you own grad-like operator.
> 
> Garth
> 
>> Alessio
>> 
>> 
>>> (psi, phi) = TestFunctions(Mh) (b, w) = TrialFunctions(Mh) k =
>>> dot( grad(phi),grad(b) )
>>> 
>>> 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
>>> )
>>> 
>> VectorElement("Lagrange", "triangle", q - 1, 2) Mh = EL1 + EL2
>> _______________________________________________ FFC-dev mailing
>> list [email-masked]