ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Mar 06, 2009 at 05:18:24PM +0000, Garth N. Wells wrote:
> 
> 
> Anders Logg wrote:
> > On Fri, Mar 06, 2009 at 05:49:45PM +0100, Kristian Oelgaard wrote:
> >> Quoting "Garth N. Wells" <gnw20@xxxxxxxxx>:
> >>
> >>>
> >>> Kristian Oelgaard wrote:
> >>>> Quoting "Garth N. Wells" <gnw20@xxxxxxxxx>:
> >>>>
> >>>>> 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?
> >>>> As far as I remember the UFL algorithm currently just looks at order of the
> >>>> basis functions:
> >>>>
> >>> So it only looks at the test and trial functions, and not the other
> >>> functions? And it doesn't take into account derivatives?
> >> Yes.
> >>
> >>> What, if anything, specific will UFL provide to the compiler to help
> >>> 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
> > 
> > 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).
> > 
> 
> 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

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

> Exact integration requires 2 points for L which is what I presume that 
> the old code uses.
> 
> Garth
> 
> > 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?

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

Attachment: signature.asc
Description: Digital signature