On 15 August 2014 09:45, Jan Blechta <[email protected]> wrote:
>
> On Fri, 15 Aug 2014 09:07:30 +0200
> Kristian Ølgaard <[email protected]> wrote:
>
> > On 14 August 2014 11:09, Martin Sandve Alnæs <[email protected]>
> > wrote:
> >
> > > On 14 August 2014 10:38, Garth N. Wells <[email protected]> wrote:
> > >
> > >> On Tue, 12 Aug, 2014 at 7:51 PM, Nico Schlömer
> > >> <[email protected]> wrote:
> > >>
> > >>> Hi all,
> > >>>
> > >>> recently, a discussion [1] came up about a long-standing weirdness
> > >>> with Dolfin `Expression`s:
> > >>> ```
> > >>> from dolfin import *
> > >>> mesh = UnitIntervalMesh(1)
> > >>> e1 = Expression("x[0] * x[0]")
> > >>> print assemble(e1 * dx(mesh))
> > >>> ```
> > >>> This looks like it's integrating the function x^2 over the
> > >>> interval [0, 1], but the result – surprisingly to some –, is 0.5
> > >>> instead of 0.33333. The reason for this behavior is that `e1`
> > >>> gets interpreted as a linear function by the quadrature. (The
> > >>> optional argument `degree` of `Expression` defaults to 1.)
> > >>>
> > >>
> > >>
> > >> This is a good example of where the current approach breaks down
> > >> and needs to be fixed.
> > >>
> > >> At present, FFC uses some heuristics to interpolate an Expression
> > >> as to not mess up the order of the method. It does this by looking
> > >> at the polynomial order of the test and trial functions. In the
> > >> above case there are no test or trial functions, so it should
> > >> throw an error rather than do something that is completely
> > >> arbitrary.
> > >
> > >
> > >
> > > Just to inform the discussion:
> > > For technical implementation reasons these heuristics are now
> > > located in UFL preprocessing, more specifically in
> > > ufl/algorithms/compute_form_data.py:_auto_select_degree, which
> > > looks like this:
> > >
> > > # Use max degree of all elements, at least 1 (to work with
> > > Lagrange elements)
> > > return max({ e.degree() for e in elements } - { None } | { 1 })
> > >
> > > here elements include all elements (including subelements of mixed
> > > elements) of both test/trial and coefficient functions.
> > >
> > >
> > >
> > > While most would agree that this is surprising and that probably
> > >>> numerous FEniCS codes suffer from bugs related to this, it is
> > >>> unclear what we should best do about this.
> > >>>
> > >>> The options include:
> > >>>
> > >>> * improve the documentation
> > >>> * add a log message/warning when `Expression` interpolation
> > >>> happens
> > >>> * add an extra mandatory argument to `Expression`, specifying
> > >>> - the element or
> > >>> - the function space;
> > >>>
> > >>
> > >> This used to be the case. A concern fed back to us was that this
> > >> complicates code. My opinion at this moment is that the argument
> > >> should be mandatory. I like explicit, and it would a allow cleaner
> > >> design of some components, e.g. UFL.
> > >>
> > >>
> > >> * introduce functions that evaluate pointwise ("Expressions") to
> > >>> UFL/UFC/FCC and let dolfin Expression inherit from those;
> > >>>
> > >>
> > >> Not sure what you mean above.
> > >>
> > >
> > > That was my wording from the issue discussion: basically
> > > introducing to UFL the concept of functions without a finite
> > > element, that can be evaluated at arbitrary points. The "quadrature
> > > element" approach is flawed and something else is needed (see
> > > below).
> > >
> > >
> > >
> > >> In the absence of an argument indicating the finite element type,
> > >> I think most users would expect an Expression to be evaluated at
> > >> quadrature points. I'd have to look back at logs (which I don't
> > >> want to do) to see how developed FFC quadrature representation was
> > >> at the time the design decisions w.r.t. Expressions was made. If
> > >> Expressions are evaluated at quadrature points by default, then
> > >> most DOLFIN demos would not be able to use the FFC tensor
> > >> contraction representation. This isn't an argument either way,
> > >> just an observation.
> > >
> > >
> > > History is past, no reason to dig that deep. I would argue it's
> > > sensible to require expressions with a degree if a user wants to
> > > make tensor representation possible. In auto-representation, we can
> > > look at all functions and skip tensor representation if a function
> > > without a finite element is found. We can document this in
> > > Expression docstring: providing the degree allows for optimizations
> > > in assembly.
> > >
> > >
> > > If Expressions are evaluated at quadrature points by default, FFC
> > >> could/should throw an error if it has no information to determine
> > >> automatically a quadrature scheme, e.g. for a functional.
> > >
> > >
> > > With this notation for integration degree (there's an issue for
> > > this in ULF):
> > >
> > > M = expression*dx(degree=3)
> > >
> > > that will still be easy to use.
> > >
> > >
> > >
> > >>> * ...
> > >>>
> > >>> Arguments in favor of adding a mandatory argument include
> > >>> * The Python Zen: "Explicit is better than implicit." -- It's
> > >>> just too easy to make a mistake here.
> > >>>
> > >>> Arguments for leaving the API as it is include:
> > >>> * The user cannot expect to have arbitrary Expressions integrate
> > >>> exactly. Why should he/she expect it for any function?
> > >>>
> > >>> What do you think about this? Does the above `assemble()` result
> > >>> surprise you at all? Do you think we should force user code to be
> > >>> explicit here or would you think that adding documentation is
> > >>> sufficient?
> > >>>
> > >>
> > >> I favour (a) explicit provision of the element/function space; or
> > >> (b) evaluation at quadrature points with errors for cases where no
> > >> data is available for deciding on a sensible quadrature scheme.
> > >> Using quadrature points would fix some other awkward issues, like
> > >> specifying boundary conditions on polygon faces which 'creep'
> > >> around corners is subdomains are not marked.
> > >>
> > >
> > >
> > > I agree with both (a) and (b).
> > >
> >
> > I vote for a), there's no easy way to determine that "x[0]*x[0]"
> > should be represented by a quadratic function.
>
> Yes, I think that only sensible way using b) is that polynomial degree
> of this 'pointwise expression' is still prescribed.
If a degree is to be provided, I'd rather set the integration degree.
Besides, I want the 'degree' of a function in a 'pointwise' space to be the
number of times it can be differentiated, such that we can have pointwise
expressions with manual derivatives.
It's common to select integration degree equal to the sum of test and trial
function degrees.
Take the form (with f = Expression("..."), u,v trial and test functions)
a = f*u*v*dx
Letting the degree of f be 0 in the estimation and thus choosing
integration degree = u.degree + v.degree would be that common behaviour.
Letting the degree of f be 1 instead would be more conservative. It would
still not be the same as todays behaviour, as today the expression is first
interpolated to an element, and then that interpolation is evaluated at
quadrature points.
Martin
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
