I am confused by second derivatives with fenics. I would like to compute the tensor field associated to the hessian of a vector field.
I started with the matrix field corresponding to the hessian of a scalar field and tried 2 different approaches (see code below):
- by projecting grad(grad(u)),
- by projecting u.dx(i).dx(j),
but none of them returns a good result.
Is there any way to get the proper field in the scalar case?
Thanks!
from __future__ import (division, print_function)
import dolfin as df
import sympy as sy
import sympy.printing as syp
pt = 0.5, 0.5
x, y = sy.symbols('x y')
uexpr = x ** 3 + y ** 3
jacexpr = [sy.diff(uexpr, s) for s in [x, y]]
hessexpr = [sy.diff(ex, s) for ex in jacexpr for s in [x, y]]
dfexpr = syp.ccode(uexpr).replace('x', 'x[0]').replace('y', 'x[1]')
print('Sympy expressions:')
print(uexpr)
print(jacexpr)
print(hessexpr)
print('\nDolfin expression:')
print(dfexpr)
print('\nANALYTIC SOLUTION:')
print('u(pt) =', uexpr.subs([(x, pt[0]), (y, pt[1])]))
print('jac(pt) =', [ex.subs([(x, pt[0]), (y, pt[1])]) for ex in jacexpr])
print('hess(pt) =', [ex.subs([(x, pt[0]), (y, pt[1])]) for ex in hessexpr])
expr = df.Expression(dfexpr)
n = 20
mesh = df.UnitSquareMesh(n, n)
deg = 3
Vu = df.FunctionSpace(mesh, 'Lagrange', deg)
umesh = df.project(expr, Vu)
Vg = df.VectorFunctionSpace(mesh, 'Lagrange', deg - 1)
jacmesh = df.project(df.grad(umesh), Vg)
Vh = df.TensorFunctionSpace(mesh, 'Lagrange', deg - 2)
hessmesh = df.project(df.grad(df.grad(umesh)), Vh)
print('\nPROJECTION:')
print('u(pt) =', umesh(pt))
print('jac(pt) =', jacmesh(pt))
print('hess(pt) =', hessmesh(pt))
Vgc = df.FunctionSpace(mesh, 'Lagrange', deg - 1)
Vhc = df.FunctionSpace(mesh, 'Lagrange', deg - 2)
jaccomp = [df.project(umesh.dx(i), Vgc) for i in range(2)]
hesscomp = [df.project(umesh.dx(i).dx(j), Vhc)
for i in range(2) for j in range(2)]
print('\nCOMPONENTS:')
print('jac(pt) =', [j(pt) for j in jaccomp])
print('hess(pt) =', [h(pt) for h in hesscomp])
