Derivatives of quadrature elements

[Marco Morandini]

Derivatives of quadrature elements are not supported:
with FFC the error is trapped with

ffc_assert(not (derivatives and is_quad_element), \
"Derivatives of Quadrature elements are not supported: " + repr(o))

in quadraturetransformerbase.py ;

with uflacs one gets a less clear error:

ufl.log.UFLException: Invalid cell None.

I have to compute an expression defined on quadrature elements:
it can assume only integer values and I fear that the interpolation
error of an interpolated expression would screw everything.

Unfortunately, the expression enters into a term for which I need to
compute the spatial gradient. The correct solution (for my particular
case) is to have the gradient of the expression equal to zero, i.e. to
treat an expression defined on quadrature elements as a locally constant
quantity.

Do you think this could be a sensible default, maybe with a warning
explaining what is going on to the user?

Thanks,

Marco

[Jan Blechta]

On Mon, 1 Feb 2016 16:46:40 +0100
Marco Morandini <marco.m...@polimi.it> wrote:

> Derivatives of quadrature elements are not supported:
> with FFC the error is trapped with
>
> ffc_assert(not (derivatives and is_quad_element), \
> "Derivatives of Quadrature elements are not supported: " +
> repr(o))
>
> in quadraturetransformerbase.py ;
>
> with uflacs one gets a less clear error:
>
> ufl.log.UFLException: Invalid cell None.
>
> I have to compute an expression defined on quadrature elements:
> it can assume only integer values and I fear that the interpolation
> error of an interpolated expression would screw everything.
>
> Unfortunately, the expression enters into a term for which I need to
> compute the spatial gradient. The correct solution (for my particular
> case) is to have the gradient of the expression equal to zero, i.e.
> to treat an expression defined on quadrature elements as a locally
> constant quantity.

You can overload Coefficient.is_cellwise_constant().

===========================================================
from dolfin import *

class E(Expression):
def eval(self, val, x):
val[:] = x[0]
def is_cellwise_constant(self):
return True

mesh = UnitSquareMesh(3, 3)
e = E(domain=mesh, element=FiniteElement('Quadrature', mesh.ufl_cell(), 1))
ex = e.dx(0) + Constant(0.0)

plot(e, mesh=mesh)
plot(ex, mesh=mesh)
===========================================================

>
> Do you think this could be a sensible default, maybe with a warning
> explaining what is going on to the user?

I don't think so.

Jan

>
> Thanks,
>
> Marco
>

[Jan Blechta]

On Wed, 3 Feb 2016 16:56:13 +0100
Jan Blechta <ble...@karlin.mff.cuni.cz> wrote:

> On Mon, 1 Feb 2016 16:46:40 +0100
> Marco Morandini <marco.m...@polimi.it> wrote:
>
> > Derivatives of quadrature elements are not supported:
> > with FFC the error is trapped with
> >
> > ffc_assert(not (derivatives and is_quad_element), \
> > "Derivatives of Quadrature elements are not supported: " +
> > repr(o))
> >
> > in quadraturetransformerbase.py ;
> >
> > with uflacs one gets a less clear error:
> >
> > ufl.log.UFLException: Invalid cell None.
> >
> > I have to compute an expression defined on quadrature elements:
> > it can assume only integer values and I fear that the interpolation
> > error of an interpolated expression would screw everything.
> >
> > Unfortunately, the expression enters into a term for which I need
> > to compute the spatial gradient. The correct solution (for my
> > particular case) is to have the gradient of the expression equal to
> > zero, i.e. to treat an expression defined on quadrature elements as
> > a locally constant quantity.
>
> You can overload Coefficient.is_cellwise_constant().

But note that I'm not sure that doing this does not break something
else. You need to test it carefully. Maybe Martin would know if that's
correct solution. There is probably cleaner one employing form
manipulations infrastructure of UFL.

Jan

[Marco Morandini]

