Compute Gradient and Hessian of power functional?

Question

I need to minimize following energy involving power functional,
$F=\left(\int f(u,\nabla u)d\Omega\right)^2+\int g(u,\nabla u)d\Omega$.
My question is that can the form of Gradient and Hessian be calculated automatically by taking the Gateaux derivative through the function derivative ?

Answer 1

You can calculate derivatives as long as the functional form is linear in the test functions you are using (see https://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html#ad). You can try the following:

gradient = derivative(F, u)
hessian = derivative(gradient,u)

Please let me know if that solves your problem.

Comments

thanks! But I think you do not notice the first term of the functional
$\left(\int f(u,\nabla u)d\Omega\right)^2$, which is the power form.

---

I see the ^2, but I do not see the test functions. You can try the statements above.

---

Just writing the functional (a specific example) as

F = pow(w(alpha) + dot(grad(alpha), grad(alpha)))*dx, 2.0)

ending up with

TypeError: unsupported operand type(s) for ** or pow(): 'Form' and 'float'.

---

Two comments:

• May you please post your code until the expression for F?
• For pow(f,g), the UFL documentation states (at https://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html#basic-nonlinear-functions):

These functions do not accept non-scalar operands or operands with
free indices or Argument dependencies.

The data type Argument represents a basis function on a given finite element (see UFL documentation).

---

from dolfin import *

parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

mesh = UnitSquareMesh(20, 20)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

c = Constant((0.0, 0.0))
r = Constant((0.3, 0.0))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]

du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

d = len(u)
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

Ic = tr(C)
J  = det(F)

E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

Pi = psi*dx

F = derivative(Pi, u, v)

J = derivative(F, u, du)

solve(F == 0, u, bcs, J=J,
    form_compiler_parameters=ffc_options)

file = File("displacement.pvd");
file << u;

If the functional has the form

$\Pi = \int( \psi \, d\Omega)^2$,

how I change the code ?

---

Hi bin.li. Just ran the code you posted above without modifications. It finished in 4 iterations of the Newton solver and produced a result. I am using FEniCS 2016.2.0.