Subject: [UFL-dev] Quadrature degree estimation

From: Martin Sandve Alnæs <[email-masked]>
Date: Fri, 17 Apr 2009 12:02:04 +0200

Take a look at
  ufl/sandbox/algorithmtests/degree_estimation.py
just added.

Is this good enough for the initial version?
It handles mixed elements by using the max
degree of subelements, and nonlinear operators
by just using the degree of the operands.
Should work fine for most forms where accurate
estimation is possible.

I'll add it into UFL soon.

Martin

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Fri, 17 Apr 2009 12:27:02 +0200

However, the rule

    def spatial_derivative(self, v, f, i):
        return max(f-1, 0)

isn't really correct, e.g. a basis function like "x*y*z"
differentiated in one direction would still be

Martin

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Fri, 17 Apr 2009 12:32:51 +0200

Oops, hit the send button. Forget this, I think it's ok.

Martin

---

From: Anders Logg <[email-masked]>
Date: Sat, 18 Apr 2009 21:01:01 +0200

> I'll add it into UFL soon.

It looks good to me, but I haven't tried it.

UFL now seems to have at least three ways to estimate the quadrature
degree:

  extract_quadrature_order
  estimate_quadrature_order
  estimate_max_quad_degree
  estimate_max_quadrature_order

Which one should be used and what is the difference between these?

-- 
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Sun, 19 Apr 2009 14:11:12 +0800

>> Should work fine for most forms where accurate

Should it be called 'estimate_max_poly_order' rather than 
'estimate_max_quadrature_order'?

> UFL now seems to have at least three ways to estimate the quadrature

This gets the order of a quadrature element.

>   estimate_quadrature_order

I don't follow exactly what this function is doing.

>   estimate_max_quad_degree

I can't find this function in the code.

>   estimate_max_quadrature_order

This computes the polynomial order of the integrand.

Garth

> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev

---

From: Anders Logg <[email-masked]>
Date: Sun, 19 Apr 2009 08:46:08 +0200

> 'estimate_max_quadrature_order'?

My suggestion would be

  estimate_max_polynomial_degree

Then it's very clear what it is.

> This gets the order of a quadrature element.

ok. I find the naming a bit confusing. How about we change the name of
this as well. Then we have two functions:

  estimate_max_polynomial_degree
  extract_max_quadrature_element_degree

Maybe the first one should also be named 'extract'. Or is it
'estimate' because it might only make a qualified guess if the
integrand is not a polynomial? So we could have

  extract_max_polynomial_degree
  extract_max_quadrature_element_degree

Then it would be very clear what these two functions do.

Martin?

> I don't follow exactly what this function is doing.

This is still in the sandbox.

Do you have a good test example that can be added to the demos so we
verify things work as expected?

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Sun, 19 Apr 2009 14:02:57 +0200

With some renaming, we now have:

extract_max_quadrature_element_degree:
Returns max degree of found quadrature elements, None with no
quadrature elements

estimate_quadrature_degree:
Uses the formula 2*v.element().degree() for linear forms,
v.element().degree()+u.element().degree() for bilinear forms,
don't know if we need this? Returns None for functionals now.
The form compilers could have an option to use this variant.

estimate_max_polynomial_degree:
Estimates the highest polynomial degree of the expression.
It is an estimate since not all operators are polynomials.

I've added a unit test, but only for some basic cases.

Martin

On Sun, Apr 19, 2009 at 8:46 AM, Anders Logg <[email-masked]> wrote:
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAknqyLAACgkQTuwUCDsYZdEPsACfbtx/920zZNHQOBUUyBIsDbTZ
> KLcAmgK6iXnkuYTCYoMN+tGbykPeTEy9
> =KkbD
> -----END PGP SIGNATURE-----
>

---

From: Anders Logg <[email-masked]>
Date: Sun, 19 Apr 2009 14:48:02 +0200

> Returns max degree of found quadrature elements, None with no

Sounds good.

> The form compilers could have an option to use this variant.

Is it to be used as a fallback if the others fail?

> It is an estimate since not all operators are polynomials.

Does it return something sensible even if the expression is not a
polynomial?

If so, then estimate_quadrature_degree may not be needed.

It would be enough for a form compiler to do this:

  quadrature_degree = extract_max_quadrature_element_degree(integral)
  if quadrature_degree is None:
      quadrature_degree = estimate_max_polynomial_degree(integral)

Perhaps this should be in UFL?

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Sun, 19 Apr 2009 15:03:44 +0200

> Is it to be used as a fallback if the others fail?

No, from what I understand this is the default in e.g. Diffpack,
and may often be sufficient? I don't know the conditions for this.
I was thinking that we can give the choice between the two
approaches to the user through options to the form compilers.

Kent, what do you say?

> Does it return something sensible even if the expression is not a

Yes, the default rule for unhandled types are to return the max
degree of the operands, so the degree of exp(f**2) will be the
same as the degree of f**2. We could define heuristic rules like
degree(exp(f)) -> 2*degree(f), if that makes any sense.

> Perhaps this should be in UFL?

As I said above, I think this should be controllable by form compiler options,
so in that case no. Also, Kristian wants splitting of expression terms
by polynomial degree in the future, which will be yet another way to
handle this.

Martin

---

From: Anders Logg <[email-masked]>
Date: Sun, 19 Apr 2009 15:21:30 +0200

> approaches to the user through options to the form compilers.

ok. That may be useful to have.

We can decide later what to name that option.

> degree(exp(f)) -> 2*degree(f), if that makes any sense.

ok.

> by polynomial degree in the future, which will be yet another way to

ok.

-- 
Anders

---

From: [email-masked]
Date: Sun, 19 Apr 2009 17:06:38 +0200 (CEST)

> approaches to the user through options to the form compilers.

Yes, this was the default in Diffpack. I think
this is common for elliptic equations even with coefficients.

As long as the default option are well documented and the order can be
changed easily, I think this is a reasonable default option. But I
don't have any strong opinion. I do however think it might be difficult to
find a general formula that is good in the sense that it will not
overestimate or underestimate in various situations.

Kent

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 20 Apr 2009 18:14:57 +0100

> this is common for elliptic equations even with coefficients.

I think that it's more common to use a scheme that integrates the 
product of the test and trial functions  exactly (coefficients are ignored).

Garth

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 20 Apr 2009 18:18:37 +0100

> The form compilers could have an option to use this variant.

I find the use of 'degree' in the above contexts confusing. Is 
'max_quadrature_element_degree' the number of points? Is it the order of 
accuracy?

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Apr 2009 19:29:57 +0200

> 'max_quadrature_element_degree' the number of points? Is it the order of  

It's the maximum degree of all quadrature elements in the form, which
is the degree argument given as input:

  element = FiniteElement("Quadrature", shape, degree)

Quadrature element shares this interface with other element families so
the name of the variable is "degree". UFL does not know anything about
what degree means for a quadrature element (nor does it know what it
means for RT elements).

The form compiler extracts this degree and in the case of FFC it sets
quadrature_order to be equal to this degree unless quadrature_order
has been specified by a user.

-- 
Anders