Sorry, for my ignorance but it seems to me that the degree can't be
guessed from Expression::eval code itself when not specified:
e = Expression("x[0]*x[0]")
print ufl.algorithms.estimate_total_polynomial_degree(e) # 1
e = Expression("x[0]*x[0]", degree=2)
print ufl.algorithms.estimate_total_polynomial_degree(e) # 2
Jan
On Wed, 27 Aug 2014 18:01:33 +0200
Martin Sandve Alnæs <[email protected]> wrote:
> 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