>> You can overload Coefficient.is_cellwise_constant().
>
> But note that I'm not sure that doing this does not break something
> else. You need to test it carefully. Maybe Martin would know if that's
> correct solution. There is probably cleaner one employing form
> manipulations infrastructure of UFL.

Seems to work, at least for the tests I've been able to perform up to now.

Thank you very much!

>>> Do you think this could be a sensible default, maybe with a warning
>>> explaining what is going on to the user?
>>
>> I don't think so.

Now that I know how to deal with this I agree with you.

Thank you again,

Marco

[Martin Sandve Alnæs]

That's not really intended for overloading, and it might not work in all scenarios or after future changes. Your overload could be "lost" during symbolic transformations within the form compiler.

Also, are you sure you want to just set grad(f) = 0?

Martin

[Martin Sandve Alnæs]

In particular, you're not justifying the use of quadrature elements instead of DG0 in this email.

Martin

[Marco Morandini]

On 02/04/2016 02:29 PM, Martin Sandve Alnæs wrote:
> That's not really intended for overloading, and it might not work in all
> scenarios or after future changes. Your overload could be "lost" during
> symbolic transformations within the form compiler.

Ok. I'll use DG0.
See the verbose explanation below for details.

> Also, are you sure you want to just set grad(f) = 0?

Yes, definitely.

To explain why I need to give you some (tedious) background: I'm trying
to write a structural finite element beam with finite rotations is
dolfin, and I'm using the rotation vector for the parametrization.

One of the tensor that must be computed is equal to the identity for a
null rotation, but is singular for a rotation of 2\pi.
This should be not a problem provided the rotation vector is
re-normalized into a rotation with $-\pi \leq
magnitude_of_rotation_vector \leq \pi$ _after_ the interpolation.

So (some code):
-----------------------------------------
parameters["form_compiler"]["quadrature_degree"] = 3

QUADRATURE_NPHI_UNWRAP_ELEMENT_VOLUME = FiniteElement("Quadrature",
mesh.ufl_cell(), 3)

cppcode_nphi_unwrap = """
class nphi_unwrap : public Expression
{
public:

nphi_unwrap() : Expression() {}

std::shared_ptr<Function> phi_function;

void eval(Array<double>& values, const Array<double>& x, const
ufc::cell& ufc_cell) const {
Array<double> phi(3);
Cell dolfin_cell(*phi_function->function_space()->mesh(),
ufc_cell.index);
phi_function->eval(phi, x, dolfin_cell, ufc_cell);
double phi1 = sqrt(phi[0] * phi[0] + phi[1] * phi[1] + phi[2] *
phi[2]);
if (phi1 > pi) {
double pi2 = 2. * pi;
int nphi = trunc(phi1 / pi2);
if ((phi1 - nphi * pi2) > pi ) nphi += 1;
values[0] = nphi * pi2;
} else {
values[0] = 0.;
}
}
};"""


NphiUnrwapVolume = Expression(cppcode_nphi_unwrap,
element=QUADRATURE_NPHI_UNWRAP_ELEMENT_VOLUME, domain=mesh)

# phi is the rotation vector
NphiUnrwapVolume.phi_function = phi

def phi_unwrapped(phi, nphi):
phi2 = dot(phi, phi)
phi1 = sqrt(phi2)
dphi = conditional(nphi > .5, nphi / phi1, 0.)
return phi * (1. - dphi)

phi_unwrapped_volume = phi_unwrapped(phi, NphiUnrwapVolume)

Gamma = function_that_compute_the_tensor(phi_unwrapped_volume)

-----------------------------------------

The functional is a function of Gamma and of grad(phi_unwrapped_volume)
and while computing grad(phi_unwrapped_volume) ufl tries to compute
grad(nphi) (if I'm understanding correctly what is going on).

Everything should work because

1) you've fixed the derivative of conditionals in uflacs
(at least for scalars)

2) when computing the derivative nphi is just a number (switch or not
switch), and as a consequence the derivative of phi_unwrapped_volume
with respect to the set on unknowns should be correct at each
integration point.

