dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #16538
Re: Expression class
Hans Petter Langtangen wrote:
I think users will find it strange that a Constant and Expression must
be associated with a mesh. Why not allow
g = Constant(6.0)
f = Expression('sin(x[0])')
For quadrature one only needs to evaluate a Constant or Expression so
I don't see any need for a mesh or space if there are other components
in the form that have a space or mesh to be used for integration.
We need to re-assess the design around Constant/Expression/Function.
We've lost our way a little with the re-designs around the function
classes. A first step is to clean up some issues on the C++ side (like
why we need the geometric dimension) which will trickle down to the
Python side. The following step would be to address some of the
Python-specific issues.
Now that we have Expression, perhaps we should get rid of Constant since
g = Constant(6.0)
h = Expression('6.0')
serve essentially the same purpose (a small difference, at least on the
C++ side, is that a Constant can be treated as a double, e.g. it can be
incremented).
For an Expression, we may be able to default to a QuadratureFunction to
avoid providing the mesh/function space at initialisation. The possible
issue at present is that UFL may insist on the order of the
QuadratureFunction being known, but this could be changed in UFL.
(The mesh determines the no of space dimensions and the length of the
x argument when an Expression is evaluated.)
For non-quadrature approaches, f must be interpolated onto a space
and an error message could tell the user to explicitly do so.
In case one does
M = f*dx # or g*dx
assemble(M)
the form or assemble lacks a mesh to do the integration and then an
error message can be issued, telling the user to insert a mesh or
space argument in assemble (I think this is more logical than
attaching a mesh or space to f). I would guess that one needs a space
and not only a mesh since it's otherwise hard to figure out the
appropriate order of the quadrature rule?
In situations where grad(f) is needed, one must follow Doug's suggestion
and explicitly interpolate:
f_approx = interpolate(f, V)
dot(grad(f_approx), grad(u))*dx
grad(Expression) should just give an error message saying that the
operation is impossible, which is intuitive: it's a formula and there is no
information about the derivative. (The optimal solution would be to
accept an eval_derivative function or an optional string for the
derivative: Expression('sin(x[0])', derivative=['cos(x[0])', 0, 0])).
This could be added easily.
In most cases, an Expression enters a form together with other quantities
that have associated function spaces. Isn't that sufficient?
If one is happy to evaluate the Expression at quadrature points but let
the form compiler determine the number of points, this should be OK.
As a tutorial writer, it is quite a problem to explain why the very
simplest examples have non-intuitive constructions when everything else
in FEniCS is so intuitive.
It's very useful to have the perspective of the tutorial writer!
Garth
These suggestions come from a naive user's point of view - I realize
that some of the ideas may be impossible because of inner workings in
UFL/DOLFIN/PyDOLFIN. My main point is to start with what's intuitive
and then issue error messages when a mesh or space needs to be
associated with a Constant or Expression.
Hans Petter
Fri, 13 Nov Anders Logg wrote:
On Fri, Nov 13, 2009 at 09:46:21AM +0000, Garth N. Wells wrote:
Anders Logg wrote:
On Fri, Nov 13, 2009 at 08:18:38AM +0000, Garth N. Wells wrote:
Anders Logg wrote:
On Fri, Nov 13, 2009 at 07:28:07AM +0100, Johan Hake wrote:
On Thursday 12 November 2009 21:15:51 Anders Logg wrote:
I have received some complaints on the new Expression class. It works
reasonably well from C++ but is still confusing from Python:
1. Why does a function space need to be specified in the constructor?
f = Expression("sin(x[0])", V=V)
Does this mean f is a function on V? (No, it doesn't.)
2. The keyword argument V=foo is confusing:
f = Expression(("sin(x[0])", "cos(x[1])"), V=V)
g = Expression("1 - x[0]", V=Q)
The reason for the function space argument V is that we need to know
how to approximate an expression when it appears in a form. For
example, when we do
L = dot(v, grad(f))*dx
we need to interpolate f (locally) into a finite element space so we
can compute the gradient.
Sometimes we also need to know the mesh on which the expression is defined:
M = f*dx
This integral can't be computed unless we know the mesh which is why
one needs to do
assemble(M, mesh=mesh)
My suggestion for fixing these issues is the following:
1. We require a mesh argument when constructing an expression.
f = Expression("sin(x[0])", mesh)
2. We allow an optional argument for specifying the finite element.
f = Expression("sin(x[0])", mesh, element)
We could also have
f = Expression("sin(x[0])", mesh, k)
where k is the order of the continuous Lagrange basis since that's the
most commonly used.
3. When the element is not specified, we choose a simple default
element which is piecewise linear approximation. We can derive the
geometric dimension from the mesh and we can derive the value shape
from the expression (scalar, vector or tensor).
This is bad. If a user increases the polynomial order of the test/trial
functions and f remains P1, the convergence rate will not be optimal.
A better solution would be to define it on a QuadratureElement by
default. This, I think, is the behaviour that most people would expect.
This would take care of higher-order cases.
Yes, I've thought about this. That would perhaps be the best solution.
Does a quadrature element have a fixed degree or can the form compiler
adjust it later to match the degree of for example the basis functions
in the form?
Kristian?
If one should happen to differentiate a coefficient in a form, we just
need to give an informative error message and advise that one needs to
specify a finite element for the approximation of the coefficient.
This will remove the need for the V argument and the confusion about
whether an expression is defined on some function space (which it is
not).
But it is when it's used in a form since it's interpolated in the given
space.
The important point is that we only need to be able to interpolate it
locally to any given cell in the mesh, so we just need a mesh and a
cell. The dofmap is never used so it's different from a full function
space.
I think you mean element rather than cell?
I actually meant cell. But the restriction operator (extracting the
information we need to integrate over the cell) requires an element
(which implements evaluate_dof).
--
Anders
Garth
It also removes the need for an additional mesh=mesh argument
when assembling functionals and computing norms.
This change will make the Constant and Expression constructors very
similar in that they require a value and a mesh (and some optional
stuff). Therefore it would be good to change the order of the
arguments in Constant so they are the same as in Expression:
f = Constant(0, mesh)
f = Expression("0", mesh)
Yes.
Good.
And we should change Constant rather than Expression since Expression
might have an optional element argument:
f = Expression("0", mesh, element)
Does this sound good?
Yes, I think so. I suppose you mean
f = CompiledExpression("0", mesh)
Just referring to the Blueprint about the simplification of the Expression
class in PyDOLFIN.
I'm not so sure anymore. Calling it Expression looks simpler.
Agree.
Good.
Garth
What
were the reasons for splitting it into Expression and
CompiledExpression? Is it the problem with the non-standard
constructor when implementing subclasses?
It's just that the most common use of expressions in simple demos will
be stuff like
f = Expression("sin(x([0])", mesh)
so one could argue that this should be as simple as possible (and just
be named Expression).
Should an Expression in PyDOLFIN then always have a mesh? This will make an
Expression in PyDOLFIN and DOLFIN different, which is fine with me.
Yes, to avoid needing to pass the mesh to assemble() and norm() in
some cases and to automatically get the geometric dimension.
If others agree, can you add it to the Blueprint, mentioned above, and I can
do the change some time next week, or after a release (?).
Let's hear some more comments first.
_______________________________________________
DOLFIN-dev mailing list
DOLFIN-dev@xxxxxxxxxx
http://www.fenics.org/mailman/listinfo/dolfin-dev
_______________________________________________
DOLFIN-dev mailing list
DOLFIN-dev@xxxxxxxxxx
http://www.fenics.org/mailman/listinfo/dolfin-dev
Follow ups
References
