Extracting nonlinear stiffness tensors

The stiffness matrix of a system can be seen as the second derivative of the energy with respect to the node displacements. Once the stiffness matrix has been calculated, one can easily formulate time-dependent problems as

M\ddot{x} + Kx = F

Computing the stiffness matrix K is very easy in FEniCSx (see MWE below).

To speed up nonlinear calculations, I would like to pre-compute a higher derivative of the energy, i.e. a nonlinear tensor (Q) that is the fourth derivative of the energy wrt the node displacements, to make a problem of the form:

M\ddot{x} + Kx + Qxxx = F

However, I’m not able to find any examples of such calculation. Is this even possible in FEniCS?

F = ufl.variable(I + ufl.grad(u))
E = ufl.variable(0.5*(ufl.dot(F.T, F) - I))
TrE = ufl.variable(ufl.tr(E))
TrE2 = ufl.variable(ufl.tr(ufl.dot(E,E)))
psi = 0.5*lmbda*TrE*TrE + mu*TrE2
P = ufl.diff(psi, F)
metadata = {"quadrature_degree": 4}
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
variationalForm = ufl.inner(P, ufl.grad(v))*dx 
problem = fem.petsc.NonlinearProblem(variationalForm, u, bcs)
stiffnessMatrix = fem.petsc.assemble_matrix(problem.a, bcs=bcs, diagonal=1e18)