Subject: [Ffc] Bug in quadrature code?

From: Anders Logg <[email-masked]>
Date: Fri, 29 Jan 2010 21:42:17 +0100

I've tracked down the error in the DOLFIN fem unit test to the
following simple test:

  from dolfin import *

  mesh = UnitSquare(1, 1)
  V = FunctionSpace(mesh, "DG", 1)

  v = TestFunction(V)
  u = TrialFunction(V)
  h = CellSize(mesh)

  a = (1/h)*v*u*ds

  A = assemble(a)
  info(A, True)

With the old FFC, the result is

row 0: (0, 0.235702)  (1, 0.117851)  (2, 0)
row 1: (0, 0.117851)  (1, 0.707107)  (2, 0.235702)
row 2: (0, 0)  (1, 0.235702)  (2, 0.471405)
row 3: (3, 0.707107)  (4, 0.353553)  (5, 0)
row 4: (3, 0.353553)  (4, 1.64992)  (5, 0.471405)
row 5: (3, 0)  (4, 0.471405)  (5, 0.942809)

With the new FFC, the result is instead

row 0: (0, 0.235702)  (1, 0.117851)  (2, 0)
row 1: (0, 0.117851)  (1, 0.471405)  (2, 0.117851)
row 2: (0, 0)  (1, 0.117851)  (2, 0.235702)
row 3: (3, 0.235702)  (4, 0.117851)  (5, 0)
row 4: (3, 0.117851)  (4, 0.471405)  (5, 0.117851)
row 5: (3, 0)  (4, 0.117851)  (5, 0.235702)

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 29 Jan 2010 22:25:52 +0100 (CET)

I don't know if it is a bug in the new quadrature code, it depends on how you look at it. :)
The generated code in tabulate_tensor is exactly identical (apart from formatting of float values and some comments).
However, this means that I no longer reset the values of A before going crazy with the += operator. This might explain the different results. I guess I should add back the reset A code since we decided to move away form allowing different representations on the same subdomain.

Kristian

On 29 January 2010 21:42, Anders Logg <[email-masked]> wrote:
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAktjSCkACgkQTuwUCDsYZdFBmwCdFX69dDe+T61eY1XePeKmGVaz
> wMoAn0nyAR4BSNj6MyplYCOsNq88ZTV/
> =xJMD
> -----END PGP SIGNATURE-----
>
> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp
>
>

---

From: Anders Logg <[email-masked]>
Date: Sat, 30 Jan 2010 00:24:49 +0100

> However, this means that I no longer reset the values of A before going crazy with the += operator. This might explain the different results. I guess I should add back the reset A code since we dec…

I think that might be the problem. You need to set all entries to zero
in case you happen not to assign to all of them in the code. DOLFIN
does not reset the element tensor when it assembles so you need to do
that as part of the generated code.

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Sat, 30 Jan 2010 20:32:50 +0100 (CET)

> that as part of the generated code.

Did you check this again? I have pushed a changeset that should fix this issue.

Kristian

> iEYEARECAAYFAktjbkEACgkQTuwUCDsYZdFKGQCfdU1VUCaPDmiYJEFD0C43L/lh
> 62oAnAnwoz7mz3Q8gGii8ctdjfqyZd7b
> =e4a3

---

From: Anders Logg <[email-masked]>
Date: Sat, 30 Jan 2010 20:38:39 +0100

> Did you check this again? I have pushed a changeset that should fix this issue.

Yes, it works perfect. The DOLFIN unit tests all pass again.

--
Anders