Correct. And nobody has claimed that so I don't get your confusion.
Martin
27. aug. 2014 18:44 skrev "Jan Blechta" <[email protected]>
følgende:
> 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
