confusion with grad(f), f.dx(0), f.dx(1), f.dx(0,1)

Asked by Chaffra Affouda

In ufl, what is the difference between grad(f) and f.dx(0)? If f is a function what does it mean to write f.dx(1) instead of f.dx(0)?

You can't extract vector component from f.dx(0) for example but I thought it was the same as grad(f): f.dx(0)[0] breaks.

Also how do you extract the vector component from graf(f) ayt a given position

mesh = UnitSquare(10,10)
Q = FunctionSpace(mesh,'CG',1)

expr = Expression('exp(x[0])')
f = project(expr,Q)

df = grad(f)
df[0]([0.5,0.5]) = 1.5689657305885802
df[1]([0.5,0.5]) = 1.5689657305885802 #SHOULD BE ZERO!

df([0.5,0.5]) #breaks with Keyerror

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Martin Sandve Alnæs
Solved:
Last query:
Last reply:
Revision history for this message
Chaffra Affouda (chaffra) said :
#1

dfx = project(df[0])
dfy = project(df[1]) seem to give the right answers.

Revision history for this message
Anders Logg (logg) said :
#2

On Tue, Jun 26, 2012 at 04:01:10PM -0000, Chaffra Affouda wrote:
> New question #201507 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/201507
>
> In ufl, what is the difference between grad(f) and f.dx(0)? If f is a function what does it mean to write f.dx(1) instead of f.dx(0)?

grad(f) = (df/dx, df/dy, df/dz)
f.dx(0) = df/dx
f.dx(1) = df/dy
f.dx(2) = df/dz

--
Anders

> You can't extract vector component from f.dx(0) for example but I thought it was the same as grad(f): f.dx(0)[0] breaks.
>
> Also how do you extract the vector component from graf(f) ayt a given position
>
> mesh = UnitSquare(10,10)
> Q = FunctionSpace(mesh,'CG',1)
>
> expr = Expression('exp(x[0])')
> f = project(expr,Q)
>
> df = grad(f)
> df[0]([0.5,0.5]) = 1.5689657305885802
> df[1]([0.5,0.5]) = 1.5689657305885802 #SHOULD BE ZERO!
>
> df([0.5,0.5]) #breaks with Keyerror
>
>
>

Revision history for this message
Chaffra Affouda (chaffra) said :
#3

Thanks for the reply. So why is it that df/dx = df/dy in my example above. df/dy should be zero?

Revision history for this message
Chaffra Affouda (chaffra) said :
#4

I have been trying to evaluate derivatives directly without relying on project. I have been using the modified version of eval below to compute the n-derivative of a function. It works well in 1D but not in 2D where grad is a vector. Do you have any hint on how I can change this for 2D?

-----modified version of eval where n is the derivative order----
void Function::eval(Array<double>& values,
                    const Array<double>& x,
                    const Cell& dolfin_cell,
                    const ufc::cell& ufc_cell, uint n) const
{
  // Developer note: work arrays/vectors are re-created each time this function
  // is called for thread-safety

  dolfin_assert(_function_space->element());
  const FiniteElement& element = *_function_space->element();

  // Compute in tensor (one for scalar function, . . .)
  const uint value_size_loc = value_size();

  dolfin_assert(values.size() == value_size_loc);

  // Create work vector for expansion coefficients
  std::vector<double> coefficients(element.space_dimension());

  // Restrict function to cell
  restrict(&coefficients[0], element, dolfin_cell, ufc_cell);

  // Create work vector for basis
  std::vector<double> basis(value_size_loc);

  // Initialise values
  for (uint j = 0; j < value_size_loc; ++j)
    values[j] = 0.0;

  // Compute linear combination
  for (uint i = 0; i < element.space_dimension(); ++i)
  {
    element.evaluate_basis_derivatives(i, n, &basis[0], &x[0], ufc_cell);
    for (uint j = 0; j < value_size_loc; ++j)
      values[j] += coefficients[i]*basis[j];
  }
}

Revision history for this message
Chaffra Affouda (chaffra) said :
#5

Bassically i would like f.dx(0)(x0) to give (df/dx)(x=x0)

Revision history for this message
Best Martin Sandve Alnæs (martinal) said :
#6

On 26 June 2012 22:55, Chaffra Affouda
<email address hidden> wrote:
> Question #201507 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/201507
>
>    Status: Answered => Open
>
> Chaffra Affouda is still having a problem:
> Thanks for the reply. So why is it that df/dx = df/dy in my example
> above. df/dy should be zero?

This seems to be behaviour you have implemented?
Evaluating df[0]([0.5,0.5]), I get this:

/home/martinal/opt/fenics/dorsal-unstable/lib/python2.7/site-packages/dolfin/functions/function.pyc
in ufl_evaluate(self, x, component, derivatives)
    227 import numpy
    228 import ufl
--> 229 assert derivatives == () # TODO: Handle derivatives

ufl_evaluate is the function that UFL calls to pass evaluation on to DOLFIN.
UFL can do symbolic evaluation of expressions, but DOLFIN is responsible
for evaluating its own Functions. The derivatives tuple contains any number
of derivative directions, e.g.:

Evaluation of a non-differentiated function gives derivatives=().
Evaluation of f.dx(0) gives derivatives=(0,).
Evaluation of f.dx(1) gives derivatives=(1,).
Evaluation of f.dx(0).dx(1) gives derivatives=(0,1).
etc.

The symbolic evaluation in UFL is slow and mainly used for testing,
but since it's used for testing it is well tested :)
With large meshes it is of course still much faster than project.

With your C++ Function::eval modification you are on the right track!
I suggest you make eval with n=1 compute the whole gradient.
That means at least this line needs to be modified to hold the entire gradient:
  const uint value_size_loc = value_size();
f.ex. (untested):
  uint value_size_loc = value_size();
  for (uint j = 0; j < n; ++j)
    value_size_loc *= x.size();
That should give the right size to the array basis, as expected by this line:
  element.evaluate_basis_derivatives(i, n, &basis[0], &x[0], ufc_cell);

If you post your results as a branch on launchpad
it would be nice to get this into dolfin officially.

Note that UFL will call ufl_evaluate for each derivative direction,
while Function::eval(..., n) should compute all derivatives of
order n since that is what evaluate_basis_derivatives does.
This can be improved in UFL at some later point.

Martin

Revision history for this message
Chaffra Affouda (chaffra) said :
#7

Thanks Martin! It worked like a charm. I pushed my changes to lp:~chaffra/dolfin/main. I also have to make some modifications to dolfin.functions.function.ufl_evaluate that might need more testing but I can now do:

f = Function(Q)
df = grad(f)
df[0]([0.0,0.0]
df[1](0.0,0.0])

Chaffra

Revision history for this message
Chaffra Affouda (chaffra) said :
#8

Thanks Martin Sandve Alnæs, that solved my question.