Subject: [FFC-dev] bug in newer versions

From: "Alessio Quaglino" <[email-masked]>
Date: Tue, 30 Oct 2007 19:45:56 +0100 (CET)

I've been debugging my code for several days without understanding why I
got I wrong solution at home on my laptop and at the office. I finally
discovered that the problem is not present in ffc 0.4.1 but only appears
in 0.4.3 (and maybe 0.4.2). I've just cloned the last repository and the
issue is still present.

I actually don't know what the bug is, but it's the kind of error I would
expect if the DiscreteFunction (i.e. the Function Foo I modify using
Foo.vector().vec() in my code) of my form is not evaluated correcly. I can
tell this because is the kind of expected wrong solution obtained without
the corrective terms that my DiscreteFunction represents.

Hope it helps,
Alessio Quaglino

---

From: Anders Logg <[email-masked]>
Date: Tue, 30 Oct 2007 22:32:51 +0100

> the corrective terms that my DiscreteFunction represents.

We need something more informative than that. Try to create a minimal
example that can reproduce this error and post it to the list.

Thanks
/Anders

---

From: "Alessio Quaglino" <[email-masked]>
Date: Wed, 31 Oct 2007 10:00:43 +0100 (CET)

> example that can reproduce this error and post it to the list.

Don't know about the minimal, but this is mine which is still quite simple:

element1 = VectorElement("Lagrange", "tetrahedron", 1)
element2 = VectorElement("Discontinuous Lagrange", "tetrahedron", 0, 6)
element4 = FiniteElement("Discontinuous Lagrange", "tetrahedron", 0)

xvn		= 	TrialFunction(element1)
psi	   	= 	TestFunction(element2)
sigma0      	= 	Function(element2)
lmbda     	= 	Function(element4) # Lame coefficient
mu        	= 	Function(element4) # Lame coefficient
dt_el4        	= 	Function(element4)

d = len(v)

def epsilon(u):
    return 0.5 * (grad(u) + transp(grad(u)))

def E(e, lmbda, mu):
    Ee = 2.0 * mult(mu, e) + mult(lmbda, mult(trace(e), Identity(d)))
    return Ee

def tomatrix(q):
    return [ [q[i + j + i * j - 4 * (j == 2) * (i == 2)] for i in
range(3)] for j in range(3) ]

psimatrix = tomatrix(psi)
sigma0matrix = tomatrix(sigma0)

aelast 	= dt_el4 * dot( E(epsilon(xvn), lmbda, mu) , psimatrix )
arot = dt_el4 * dot( mult(grad(xvn), sigma0matrix) + mult(sigma0matrix,
transp(grad(xvn))) , psimatrix )

a = ( arot - aelast ) * dx

It is correct with 0.4.1 and not with 0.4.2 or later.

Alessio

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

---

From: Anders Logg <[email-masked]>
Date: Wed, 31 Oct 2007 10:31:44 +0100

> It is correct with 0.4.1 and not with 0.4.2 or later.

It's still too much (if you want me to debug it). If there is an error
somewhere, there must be at least one term that is wrong. Try to get
it down to 2 or 3 (short) lines of code and then I can take a look.

/Anders

---

From: "Alessio Quaglino" <[email-masked]>
Date: Wed, 31 Oct 2007 10:54:20 +0100 (CET)

>> >> Foo.vector().vec() in my code) of my form is not evaluated correcly.
>> I
>> >> can

Well the thing is I don't know how to isolate it. I have 5 other forms in
my code, all working fine even if they do similar operations. The error
shows up only here and the only additional term is:

the issue cannot come from tomatrix() because I use it in other forms. I
use also grad(), transp(), mult() and dot() somewhere else, but I use a
Function which is not constant on the mesh only here, sigma0. In fact, in
accordance to the theory for some cases this additional term does not
contribute to the solution and the bug does not even show up. I'll try to
think more about how to make a simpler example.

Alessio