Lets keep the guesswork out of it, here's the code for automatic element
selection, currently residing in ufl/algorithms/compute_form_data.py:
def _auto_select_degree(elements):
"""
Automatically select degree for all elements of the form in cases
where this has not been specified by the user. This feature is
used by DOLFIN to allow the specification of Expressions with
undefined degrees.
"""
# Use max degree of all elements, at least 1 (to work with Lagrange
elements)
return max({ e.degree() for e in elements } - { None } | { 1 })
def _compute_element_mapping(form):
"Compute element mapping for element replacement"
# Extract all elements and include subelements of mixed elements
elements = [obj.element() for obj in chain(form.arguments(),
form.coefficients())]
elements = extract_sub_elements(elements)
# Try to find a common degree for elements
common_degree = _auto_select_degree(elements)
# Compute element map
element_mapping = {}
for element in elements:
# Flag for whether element needs to be reconstructed
reconstruct = False
# Set domain/cell
domain = element.domain()
if domain is None:
domains = form.domains()
ufl_assert(len(domains) == 1,
"Cannot replace unknown element domain without
unique common domain in form.")
domain, = domains
info("Adjusting missing element domain to %s." % (domain,))
reconstruct = True
# Set degree
degree = element.degree()
if degree is None:
info("Adjusting missing element degree to %d" %
(common_degree,))
degree = common_degree
reconstruct = True
# Reconstruct element and add to map
if reconstruct:
element_mapping[element] = element.reconstruct(domain=domain,
degree=degree)
else:
element_mapping[element] = element
return element_mapping
On 27 August 2014 17:48, Jan Blechta <[email protected]> wrote:
> On Wed, 27 Aug 2014 12:13:11 +0200
> Kristian Ølgaard <[email protected]> wrote:
>
> > On 27 August 2014 11:36, Garth N. Wells <[email protected]> wrote:
> >
> > >
> > >
> > > On Wed, 27 Aug, 2014 at 10:21 AM, Kristian Ølgaard
> > > <[email protected]> wrote:
> > >
> > >> On 27 August 2014 10:44, Garth N. Wells <[email protected]> wrote:
> > >>
> > >>>
> > >>>
> > >>> On Wed, 27 Aug, 2014 at 9:15 AM, Martin Sandve Alnæs
> > >>> <[email protected]> wrote:
> > >>>
> > >>>> On 27 August 2014 08:25, Kristian Ølgaard
> > >>>> <[email protected]> wrote:
> > >>>>
> > >>>>> On 26 August 2014 15:59, Martin Sandve Alnæs
> > >>>>> <[email protected]> wrote:
> > >>>>>
> > >>>>>> On 26 August 2014 15:21, Kristian Ølgaard
> > >>>>>> <[email protected]> wrote:
> > >>>>>>
> > >>>>>>>
> > >>>>>>>
> > >>>>>>> ---------- Forwarded message ----------
> > >>>>>>> From: Kristian Ølgaard <[email protected]>
> > >>>>>>> Date: 26 August 2014 15:20
> > >>>>>>> Subject: Re: [FEniCS] `Expression`s and their silent
> > >>>>>>> interpolation To: Jan Blechta <[email protected]>
> > >>>>>>>
> > >>>>>>>
> > >>>>>>> On 26 August 2014 14:18, Jan Blechta
> > >>>>>>> <[email protected]> wrote:
> > >>>>>>>
> > >>>>>>>> On Tue, 26 Aug 2014 09:50:23 +0100
> > >>>>>>>> "Garth N. Wells" <[email protected]> wrote:
> > >>>>>>>>
> > >>>>>>>> > To summarise this thread, it seems we need to introduce the
> > >>>>>>>> concept
> > >>>>>>>> > of an 'Expression' that can be evaluated at arbitrary
> > >>>>>>>> > points. It should not be a Quadrature{Element/Function}
> > >>>>>>>> > because the proposed object could be used in different
> > >>>>>>>> > forms with different evaluation
> > >>>>>>>>
> > >>>>>>>
> > >>>>>> Agree. Also it should not have any notion of degree. Pointwise
> > >>>>>> is pointwise.
> > >>>>>>
> > >>>>>>
> > >>>>>> > points. The follow-on on issue is then how a 'point-wise'
> > >>>>>>>> expression
> > >>>>>>>> > should be treated in forms. We could estimate the
> > >>>>>>>> > quadrature
> > >>>>>>>> scheme
> > >>>>>>>> > when test/trial functions are present, and in the case of
> > >>>>>>>> functionals
> > >>>>>>>> > throw an error if the user doesn't supply the quadrature
> > >>>>>>>> > degree.
> > >>>>>>>>
> > >>>>>>>> There's no principal difference regarding rank of the form.
> > >>>>>>>> Consider
> > >>>>>>>>
> > >>>>>>>> f = PointwiseExpression(eval_formula)
> > >>>>>>>> u, v = TrialFunction(V), TestFunction(V)
> > >>>>>>>> a = f*u*v*dx
> > >>>>>>>> L = f*v*dx
> > >>>>>>>> F = f*dx
> > >>>>>>>>
> > >>>>>>>> Still, you need to know what is the polynomial degree of f
> > >>>>>>>> to have exact quadrature of any of these forms. Ignoring
> > >>>>>>>> non-zero degree of f
> > >>>>>>>> (which seems to me you do suggest for a and L) means that
> > >>>>>>>> you're underintegrating any of those three forms. This is
> > >>>>>>>> analogical to integrating F with scheme of order zero. I
> > >>>>>>>> don't see any good reason why having distinct behaviour
> > >>>>>>>> based on rank of the respective form.
> > >>>>>>>>
> > >>>>>>>
> > >>>>>> Agree.
> > >>>>>>> For PointwiseExpression, one should define EITHER the
> > >>>>>>> polynomial degree that the user would like the use for the
> > >>>>>>> approximation (of e.g., 'sin(x)') OR the (degree of)
> > >>>>>>> quadrature rule for the measure. The latter should take
> > >>>>>>> precedence if both are defined, just as it does currently.
> > >>>>>>>
> > >>>>>>
> > >>>>>> Please, no. Isn't that basically the situation we're trying to
> > >>>>>> get away from? A pointwise expression doesn't have a degree
> > >>>>>> and it's not a good abstraction to assign one to it. The rules
> > >>>>>> become complex which makes the source code hard to follow, the
> > >>>>>> documentation poor, and confuses the users and developers
> > >>>>>> alike.
> > >>>>>>
> > >>>>>
> > >>>>> Perhaps I got a little confused. Is PointwiseExpression
> > >>>>> supposed to replace Expression or will they co-exist?
> > >>>>>
> > >>>>> If PointwiseExpression will replace Expression, my suggestion to
> > >>>>> forcing the user to supply the degree (or element) was intended
> > >>>>> to solve Nico's original problem.
> > >>>>>
> > >>>>
> > >>>> Good question. There are two sides to that:
> > >>>>
> > >>>> From the UFL/FFC/UFC/Assembler side there needs to be two
> > >>>> distinct concepts: If a function is pointwise evaluated FFC must
> > >>>> generate function calls instead of using basis tables and dofs.
> > >>>> Whether this is implemented by adding a PointwiseExpression or a
> > >>>> "Pointwise" family to UFL is an implementation detail that
> > >>>> affects other work in progress so lets not go there in this
> > >>>> thread at least. This is basically the new functionality we're
> > >>>> discussing here, and I don't think anyone disagrees with this,
> > >>>> apart from annotating the PointwiseExpression with a "virtual
> > >>>> degree" for integration degree purposes.
> > >>>>
> > >>>>
> > >> Agree, let's discuss this later.
> > >>
> > >> From the DOLFIN user interface side, there are several ways to
> > >> go, and
> > >>>> this is where most opinions in this thread differ.
> > >>>>
> > >>>> The main interface point is whether to introduce new notation
> > >>>> PointwiseExpression("x[0]") and deprecate Expression("x[0]"), or
> > >>>> to change Expression("x[0]") to mean pointwise and remove the
> > >>>> implicit interpolation.
> > >>>>
> > >>>> Deprecating Expression("x[0]") breaks all demos and lots of user
> > >>>> programs.
> > >>>>
> > >>>> Changing the behaviour of Expression("x[0]") changes the
> > >>>> numerical results of lots of programs.
> > >>>>
> > >>>
> > >>> I'd be reluctant to introduce yet another object in the form of
> > >>> PointwiseExpression. How about passing a string or element
> > >>> argument to Expression? For now, the default could be the present
> > >>> behaviour (interpolation) with a warning that the default will
> > >>> change in release 1.foo to pointwise?
> > >>>
> > >>
> > >> Do we need two arguments to do that?
> > >>
> > >> Expression('x[0]*x[0]', pointwise=False, element=None)
> > >>
> > >
> > > We need to allow passing of different types, but it could be just
> > > one argument and the Expression constructor can do whatever it
> > > needs to do depending on the type.
> >
> >
> > OK.
> >
> >
> > >
> > >
> > > This should then default to the current behaviour and interpolate
> > > the
> > >> expression using a linear element.
> > >>
> > >
> > >
> > > For now, yes. The default isn't a linear element. I think there is
> > > some magic to pick a suitable order.
> >
> >
> > Could be, I don't recall it exactly.
>
> I don't think there is any parser of C++/Python code guessing the
> degree. It is not anywhere in DOLFIN and UFL does not see
> Expression::eval code.
>
> Jan
>
> >
> >
> > >
> > >
> > > One can then supply the element to get the 'exact' integration
> > >>
> > >> Expression('x[0]*x[0]', pointwise=False, element=quadratic_element)
> > >>
> > >
> > > No need to pass pointwise. Just an element would do.
> >
> >
> > Sure, that was just to show the complete argument list, but passing
> > one argument should be enough.
> >
> >
> > >
> > >
> > > and test the implementation of the pointwise feature by
> > >>
> > >> Expression('x[0]*x[0]', pointwise=True)
> > >>
> > >> which, depending on implementation details, will have a virtual
> > >> degree equal to zero.
> > >>
> > >
> > > If one did
> > >
> > > f = Expression('x[0]*x[0]', pointwise=True)
> > > assemeble(f*dx)
> > >
> > > it should throw an error. Integrating f would require something like
> > >
> > > assemeble(f*dx(degree=2))
> >
> >
> > Yes.
> >
> > Kristian
> >
> >
> > >
> > >
> > >
> > > In the future we can either set pointwise=True by default, or
> > > remove this
> > >> argument such that failing to provide the element implies pointwise
> > >> evaluation.
> > >>
> > >
> > > Yes.
> > >
> > > Garth
> > >
> > >
> > > Kristian
> > >>
> > >>
> > >>
> > >>
> > >>
> > >>> Garth
> > >>>
> > >>>
> > >>>
> > >>> Martin
> > >>>>
> > >>>>
> > >>>> These are two distinct issues:
> > >>>>>> 1) We need a "PointwiseExpression" with no degree and no hidden
> > >>>>>> interpolation under the hood. This expression is evaluated in
> > >>>>>> quadrature points - this is a clean concept and easy to
> > >>>>>> understand.
> > >>>>>>
> > >>>>>
> > >>>>> If PointwiseExpression and Expression will co-exist, I agree to
> > >>>>> this definition. If it makes implementation cleaner for
> > >>>>> quadrature degree estimation, we can set the degree equal to
> > >>>>> zero, but hidden from users.
> > >>>>>
> > >>>>>
> > >>>>> 2) Degree estimation is not exact and some people are confused
> > >>>>> by
> > >>>>>> that. But it is not exact today, never was claimed to be, and
> > >>>>>> never will be. If that's not acceptable, we can just as well
> > >>>>>> disable it completely. Disabling it where it isn't exact will
> > >>>>>> break a _lot_ of programs. What we _can_ do without breaking
> > >>>>>> programs or making the interface more cumbersome than today,
> > >>>>>> is to make it more obvious how to control the integration
> > >>>>>> degree, and to document it better.
> > >>>>>>
> > >>>>>
> > >>>>> I don't think anyone can disagree with this.
> > >>>>>
> > >>>>> Kristian
> > >>>>>
> > >>>>>
> > >>>>>> Martin
> > >>>>>>
> > >>>>>
> > >>>>>
> > >>>>
> > >>>
> > >>
> > >
>
> _______________________________________________
> fenics mailing list
> [email protected]
> http://fenicsproject.org/mailman/listinfo/fenics
>
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
