ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

/Anders