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