1. FEniCS Project Avatar FEniCS Project
  2. Core
  3. DOLFIN
  4. Issues

`Expression`s interpolate to linear functions *silently*

Issue #355 resolved
Nico Schlömer created an issue

As @blechta pointed out to me several times, Expressions interpolate to linear function if the degree is not specified explicitly. 

from dolfin import *
mesh = UnitIntervalMesh(1)
e1 = Expression("x[0]*x[0]")
e2 = Expression("x[0]*x[0]", degree=2)
print assemble(e1*dx(mesh))
print assemble(e2*dx(mesh))

Yikes! I bet this is a source of error for many Dolfin codes. (At least it is for me and one of my colleagues.)

How could we improve the situation?

Suggestions:

  • Make the default degree infinity, and warn or error out every time that integrals cannot be computed exactly.
  • Make the degree argument mandatory.
  • Warn if the degree argument is not specified.

Comments (36)

  2. Jan Blechta

    I'm for mandatory specification of element - not just a degree - which defaults to Lagrange. Explicit is better than implicit. Note that this is the approach done in UFL/C++DOLFIN codes. UFL does not see nothing but element of every GenericFunction/Coefficient.

    Infinite degree approach is not feasible. DOLFIN/UFL do not have tools to analyze Expressions for appropriate degree; mostly because Expression values are computated by imperative code. If user wants to take advantage of UFL algorithms to do this he may express his expression as UFL expression whenever possible.

  3. Martin Sandve Alnæs

    Without going in detail, I see two good solutions to this old issue: mandatory function space (not just element), and/or calling expressions from generated code for each quadrature point. The former idea has been shot down repeatedly over the years. There are multiple ways to solve the latter. There is no way to 'fix' this without changing the numerical result of assembly involving Expressions.

  4. Jan Blechta

    What would be the difference between functions space and element specification? It seems to me now, that Expressions don't need any additional information contained in FunctionSpace except of geometrical dimension (which is now AFAIK provided by form using the expression).

    The latter still requires that quadrature degree will be determined somehow. The 'numerical result' up to rounding error would stay the same for polynomial Expressions (defined using lossless element).

  5. Martin Sandve Alnæs

    The mesh. Consistency with Functions. Why define both a FunctionSpace and a FiniteElement with the same properties in a dolfin program? UFL is moving toward having function spaces and more domain/mesh knowledge.

    Expression("sin(x[0])") will definitely give a different result when evaluating in quadrature points with or without prior interpolation to linear elements.

  6. Nico Schlömer reporter

    Any solution that make the element specification explicit is an improvement. If it make sense in the context of UFL to associate a FunctionSpace with every Expression, as @martinal suggests, then be it. For one, this is certainly a paradigm that makes clear what actually happens when integrating.

  7. Martin Sandve Alnæs

    "If it make sense in the context of UFL to treat Expressions like Functions, as Martin Alnæs suggests,"

    I have no idea what you mean by that. I suspect I didn't explain myself well enough, but it doesn't matter.

  8. Jan Blechta

    I believe that @nschloe meant that Expressions are represented by Coefficients (with element specified) on UFL sides which is explicit when using UFL/C++DOLFIN but hidden behind the scene when using PyDOLFIN.

  9. Jan Blechta

    @martinal BTW, wouldn't the approach of switching to Expressions evaluated directly at quadrature points without 'interpolating' to element, you suggested, be equivalant to using Expression on Quadrature element?

  11. Jan Blechta

    @martinal There's one good reason why not to require function space for Expression construction but just element (with domain specified as you wish): the construction of the DOFmap is costly and may be unnecessary in this case.

  12. Martin Sandve Alnæs

    @blechta There's also no chance at all that you'll be able to convince certain other core FEniCS developers to add required arguments to Expression.

  13. Nico Schlömer reporter

    @martinal There are strong arguments for adding a required argument to Expression. From what I understand, the three of us more or less agree upon this. We can certainly discuss this with "certain other core FEniCS developers" (?). Would you add them to this bug thread so we can discuss the options?

  14. Martin Sandve Alnæs

    Let me clarify something else: By introducing functions with pointwise evaluation to UFL/UFC/FFC, there's no need for a finite element or function space for Expression.

    However the amount of design issues that pops up around changes to Expression makes this a very little tempting issue to fix. My involvement in this discussion is not volunteering to do anything.

  15. Jan Blechta

    @martinal Do you mean @logg? @garth-wells is able to abandon wrong design choices AFAIK. Also check votes - 3 votes in 2 days...

    If you're correct last option is to improve documentation and demos. I bet that there are few demos using quadratic polynomial expression defined on Lagrange degree 1 element.

  16. Martin Sandve Alnæs

    @nschloe There's nothing new in this thread. If you want to raise it to a general discussion, try the email list.

  17. Jan Blechta

    @martinal

    introducing functions with pointwise evaluation

    Expression on Quadrature element behaves likes this. Is there there some other reason for introducing them separately?

  18. Martin Sandve Alnæs

    Instead of fighting windmills, the constructive approach would be to introduce functions with pointwise evaluation (i.e. Expressions) to UFL/UFC/FFC and let dolfin Expression inherit from those if no element is provided using the metaclass magic framework.

  19. Nico Schlömer reporter

    @martinal How would the quadrature method be determined when assemble()ing those Expressions?

  20. Martin Sandve Alnæs

    @blechta I don't know. You could try making Expression use Quadrature element if none is provided in the dolfin interface. Maybe that's sufficient for this particular issue.

  21. Martin Sandve Alnæs

    @nschloe Through heuristics, e.g. let it count as degree 1 or 2 by default. The automatic quadrature degree is a somewhat flawed concept anyway, expressions as simple as 1/(1+x) or exp(x) cannot be integrated exactly with fixed quadrature. I've added an issue on making it easier to specify manually like f*dx(degree=3).

  22. Nico Schlömer reporter

    [mid-air collision post] What I think is important from the user perspective is explicitness. The code

    e1 = Expression("x[0]*x[0]")
