Subject: [FFC-dev] bug with tensor fields

From: [email-masked]
Date: Tue, 22 Sep 2009 09:49:38 +0200

The following variational forms should produce the same matrices,

Ai = Function(T, ( ("1.0", "0.0"), ("0.0", "1.0")) )
a1 = dot(dot(Ai, curl(v)), curl(u))*dx
a2 = dot(curl(v), curl(u))*dx

But it doesn't.

I am not sure whether the bug is in Dolfin, FFC or UFL.

The following test script can be used to test the problem.

from dolfin import *

N = 2
mesh = UnitCube(N,N,N)
V = FunctionSpace(mesh, "N1curl", 1)
T = TensorFunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)

Ai = Function(T, ( ("1.0", "0.0"), ("0.0", "1.0")) )
a1 = dot(dot(Ai, curl(v)), curl(u))*dx
a2 = dot(curl(v), curl(u))*dx

A1 = assemble(a1)
A2 = assemble(a2)

file = File("A1.m")
file <<A1
file = File("A2.m")
file <<A2

file_str = """
A1;
A1 = A;
clear A;
A2;
A2 = A;
clear A;

max(max(A1-A2))
"""

f = open("test_tensor.m", 'w')
f.write(file_str)
f.close()

import os
os.system("octave test_tensor.m")

---

From: [email-masked]
Date: Tue, 22 Sep 2009 09:54:36 +0200

Sorry about the early posting. The bug is in compile_function in  Dolfin
(I guess). The generated  tensor is 2x2 instead of 3x3.

Kent

> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 22 Sep 2009 08:57:55 +0100

Tenor fields have not been tested in DOLFIN. See

     https://blueprints.launchpad.net/dolfin/+spec/tensor-functions

Garth

---

From: Johan Hake <[email-masked]>
Date: Tue, 22 Sep 2009 09:58:58 +0200

> Ai = Function(T, ( ("1.0", "0.0"), ("0.0", "1.0")) )

This creates a 2x2 tensor function, do you not want a 3x3 tensor?

I have plans on fixing rigorous tests for value shapes for function 
expressions. Which are automatically checked with the value shape of the 
FunctionSpace. This could prevent your error and some nasty segfaults...

Johan

---

From: [email-masked]
Date: Tue, 22 Sep 2009 10:03:12 +0200

> FunctionSpace. This could prevent your error and some nasty segfaults...

That would be great. I think a simple check on the value shape and
the tuples is sufficient.

Kent

---

From: Johan Hake <[email-masked]>
Date: Tue, 22 Sep 2009 10:22:45 +0200

> the tuples is sufficient.

Yes something like that. What makes it not straightforwardly is that a user 
can define his own Function class:

  class MyFunc(Function):
     cpparg = (("2","3"),("3","4"))

Then I need to create some value_shape function that is used when this 
function is instantiated with a FunctionSpace.

It is not nessesary to have this for pure Python expressions, like:

  class MyFunc(Function):
       def eval(values,x):
          .....

Here will the callback typemap produce the correct value shape on the NumPy 
array that is passed to the user, with all its index checks.

Johan