Assembly of product of 3 basis functions

Question

Hello everyone,

I am interested in assembling a product of 3 FEM basis functions, namely

$$\int_{\Omega} \phi_i * \phi_j * \frac{\partial \phi_k}{\partial x} dX $$

for all possible combinations of i, j and k. Is there a simple way to do this?

I have tried several approaches like

test1 = assemble(inner(inner(phi_j, phi_k.dx(1)), phi_i))

or

test2 = assemble(inner(phi_j * phi_k.dx(1), phi_i))

but I always get the error

*** Error:   Unable to creating dolfin.Form.
*** Reason:  Expected a ufl.Form, ufc.form or a dolfin.Form.

Best regards,
Sebastian

Answer 1

I'm not sure what you are trying to achieve here (build the tensor associated to the trilinear form?), but you're missing the measure, FEniCS needs to know over which domain you're integrating over. Are you meaning this ?

 assemble(inner(inner(phi_j, phi_k.dx(1)), phi_i) * dx)

Also, this looks like a duplicate of this question.

Comments

Hi,

first of all, thanks for the quick answer. What I am trying to achieve ist to build a Matrix K \in \mathbb{R}^{N,N} with the entries looking like

K[k][l] = "sum over i: inner(phi_k * phi_l.dx(1) * phi_i)" 
(phi_k being the k-th basis function)

for Model Order Reduction purposes. My question now is how to access the individual basis functions and then assemble the entries of my Matrix K. Right now I am thinking about inserting the k-th unit vector to the nodes and get acces to the k-th basis function this way. But I am not sure wether this is the right way to do and also it seems really costly to me. I have been working with finite volume codes until just recently and am quite new to the finite element method.

---

Hi again,

I think I found a solution, but it is terribly costly.

mesh = Mesh(meshfile)
N = mesh.num_vertices()

Q = FunctionSpace(mesh, "Lagrange", 1)

phi_i = Function(Q)
phi_k = Function(Q)
phi_l = Function(Q)

Dphi_lDx = project(phi_l.dx(0), Q)

e_vector = np.zeros(shape=(N))
K = np.zeros(shape=(N, N))

for k in range(N):
    e_vector[k] = 1.0
    phi_k.vector().set_local(e_vector)
    e_vector[k] = 0.0

    for l in range(N):
      e_vector[l] = 1.0
      Dphi_lDx.vector().set_local(e_vector)
      e_vector[l] = 0.0

      for i in range(N):
        e_vector[i] = 1.0
        phi_i.vector().set_local(e_vector)
        e_vector[i] = 0.0

        K[k][l] += assemble(inner(inner(phi_k, Dphi_lDx), phi_i) * dx)

Can someone tell me if my idea is correct and if there is possibly a better way to implement this?

Best regards,
Sebastian

---

The tensor you're trying to assemble is $i$-$j$-symmetric, so no need to do the full loop over $j$. Then, FEniCS can assemble matrices (for fixed $k$) for you, so let's use that, and let's note $\mathcal I_{i,j}[k] = \int_\Omega \phi_i \phi_j \frac{\partial \phi_k}{\partial x_1}$. This is the derivative w.r.t. to second variable, as originally written in your question.

phi_i = TrialFunction(Q)
phi_l = TestFunction(Q)

phi_k = Function(Q)
e_vector = np.zeros(shape=(N))

K = []

for k in range(N):
    e_vector[k] = 1.0
    if k > 0:
        e_vector[k] = 0.0

    # not sure if really needed after first iteration
    phi_k.vector()[:] = e_vector

    K.append(assemble(inner(inner(phi_j, phi_k.dx(0)), phi_i) * dx))

Then K will be a list of FEniCS matrices corresponding to the $\mathcal I_{i,j}[k]$ for fixed $k$. Assembling the tensor in a 3-dimenional numpy object may or may not possible depending on the backend you're using.

---

Hello gjankowiak,

thanks alot for your help, that worked. Since I'm only looking for nonzero entries, I simply wrote all i j k + value into a table.

Thanks again!

Best regards,
Sebastian

Answer 2

Hi again,

I'm still struggling with my problem but after putting more thought into it I can now hopefully clarify what I am looking for.

Basically, I want to assemble the rank-3 tensor C_{ijk}, i,j,k \in \mathbb{N}_{+}^N (N = dof),
whose elements are computed by

\int_{Omega} \phi_i * \phi_j * \frac{\partial \phi_k}{\partial x} dX.

The approach I posted as a comment comments seems way too complicated on the one hand I am also not sure wether it delivers the correct solution, so I am hoping that someone can help me out.

I started looking into the FFC documentation but it seems to be outdated. I would be grateful for any kind of advice.

Best regards,
Sebastian

Comments

What do you need this tensor for? There is a reason why FEniCS does not support rank-3 assembly (at least to my knowledge)... it is rarely needed. If you really need this, then look at gjankowiak's answer. But I cannot see any sense in what you wrote above, i.e., that you need the sum over i of inner(phi_k * phi_l.dx(1) * phi_i). You know that the sum of all phi_i is exactly one?

---

I need this tensor for a reduced order model of a Navier-Stokes problem I'm solving. By inserting a galerkin ansatz into the NSE, I get an ODE model for the time dependent coefficients of my basis functions which I created via POD. Also, I am not summing over all i but over a product of phi_i times an element of my POD basis vector, so it's not always 1.