print assemble(e1*dx(mesh))

    suggests to exactly compute the integral of $x^2$ over mesh, but doesn't. @martinal , would your suggestion change that? Maybe yield an error in assemble?

  23. Nico Schlömer reporter

    @martinal

    expressions as simple as 1/(1+x) or exp(x) cannot be integrated exactly with fixed quadrature.

    Indeed. Those things must be made obvious to the user by, e.g., throwing an according error message ("Cannot integrate! Please specify a degree. Note that the integration may not be exact then."). Instead of

    e1 = Expression("1.0 / (1+x[0]*x[0])")
print assemble(e1*dx(mesh))

    we should force the user to write

    e1 = Expression("1.0 / (1+x[0]*x[0])", degree=3)
print assemble(e1*dx(mesh))

    or

    e1 = Expression("1.0 / (1+x[0]*x[0])")
print assemble(e1*dx(mesh, degree=3))

    This makes clear what happens.

  24. Martin Sandve Alnæs

    We'll have to assume users are aware that numerical integration is not exact.

    Thus your example is mainly an issue of having good enough documentation.

    We can't throw errors for everything that needs explaining.

    To do what you wanted to do there, you should write

    x = SpatialCoordinate(mesh)
e1 = x[0]*x[0]
print assemble(e1*dx(mesh))

    which allows degree estimation to work.

  25. Nico Schlömer reporter

    I can only speak for myself here of course, but I've always assumed that

    e1 = Expression("x[0]*x[0]")
print assemble(e1*dx(mesh))

    does the correct thing. My notion is that many make the same mistake. Would you agree?

    If one is not careful, or if the definition of e1 happens in a different part of a large code, it's easy to miss that 

    e1 = Expression("1.0/(1+x[0]*x[0])")
[...]
print assemble(e1*dx(mesh))

    doesn't to the right thing either.

    The question is how to enforce this, and there is indeed some room for discussion. Your suggestion of documenting this more clearly is certainly valid. I would argue that this mistake is so easy to make that it needs stronger enforcement though. To get a more objective view, we could possibly start a poll on the mailing list asking people what they would do to integrate $x^2$ on $[0,1]$ and see what comes back to us.

  27. Martin Sandve Alnæs

    I disagree wholeheartedly with the silent interpolation that Expression does.

    I do not agree that the expectation of exact integration is reasonable.

    I agree that good enough documentation is reasonable to expect.

    I do not agree that the interface should or can protect users from making mistakes if they do not read the documentation and test the behaviour of the library themselves.

  28. Christian Clason

    Sorry for butting in, but having cast one of the three votes, I thought I'd add my thoughts from the perspective of a user (purely of the Python frontend). Like @nschloe, I was a bit surprised by what Expression was doing under the hood -- in fact, I had assumed that evaluation at quadrature points was exactly what was happening (based on the phrase "Expressions can be used as coefficients in variational forms" in the first line of the full reference, because this is how variable coefficients are usually treated in assembly). Even checking the documentation (Tutorial, brief and full reference) now, this behavior -- and default -- is apparently nowhere mentioned. Since Expressions are useful in a variety of contexts (right-hand sides, coefficients, boundary conditions, exact solutions), each with a different sensible default, I think the current situation is fine as long as

    a) it is possible to control what happens (such as decide whether interpolation or quadrature is used) and

    b) it (i.e., its possible uses and default) is clearly documented (preferably already in the tutorial).

    Regarding the silent interpolation, maybe a log/console message (not a warning, much less an error) such as "Interpolating Expression" (possibly "with degree/element") would be sufficient to alert the user that something is going on they might wish to investigate?

  29. Prof Garth Wells

    I suggest raising this issue on the mailing list for discussion. Maybe @nschloe can kick it off with a summary message. Some things have changed in FEniCS over the years, so it doesn't hurt to revisit this issue.

  34. Martin Sandve Alnæs

    Working with the code surrounding Expression is a painful experience, I would suggest adding a new experimental type on the side with a different name in ufl/ffc/ufc/dolfin, allow some experimentation period to make it work the way we want, then deprecate Expression and remove it or replace it with the new implementation later.

  35. Prof Garth Wells

    Fixed in 58f6362 by requiring at least the polynomial degree for Expressions.

    Longer term fix is to separate out FE functions and point-wise expressions.

  37. Log in to comment
Assignee
–
Type
bug
Priority
major
Status
resolved
Component
–
Milestone
–
Version
dev
Votes
5
Watchers
8
Jira: the preferred issue tracker for Bitbucket. Join the team!