This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Second derivative of a form

0 votes

I would like access to the second derivative of a form. That results in a 3-tensor, which dolfin doesn't want to assemble. Instead it would be enough to access the action of the 3-tensor, i.e., to contract it to a matrix. Unfortunately, action() fails as follows:

from dolfin import *
mesh = UnitSquareMesh(2, 2)
U = FunctionSpace(mesh, "CG", 2)
V = U * U
uf = Function(V)
u, f = split(uf)
v, w = split(TestFunction(V))

a = dot(grad(v), grad(u)) * dx - v * f * dx  # assembly results in a vector
da = derivative(a, uf)  # assembly results in a matrix
d2a = derivative(da, uf)  # assembly fails with: RuntimeError: Unable to create tensors of rank 3.

p = Function(V)
uf.vector()[:] = 1
p.vector()[:] = 1
H = action(d2a, p)

The error message is

/usr/local/lib/python2.7/site-packages/ufl/formoperators.pyc in action(form, coefficient)
    100     form = as_form(form)
    101     form = expand_derivatives(form)
--> 102     return compute_form_action(form, coefficient)
    103 
    104 def energy_norm(form, coefficient=None):

/usr/local/lib/python2.7/site-packages/ufl/algorithms/formtransformations.pyc in compute_form_action(form, coefficient)
    409 
    410     # Pick last argument (will be replaced)
--> 411     u = arguments[-1]
    412 
    413     e = u.element()

IndexError: tuple index out of range

The reason is that the call to expand_derivatives() returns Form([]).

I tried computing d2a = derivative(da, uf, Argument(V, k)) for various values of k without success.

What am I doing wrong? I'm using version 1.6.0.

Thanks!

asked Mar 17, 2016 by dpo FEniCS Novice (150 points)
edited Mar 17, 2016 by dpo

1 Answer

+1 vote

I get the same error message on

H = action(d2a, p)

Replacing it with the following appears to work.

H = derivative(da, uf, p)
answered Mar 18, 2016 by Magne Nordaas FEniCS Expert (13,820 points)

Of course! A directional derivative. Thanks!!

Could you please tell me what exactly the meaning of 'action()' statement is ?
I am a beginner of fenincs. And I would appreciate your help.

...