ffc team mailing list archive

["Alessio Quaglino" <alessio@xxxxxx>]

> On Wed, Oct 31, 2007 at 10:00:43AM +0100, Alessio Quaglino wrote:
>> > 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
>
> 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.

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:

arot = dt_el4 * dot( mult(grad(xvn), sigma0matrix) + mult(sigma0matrix,
transp(grad(xvn))) , psimatrix )

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


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