dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Nov 13, 2009 at 06:46:49PM +0000, Garth N. Wells wrote:
>
>
> 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.

Yes. Let's try hard to get it right this time. The new Function class
design (as of 0.9.4) is a step in the right direction. We just need to
figure out the right way to handle expressions/user-defined
coefficients.

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

In Python, it would be fairly easy (for Johan Hake) to have just one
class named Expression that handled constant values like

  k = Expression(0.1)

or

  dt = 0.1
  k = Expression(dt)

We could let Expression handle both strings and floats, but I believe
it makes for a cleaner design to have two classes. Mixing classes that
do different things have gotten us into trouble before (like the old
Function class).

The important difference between

  k = Constant(0.1)

and

  k = Expression(0.1)

is that the first allows the JIT compiler to reuse old code when for
example changing the size of a time step. The second does not since
Expression(0.1) and Expression(0.05) will generate different code.

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

Yes, we could even let it default to None or "auto". Nothing in the
abstract handling of an UFL expression should require knowledge of the
number of quadrature points.

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

The error message needs to say that f needs to be interpolated, either
*explicitly* by

  f_approx = interpolate(f, V)

or *implicitly* by attaching a finite element to the expression. In
the latter case, the assembler will interpolate f locally on each cell
of the mesh which can be more efficient (in terms of memory usage)
since we never need to form the global vector of dofs. And this is
already part of the machinery (in UFL/FFC/UFC/DOLFIN).

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

I think this is a good option, if we also do the following:

1. Allow setting the element explicitly at the time of construction
(could be a specific quadrature rule or a regular finite element).

2. An informative error message explaining how to interpolate an
expression, either explicitly by interpolate() or implicitly by
setting the element (1).

--
Anders


> > 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).
> >>
> >>
> >>
> >>> 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
> >
> >
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev

Attachment: signature.asc
Description: Digital signature