Hi everyone,
I am trying to solve NS equation and need to get a 3-tensor induced by the nonlinear term u\cdot \nabla u. In other words, I want to get a tensor A defined by A_{i j k} = \int_\Omega (\phi_i \cdot \nabla) \phi_j \cdot \phi_k d x, where \phi_i-s are the basis functions.
There are some similar questions before but it seems that there are never a solution to this problem. But those questions are mostly 8 or 10 years ago. So I just want to take a chance to see whether there are some new solutions.
I have tried to import āArgumentā but failed. Here is the code, for a simple trilinear form:

=======================================

from dolfin import *
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, āPā, 1)
u = Argument(V, 0) # Trial function argument
v = Argument(V, 1) # Test function argument
w = Argument(V, 2) # Another test function argument
m = u * v * w * dx # Define a variational problem
m_mat = assemble(m)

=======================================
This does not work because the function āassembleā seems unable to deal with trilinear form. However, in the documentation, it says:

It would be good if @yangchanghe gave us more context about what they are trying to do with this trilinear form. In principle, you could have w be a dolfin.Function which only has a single DOF equal to one and assemble that bilinear form, then repeat that for all possible w. However, do not try this unless you are sure that it really is what you need, because it will likely eat up all your RAM

I think scikit fem can assemble trilinear form. But I met another problem with scikit fem. So I hope fenics can do that because I am more familiar with fenics.

I want to export this tensor and matrices induced by other terms in NS equation. Then I need to solve an optimization problem, which I have to use some tools in matlab.
Of course I can also use matlab.engine in every iteration, but I believe it would be very difficult to debug so this is my last choice.

Thank you for the helpful code example. I have a similar problem as @yangchanghe for assembling a trilinear form, but a little bit different. Instead of having a written trilinear variational form, I would like to take second derivative of an energy form for linear elasticity to get the sensitivity of the stiffness matrix with respect to the material property. I know it would be crazy costly, but this is required to couple my solver with external optimization paradigm.

Here is a minimal example of the problem, which is adapted from your hyperelasticity tutorial (Hyperelasticity ā FEniCSx tutorial). Iām having a hard time trying to adapt this to the method you showed in the code example.

from mpi4py import MPI
import dolfinx
import ufl
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 1, 1, 1)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, )))
u = dolfinx.fem.Function(V)
v = ufl.TestFunction(V)
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# Elasticity parameters
VE = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
E = dolfinx.fem.Function(VE)
E.x.array[:] = 1.0E4
nu = 0.3
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
P = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I
F = ufl.inner(ufl.grad(v), P) * ufl.dx
K = ufl.derivative(F, u)
dKdt = ufl.derivative(K, E)
# Assemble stiffness matrix
K_compiled = dolfinx.fem.form(K)
K_mat = dolfinx.fem.assemble_matrix(K_compiled)
# Assemble sensitivity of stiffness matrix wrt material property
#******************** Want to have ***********************
dKdt_compiled = dolfinx.fem.form(dKdt)
dKdt_mat = dolfinx.fem.assemble_matrix(dKdt_compiled)
#*********************************************************

As expected, I got an error in the last line of code saying āRuntimeError: Cannot create sparsity pattern. Form is not a bilinear formā. Any ideas will be appreciated!

If you look at help(ufl.derivative), youāll see that the signature is derivative(form, coefficient, argument=None, coefficient_derivatives=None). You could:

create an argument (in the space VE) to be passed as third input

once you get dKdt, use ufl.replace to replace that argument with dolfinx.fem.Function

proceed as in the previous post.

Of course, this only makes sense if dKdt is actually multilinear. Continuum mechanics is not my field, but I would be worried about the ln(J) term