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).
The problem with the quadrature element approach is that a quadrature
element has a fixed set of quadrature points, and cannot be used both in a
cell and on a facet, where a "boundary quadrature element" is needed. It's
possible that previous quadrature point evaluation efforts have stranded on
these complications. What's needed for evaluation of expressions at
quadrature points in generated code is one of two approaches:
1) Let the integral object return quadrature points at run time such that
the expression can be evaluated in those and the values passed into
tabulate_tensor as dofs of a "quadrature element".
or
2) Let tabulate_tensor take a vector of ufc::function pointers (or a new
ufc::expression interface) and make calls through those. For efficiency,
ufc::expression can have a vectorized eval version that gets all quadrature
points.
The latter will pave the road for an apply operator in UFL, e.g.
f = Expression(...)
M = apply(f, u**2)*dx
meaning
\int_Omega f(u^2) dx
Martin
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
