Sorry, I've misunderstood Garth and Kristian. Jan
On Wed, 27 Aug 2014 18:49:52 +0200 Martin Sandve Alnæs <[email protected]> wrote: > 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
