← Back to team overview

ffc team mailing list archive

Re: bug in newer versions

 

> On Tue, Oct 30, 2007 at 07:45:56PM +0100, Alessio Quaglino wrote:
>>
>> 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
>
> We need something more informative than that. Try to create a minimal
> 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
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/ffc-dev
>




Follow ups

References