On Thu, 14 Aug 2014 09:38:48 +0100
"Garth N. Wells" <gn...@cam.ac.uk> wrote:

> On Tue, 12 Aug, 2014 at 7:51 PM, Nico Schlömer 
> <nico.schloe...@gmail.com> 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.

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
> > fenics@fenicsproject.org
> > http://fenicsproject.org/mailman/listinfo/fenics
> 
> _______________________________________________
> fenics mailing list
> fenics@fenicsproject.org
> http://fenicsproject.org/mailman/listinfo/fenics

_______________________________________________
fenics mailing list
fenics@fenicsproject.org
http://fenicsproject.org/mailman/listinfo/fenics

Reply via email to