# `Expression`s interpolate to linear functions *silently*

Issue #355 resolved
Nico Schlömer created an issue

As @blechta pointed out to me several times, `Expression`s 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.

1. reporter
• edited description

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.

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.

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).

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.

2. 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.

"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.

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

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

@blechta at least mathematically equivalent, if not technically.

@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.

@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.

3. 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?

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.

@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.

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

@martinal

introducing functions with pointwise evaluation

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

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.

4. reporter

@martinal How would the quadrature method be determined when `assemble()`ing those `Expression`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.

@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).

5. 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`?

6. 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.

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.

7. 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.

8. reporter

@martinal BTW, great tip with the `SpatialCoordinate`s. Didn't know about this!

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.

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 `Expression`s 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?

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.

• changed milestone to 1.7

@blechta Nice digging. Looks like we were doing it right in 2009, and then changed.

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.

9. 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.

• removed milestone

Removing milestone: 1.7 (automated comment)