ffc team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

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

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?
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> UFL-dev mailing list
> UFL-dev@xxxxxxxxxx
> http://fenics.org/mailman/listinfo/ufl-dev