On 15 Aug 2014, at 08:50, Jan Blechta <[email protected]> wrote:

> On Thu, 14 Aug 2014 09:38:48 +0100
> "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.
>> 
>>> 
>>> 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.
>> 
>> 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.
>> 
>> 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.
>> 
>>> 
>>> * ...
>>> 
>>> 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.
> 
> Note that a situation in DirichletBC is different. Expression valued
> DirichetBCs are restricted to bc->_function_space->element()
> rather than to bc->_g->element() in order to compute Dirichlet values at
> respective DOFs of bc->_function_space.
> 

The issue I describe is for Neumann conditions.

Garth

> So degree/element specified within Expression plays absolutely no role
> when evaluating it within DirichletBC.
> 
> Jan
> 
>> 
>> Garth
>> 
>>> 
>>> Cheers,
>>> Nico
>>> 
>>> 
>>> [1] 
>>> https://bitbucket.org/fenics-project/dolfin/issue/355/expression-s-interpolate-to-linear
>>> _______________________________________________
>>> fenics mailing list
>>> [email protected]
>>> http://fenicsproject.org/mailman/listinfo/fenics
>> 
>> _______________________________________________
>> fenics mailing list
>> [email protected]
>> http://fenicsproject.org/mailman/listinfo/fenics
> 

_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics

Reply via email to