Buggy degree estimation of non-integer powers

[Jan Blechta]

For

e = FiniteElement("P", triangle, 6)
u = Coefficient(e)
M = u**0.5*dx

effective quadrature degree is 0/1 (1 quadrature point). Test by running FFC using the UFL above and check the comments in the generated tabulate_tensor.

The behaviour is not intended and does not match the docstring.

Thanks @amcrae for reporting.

[Andrew McRae]

This is because the code is [stripping out TPE stuff to avoid confusion]

    def power(self, v, a, b):
        """If b is an integer:
        degree(a**b) == degree(a)*b
        otherwise use the heuristic
        degree(a**b) == degree(a)*2"""
        f, g = v.ufl_operands
        try:
            gi = abs(int(g))
            return a*gi
        except:
            pass
        # Something to a non-integer power, this is just a heuristic
        # with no background
        return a*2

I.e., degree(a) is multipled by abs(int(g)) instead of 2.

Instead of try/except, I guess there should be a proper check on whether g is an int. (Remember to check cases such as 0.5**f and f**f, not just f**0.5).

[Miklós Homolya]

Perhaps we should attempt to cast it to a float instead.

[Lawrence Mitchell]

I think you want isinstance(g, numbers.Integral) no? Casting to a float doesn't dtrt if you want to implement the documented behaviour as is.

[Miklós Homolya]

I think what happens is that the exponent is a UFL expression, so it will never be an instance of numbers.Integral. But it may be castable to int or float.

[Martin Sandve Alnæs]

a, b are the results of degree estimation of f, g.

In my original implementation they were always integers, but someone from the firedrake team (Andrew?) added the possibility of tuples to have separate degrees per dimension for TPE.

Instead of try: I guess something like isinstance(g, ufl.IntValue) and g >= 0 is more correct.

[Miklós Homolya]

Fixed in pull request #75.