> In particular, you're not justifying the use of quadrature elements
> instead of DG0 in this email.

You are right. It occurred to me _after_ sending the email
(after and not before, as always, silly me)
that I can use DG0: what should be important is to stay away from 2\pi
at each integration point. However,
if the elements are small enough this should be guaranteed with DG0
because the switch occurs if the rotation magnitude is $\ge \pi$, well
away from the critical value of 2\pi.

Thanks again,

Marco

[Martin Sandve Alnæs]

Thanks for explaining. I think your original thought of defining

grad(quadrature_function) to be 0 by default with a warning

is _better_ than today's error, so I'm fine with adding that solution to ufl.


A further improvement over that would be user-specified

quadrature functions for the gradient of quadrature function,

but that's more work.


Martin

[Jan Blechta]

On Thu, 4 Feb 2016 16:45:27 +0100
Marco Morandini <marco.m...@polimi.it> wrote:

> dphi = conditional(nphi > .5, nphi / phi1, 0.)

^^^^^^^^^
This will not work as expected.

Jan

[Jan Blechta]

On Fri, 5 Feb 2016 10:13:10 +0100
Martin Sandve Alnæs <mart...@simula.no> wrote:

> Thanks for explaining. I think your original thought of defining
>
> grad(quadrature_function) to be 0 by default with a warning
>
> is _better_ than today's error, so I'm fine with adding that solution
> to ufl.

Martin, could you give some reasoning why grad(quadrature_f) == 0 is
sensible choice.

I just don't see from Marco's complicated example. (There was even no
reasoning why it should be quadrature element. Maybe that it is
volume, so that NphiUnrwapVolume*dx naturally describes a physical
volume measure...? But why should derivatives of it appear anywhere?!)

Anyway I'd expected somewhat more general reasoning for this decision.

Jan

[Martin Sandve Alnæs]

Why not?

Martin

[Jan Blechta]

On Fri, 5 Feb 2016 11:16:35 +0100

Martin Sandve Alnæs <mart...@simula.no> wrote:

Citing from relevant section of UFL manual:

Because of details in the way Python behaves, we cannot over-
load the builtin comparison operators for this purpose, hence these
named operators [eq, ne, le, ge, lt, gt].

I haven't noticed that this language-design limitation had been
overcome :)

Jan

>
> Martin
>

[Martin Sandve Alnæs]

The manual needs more updating...

It has for < > <= >= but not for ==, !=.

>>> import ufl

>>> from ufl import *

>>> x = SpatialCoordinate(triangle)

>>> print x[0] < x[1]

(x[0]) < (x[1])

Martin

[Martin Sandve Alnæs]

Jan, I take it back, it's not a good idea :)

Setting grad(quadrature_f) == 0 will lead to silently

wrong computations and is therefore not a good idea.

The details of his example are not important for this decision.

We will at some point have Expression evaluated in quadrature points

instead of silently interpolated. Then this will become a more widespread

issue. Solutions involve either interpolation similar to today,

or manually specifying gradients as another expression in quadrature points.

Lets leave the details until someone gets time to do this work.

Martin

[Jan Blechta]

On Fri, 5 Feb 2016 11:31:59 +0100
Martin Sandve Alnæs <mart...@simula.no> wrote:

> The manual needs more updating...
> It has for < > <= >= but not for ==, !=.

The current one (in ufl source tree) is ok.

I was citing from the old one. There does not seem to be any up-to-date
manual available on the webpage,
http://fenicsproject.org/pub/documents/ufl/ufl-user-manual/ufl-user-manual.pdf
is terribly old. There is no reasonable link to the manual on
http://fenicsproject.org/documentation/,
only to the book excerpts
http://launchpad.net/fenics-book/trunk/final/+download/fenics-manual-2011-10-31.pdf

This will be resolved by RTD,
https://fenics.readthedocs.org/projects/ufl/en/latest/user/form_language.html#conditions
is up-to-date and correct on the conditions.

Jan