How to assemble a trilinear form?

General
1

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:

image
image753Ɨ268 73.1 KB

So I guess it should be able to assemble a trilinear form. I will appreciate it if anyone could help me with this.

2

Trilinear forms are not supported. Dolfin supports assembling tensors of order 0,1,2.

I am not aware of any linear algebra package supporting sparse tensors of third order (hopefully someone might know)

3

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 :wink:

1 Like
4

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.

5

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.

6

You can always create these in the way Francesco suggested, i.e.

from mpi4py import MPI
import dolfinx
import ufl

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

phi_i = ufl.TrialFunction(V)
phi_j = ufl.TestFunction(V)
phi_k = dolfinx.fem.Function(V)


a = ufl.inner(phi_k*ufl.grad(phi_i), ufl.grad(phi_j)) * ufl.dx

# general A to assemble into
a_compiled = dolfinx.fem.form(a)
A = dolfinx.fem.create_matrix(a_compiled)
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
local_dof_range = V.dofmap.index_map.local_range * V.dofmap.index_map_bs
A_kij = []
print(num_dofs_global)
for i in range(num_dofs_global):
    phi_k.x.array[:] = 0.0
    if i >= local_dof_range[0] and i < local_dof_range[1]:
        phi_k.x.array[i-local_dof_range[0]] = 1.0
    phi_k.x.scatter_forward()
    A.set_value(0)
    A_kij.append(dolfinx.fem.assemble_matrix(A, a_compiled).to_scipy())
print(A_kij)

But this will eat up your RAM for any sufficiently large problem.

7

Thanks a lot! I will try it.

8

Hi @dokken,

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!

9

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

  1. create an argument (in the space VE) to be passed as third input
  2. once you get dKdt, use ufl.replace to replace that argument with dolfinx.fem.Function
  3. 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 :wink:

1 Like
10

It works perfectly for my case! Thanks a lot :slight_smile: