Subject: [FFC-dev] Quadrature order

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 06 Mar 2009 16:28:48 +0000

I just tested the new FFC generated code (with quadrature) from the UFL 
input for Poisson. From the input code

     element = FiniteElement("Lagrange", "triangle", 1)
     v = TestFunction(element)
     u = TrialFunction(element)
     f = Function(element)

     a = dot(grad(v), grad(u))*dx(0, {"quadrature_order":1})
     L = v*f*dx(0, {"quadrature_order":2})

I computed exactly the same solution as with the FFC .form code, but with

     element = FiniteElement("Lagrange", "triangle", 1)
     v = TestFunction(element)
     u = TrialFunction(element)
     f = Function(element)

     a = dot(grad(v), grad(u))*dx
     L = v*f*dx

I see differences. How is the order of integration being selected?

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 06 Mar 2009 17:37:22 +0100

> I see differences. How is the order of integration being selected?

As far as I remember the UFL algorithm currently just looks at order of the
basis functions:

a: 1*1 = 2
L: 1   = 1

So the difference is expected.

Kristian

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

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 06 Mar 2009 16:41:53 +0000

> As far as I remember the UFL algorithm currently just looks at order of the

So it only looks at the test and trial functions, and not the other 
functions? And it doesn't take into account derivatives?

What, if anything, specific will UFL provide to the compiler to help 
decide on the integration order?

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 06 Mar 2009 17:49:45 +0100

> functions? And it doesn't take into account derivatives?

Yes.

> decide on the integration order?

The algorithm to compute this is present in UFL now it *just* needs to be
finished, but I don't know the priority of this.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Fri, 6 Mar 2009 18:08:01 +0100

> finished, but I don't know the priority of this.

But it's strange that there should be a difference. The old code
presumably used quadrature of high enough order to integrate
exactly. Even if the new code uses quadrature of even higher order,
there shouldn't be any difference (to within machine precision).

By the way, is quadrature_order the polynomial degree of the
polynomials that should be integrated exactly, or is it the order of
convergence. In case it is the former, maybe we should rename it to
quadrature_degree to avoid confusion?

-- 
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 06 Mar 2009 17:18:24 +0000

> there shouldn't be any difference (to within machine precision).

The new generated code is not using enough points for L with the default 
options. Below is part of Kristian's response which has somehow been 
snipped from the dialogue,

     a: 1*1 = 2
     L: 1   = 1

Exact integration requires 2 points for L which is what I presume that 
the old code uses.

Garth

> ------------------------------------------------------------------------
> 
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 06 Mar 2009 18:18:34 +0100

> there shouldn't be any difference (to within machine precision).

old code:
order of a = 0 --> 1 point
order of L = 2 --> 2 points

new code:
order of a = 2 --> 2 points
order of L = 1 --> 1 points, which is the problem.

> quadrature_degree to avoid confusion?

As I understood it quadrature_order is the polynomial degree that the quadrature
rule must be able to integrate.

Kristian

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 06 Mar 2009 18:21:04 +0100

> quadrature_degree to avoid confusion?

I actually put an explanation on the man page
man ffc
where both -f quadrature_order and -f quadrature_points are explained.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Fri, 6 Mar 2009 18:46:15 +0100

> snipped from the dialogue,

ok, I missed the right-hand side. Yes, this is obviously the cause.

> > quadrature_degree to avoid confusion?

It seems that quadrature order is also sometimes defined as the number
of points. But let's keep our current definition since it's already in
the manual/man page.

-- 
Anders